Post Reply 
HP Prime complex division problem in CAS
11-27-2020, 11:05 PM
Post: #21
RE: HP Prime complex division problem in CAS
(11-27-2020 08:30 AM)lyuka Wrote:  I found relatively simple workaround for complex division in CAS,
that is to apply scaling process to Smith's method (1962).

Unlike IEEE double, HP Prime does not support gradual underflow.
It does not do flush-to-zero either.

CAS> MINREAL/2             → 2.2250738585e−308
CAS> MINREAL/2^47       → 2.2250738585e−308
CAS> MINREAL/2^48       → 0.

With this underflow behavior, it may be best to scale up, to some safe maximum value.

Smith's method, with u = IM(y)/RE(y), |u| ≤ 1 : x/y = x * (1 - u*i) / (RE(y) + u*IM(y))

As long as x, y parts at or below MAXREAL/2, RHS numerator and denominator will not overflow.

x/y = cdiv_smith(x*s, y*s), where s = (MAXREAL/2) / max(|RE(x)|,|IM(x)|,|RE(y)|,|IM(y)|)

Example, with g = MAXREAL/2:

(2g + 2g*i) / (g + 2g*i)
= (g + g*i) / (0.5*g + g*i)                                        // scaled by s = g / (2g) = 0.5
= conj((g + g*i) / (g + 0.5*g*i))                               // flip to force u=0.5, |u| ≤ 1
= conj((g + g*i) * (1 - 0.5*i)) / (g + 0.5*(0.5*g))      // smith method
= (1.5*g - 0.5*g*i) / (1.25*g)
= 1.2 - 0.4*i
Find all posts by this user
Quote this message in a reply
11-28-2020, 02:00 AM
Post: #22
RE: HP Prime complex division problem in CAS
For IEEE double, we cannot scale to DBL_MAX/2 = 0x1.fffffffffffffp+1022
With round-to-nearest, halfway-to-even, it is possible scaled parts is slightly bigger.
To ensure no overflow, we reduced it a bit, to say, 0x1p+1022

Also, to flip the parts, we can use x/y = (x*i)/(y*i)
Here is compdiv_improved, with scaling.

Code:
function cdiv(xr, xi, yr, yi)
    local s = 0x1p1022 / max(abs(xr),abs(xi),abs(yr),abs(yi))
    xr, xi, yr, yi = xr*s, xi*s, yr*s, yi*s -- scaling
    if abs(yr) < abs(yi) then xr,xi,yr,yi = -xi,xr,-yi,yr end
    local u = yi / yr
    local v = 1/(yr + u*yi)
    if u ~= 0 then return (xr + xi*u)*v, (xi - xr*u)*v end
    return (xr + xi/yr*yi)*v, (xi - xr/yr*yi)*v
end

lua> cdiv(1e154, 2e154, 2e154, 1e154)
0.8             0.6000000000000001
lua> cdiv(1e307, 1e-307, 1e204, 1e-204)
1e+103      -1.0000000000000001e-305

Example from improved_cdiv_v0.3.pdf, "4.7 A failure of the robust algorithm", page 35.

lua> cdiv(0x1p-912, 0x1p-1029, 0x1p-122, 0x1p46)
5e-324      -4.1045368012983762e-289

With scaling, we get the correctly rounded parts. Smile
Without scaling, we lost the real part, same as compdiv_robust
Find all posts by this user
Quote this message in a reply
11-28-2020, 01:13 PM (This post was last modified: 11-28-2020 02:08 PM by Albert Chan.)
Post: #23
RE: HP Prime complex division problem in CAS
(11-27-2020 11:05 PM)Albert Chan Wrote:  x/y = cdiv_smith(x*s, y*s), where s = (MAXREAL/2) / max(|RE(x)|,|IM(x)|,|RE(y)|,|IM(y)|)

I just noticed s itself might overflowed, if denominator is small.
We should do scaling in steps:

m = MAXREAL * MINREAL ≈ 2^(1024-1022) = 4.0
s1 = m / max(|RE(x)|,|IM(x)|,|RE(y)|,|IM(y)|)
s2 = (MAXREAL/2)/m ≈ MAXREAL/8

x/y = cdiv_smith(x*s1*s2, y*s1*s2)

Same issue with previous post Lua code.
To avoid rounding errors, I scaled parts with pow-of-2.

I also stored d = 1/v instead, to avoid v possibly losing precision.
(v = 1/d turned subnormal if 2^1022 < |d| < 2^1024)

Code:
logb = require'mathx'.logb

function cdiv(xr, xi, yr, yi)
    local e = 1022 - max(logb(xr),logb(xi),logb(yr),logb(yi))
    xr,xi,yr,yi = ldexp(xr,e),ldexp(xi,e),ldexp(yr,e),ldexp(yi,e)
    if abs(yr) < abs(yi) then xr,xi,yr,yi = -xi,xr,-yi,yr end    
    local u = yi / yr
    local d = yr + u*yi
    if u ~= 0 then return (xr + xi*u)/d, (xi - xr*u)/d end
    return (xr + xi/yr*yi)/d, (xi - xr/yr*yi)/d
end


Example from improved_cdiv_v0.3.pdf, "Figure 3: A collection of difficult complex divisions", page 23.
Note: examples below have s ≥ 2^1024; shifting binary exponents avoided the problems.

lua> cdiv(0x1p-347, 0x1p-54, 0x1p-1037, 0x1p-1058) -- case 7
3.8981256045591133e+289      8.174961907852354e+295
lua> cdiv(0x1p-1074, 0x1p-1074, 0x1p-1073, 0x1p-1074) -- case 8
0.6                                          0.2
lua> cdiv(0x1p-622, 0x1p-1071, 0x1p-343, 0x1p-798) -- case 10
1.0295115178936058e-084       6.971459875150762e-220
Find all posts by this user
Quote this message in a reply
11-30-2020, 04:52 PM
Post: #24
RE: HP Prime complex division problem in CAS
(11-28-2020 01:13 PM)Albert Chan Wrote:  I just noticed s itself might overflowed, if denominator is small.
We should do scaling in steps:
...

Perhaps multi-step scaling is overkill. Limiting s ≤ MAXREAL work too.

s = (MAXREAL/2) / max(0.5, |RE(x)|,|IM(x)|,|RE(y)|,|IM(y)|)
x/y = cdiv_smith(x*s, y*s)

Patching similar idea to previous cdiv Lua code:
Code:
< xr,xi,yr,yi = ldexp(xr,e),ldexp(xi,e),ldexp(yr,e),ldexp(yi,e)

> e = 2 ^ min(e,1023)
> xr,xi,yr,yi = xr*e, xi*e, yr*e, yi*e
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)