HP Prime complex division problem in CAS

[lyuka]

May future firmware at least pass the next simple division test.

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); 
    z := (x * g) * inv(y * g); PRINT(k+":*inv-> "+z); 
    z := (x * g) / (y * g); //PRINT(k+":#/# -> "+z); // stall
  END;
END;
#END

This is the execution result. The calculation result of (x * g) / (y * g) cannot be printed because the calculator terminates abnormally when trying to print undef-undef * i.


Anyway, I ported Robert L. Smith and Michael Baudin's "Robust Complex Division".

Code:

#CAS
// 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).
// An implimentation of the following algorithm by Takayuki HOSODA (aka Lyuka).
// "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) returns complex division x / y.
cdiv(x, y):=
BEGIN
  LOCAL z;
  LOCAL r, t, AB, CD, B, S;
  LOCAL OV2, eps, UNBeps, Be;
  AB := fmax(fabs(RE(x)), fabs(IM(x)));
  CD := fmax(fabs(RE(y)), fabs(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 cabs(IM(y)) <= cabs(RE(y)) THEN
    r := IM(y) / RE(y); 
    t := 1 / (RE(y) + IM(y) * r);
    IF r == 0 THEN
      z :=  (RE(x) + IM(y) * (IM(x) / RE(y))) * t + (IM(x) - IM(y) * (RE(x) / RE(y))) * t * i;
    ELSE
      z :=  (RE(x) + IM(x) * r) * t + (IM(x) - RE(x) * r) * t * i;
    END;
  ELSE
    r := RE(y) / IM(y); 
    t := 1 / (RE(y) * r + IM(y));
    IF r == 0 THEN
      z :=   (RE(y) * (RE(x) / IM(y)) + IM(x)) * t + (RE(y) * (IM(x) / IM(y)) - RE(x)) * t * i;
    ELSE
      z :=   (RE(x) * r + IM(x)) * t + (IM(x) * r - RE(x)) * t * i;
    END;
  END;
  return S * z;
END;
//fmax(a, b) -- returns the larger value from a and b.
fmax(a, b):=
BEGIN
  IF a < b THEN
    return b;
  ELSE
    return a;
  END;
END;
//fabs(a) -- returns the absolute value of 'a'.
fabs(a):=
BEGIN
  IF a < 0 THEN
    return -a;
  ELSE
    return a;
  END;
END;
// This is a workaround for the dynamic range problem of complex division and abs() in CAS.
// Conforms to HP Prime software version 2.1.14433 (2020 01 21).
// cabs(z) -- returns the magnitude of the complex number z.
cabs(z):=
BEGIN
  LOCAL x, y, t, s;
  x := RE(z);
  y := IM(z);
  IF x < 0 THEN x := -x; END;
  IF y < 0 THEN y := -y; END;
  IF x == 0 THEN return y; END;
  IF y == 0 THEN return x; END;
  IF (y > x) THEN
    t := x; x := y; y := t;
  END;
  t := x - y;
  IF t > y THEN
    t := x / y;
    t := t + sqrt(1 + t * t);
  ELSE
    t := t / y;
    s := t * (t + 2);
    t := t + s / (sqrt(2) + sqrt(2 + s));
    t := t + sqrt(2) + 1;
  END;
  return(x + y / t);
END;
#END

Have fun.
Lyuka