Re: Apple atan2f

[Rik]

> 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