Re: Apple atan2f

[Jarno Rajahalme]

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