HP Prime complex division problem in CAS

[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