(/ 1.0 0+inf.0i) produces NaN, which also affected result of atan and atanh

[LdBeth]

I was investigating a bug in arctan implementation used in J

And I used Chez Scheme and Clozure CL for comparison, and find there is also a bug (but a different one) in Chez.

> (atan +inf.0+0.0i)
1.5707963267948966+nan.0i
> (atan +inf.0+1.2i)
1.5707963267948966+nan.0i
> (atanh 0+inf.0i)
+nan.0+1.5707963267948966i
> (atanh 0+1.0e18i)
0.0+1.5707963267948966i

That Clozure CL gives correct result

CL-USER> (atanh #C(0D0 1D++0))
#C(0.0D0 1.5707963267948966D0)

Where the issue itself is in https://github.com/cisco/ChezScheme/blob/8c5f5cc27cf8f2b432aafc4f2fc640435294b74a/s/5_3.ss#L556

since clfatan calls clfatanh to do the actual job.

I know Chez copied algorithm from Kahan's paper, I checked and seems the bug is not from the original algorithm.

[LdBeth]

The real cause is

> (/ 1.0 0+inf.0i)
0.0+nan.0i

Which should be just 0.0.

[LdBeth]

https://github.com/cisco/ChezScheme/blob/8c5f5cc27cf8f2b432aafc4f2fc640435294b74a/s/library.ss#L217

(let ([c ($inexactnum-real-part y)] [d ($inexactnum-imag-part y)])
          (let ([t (fl/ x (fl+ (fl* c c) (fl* d d)))])
             (fl-make-rectangular (fl* c t) (fl- (fl* d t)))))]

with d is +inf.0 this would cause (fl* d t) becomes (fl* +inf.0 0.0) thus nan.

By symmetry of the algorithm,

> (/ 1.0 +inf.0+1.0i)
+nan.0-0.0i

This is also wrong.

The else branch is also using a similar equation and if inf is involved is would also give wrong result.

> (/ 1.0+0.0i +inf.0+1.0i)
+nan.0+nan.0i

[LdBeth]

The Algorithm 116 is the more robust division algorithm and is the one used in Clozure CL: https://dl.acm.org/doi/pdf/10.1145/368637.368661

[LdBeth]

Unfortunately, cfltanh also has a bug!

This part from cfltanh corresponding to Re[1/(x+y*i)], would also produce nan if y is +inf.0

(cond
 [(fl> x ay) (fl/ (fl+ x (fl* (fl/ y x) y)))]
 [(fl< x ay) (let ([r (fl/ y x)])
               (fl/ r (fl+ (fl* x r) y)))]
 [else (fl/ (fl+ x ay))])