Define value for `_3&o.` (arctan) when argument is a complex number with infinity components

[LdBeth]

Currently _3 o. _ (arctan(infinity)) gives 1.5708 and _3 o. __ gives _1.5708, however _3 o. 0j_ or _3 o. 0j__ would give NaN error (until J9.4).

Experimenting with other programming languages that has complex float and atan function such as Clozure Common Lisp or Chez Scheme gives these results and they are pretty much agree with each other if the imaginary part has infinity:

CL-USER> (ccl:set-fpu-mode :division-by-zero nil)
CL-USER> (atan #C(-432.0 1E++0))
#C(-1.5707964 0.0)
CL-USER> (atan #C(+432.0 1E++0))
#C(1.5707964 0.0)
CL-USER> (atan #C(+432.0 -1E++0))
#C(1.5707964 -0.0)
CL-USER> (atan #C(+432.0 +1E++0))
#C(1.5707964 0.0)

> (atan -432.0+inf.0i)
-1.5707963267948966+0.0i
> (atan 432.0+inf.0i)
1.5707963267948966+0.0i

But when the real part is infinity they are indifferent on what is NaN:

CL-USER> (ccl:set-fpu-mode :invalid nil)
6528 (13 bits, #x1980)
CL-USER> (atan #C(+1E++0 1.32))
#C(1.5707964 0.0)
CL-USER> (atan #C(+1E++0 -1.32))
#C(1.5707964 -0.0)
CL-USER> (atan #C(-1E++0 -1.32))
#C(-1.5707964 -0.0)
CL-USER> (atan #C(+1E++0 +1E++0))
#C(1.5707964 1E+-0 #| not-a-number |#)
CL-USER> (atan #C(+1E++0 -1E++0))
#C(1.5707964 1E+-0 #| not-a-number |#)

> (atan +inf.0-1.2i)
1.5707963267948966+nan.0i
> (atan +inf.0+1.2i)
1.5707963267948966+nan.0i
> (atan +inf.0+inf.0i)
1.5707963267948966+0.0i
> (atan +inf.0-inf.0i)
1.5707963267948966-0.0i
> (atan -inf.0-inf.0i)
-1.5707963267948966-0.0i

[LdBeth]

Version 9.4 gives an error

   _3 o. 0j_
|NaN error, executing dyad o.
|you have calculated the equivalent of _-_ or _%_
|   _3     o.0j_

Which makes curious about the internal algorithm used, since Clozure CL can compute (atan #C(0 -1E++0)) or (atan #C(0 +1E++0)) without raise floating point exceptions in default setting while (- +1E++0 +1E++0) or (/ +1E++0 +1E++0) would still raise a FLOATING-POINT-INVALID-OPERATION exception.

[moon-chilled]

Will take a look, thanks. Rubberstamp is in c99/c11§G.6 (I could have sworn there was a more authoritative source for these but can't find anything in ieee 754). The j interpreter uses some self-implemented algorithms which simply happened not to consider this case, as it predates the introduction of complex math routines to c; clozure and chez are almost certainly punting to the c standard library, and it might be a good idea for us to do so too.

[LdBeth]

IEEE 754 does not tell much thing about complex numbers indeed.

Chez defines arctan(z) to be arctanh(z*i)/i when z is a complex number, where arctanh(z)=(log(1+z) - log(1-z))/2 (there is a special case of z has infinity).

This means that, _7 o. _ and _7 o. __ should both be defined to be equal to j. o. 0.5 but they would also raise NaN error in J.

The source code of Chez Scheme referenced W. Kahan's "Branch Cuts for Complex Elementary Functions" as the source of definitions used to compute them. This is also referenced by Common Lisp the Language 2nd Ed for the formulas.

[LdBeth]

A interesting discovery: Kahan's formulas are based the work done when complex values are introduced to APL: https://dl.acm.org/doi/10.1145/390007.805368

[moon-chilled]

Indeed; penfield himself weighs in on the issue, a couple of years ago.

[LdBeth]

It seems J uses the straight up definition for atan and atanh that has the same branch cut to Kahan’s

jsource/jsrc/vz.c

Line 208 in dee60e7

static ZF1(jtzatanh){R ztymes(zrj0((D)0.5),zlog(zdiv(zplus(z1,v),zminus(z1,v))));}

either we could adopting Kahan’s implementation (which rises the issue of checking of other elementary functions) or simply add the case for infinities like what already been done to handle _ and __

Or maybe as you said earlier, adhere to what C and C++ provides so can benefit from modern computer arch optimizations.

[LdBeth]

Actually, the current J atan/atanh implementation ignores the sign of real part when the complex part goes infinity. While in Chez or CCL the sign of real component is use for the result:

> (atan 0.0+1e18i)
1.5707963267948966+0.0i
> (atan -0.0+1e18i)
-1.5707963267948966+0.0i

   _3 o. _1.0j_1e18
1.5708
   _3 o. 1.0j_1e18
1.5708

This also comes to the notorious signed zero problem. I know J does not differentiate pos zero and neg zero. There is really a design choice need to be made for the problem.

[moon-chilled]

There is no reason not to simply follow c99.

differentiate pos zero and neg zero

   _ -: % _1 * %_
0

[LdBeth]

Ah, that's neat, good to know about that. My bad on presuming the repl would display all the necessary details.

[LdBeth]

Closing as been fixed in j9.5.0-beta8