HP Prime complex division problem in CAS

[lyuka]

I found relatively simple workaround for complex division in CAS,
that is to apply scaling process to Smith's method (1962).

Code:

#CAS
// cdiv(x, y) returns complex division x / y. 
// written by Takayuki HOSODA (aka Lyuka) Rev. 0.3.2 (Nov. 25 2020)
// This is a workaround for the dynamic range problem of complex division in CAS.
// Conforms to HP Prime software version 2.1.14433 (2020 01 21).
// This cdiv() uses Smith's method (1962) with scaling process in the following paper.
// "Improved Complex Division in Scilab", Michael Baudin, Robert L. Smith, Jun. 2011
// http://forge.scilab.org/index.php/p/compdiv/downloads/get/improved_cdiv_v0.3.pdf
cdiv(x, y):=
BEGIN
  LOCAL z;
  LOCAL r, t, AB, CD, B, S;
  LOCAL OV2, eps, UNBeps, Be;
  AB := MAX(abs(RE(x)), abs(IM(x)));
  CD := MAX(abs(RE(y)), abs(IM(y)));
  B := 2;
  S := 1;
  // Machine dependent constants (for HP Prime)
  OV2 :=  8.988465674311581E307;  // MAXREAL / 2; 
  eps := 7.105427357601002E-15;    // 1 / 2^47; 
  UNBeps := 6.263026125028041E-294; // MINREAL * B / eps;
  Be := B / (eps * eps);
  IF (AB >= OV2)   THEN x = x / 2;  S = S * 2;  END; // Scale down x
  IF (CD >= OV2)   THEN y = y / 2;  S = S / 2;  END; // Scale down y
  IF (AB < UNBeps) THEN x = x * Be; S = S / Be; END; // Scale up   x
  IF (CD < UNBeps) THEN y = y * Be; S = S * Be; END; // Scale up   y
  IF abs(RE(y)) < abs(IM(y)) THEN
    r := RE(y) / IM(y); 
    t := (RE(y) * r  + IM(y));
    z := (RE(x) * r + IM(x)) / t  + (IM(x) * r - RE(x)) / t * i;
  ELSE
    r := IM(y) / RE(y); 
    t := (RE(y) + IM(y) * r);
    z := (RE(x) + IM(x) * r) / t  + (IM(x) - RE(x) * r) / t * i;
  END;
  return S * z;
END;
#END

This cdiv shown above has passed McLaren's difficult complex division test shown below.

Code:

#CAS
// from Mc Laren's difficult complex divisions test
cdiv_test():=
BEGIN
  LOCAL c, x, y, z;
  LOCAL g, k;
  PRINT();
  g := MAXREAL / 2;
  FOR k FROM 0.0 TO 2.0 STEP 0.5 DO
    x := 1 + i;
    y := 1 + i * k;
    c := x / y;
    z := (x * g) / (y * g); PRINT(k + ":Ref.-> "+c); 
    z := cdiv(x *g, y * g); PRINT(k+":cdiv-> "+z); 
  END;
END;
#END

For more information about complex division and implimentation to C99,
The "compdiv_robust" algorithm rewritten for C99 And a comparison with Smith's method (1962) using the scaling process used in it.