HP Prime complex division problem in CAS

[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