Re: [LAD] Is -ffast-math safe for audio?

On Thu, 22 Nov 2018 23:29:11 +0100
Robin Gareus <robin@gareus.org> wrote:

>Hi Will,
>
>I just ran your code and -ffast-math does not make any difference.
>
>With or without --ffast-math I get "int: 5 rem: 0.049994"
>
>However optimizing the code with `-O2 --ffast-math` does make a
>difference because SSE is used.
>
>Do you also use -O2, or -O3 along with --fast-math?
>
>On 11/22/2018 06:30 PM, Will Godfrey wrote:
>[...]
>> My suspicion is that the difference is due to accumulated rounding
>> errors.
>>
>> Curiously without the decrements the behavior with and without
>> -ffast-math seems to be identical well into the millions.
>Yep. you do loose precision. In IEEE 32bit float, there is no 0.1.
>
>0.1 would be approximated as (13421773) * 2^(-27), and 0.05 as
>(13421773) * 2^(-28) and truncated.
>
>Note that this is an odd number. The last bit of the mantissa
>(0x004ccccd) is lost when the exponent changes.
>
>A simpler example to show this is
>
>#include <stdio.h>
>#include <math.h>
>int main (int argc, char** argv) {
> float a = 0;
> for (int i = 0; i < 100; ++i) {
> a += 0.1f;
> a -= 0.05f;
> a = fmodf (a, 1.f);
> }
> printf ("%f\n", a);
> return 0;
>}
>
>using gcc 6.3.0-18+deb9u1, x86_64, this
>prints 1.000000 (when compiled with -O0)
>and 0.000001 (when compiled with -O2 --fast-math)
>
>
>The difference here is that the former calls fmodf(), while in the
>latter, optimized, case the compiler uses cvtss2sd SSE instead. The
>internal limits of float precision differ for those cases.
>
>
>--ffast-math issues in audio-dsp are usually due to re-ordering and
>reciprocal approximations. e.g. instead of (x / 2) the compiler calls
>(0.5 * x), which can lead to different results if the inverse value does
>not a precise floating point representation. e.g. 1.0 / 0.1 != 10.0
>
>A good read on the subject is
>http://www.itu.dk/~sestoft/bachelor/IEEE754_article.pdf (Goldberg,
>David: "What Every Computer Scientist Should Know About Floating-Point
>Arithmetic").
>
>If you were change the two constants to
> a += 0.5f;
> a -= 0.25f;
>the issue vanishes.
>
>Cheers!
>robin

Thanks for going into this in such detail Robin. I never realised fp stuff
could be *quite* so, umm, approximate!

I was using O3. Following on from this I found that if I removed either O3 or
ffast-math the problem disappeared.

It's quite possible that I'm over-thinking this as there are actually only a
few iterations before the initial figures are re-calculated, but I don't like
mysteries :(

-- 
Will J Godfrey
http://www.musically.me.uk
Say you have a poem and I have a tune.
Exchange them and we can both have a poem, a tune, and a song.
_______________________________________________
Linux-audio-dev mailing list
Linux-audio-dev@lists.linuxaudio.org
https://lists.linuxaudio.org/listinfo/linux-audio-dev