Re: [PATCH] d3dx9: Reduce rounding error in D3DXQuaternionToAxisAngle

[Henri Verbeet]

On 9 November 2017 at 16:45, Alex Henrie <alexhenrie24(a)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)"?

[PATCH] d3dx9: Reduce rounding error in D3DXQuaternionToAxisAngle

[Alex Henrie]

2017-11-09 6:33 GMT-07:00 Henri Verbeet <hverbeet(a)gmail.com>:
> On 9 November 2017 at 16:45, Alex Henrie <alexhenrie24(a)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

[PATCH] d3dx9: Reduce rounding error in D3DXQuaternionToAxisAngle

[Chris Robinson]

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).

[PATCH] d3dx9: Reduce rounding error in D3DXQuaternionToAxisAngle

[Alex Henrie]

2017-11-09 8:05 GMT-07:00 Chris Robinson <chris.kcat(a)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