HP Prime complex division problem in CAS

[Albert Chan]

(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

[Albert Chan]

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.
Without scaling, we lost the real part, same as compdiv_robust

[Albert Chan]

(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

[Albert Chan]

(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