On 9 November 2017 at 16:45, Alex Henrie alexhenrie24@gmail.com wrote:
*pangle = 2.0f * acosf(pq->w);
*pangle = 2.0f * acos(pq->w); /* acosf has too much rounding error */
Does that rounding error by any chance get better if you replace "acosf(x)" with "atan2f(sqrtf((1.0f - x) * (1.0f + x)), x)"?
2017-11-09 6:33 GMT-07:00 Henri Verbeet hverbeet@gmail.com:
On 9 November 2017 at 16:45, Alex Henrie alexhenrie24@gmail.com wrote:
*pangle = 2.0f * acosf(pq->w);
*pangle = 2.0f * acos(pq->w); /* acosf has too much rounding error */Does that rounding error by any chance get better if you replace "acosf(x)" with "atan2f(sqrtf((1.0f - x) * (1.0f + x)), x)"?
That actually makes the rounding error worse on my machine:
printf("acos: %.9f\n", 2*acos(10.0f/22.0f)); printf("acosf: %.9f\n", 2*acosf(10.0f/22.0f)); printf("atan: %.9f\n", 2*atan2f(sqrtf((1.0f - 10.0f/22.0f) * (1.0f + 10.0f/22.0f)), 10.0f/22.0f));
acos: 2.197868979 acosf: 2.197869062 atan: 2.197868824
-Alex
On 11/09/2017 05:41 AM, Alex Henrie wrote:
That actually makes the rounding error worse on my machine:
printf("acos: %.9f\n", 2*acos(10.0f/22.0f)); printf("acosf: %.9f\n", 2*acosf(10.0f/22.0f)); printf("atan: %.9f\n", 2*atan2f(sqrtf((1.0f - 10.0f/22.0f) * (1.0f + 10.0f/22.0f)), 10.0f/22.0f));
acos: 2.197868979 acosf: 2.197869062 atan: 2.197868824
Should be noted that the difference between 2.197869062 and 2.197868824 is a single bit on floats (0x1.1953c6p+1 and 0x1.1953c4p+1 respectively, in hex float format). 2.197868979 is a double-precision result (0x1.1953c54cea1bp+1), which will get truncated to one of the other two values when assigned to a float depending on how it rounds.
Interestingly, when I do:
volatile float f = 22.0f; printf("%.9f", 2*acosf(10.0f/f));
I get 2.197868824, the same as your atan2f result. I only get 2.197869062 when I do:
printf("%.9f", (float)(2*acos(10.0/f)));
Note that 10 is a double. If I give a single-precision input:
volatile float f = 22.0f; printf("%.9f", 2*acos(10.0f/f));
I again get 2.197868824 (0x1.1953c4p+1).
2017-11-09 8:05 GMT-07:00 Chris Robinson chris.kcat@gmail.com:
On 11/09/2017 05:41 AM, Alex Henrie wrote:
That actually makes the rounding error worse on my machine:
printf("acos: %.9f\n", 2*acos(10.0f/22.0f)); printf("acosf: %.9f\n", 2*acosf(10.0f/22.0f)); printf("atan: %.9f\n", 2*atan2f(sqrtf((1.0f - 10.0f/22.0f) * (1.0f + 10.0f/22.0f)), 10.0f/22.0f));
acos: 2.197868979 acosf: 2.197869062 atan: 2.197868824
Should be noted that the difference between 2.197869062 and 2.197868824 is a single bit on floats (0x1.1953c6p+1 and 0x1.1953c4p+1 respectively, in hex float format). 2.197868979 is a double-precision result (0x1.1953c54cea1bp+1), which will get truncated to one of the other two values when assigned to a float depending on how it rounds.
Interestingly, when I do:
volatile float f = 22.0f; printf("%.9f", 2*acosf(10.0f/f));
I get 2.197868824, the same as your atan2f result. I only get 2.197869062 when I do:
printf("%.9f", (float)(2*acos(10.0/f)));
Note that 10 is a double. If I give a single-precision input:
volatile float f = 22.0f; printf("%.9f", 2*acos(10.0f/f));
I again get 2.197868824 (0x1.1953c4p+1).
Interesting, it seems that GCC precomputes the value when it's not stored in a variable. A better test program would be:
float x = 10.0f/22.0f; printf("perfect: %.9f\n", 2.197868979f); printf("acos: %.9f\n", (float)(2*acos(x))); printf("acosf: %.9f\n", 2*acosf(x)); printf("atan2f: %.9f\n", 2*atan2f(sqrtf((1.0f - x) * (1.0f + x)), x));
Which gives:
perfect: 2.197869062 acos: 2.197869062 acosf: 2.197868824 atan2f: 2.197869062
So Henri was right, the atan2f solution does work.
-Alex
-
Alex Henrie -
Chris Robinson -
Henri Verbeet