Re: Apple atan2f

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Rik-5
> Looking at the code itself shows that the intentional rounding towards zero is done ONLY when returning one for the prescribed pi/-pi values. It appears that when in the Polynomial, it will automatically stay within the prescribed range. Also,  there are no multiplications or divisions (as documented) when returning pi or -pi. And even if there were, calculations are performed in double precision, so any roundoffs are insignificant when casting to single.
>
> The code that returns atan2f(0.0f, -1.0f) is:
>
> NearNegativeXAxis:
> // Return -pi or +pi, matching the sign of y, rounded toward zero.
>
> movsd Address(AlmostpPi), p0
> xorpd Base, p0
> jmp ReturnDouble
>
> AlmostpPi was defined earlier as:
>
> /* Define values near +/- pi that yield +/- pi rounded toward zero when
> converted to single precision.  This allows us to generate inexact and
> return the desired values for atan2f on and near the negative side of the x
> axis.
> */
> AlmostpPi: .double +3.1415924
>
> Also, AlmostpPi is never used in the Polynomial itself, but used exclusively and always for returning either pi or -pi (Base contains the negative sign above when needed).
>
> So, obviously the implementer of the code read the requirement for the range [-pi,pi] to mean the whole function. And the point is that to be in that range, the range is actually (-PI,PI) (PI here representing the actual (symbolic) value of pi), as no numeric precision is enough to contain the exact value of pi. Maybe it happens that the double precision pi rounds naturally towards zero?
Good digging into the code.  I just read the comments and didn't bother to
check the assembly language.  You're right that the implementer obviously
intended this, although I still think they are wrong.

The C language specification says that the code must return pi.  However,
it doesn't specify the actual representation.  According to my quick web
search
(http://developer.apple.com/library/mac/documentation/Darwin/Reference/ManPages/man3/math.3.html)
the libm functions on Mac OS X use IEEE-754 representations.  So how should
the irrational number pi be represented according to IEEE-754?  The default
rounding rule for the binary representation is to round to the nearest
floating point value (See, for example,
http://en.wikipedia.org/wiki/IEEE_floating_point).  The atan2f author
violates the IEEE-754 specification when he uses a rounding function which
rounds towards zero as is noted in the code comments.

Imagine representing pi in base10 with only four decimal digits.
pi = 3.14159265...
pi(nearest value) = 3.1416
pi(round->0) = 3.1415

Interestingly, the 754 specification does allow 4 other types of rounding
rules, including rounding towards zero.  But, if Apple wants to use
rounding towards zero it needs to do that uniformly across all their math
functions.

Here are results from Octave on a Linux machine:
octave:1> almostpi = +3.1415924;
octave:2> format hex
octave:3> single(almostpi)
ans = 400921fb40000000
octave:4> single(pi)
ans = 400921fb60000000
octave:5> pi
ans = 400921fb54442d18

The value of pi from double precision shows that it needs a hex '5' in the
last position.  Unfortunately, the single precision does not have the ones
bit so the final hex digit must be either 4 or 6.  Linux seems to follow
the IEEE specification and rounds to the nearest value which is 6.

--Rik
Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Michael Godfrey
On 11/10/2010 08:27 PM, Rik wrote:
Imagine representing pi in base10 with only four decimal digits.
pi = 3.14159265...
pi(nearest value) = 3.1416
pi(round->0) = 3.1415

Interestingly, the 754 specification does allow 4 other types of rounding
rules, including rounding towards zero.  But, if Apple wants to use
rounding towards zero it needs to do that uniformly across all their math
functions.

Here are results from Octave on a Linux machine:
octave:1> almostpi = +3.1415924;
octave:2> format hex
octave:3> single(almostpi)
ans = 400921fb40000000
octave:4> single(pi)
ans = 400921fb60000000
octave:5> pi
ans = 400921fb54442d18

The value of pi from double precision shows that it needs a hex '5' in the
last position.  Unfortunately, the single precision does not have the ones
bit so the final hex digit must be either 4 or 6.  Linux seems to follow
the IEEE specification and rounds to the nearest value which is 6.

--Rik
It is worth asking: for what value of n does pi - n*eps become different from pi
using round to zero against round to nearest, or other choices.  Obviously in
this case n for round to zero is a lot bigger.  Iterative algorithms which depend
on differences need to consider when they are within a loss of significance
due to rounding region.  I am really surprised that Intel has made this choice.

Michael
 

Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

bpabbott
Administrator
On Nov 11, 2010, at 1:45 PM, Michael D Godfrey wrote:

> On 11/10/2010 08:27 PM, Rik wrote:
>> Imagine representing pi in base10 with only four decimal digits.
>> pi = 3.14159265...
>> pi(nearest value) = 3.1416
>> pi(round->0) = 3.1415
>>
>> Interestingly, the 754 specification does allow 4 other types of rounding
>> rules, including rounding towards zero.  But, if Apple wants to use
>> rounding towards zero it needs to do that uniformly across all their math
>> functions.
>>
>> Here are results from Octave on a Linux machine:
>> octave:1> almostpi = +3.1415924;
>> octave:2> format hex
>> octave:3> single(almostpi)
>> ans = 400921fb40000000
>> octave:4> single(pi)
>> ans = 400921fb60000000
>> octave:5> pi
>> ans = 400921fb54442d18
>>
>> The value of pi from double precision shows that it needs a hex '5' in the
>> last position.  Unfortunately, the single precision does not have the ones
>> bit so the final hex digit must be either 4 or 6.  Linux seems to follow
>> the IEEE specification and rounds to the nearest value which is 6.
>>
>> --Rik
>>
> It is worth asking: for what value of n does pi - n*eps become different from pi
> using round to zero against round to nearest, or other choices.  Obviously in
> this case n for round to zero is a lot bigger.  Iterative algorithms which depend
> on differences need to consider when they are within a loss of significance
> due to rounding region.  I am really surprised that Intel has made this choice.
>
> Michael

Rather than "Intel", I assume you mean "Apple".

I don't think this was intentional, but if it was, I'd qualify it as a bug since it is a violation of the IEEE standard.

Ben
Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Michael Godfrey
On 11/10/2010 11:17 PM, Ben Abbott wrote:
Rather than "Intel", I assume you mean "Apple".

Right. I meant Apple/Intel.
I don't think this was intentional, but if it was, I'd qualify it as a bug since it is a violation of the IEEE standard.

Well, not quite a violation of IEEE since there are rounding choices, but almost surely not
a good choice.  I think there is a reference implementation of the elementary functions,
but right now I cannot remember where it is.  It may include specification about rounding
as well as the algorithms.

Michael


Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Jarno Rajahalme

On Nov 11, 2010, at 9:25 , ext Michael D Godfrey wrote:

On 11/10/2010 11:17 PM, Ben Abbott wrote:
Rather than "Intel", I assume you mean "Apple".

Right. I meant Apple/Intel.
I don't think this was intentional, but if it was, I'd qualify it as a bug since it is a violation of the IEEE standard.

Well, not quite a violation of IEEE since there are rounding choices, but almost surely not
a good choice.  I think there is a reference implementation of the elementary functions,
but right now I cannot remember where it is.  It may include specification about rounding
as well as the algorithms.

Michael



See below the related C standard text (as found from http://std.dkuug.dk/JTC1/SC22/WG14/www/docs/n843.htm). Note the paragraph #3, which specifies the range of returned values. Returning single(pi) using the normal rounding rules is simply wrong, as the result is outside of the prescibed range, since single(pi) > pi. Note that it happens that the double representation of PI happens to round towards zero, so double(pi) is in the right range.

Returning a value on the wrong side of PI puts the returned angle to the wrong quadrant of the circle, flipping the sign of the tan function, for example:

single(pi) > PI:

octave:76> tan(double(single(pi)))
ans =  8.74227800037251e-08

rounding towards zero:

octave:77> tan(double(single(pi-1e-7)))
ans = -1.50995799097839e-07
octave:78> 


  242          Committee Draft  --  August 3, 1998   WG14/N843


       7.12.4.4  The atan2 functions

       Synopsis

       [#1]

               #include <math.h>
               double atan2(double y, double x);
               float atan2f(float y, float x);
               long double atan2l(long double y, long double x);

       Description

       [#2] The atan2 functions compute the principal value of  the
       arc  tangent  of  y/x,  using the signs of both arguments to
       determine the quadrant of the return value.  A domain  error
       may occur if both arguments are zero.

       Returns

       [#3]  The  atan2 functions return the arc tangent of y/x, in
       the range [-pi, +pi] radians.

Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Michael Godfrey
On 11/11/10 12:19 AM, Jarno Rajahalme wrote:
See below the related C standard text (as found from http://std.dkuug.dk/JTC1/SC22/WG14/www/docs/n843.htm). Note the paragraph #3, which specifies the range of returned values. Returning single(pi) using the normal rounding rules is simply wrong, as the result is outside of the prescibed range, since single(pi) > pi. Note that it happens that the double representation of PI happens to round towards zero, so double(pi) is in the right range.

Returning a value on the wrong side of PI puts the returned angle to the wrong quadrant of the circle, flipping the sign of the tan function, for example:

single(pi) > PI:

octave:76> tan(double(single(pi)))
ans =  8.74227800037251e-08

rounding towards zero:

octave:77> tan(double(single(pi-1e-7)))
ans = -1.50995799097839e-07
octave:78> 


  242          Committee Draft  --  August 3, 1998   WG14/N843


       7.12.4.4  The atan2 functions

       Synopsis

       [#1]

               #include <math.h>
               double atan2(double y, double x);
               float atan2f(float y, float x);
               long double atan2l(long double y, long double x);

       Description

       [#2] The atan2 functions compute the principal value of  the
       arc  tangent  of  y/x,  using the signs of both arguments to
       determine the quadrant of the return value.  A domain  error
       may occur if both arguments are zero.

       Returns

       [#3]  The  atan2 functions return the arc tangent of y/x, in
       the range [-pi, +pi] radians.

OK.  I think this is correct, but it is time to do some more checking
if I can find the time.  Anyhow, very good that you looked this up.
Ending up in a wrong quadrant is also bad!

Michael

Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Rik-5
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 (2K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

bpabbott
Administrator
On Nov 16, 2010, at 1:29 PM, Rik wrote:

>
> 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
>

Using Apple's compiler (gcc 4.2.1 with their patches) I get the error below.

$ g++ -O0 -frounding-math -std=c++0x    rndtst.cpp   -o rndtst
cc1plus: error: unrecognized command line option "-std=c++0x"

Ben


Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Jordi Gutiérrez Hermoso
On 16 November 2010 17:29, Ben Abbott <[hidden email]> wrote:

> Using Apple's compiler (gcc 4.2.1 with their patches) I get the error below.
>
> $ g++ -O0 -frounding-math -std=c++0x    rndtst.cpp   -o rndtst
> cc1plus: error: unrecognized command line option "-std=c++0x"

gcc 4.3 was where C++0x features (updates to the C++ standard) were
beginning to be implemented, but gcc 4.3 is also where the move to
GPLv3 was done. Apple seems intent on shunning GPLv3 as much as
possible. This is largely whence Apple's interest in clang: it's not
copylefted. They pay some lip service to technical merit, but the fact
is that Apple is trying to stifle copylefted software. Pity, because
clang doesn't implement much or any of C++0x, and I expect it won't
for some time.

Sorry for the lecture, but I actually find this worrying. It might
create a schism between developers in GNU-based systems and Mac OS X
users, where we'll now have to worry about supporting two major
compilers instead of focussing mostly on gcc. True that Octave tries
to be very compatible across compilers, but clang would be yet another
compiler to support, and a major one at that. I expect that packaing
higher versions of gcc for macports or similar will be a chore few
will attend to.

- Jordi G. H>

Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

bpabbott
Administrator
On Nov 16, 2010, at 6:42 PM, Jordi Gutiérrez Hermoso wrote:

> On 16 November 2010 17:29, Ben Abbott <[hidden email]> wrote:
>
>> Using Apple's compiler (gcc 4.2.1 with their patches) I get the error below.
>>
>> $ g++ -O0 -frounding-math -std=c++0x    rndtst.cpp   -o rndtst
>> cc1plus: error: unrecognized command line option "-std=c++0x"
>
> gcc 4.3 was where C++0x features (updates to the C++ standard) were
> beginning to be implemented, but gcc 4.3 is also where the move to
> GPLv3 was done. Apple seems intent on shunning GPLv3 as much as
> possible. This is largely whence Apple's interest in clang: it's not
> copylefted. They pay some lip service to technical merit, but the fact
> is that Apple is trying to stifle copylefted software. Pity, because
> clang doesn't implement much or any of C++0x, and I expect it won't
> for some time.
>
> Sorry for the lecture, but I actually find this worrying. It might
> create a schism between developers in GNU-based systems and Mac OS X
> users, where we'll now have to worry about supporting two major
> compilers instead of focussing mostly on gcc. True that Octave tries
> to be very compatible across compilers, but clang would be yet another
> compiler to support, and a major one at that. I expect that packaing
> higher versions of gcc for macports or similar will be a chore few
> will attend to.
>
> - Jordi G. H>

Running g++ 4.4, I get

$ ./rndtst
Please enter '0.0' and hit <Return>
0.0
--------------------
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.14159250 (0x40490fda)
atan2f (0.0f, -1.0f) - (float) M_PI = -0.00000024 (0xb4800000)
2*asinf (1.0f) = 3.14159250 (0x40490fda)
2*asinf (0.0f, -1.0f) - (float) M_PI = -0.00000024 (0xb4800000)
--------------------
Computed with TOWARDZERO
atan2f (0.0f, -1.0f) = 3.14159226 (0x40490fd9)
atan2f (0.0f, -1.0f) - (float) M_PI = -0.00000023 (0xb4800000)
2*asinf (1.0f) = 3.14159226 (0x40490fd9)
2*asinf (0.0f, -1.0f) - (float) M_PI = -0.00000023 (0xb4800000)
--------------------

Regarding Apple, it is possible to use non-Apple compilers. It is just a hassle to manage all the dependencies. Often needing two copies. On for Apple and another for gnu.

Ben





Reply | Threaded
Open this post in threaded view
|

Re: Apple atan2f

Jarno Rajahalme
In reply to this post by Rik-5

On Nov 16, 2010, at 21:29 , ext Rik wrote:

> 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.
>

This may be nitpicking, but functions by definition cannot be multi-valued. Also, the C language specification says that atan2(-0, -x) = -pi, so there are two different results for the negative x-axis, which may be necessary because the value of pi cannot be exactly represented. This also means that there is no atan2 () return value describing the negative x-axis, exactly.

> 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.

The standard does, however, specify that:

"Returns
The atan2 functions return arctan y/x in the interval [−π , +π ] radians."

This is the same regardless of the precision (double, float, long double). Since the representation is not specified, there is no reason to doubt that the range is what it needs to be mathematically, i.e. none of the returned values should be < −π or > +π.

>
> 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)".
>

But that is only natural, given the definition of the functions. Assuming the default rounding rules, float(M_PI) rounds to the nearest float representation, and atan2(0, -1) to the nearest double representation within the allowed range for atan2 function. You do not expect these representation to be the same, right?

Also, as a side note, I am not aware of a IEEE specification for pi, atan2, or IEEE interpretation of the mathematical range  [−π , +π ]. Do you have a reference for this?

> 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)].

Doing this, you are expanding the range whenever Frep() results in a bigger absolute value than it's argument's absolute value. Given your interpretation atan2f () return values would not be confined to the principal values, as Frep(-pi) and Frep(pi) are crossing to the wrong side of the negative x-axis.

> 2) atan2(0, -1) = Frep (pi).

For atan2f () this gives a return value on the negative side of the y-axis, effectively -pi, when considering the principal range of the atan2 function.

> 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.
>

If the intent of the range specification for the atan2 function is to retain the result in the principal range, then you are mistaken above. Also, there is no "the representation of pi", there are only approximations of pi, each with differing accuracy, given the number of available bits, base, and rounding rules. There is, however, the C language specification requirement that the atan2 function return a value that stays within the principal range, and within the right quadrant (-pi or pi, for example). Given your interpretation the requirement for returning pi, for example, makes no sense as it's rounded-to-nearest float representation is actually "-pi degrees" (since it is on the wrong side of the y-axis).

> 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 point is not that they differ (which is only natural, given the different representations), but that the returned values, given their representation, stay within the range specified for the function in question. In your example, since floor () is rounding towards zero, the range requirement for the atan2 () is not violated, assuming someone would be interested in integer only tan2 (tan2i ?).

Continuing on your example, and *imagining, for the purposes of this example*, that pi would naturally round to 4, when rounded to integer using Frep(), which of the following would be correct set of return values for atan2i ():
[-4 -3 -2 -1 0 1 2 3 4], or
[-3 -2 -1 0 1 2 3],
given the range requirement of [-pi, -pi] (in a standard that does not define a representation for pi), and the fact that -0 has no integer representation?

How is the above example different from the case of atan2f (), other than the different accuracy of numeric representation between integers and floats?

> 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?
>

This point is moot, as the negative x-axis has no exact radian representation at all. Any representation of the negative x-axis in radians is necessarily in either of the respective quadrants. That is why the difference between pi and -pi is significant.

Finally, it seems that the typical mathematical range for the atan2 function is actually (-pi,pi], excluding the symbolic -pi, since it is equivalent to radian (symbolic) pi. From this viewpoint it can be expected that in the C-language specification -pi and pi map to different, and obviously non-overlapping radian values.

  Jarno