Ellipsoid surface area

[Albert Chan]

Extended RFD(x,y,z) for complex numbers, replaced comparison test with equality (with fuzzy factor)

Code:

RFD(x,y,z)  // return RF(x,y,z), RD(x,y,z)
BEGIN
    LOCAL rx,ry,rz,k;
    rx := sqrt(approx(x));
    ry := sqrt(approx(y));
    rz := sqrt(approx(z));
    k := rx*(ry+rz) + ry*rz;
    rx := (x+y+z-k)*.01; // eps
    ry := (rx+k-k == 0); // x,y,z close?
//  print(k, k+20*rx, k+100*rx);
    IF ry THEN k:=sqrt(3/(k+20*rx)); return [k, k^3] END 
    rx, ry := RFD((x+k)/4, (y+k)/4, (z+k)/4);
    return [rx, ry/4 + 3/(rz*(z+k))];
END;

Previously, recursed arguments quadrupled in size.
Arguments (x,y,z) now scaled to be similar in size.

If we uncomment print(...) code, we get (x,y,z), convergence trend.
If x==y==z, then (k = √(xy)+√(yz)+√(xz)) == (x+y+z)

CAS> ellipsoid_area(2,4,9)      → 283.427384268

Code:

k              k+(x+y+z-k)/5  (x+y+z)
3.33333333333  3.70617283951  5.1975308642
3.65931727734  3.68733036508  3.79938271605
3.68498213244  3.68685243336  3.69433363702
3.68672507777  3.68684406393  3.68732000859
3.68683645911  3.68684392938  3.68687381047
3.68684345985  3.68684392727  3.68684579695
3.68684389801  3.68684392723  3.68684404412
3.68684392541  3.68684392723  3.68684393454
3.68684392712  3.68684392723  3.68684392769
3.68684392723  3.68684392723  3.68684392726
3.68684392723  3.68684392723  3.68684392723
3.68684392723  3.68684392723  3.68684392723

k seems to converge 4x faster, than (x+y+z), but from opposite side.
Extrapolate this, we have (k + (x+y+z-k)/5) converge at doubled rate.

Trivia: HP Prime float equality is really approximately equal (about 12 digits)

CAS> 9*(1+1e-12)            → 9.00000000001
CAS> Ans>9., Ans==9.      → [1, 1]

If our convergence test were "rx+k == k", we may lose needed recursions.
Test is written as "rx+k-k == 0", to ensure what we meant, bits equality.