Re: test failure for mappers.cc

[Jarno Rajahalme]

On Nov 10, 2010, at 20:20 , ext Rik wrote:

>> 
>>> I get the following
>>> 
>>> M_PI = 3.141592741012573
>>> (float) M_PI = 3.1415927 (0x40490fdb)
>>> atan2f (0.0f, -1.0f) = 3.14159250 (0x40490fda)
>>> atan2f (0.0f, -1.0f) - (float) M_PI = -0.00000024 (0xb4800000)
>>> 
>>> Looking through Apple's sources, I don't see a atan2f.c. But I do see 
>>> atan2f.s ...
>>> 
>>>     
>>> http://www.opensource.apple.com/source/Libm/Libm-315/Source/Intel/atan2f.s
>> 
>> Looking at that code, it seems it is working as specified. 
> I disagree.  The notes in the code begin with "Return value for atan2f(y,
> x) (C F.9.1 12 and C F.9.1.4):" and then outline a series of prescribed
> return values.  For the case in question they list the return value as
> 4*pi/4 (incidentally, C F.9.1.4 lists the return value as just pi.  They
> may be incorrectly implementing multiplications/divisions such that 4*pi/4
> =! pi when roundoff errors are considered).  The documentation continues
> with "Otherwise:", i.e., in all other cases it will return a value in the
> range [-pi,pi].  The note that pi will be rounded to within the interval
> [-pi,pi] appears to apply only when the algorithm is used, and not when
> returning one of the mandated values.
> 

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?

  Jarno

> --Rik
> Note:
>> 
>> "Return a value in [-pi, +pi] (C 7.12.4.4 3).  Note that this
>> prohibits returning correctly rounded values for -pi and +pi, since
>> pi rounded to a float lies outside that interval."
>> 
>> On octave, I get:
>> 
>> octave:5> pi - double(single(pi))
>> ans = -8.7423e-08
>> 
>> which shows that pi rounded to a single actually is larger than pi as a 
>> double, and therefore outside of the range [-pi, pi].
>>