Re: Apple atan2f

[Rik]

11/16/10

I'm late in getting back to this thread, but I've been short on time to
code some necessary tests.  My executive summary, for those short on time,
is that this is indeed an Apple library issue and that Octave should not
change its code.

First, the C language specification says that 1) the return value for atan2
must be in the range [-pi, pi] and 2) that atan2 (0, -x) = pi.  The atan2
function is multi-valued with a period of 2*pi.  atan2(1,1) = {...,
-7*pi/4, pi/4, 9*pi/4, ...}.  Requirement 1 removes the ambiguity and
multiple values by specifying that only the principal value is to be
returned.  There is, however, one leftover bit of ambiguity.  The negative
x-axis could be described equally by -pi or +pi and both of those values
are within the range of requirement 1.  Therefore, requirement 2 becomes
necessary and picks the positive value of pi for the negative x-axis.

The C language specification, wisely, says nothing about the actual machine
representation of pi.  Given that C is compiled on 8-bit microcontrollers
and vector supercomputers, on little-endian machines and on big-endian
machines, it really would be difficult to set a standard that worked on all
such diverse systems.

However, CPUs don't operate on symbolic constants.  Every number needs a
representation and this is where the IEEE-754 specification takes over from
the C language specification.  I will use Frep in the following discussion
to be the representation function which translates a symbolic number to a
machine representation.  The Frep function is many-to-one and
deterministic.  It must necessarily be many-to-one because there are an
infinite number of real numbers and a finite number of representations.
Ergo, many real numbers will be mapped to the same representation.  Frep is
also deterministic and will yield the same representation for the same
input every time.  As a thought experiment, imagine that this were not true
and that Frep(1) yielded either 1 or 0.  The mathematical expression "1 -
1" would be translated to "Frep(1) - Frep(1)" === "{1,0} - {1,0}" = "{-1,
0, 1}".  Having 3 possible programming outcomes, where the symbolic math
calculation yields just one, is a disaster.  This is, in effect, what Apple
has done by breaking the IEEE specification.  Frep(pi) has two values
depending on whether it is calculated through a constructor "float (M_PI)"
or through a trig function "atan2 (0, -1)".

Continuing with this line of thought, I translate the C requirements as:
1) output of atan2 in the range Frep( [-pi,pi] ) === [Frep(-pi), Frep(pi)].
2) atan2(0, -1) = Frep (pi).
There is no requirement for the mixed conditional Frep (pi) <
symbolic_value (pi).  This is what the the Apple developer mistakenly
implemented when he forced a return value that was not the representation
of pi used elsewhere in the math libraries.

Jarno brought up an issue about the quadrant differing between
"tan(double(pi))" and "tan(double(single(pi)))".  I agree that they differ,
but I don't think this is an issue.  The functions double and single are
many-to-one and therefore do not have proper inverses.  Thus,
"Finverse(F(x)) = x" is not true as it would be for a proper function.  As
a concrete example, should "tan(double(floor(pi)))" be the same as
"tan(double(pi))"?  After all, floor is just a different representation
function.

The second reason I don't believe this to be an issue is that I can find no
concensus that the axes themselves belong to any particular quadrant.  See
this reference on Wolfram Mathworld
(http://mathworld.wolfram.com/Quadrant.html) which specifically says the
axes don't belong to any quadrant.  It is the programmer's choice,
therefore, as to whether axes get mapped into 0, 1, or 2 quadrants.  For
example, when producing a map of the globe, cities directly on the equator
should appear in both a map of the southern hemisphere and the northern
hemisphere.  This application requires axes mapped to two quadrants.  Even
if you took the position that the positive x-axis belonged to quadrant 1,
the positive y-axis to quadrant 2, etc., where would you put the point
(0,0) and why?

Cheers,
Rik

Companion Program and Testing:
C99 allows the rounding rules for the IEEE-754 representation function to
be queried and set by the programmer.  I used this to examine the behavior
of glibc under different rounding rules.  Since Apple's atan2f function
uses truncation rounding, it should be possible to set all of the math
library functions to use truncation rounding in which case there will be
only one representation of pi and their bug will disappear.  It is
interesting to note that glibc does not get atan2f correct either.  They
always use round-to-nearest which works only when the default rounding
rules are in effect.

Optimizing compilers will interfere with this test since they can detect
constant expressions and substitute return values at compile time.  I work
around this by having the user enter the number "0.0" and then do all of my
calculations as "number + user_input".  This forces the evaluation to be
done at run-time by the math libraries.  For good measure I also compiled
without optimizations and forced g++ not to use substitutions.  My
compilation line was:

g++ -O0 -frounding-math -std=c++0x    rndtst.cpp   -o rndtst

I first check that the default rounding rules are in fact round-to-nearest.
 Next, I look at the representations of pi and Apple's almostpi under
different rounding rules.  Finally, I examine 3 different methods of
obtaining the value pi and their representations (which should all be
equal).  The three methods are: 1) direct constructor, 2) atan2(0, -1), 3)
2*asin(1).

Results for glibc 2.10.1

--------------------
Default round mode is 0
TOWARDZERO round mode is 3072
TONEAREST round mode is 0
--------------------
TONEAREST : (float) almostpi = 3.1415925 (0x40490fda)
TONEAREST : (float) M_PI = 3.1415927 (0x40490fdb)
TOWARDZERO: (float) M_PI = 3.1415925 (0x40490fda)
TOWARDZERO: (float) almostpi = 3.1415925 (0x40490fda)
--------------------
Computed with TONEAREST
atan2f (0.0f, -1.0f) = 3.14159274 (0x40490fdb)
atan2f (0.0f, -1.0f) - (float) M_PI = 0.00000000 (0x00000000)
2*asinf (1.0f) = 3.14159274 (0x40490fdb)
2*asinf (0.0f, -1.0f) - (float) M_PI = 0.00000000 (0x00000000)
--------------------
Computed with TOWARDZERO
atan2f (0.0f, -1.0f) = 3.14159274 (0x40490fdb)
atan2f (0.0f, -1.0f) - (float) M_PI = 0.00000024 (0x34800000)
2*asinf (1.0f) = 3.14159250 (0x40490fda)
2*asinf (0.0f, -1.0f) - (float) M_PI = 0.00000000 (0x00000000)
--------------------

rndtst.cpp
Description: Text Data