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)

--------------------