From @tsee

This is a bug report for perl from smueller@​cpan.org,
generated with the help of perlbug 1.35 running under perl v5.8.7.

This bug is still present in 5.9.2.

[I just found out about perlbug. I hope this finally is the
preferred method of reporting bugs. This one's really easy to
fix (just apply the three line patch), but can be a real hassle
if it isn't fixed. I have tried contacting the former module
author (before Math​::Complex went into the core), but the mail
address of Mr. Lewart is no longer valid.]

The syntax of using the atan2 routine is, generally​:
atan2(y, x) = atan(y/x).

You will know that
atan(0) = 0. (Since tan(0) = 0 = sin(0) / cos(0).)

perl -MMath​::Complex -e "print atan2(0,1)"
prints '0' alright. (Which is, in fact atan(0/1).)

perl -MMath​::Complex -e "print atan2(1,0)"
prints pi/2 alright. (Which is, in fact atan(1/0)
and lim(x-->inf) atan(x) is pi/2.)

Now, in the complex plane, we get​:
perl -MMath​::Complex -e "print atan2(0,i)"
i/0​: Division by zero.
Died at c​:/perl/perl58/lib/Math/Complex.pm line 1284.

This is not correct. Obviously, 0/i is the same as 0/1
which is 0. Thus atan2(0,i) == atan2(0,1) == atan(0) == 0

At fault is the following code​:
sub atan2 {
  my ($z1, $z2, $inverted) = [at]_;
  my ($re1, $im1, $re2, $im2);
  if ($inverted) {
  ($re1, $im1) = ref $z2 ? [at]{$z2->cartesian} : ($z2, 0);
  ($re2, $im2) = [at]{$z1->cartesian};
  } else {
  ($re1, $im1) = [at]{$z1->cartesian};
  ($re2, $im2) = ref $z2 ? [at]{$z2->cartesian} : ($z2, 0);
  }
  if ($im2 == 0) {
  return CORE​::atan2($re1, $re2) if $im1 == 0;
  return ($im1<=>0) * pip2 if $re2 == 0;
  }
  my $w = atan($z1/$z2); # !!! line 1284
  my ($u, $v) = ref $w ? [at]{$w->cartesian} : ($w, 0);
  $u += pi if $re2 < 0;
  $u -= pit2 if $u > pi;
  return cplx($u, $v);
}

In the last test case above, $inverted is set. Therefore, if
$im2 = 0, the switched real and imaginary parts are used. But
if that's not the case, the code at line 1284 is executed with
the *unswitched* complex numbers $z1 and $z2.
I suppose the author omitted something like
($z1, $z2) = ($z2, $z1) if $inverted;
before line 1284. (The patch uses slightly different code.)

I attached a patch which adds another if($inverted){}else{}
block. It's a patch against the Math​::Complex in 5.8.7, but I
have verified that the version in the bleadperl is the same.

I assume you want to add a test for any bugfix, but since I
know nothing about perl's internal testing (and Complex.t in
5.9.2 is as ugly as it gets), I'll just include the simplest bit
of code that tests the patch, provided Math​::Complex was loaded​:

($eps defined in Complex.t)

my $res;
eval {$res = atan2(0,i);};
print( ($@​ or not abs($res) < $eps ) ? 'not ok' : 'ok' );

diff -u follows​:

--- Complex.pm  2005-04-01 11:43:02.000000000 +0200
+++ C:/perl/perl58/lib/Math/Complex.pm  2005-08-04 16:45:13.890625000 +0200
@@ -1281,7 +1281,13 @@
            return CORE::atan2($re1, $re2) if $im1 == 0;
            return ($im1<=>0) * pip2 if $re2 == 0;
        }
-       my $w = atan($z1/$z2);
+       my $w;
+        if ($inverted) {
+            $w = atan($z2/$z1);
+        }
+        else {
+            $w = atan($z1/$z2);
+        }
        my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
        $u += pi   if $re2 < 0;
        $u -= pit2 if $u > pi;
-----

Flags​:
  category=core
  severity=medium

Site configuration information for perl v5.8.7​:

Configured by builder at Mon Jun 6 13​:36​:05 2005.

Summary of my perl5 (revision 5 version 8 subversion 7) configuration​:
  Platform​:
  osname=MSWin32, osvers=5.0, archname=MSWin32-x86-multi-thread
  uname=''
  config_args='undef'
  hint=recommended, useposix=true, d_sigaction=undef
  usethreads=define use5005threads=undef useithreads=define
usemultiplicity=define
  useperlio=define d_sfio=undef uselargefiles=define usesocks=undef
  use64bitint=undef use64bitall=undef uselongdouble=undef
  usemymalloc=n, bincompat5005=undef
  Compiler​:
  cc='cl', ccflags ='-nologo -Gf -W3 -MD -Zi -DNDEBUG -O1 -DWIN32
-D_CONSOLE -DNO_STRICT -DHAVE_DES_FCRYPT -DBUILT_BY_ACTIVESTATE
-DNO_HASH_SEED -DUSE_SITECUSTOMIZE -DPERL_IMPLICIT_CONTEXT
-DPERL_IMPLICIT_SYS -DUSE_PERLIO -DPERL_MSVCRT_READFIX',
  optimize='-MD -Zi -DNDEBUG -O1',
  cppflags='-DWIN32'
  ccversion='12.00.8804', gccversion='', gccosandvers=''
  intsize=4, longsize=4, ptrsize=4, doublesize=8, byteorder=1234
  d_longlong=undef, longlongsize=8, d_longdbl=define, longdblsize=10
  ivtype='long', ivsize=4, nvtype='double', nvsize=8, Off_t='__int64',
lseeksize=8
  alignbytes=8, prototype=define
  Linker and Libraries​:
  ld='link', ldflags ='-nologo -nodefaultlib -debug -opt​:ref,icf
-libpath​:"c​:\perl\perl58\lib\CORE" -machine​:x86'
  libpth=\lib
  libs= oldnames.lib kernel32.lib user32.lib gdi32.lib winspool.lib
comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib
netapi32.lib uuid.lib ws2_32.lib mpr.lib winmm.lib version.lib
odbc32.lib odbccp32.lib msvcrt.lib
  perllibs= oldnames.lib kernel32.lib user32.lib gdi32.lib winspool.lib
comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib
netapi32.lib uuid.lib ws2_32.lib mpr.lib winmm.lib version.lib
odbc32.lib odbccp32.lib msvcrt.lib
  libc=msvcrt.lib, so=dll, useshrplib=yes, libperl=perl58.lib
  gnulibc_version='undef'
  Dynamic Linking​:
  dlsrc=dl_win32.xs, dlext=dll, d_dlsymun=undef, ccdlflags=' '
  cccdlflags=' ', lddlflags='-dll -nologo -nodefaultlib -debug
-opt​:ref,icf -libpath​:"c​:\perl\perl58\lib\CORE" -machine​:x86'

Locally applied patches​:
  ACTIVEPERL_LOCAL_PATCHES_ENTRY
  # if !defined(PERL_DARWIN)
  Iin_load_module moved for compatibility with build 806
  # endif
  # if defined(__hpux)
  Avoid signal flag SA_RESTART for older versions of HP-UX
  # endif
  PerlEx hacks for CGI​::Carp
  Less verbose ExtUtils​::Install and Pod​::Find
  instmodsh upgraded from ExtUtils-MakeMaker-6.25
  24699 ICMP_UNREACHABLE handling in Net​::Ping
  21540 Fix backward-compatibility issues in if.pm

@​INC for perl v5.8.7​:
  c​:/perl/perl58/lib
  c​:/perl/perl58/site/lib
  .

Environment for perl v5.8.7​:
  HOME (unset)
  LANG (unset)
  LANGUAGE (unset)
  LD_LIBRARY_PATH (unset)
  LOGDIR (unset)
  PATH=D​:\texmf\miktex\bin;c​:\perl\perl58\bin;C​:\cygwin\bin;C​:\Tcl\bin;p​:\util;p​:\util\gnuplot\bin;p​:\util\perltidy\bin;C​:\WINNT\system32;C​:\WINNT;C​:\WINNT\system32\nls;C​:\WINNT\system32\nls\ENGLISH;C​:\Programme\Network
Associates\PGPNT;C​:\Programme\Novell\ZENworks\;C​:\Programme\Gemeinsame
Dateien\GTK\2.0\bin;p​:\sdbdata;p​:\sdbdata\sdbdata_win;c​:\cygwin\usr\X11R6\bin
  PERL_BADLANG (unset)
  SHELL (unset)

From @jhi

While I certainly agree that a runtime divide-by-zero is not the right
thing to do, I cannot quite bring myself to agree on the proposed fix.

In my books atan2() is not well-defined for complex arguments at all,
and even if in some special cases it might have a reasonable definition,
I still don't think it makes sense.

See e.g. http​://mathworld.wolfram.com/InverseTangent.html

The RT System itself - Status changed from 'new' to 'open'

From @jhi

Jarkko Hietaniemi wrote​:

While I certainly agree that a runtime divide-by-zero is not the right
thing to do, I cannot quite bring myself to agree on the proposed fix.

In my books atan2() is not well-defined for complex arguments at all,
and even if in some special cases it might have a reasonable definition,
I still don't think it makes sense.

This bit​:

: This is not correct. Obviously, 0/i is the same as 0/1
: which is 0. Thus atan2(0,i) == atan2(0,1) == atan(0) == 0

seems to be incorrect... if one believes the Mathematica definition
atan2(0, i) is pi/2, not zero.

Another good reference on these matters, MATLAB ignores the complex
parts of atan2() arguments.

See e.g. http​://mathworld.wolfram.com/InverseTangent.html

From @jhi

This bit​:

: This is not correct. Obviously, 0/i is the same as 0/1
: which is 0. Thus atan2(0,i) == atan2(0,1) == atan(0) == 0

seems to be incorrect... if one believes the Mathematica definition
atan2(0, i) is pi/2, not zero.

Ahem, sorry... it seems that Mathematica's ArcTan is ArcTan[x, y] while
atan2 is atan2(y, x)... so the zero is right, after all.

Another good reference on these matters, MATLAB ignores the complex
parts of atan2() arguments.

This still is true.

From @jhi

Hence, I'd suggest patching the behaviour one way or another. There's
three things I'd consider sensible to do​:

- atan2($z1, $z2) = atan($z1/$z2) (same order as with the normal atan2)
- remove the atan2($z1, $z2) alltogether and just carp() instead.
- atan2($z1, $z2) = atan(Re($z1)/Re($z2)) (Same as MATLAB)

Either way, the behaviour should be documented. And IMHO it shouldn't be
what it is now. Maybe there are other sensible things atan2 could do. Tell
me which you consider the best and I'll write the patch.

I'll work out a patch... I have other outstanding things to do for
Math​::Complex and Math​::Trig.

Best regards,
Steffen

From st.mueller@smt.zeiss.com

Hi Jarkko,

You wrote​:

Another good reference on these matters, MATLAB ignores the complex
parts of atan2() arguments.

Well, that's fine. (Though not how I'd like to use the function.) Right
now, however, the behaviour is clearly wrong.

We have​:
  atan2($real1, $real2) == atan($real1/$real2)
and​:
  atan2($cplx1, $cplx2) == atan($cplx2/$cplx1)

(Due to overloading, it doesn't matter if any one of the $cplxX is
actually real, but you know that, of course.)

I understand your argument that atan2 of complex arguments isn't well
defined. But the current behaviour isn't even treating the parameter order
as with real numbers. That's just wrong. (Which is, in fact, well defined,
but still bad. ;-)

Hence, I'd suggest patching the behaviour one way or another. There's
three things I'd consider sensible to do​:

- atan2($z1, $z2) = atan($z1/$z2) (same order as with the normal atan2)
- remove the atan2($z1, $z2) alltogether and just carp() instead.
- atan2($z1, $z2) = atan(Re($z1)/Re($z2)) (Same as MATLAB)

Either way, the behaviour should be documented. And IMHO it shouldn't be
what it is now. Maybe there are other sensible things atan2 could do. Tell
me which you consider the best and I'll write the patch.

Best regards,
Steffen

From @jhi

The attached patch (#1) brings several pending updates to Math​::Complex
and Math​::Trig.

  - Complex​: fix for the [perl #31117]​: atan2(0, i) now works,
  as do all the (computable) complex argument cases (I adopted
  the Mathematica definition)
  - Complex​: fixes for certain bugs in the make() and emake()
  - Complex​: support returning the kth root directly
  - Complex​: support [2,-3pi/8] in emake()
  - Complex​: support 'inf' for make()/emake() (the 9**9**9 trick)
  - Complex​: document make()/emake() more visibly
  - Trig​: add more great circle routines​: great_circle_waypoint()
  and great_circle_destination()

I also made the modules scratch-proof, err, 5.00504 and 5.6.2-resistant.

Since the Complex+Trig have had some fixes/updates over the years, maybe
also them should be spun off to CPAN so that users of older Perls can
use the latest versions?

The tiny patch (#2) documents the attested fact that atan2(0, 0) is not
well-defined (i.e., or ie., or ie, the result depends on the underlying
math libraries).

Both patches should be maint-worthy.

From @jhi

cxtr.pat.bz2

From @jhi

atan2.pat.bz2

@rgs - Status changed from 'open' to 'resolved'

From @rgs

Jarkko Hietaniemi wrote​:

The attached patch (#1) brings several pending updates to Math​::Complex
and Math​::Trig.

\- Complex​: fix for the \[perl \#31117\]​: atan2\(0\, i\) now works\,
  as do all the \(computable\) complex argument cases \(I adopted
  the Mathematica definition\)
\- Complex​: fixes for certain bugs in the make\(\) and emake\(\)
\- Complex​: support returning the kth root directly
\- Complex​: support \[2\,\-3pi/8\] in emake\(\)
\- Complex​: support 'inf' for make\(\)/emake\(\) \(the 9\*\*9\*\*9 trick\)
\- Complex​: document make\(\)/emake\(\) more visibly
\- Trig​: add more great circle routines​: great\_circle\_waypoint\(\)
  and great\_circle\_destination\(\)

I also made the modules scratch-proof, err, 5.00504 and 5.6.2-resistant.

Since the Complex+Trig have had some fixes/updates over the years, maybe
also them should be spun off to CPAN so that users of older Perls can
use the latest versions?

The tiny patch (#2) documents the attested fact that atan2(0, 0) is not
well-defined (i.e., or ie., or ie, the result depends on the underlying
math libraries).

Both patches should be maint-worthy.

Both applied to blead as #25414, thanks. But may I suggest that next time,
you run tests before sending them ? Just look again at your next_test()
function for lib/Math/Complex.t (I didn't apply this part) and you'll
see that you made a trivial mistake... :)

--
Faithful to this feeble magic, he would invent, so that they might not happen,
the most atrocious particulars.
  -- Borges

From Nadeem.Douba@ec.gc.ca

There is a small technical bug in "sub great_circle_destination {...}".
An additional (optional) argument, RHO, must be passed to the
subroutine. Otherwise, the calculation will always be based on the
assumption that we are using a unit sphere. RHO will represent the
radius of the sphere. I have previously written a subroutine also named
great_circle_destination which may be of use to you :) The subroutine
assumes that the user is inputting coordinates in radians using the
following standards​:

North Western Hemisphere = (
  -pi <= Longitude <= 0,
  0 <= Latitude <= pi/2
  )

South Western Hemisphere = (
  -pi <= Longitude <= 0,
  -pi/2 <= Latitude <= 0
  )

North Eastern Hemisphere = (
  0 <= Longitude <= pi,
  0 <= Latitude <= pi/2
  )

South Eastern Hemisphere = (
  0 <= Longitude <= pi,
-pi/2 <= Latitude <= 0
  )

There is no need to subtract by pi/2 (pip2) in this case. The
destination is also required to be input in radians where​:

North is 0;
South is -pi or pi;
East is pi/2;
West is -pi/2;

Here it is (tried and tested on GIS applications)​:

=cut
Example​:

my ($dlat, $dlon) = great_circle_distance($lat, $lon, $distance,
$direction, 6378.1);

# 6378.1 is the radius of the earth in km. You could alternatively use
any
# other value for rho, for example, 3443.89849 which is the radius of
the
# earth in nautical miles;)
=cut

sub great_circle_destination {

  my (
  $lat0,
  $lon0,
  $distance,
  $direction,
  $rho
  ) = @​_;

  my (
  $lat1,
  $lon1
  );

# lat2 = asin(sin(lat1)*cos(d) + cos(lat1)*sin(d)*cos(hdng))
# lon2 = lon1 + atan2(sin(hdng)*sin(d)*cos(lat1),
cos(d)-sin(lat1)*sin(lat2))
# Reference​: http​://www.movable-type.co.uk/scripts/LatLong.html

  $rho = 1 unless defined($rho);
  $distance = $distance/$rho;
  $lat1 = asin(sin($lat0)*cos($distance) +
cos($lat0)*sin($distance)*cos($direction));
  $lon1 = $lon0 + atan2(sin($direction)*sin($distance)*cos($lat0),
cos($distance)-sin($lat0)*sin($lat1));

  return ($lat1, $lon1);

}

-----Original Message-----
From​: Jarkko Hietaniemi [mailto​:jhietaniemi@​gmail.com]
Sent​: Wednesday, September 14, 2005 2​:26 AM
To​: Mueller, Steffen
Cc​: perl5-porters@​perl.org; Douba,Nadeem [CIS]
Subject​: [PATCH] Math​::Complex and Math​::Trig updates (Re​: [perl #37117]
Math​::Complex atan2 bug)

The attached patch (#1) brings several pending updates to Math​::Complex
and Math​::Trig.

  - Complex​: fix for the [perl #31117]​: atan2(0, i) now works,
  as do all the (computable) complex argument cases (I adopted
  the Mathematica definition)
  - Complex​: fixes for certain bugs in the make() and emake()
  - Complex​: support returning the kth root directly
  - Complex​: support [2,-3pi/8] in emake()
  - Complex​: support 'inf' for make()/emake() (the 9**9**9 trick)
  - Complex​: document make()/emake() more visibly
  - Trig​: add more great circle routines​: great_circle_waypoint()
  and great_circle_destination()

I also made the modules scratch-proof, err, 5.00504 and 5.6.2-resistant.

Since the Complex+Trig have had some fixes/updates over the years, maybe
also them should be spun off to CPAN so that users of older Perls can
use the latest versions?

The tiny patch (#2) documents the attested fact that atan2(0, 0) is not
well-defined (i.e., or ie., or ie, the result depends on the underlying
math libraries).

Both patches should be maint-worthy.

From @jhi

Rafael Garcia-Suarez wrote​:

Jarkko Hietaniemi wrote​:

Both applied to blead as #25414, thanks. But may I suggest that next time,
you run tests before sending them ? Just look again at your next_test()

Ooops. In my defense I can say I did

perl -wIlib lib/Math/Complex.t|grep not

(and came up empty) but yeah, didn't run full test or use the harness...

function for lib/Math/Complex.t (I didn't apply this part) and you'll
see that you made a trivial mistake... :)

From @jhi

There is a small technical bug in "sub great_circle_destination {...}".
An additional (optional) argument, RHO, must be passed to the
subroutine. Otherwise, the calculation will always be based on the

I am confused. On a sphere the destination point d (radians)
one ends up after going from point a (radians) in direction b (radians)
for c (radians) doesn't depend on the rho (radius) ... that's the beauty
of radians, right? Or do you mean you want to specify the distance
to travel in some units other than radians?

assumption that we are using a unit sphere. RHO will represent the
radius of the sphere. I have previously written a subroutine also named
great_circle_destination which may be of use to you :) The subroutine
assumes that the user is inputting coordinates in radians using the
following standards​:

We cannot change the default coordinate system of the great_circle_*
routines. For better or worse, the pi/2 or 90 deg needed twist has been
documented to work that way for quite some time now, we cannot change
this without breaking people's code. See also the discussion in the
documentation about the theta and phi​: these is no one single "correct"
spherical coordinate system.

One possible way to incorporate your semantics would be to add
an export tag, e.g. "​:great_circle_latlon" that would affect the
great circle formulas.