9
$\begingroup$

Bug introduced in 5.2 or earlier and persisting through 12.0.0


Why ArcTan[1, 0. I ] yields -1.5708+0. I ? The result should equal to 0.

enter image description here

Is this a bug?

$\endgroup$
8
  • 1
    $\begingroup$ it does look like a bug. For a work around use exact zero ArcTan[1, 0*I] instead of non-exact zero ArcTan[1, 0.*I] $\endgroup$
    – Nasser
    Commented Aug 24, 2019 at 5:30
  • 1
    $\begingroup$ @Nasser The original expression I calculated is like ArcTan[1, p q], where p is a MachinePrecision complex number and q is a real number. Let q be 0, it yeilds -1.5708+0. I, but it should be 0. And then I found this is because the result of ArcTan[1, 0. I ] is -1.5708+0. I during the calculation. $\endgroup$ Commented Aug 24, 2019 at 10:20
  • 1
    $\begingroup$ @xzczd I have tried to plot Re[ArcTan[1, (x + I y)]] and Im[ArcTan[1, (x + I y)]] in xy plane, and there is no branch cut near x=0 and y=0. $\endgroup$ Commented Aug 24, 2019 at 10:24
  • 2
    $\begingroup$ (1) The docs "define" ArcTan[x,y] in terms of quadrants, implying real arguments, but also give a formula valid for complex numbers, -I Log[(x + I y)/Sqrt[x^2 + y^2]]. (2) The MKL seems to define atan2 only for reals. So maybe WRI botched the implementation. The second argument 0. I seems to be a special case, since any small nonzero number yields a correct value. $\endgroup$
    – Michael E2
    Commented Aug 24, 2019 at 10:40
  • 1
    $\begingroup$ I added the bugs tag. There seems to be consensus (at least no dissent). Please report it to Wolfram Research. $\endgroup$
    – Michael E2
    Commented Aug 25, 2019 at 14:53

1 Answer 1

8
$\begingroup$

It looks like a bug.

Perhaps

ArcTan[1, p q + $MinMachineNumber * I]

will be acceptable. The little noise is unlikely to be a numerical problem. The main pitfall is when p q == -$MinMachineNumber * I, which drops you back into the buggy case.

ArcTan[1. + 0. I, $MinMachineNumber*I]
(*  0. + 0. I  *)

Update: Runtime environments for existing code

You can overload ArcTan with code like this:

Internal`InheritedBlock[{ArcTan},
 Unprotect[ArcTan];
 ArcTan[x_, y_] := func[x, y];
 Protect[ArcTan];
 ArcTan[1., 0. I]
 ]

There are several issues with doing this, and the simplest thing to do is to understand your particular use-case and choose an appropriate compromise.

For instance, here is one possibility:

ClearAll[runWithNewArcTan];
SetAttributes[runWithNewArcTan, HoldAll];
runWithNewArcTan[code_] :=
  Internal`InheritedBlock[{ArcTan},
   Unprotect[ArcTan];

   (* keeps packed arrays from being unpacked *)
   ClearAttributes[ArcTan, Listable];

   (* makes all ArcTan[x,y] results Complex;
    * vectorized formula means it still works on lists
    * however, the formula can have rounding errors when y == 0
    * of around $MachineEpsilon in magnitude *)
   ArcTan[x_, y_] /; ! FreeQ[{x, y}, _Complex] :=
     -I Log[(x + I y)/Sqrt[x^2 + y^2]];

   Protect[ArcTan];

   (* run code *)
   code
   ];

Here is another:

ClearAll[runWithNewArcTan];
SetAttributes[runWithNewArcTan, HoldAll];
runWithNewArcTan[code_] :=
  Internal`InheritedBlock[{ArcTan},
   Unprotect[ArcTan];
   (* keeping Listable attribute means packed arrays will be unpacked *)

   (* fix just the buggy values; the patterns tests are not
    * vectorized, so the Listable attribute will unpack packed arrays
    * even when neither definition below is used. *)
   ArcTan[x_, y_] /; 
     Precision[{x, y}] === MachinePrecision && 
       ! FreeQ[{x, y}, _Complex] && Positive[x] && y == 0 := 0. + 0. I;
   ArcTan[x_, y_] /; 
     Precision[{x, y}] === MachinePrecision && 
       ! FreeQ[{x, y}, _Complex] && Negative[x] && y == 0 := Pi + 0. I;
   Protect[ArcTan];

   (* run code *)
   code
   ];

The last one would work well in code that does not use packed arrays. I haven't thought of a simple way that would work exactly like ArcTan[x, y] but fix the bug.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.