#4737 atan2 with bfloat arg(s) is inconsistent on negative x-axis

[Kris Katterjohn]

See also bug [#4643] where Stavros initially noted this inconsistency in passing.

atan2 with bfloat args is inconsistent (wrong) on the negative x-axis:

(%i1) assume (equal (y, 0), x < 0)$

(%i2) atan2 (y, x);
(%o2) %pi

(%i3) atan2 (0, -1);
(%o3) %pi

(%i4) atan2 (0.0, -1.0);
(%o4) 3.141592653589793

(%i5) atan2 (0.0b0, -1.0b0); /* wrong */
(%o5) -3.141592653589793b0

In particular, compare %o4 and %o5 to the atan2 result using a negative float zero -0.0 under lisp implementations that support signed zeros:

(%i6) atan2 (-0.0, -1.0);
(%o6) -3.141592653589793

The atan2 results with float args are consistent with the lisp ATAN (as expected):

(%i7) :lisp (atan 0.0 -1.0)
3.141592653589793

(%i7) :lisp (atan -0.0 -1.0)
-3.141592653589793

The underlying issue in the bfloat case is that for points on the negative x-axis, atan2(y,x) has been using the value atan(y/x)-%pi when it should use atan(y/x)+%pi.

http://functions.wolfram.com/01.15.27.0006.01
http://functions.wolfram.com/01.15.27.0004.01

(N.B. Wolfram's ArcTan[x,y] has its args swapped relative to our atan2(y,x).)

Patch with new tests below. I'll commit this sometime soon if there are no objections.

diff --git a/src/float.lisp b/src/float.lisp
index b8476abfc..ec3480c9c 100644
--- a/src/float.lisp
+++ b/src/float.lisp
@@ -920,6 +920,8 @@

 ;; atan(y/x) taking into account the quadrant.  (Also equal to
 ;; arg(x+%i*y).)
+;;
+;; See also http://functions.wolfram.com/01.15.27.0006.01
 (defun fpatan2 (y x)
   (cond ((equal (car x) 0)
     ;; atan(y/0) = atan(inf), but what sign?
@@ -934,12 +936,12 @@
    ((signp g (car x))
     ;; x > 0.  atan(y/x) is the correct value.
     (fpatan (fpquotient y x)))

-   ((signp g (car y))
-    ;; x < 0, and y > 0.  We're in quadrant II, so the angle we
-    ;; want is pi+atan(y/x).
+   ((signp ge (car y))
+    ;; x < 0, and y >= 0.  We're in quadrant II or on the
+    ;; negative real axis, so the angle we want is pi+atan(y/x).
     (fpplus (fppi) (fpatan (fpquotient y  x))))
    (t
-    ;; x <= 0 and y <= 0.  We're in quadrant III, so the angle we
+    ;; x < 0 and y < 0.  We're in quadrant III, so the angle we
     ;; want is atan(y/x)-pi.
     (fpdifference (fpatan (fpquotient y x)) (fppi)))))

diff --git a/tests/rtest_atan2.mac b/tests/rtest_atan2.mac
index e7a1cb197..988c6873b 100644
--- a/tests/rtest_atan2.mac
+++ b/tests/rtest_atan2.mac
@@ -255,6 +255,30 @@ block([trigsign : true,ans],
   ans);
 0$

+atan2 (-1.0b0, -1.0b0);
+-2.356194490192345b0;
+
+atan2 (-1.0b0, 0.0b0);
+-1.570796326794897b0;
+
+atan2 (-1.0b0, 1.0b0);
+-7.853981633974484b-1;
+
+atan2 (0.0b0, 1.0b0);
+0.0b0;
+
+atan2 (1.0b0, 1.0b0);
+7.853981633974484b-1;
+
+atan2 (1.0b0, 0.0b0);
+1.570796326794897b0;
+
+atan2 (1.0b0, -1.0b0);
+2.356194490192345b0;
+
+atan2 (0.0b0, -1.0b0);
+3.141592653589793b0;
+
 (kill(values),0);
 0$ 

@@ -262,4 +286,4 @@ facts();
 []$

 (reset(trigsign),0);
-0$
\ No newline at end of file
+0$

[Raymond Toy]

Probably a left over from when atan2 for bigfloats supported signed zeroes.

But I think we should also fix atan2 for float numbers. It would surprise many people if atan2(-0.0, -1.0) returned a different value from atan2(-0b0, -1b0)

All of the sign issues could be easily fixed by just adding +0.0 to the args of atan2 (float case).

[Kris Katterjohn]

Hi Ray. Thanks.

In [#4643] there is a more general discussion about how to handle signed zeros in Maxima. atan2 is used as an example there. (And I created this current ticket based on the discussion there.)

[Raymond Toy]

Well, to be consistent with bfloats, we should get rid of signed zeroes. We can start with the elementary special functions. Any function with a branch cut on the negative real axis needs updating.

Having said that, I actually like signed zeroes. Makes it a lot easier to reason about what the value should be because nothing is ever on the branch cut.

The alternative is to make bfloats support signed zeroes. I think that means addin an extra field to the bfloat object to hold the sign. Pretty pervasive change, though.

[Stavros Macrakis]

I'm not comfortable with adding negative zero to Maxima.
For one thing, not all CLs support it.
Would we also fix limit(1/x,x,-0.0)? If we're going to do that, maybe we'd also want limit(1/x,x,zerob) to work?
Would we also need a negative 0? I realize that the integers don't have infinitesimals, but 0 in Maxima denotes both the integer and the real zero. Should -0 parse as -0.0?
I think it's a bag of worms.

[Raymond Toy]

Agreed. We just need (for now) to look at the special functions and add 0.0 to each arg. That's easy.

And even in CL, the handling of branch cuts is quite complicated. What is the value of (atanh 2.0) Is 2 converted to #c(2d0 0d0). But that introduces a forced +0d0 for the imaginary component. I think what cmucl does (or tries to) is to go from first principles and not introduce an imaginary component until it's unvaoidable. This idea is based on how Kahan generally computed z + x. The imaginary part of z is unchanged because we're only adding the real number x to the real part of z.

I think cmucl gets different values for these cases from other lisps. It's all quite messy.

[Kris Katterjohn]

Fixed with [bcfea3]. Closing this ticket.