22 Nov
2018
22 Nov
'18
7:30 p.m.
While testing some mixed floating point and integer calculations I found a
quite surprising difference when this compiler option was set (gcc 6.x). It was
clearly different at only 100 iterations and got dramatically worse with
larger counts.
My test routine was this:
int a = 0;
float b = 0;
float c = 0;
float inc = 0.1f;
float dec = 0.05f;
int it = 100;
for (int i = 0; i < it; ++ i)
{
a = (int)truncf(b);
c = b - floorf(b);
b += inc;
a = (int)truncf(b);
c = b - floorf(b);
b -= dec;
}
cout << "int " << a << " rem " << c <<
endl;
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.
--
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.
22 Nov
22 Nov
7:57 p.m.
Am 22.11.18 um 18:30 schrieb Will Godfrey:
While testing some mixed floating point and integer
calculations I found a
quite surprising difference when this compiler option was set (gcc 6.x). It was
clearly different at only 100 iterations and got dramatically worse with
larger counts.
My test routine was this:
int a = 0;
float b = 0;
float c = 0;
float inc = 0.1f;
float dec = 0.05f;
int it = 100;
for (int i = 0; i < it; ++ i)
{
a = (int)truncf(b);
c = b - floorf(b);
b += inc;
a = (int)truncf(b);
c = b - floorf(b);
b -= dec;
}
cout << "int " << a << " rem " << c <<
endl;
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.
In the guitarix project we've disabled -ffast-math several years ago,
when I remember right it was at gcc3, as it could lead to different
un-reproduciable calculations. Last option I've disabled on gcc8 now, is
-ffinite-math-only, this one leads to nan's and inf's in several cases,
which been as well not reproducible.
8:29 p.m.
On 2018-11-22 18:57:15 (+0100), Hermann Meyer wrote:
In the guitarix project we've disabled
-ffast-math several years ago,
when I remember right it was at gcc3, as it could lead to different
un-reproduciable calculations. Last option I've disabled on gcc8 now,
is -ffinite-math-only, this one leads to nan's and inf's in several
cases, which been as well not reproducible.
Rabbit hole stuff! SuperCollider came
to a similar conclusion:
https://github.com/supercollider/supercollider/issues/4116
--
https://sleepmap.de
9:59 p.m.
On Thu, 22 Nov 2018 19:28:58 +0100
David Runge <dave(a)sleepmap.de> wrote:
On 2018-11-22 18:57:15 (+0100), Hermann Meyer wrote:
Thanks a lot people. I kept thinking I was doing something wrong!
--
It wasn't me! (Well actually, it probably was)
... the hard part is not dodging what life throws at you,
but trying to catch the good bits.
In the guitarix project we've disabled
-ffast-math several years ago,
when I remember right it was at gcc3, as it could lead to different
un-reproduciable calculations. Last option I've disabled on gcc8 now,
is -ffinite-math-only, this one leads to nan's and inf's in several
cases, which been as well not reproducible.
Rabbit hole stuff! SuperCollider came
to a similar conclusion:
https://github.com/supercollider/supercollider/issues/4116 
10:12 p.m.
On 11/22/2018 07:28 PM, David Runge wrote:
Rabbit hole stuff! SuperCollider came to a similar
conclusion:
https://github.com/supercollider/supercollider/issues/4116
This is a different issue. In SuperCollider's case they do want NaN and
not finite-math, that is not usually the case for audio.
If you have perform a calculation that can produce NaN or infinite
values or needs to distinguish +/-0 in your DSP code, there is a much
bigger issue.
ciao,
robin
10:07 p.m.
On Thu, Nov 22, 2018 at 06:57:15PM +0100, Hermann Meyer wrote:
Am 22.11.18 um 18:30 schrieb Will Godfrey:
> While testing some mixed floating point and integer calculations I found a
> quite surprising difference when this compiler option was set (gcc 6.x). It was
> clearly different at only 100 iterations and got dramatically worse with
> larger counts.
>
> ...
In the guitarix project we've disabled
-ffast-math several years ago, when
I remember right it was at gcc3, as it could lead to different
un-reproduciable calculations. Last option I've disabled on gcc8 now, is
-ffinite-math-only, this one leads to nan's and inf's in several cases,
which been as well not reproducible.
-ffast-math should be perfectly OK for any normal audio DSP code.
The only use case where I'd consider not using it is for precision
measurement applications.
If using it leads to unexpected results that means that the algorithm
used is itself numerically unstable. A classic example of this is the
'biquad' and similar structures when used for filters that operate
mainly on low frequencies. Unfortunately, that sort of code is
everywhere, cut and pasted from one bad example to the next.
Strange that this turns up right now. Earlier today I was implementing
an IIR filter originally designed in MatLab which uses doubles for
everything. When the coefficients were truncated to 32-bit floats,
the filter became unstable. This doesn't mean that such a filter
needs double precision. It just indicates bad design.
Ciao,
--
FA
10:28 p.m.
Am 22.11.18 um 21:07 schrieb Fons Adriaensen:
On Thu, Nov 22, 2018 at 06:57:15PM +0100, Hermann
Meyer wrote:
In guitarix nearly all DSP is generated by FAUST.
However, problems with NAN's and INF's when use -ffinite-math-only
occurs only when we save preset values to file, and only sometimes then.
Without -ffinite-math-only never ever.
Am 22.11.18 um 18:30 schrieb Will Godfrey:
-ffast-math should be perfectly OK for any
normal audio DSP code.
The only use case where I'd consider not using it is for precision
measurement applications.
If using it leads to unexpected results that means that the algorithm
used is itself numerically unstable. A classic example of this is the
'biquad' and similar structures when used for filters that operate
mainly on low frequencies. Unfortunately, that sort of code is
everywhere, cut and pasted from one bad example to the next.
Strange that this turns up right now. Earlier today I was implementing
an IIR filter originally designed in MatLab which uses doubles for
everything. When the coefficients were truncated to 32-bit floats,
the filter became unstable. This doesn't mean that such a filter
needs double precision. It just indicates bad design.
Ciao, While testing some mixed floating point and
integer calculations I found a
quite surprising difference when this compiler option was set (gcc 6.x). It was
clearly different at only 100 iterations and got dramatically worse with
larger counts.
...
In the guitarix project we've disabled -ffast-math several years ago,
when
I remember right it was at gcc3, as it could lead to different
un-reproduciable calculations. Last option I've disabled on gcc8 now, is
-ffinite-math-only, this one leads to nan's and inf's in several cases,
which been as well not reproducible.
24 Nov
24 Nov
9:40 a.m.
Am 22.11.18 um 22:41 schrieb Fons Adriaensen:
On Thu, Nov 22, 2018 at 09:27:58PM +0100, Hermann
Meyer wrote:
Just as a side note, we've no problem in guitarix when we use float
instead of double. I've done so for some ARM ports.
In guitarix nearly all DSP is generated by FAUST.
That doesn't make any difference for numerical stability.
This is a property of an algorithm, not of the language.
Ciao,
23 Nov
23 Nov
12:44 a.m.
On 11/22/2018 09:27 PM, Hermann Meyer wrote:
However, problems with NAN's and INF's when use -ffinite-math-only
occurs only when we save preset values to file, and only sometimes then.
A shot in the dark..
Serializing a float in most parts of the world uses comma as decimal
separator. e.g. "0,5"
Reading this in the C-locale (dot as decimal separator) with c/c++
standard library functions returns 0. This may lead to a division by
zero later.
Cheers!
robin
4:56 a.m.
Am 22.11.18 um 23:44 schrieb Robin Gareus:
On 11/22/2018 09:27 PM, Hermann Meyer wrote:
that happen during write, not during read. During read, nan's and inf's
will be replaced by standard values when detected, but it mean's that
original settings been lost.
The values been saved as 0.5 with a dot, the comma is used as separator
between the entry's.
(and no error with entry separation ever happen)
It could happen as well to values which may ain't been used/changed
during a session.
And it happen in a non-reproducible way. Very strange, and indeed only
when -ffinite-math-only
is turned on.
I've tried for weeks to find out what really happen, but could only nail
it to -ffinite-math-only.
-ffast-math, on the other hand changes the resulting sound, it's
minimal, but, notable.
The only flag from the -ffast-math macro we still use is -fno-math-errno.
regards
hermann
However, problems with NAN's and INF's when use -ffinite-math-only
occurs only when we save preset values to file, and only sometimes then.
A shot in
the dark..
Serializing a float in most parts of the world uses comma as decimal
separator. e.g. "0,5"
Reading this in the C-locale (dot as decimal separator) with c/c++
standard library functions returns 0. This may lead to a division by
zero later.
Cheers!
robin
6:09 a.m.
Am 23.11.18 um 03:56 schrieb Hermann Meyer:
Am 22.11.18 um 23:44 schrieb Robin Gareus:
As a side note, in my side project GxPlugins.lv2 I still use -ffast-math
without problems.
On 11/22/2018 09:27 PM, Hermann Meyer wrote:
that happen during write, not during read. During read, nan's and
inf's will be replaced by standard values when detected, but it mean's
that original settings been lost.
The values been saved as 0.5 with a dot, the comma is used as
separator between the entry's.
(and no error with entry separation ever happen)
It could happen as well to values which may ain't been used/changed
during a session.
And it happen in a non-reproducible way. Very strange, and indeed only
when -ffinite-math-only
is turned on.
I've tried for weeks to find out what really happen, but could only
nail it to -ffinite-math-only.
-ffast-math, on the other hand changes the resulting sound, it's
minimal, but, notable.
The only flag from the -ffast-math macro we still use is -fno-math-errno.
regards
hermann
However, problems with NAN's and INF's when use -ffinite-math-only
occurs only when we save preset values to file, and only sometimes
then.
A shot in the dark..
Serializing a float in most parts of the world uses comma as decimal
separator. e.g. "0,5"
Reading this in the C-locale (dot as decimal separator) with c/c++
standard library functions returns 0. This may lead to a division by
zero later.
Cheers!
robin
12:29 a.m.
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
2 p.m.
On Thu, 22 Nov 2018 23:29:11 +0100
Robin Gareus <robin(a)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:
[...]
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.
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
3:09 p.m.
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Thanks for going into this in such detail Robin. I
never realised fp stuff
could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
float a = (1 + 1e20) - 1e20;
float b = 1 + 1e20; b -= 1e20;
printf ("%f %f\n", a, b);
Any guesses what this code prints without compiling it? :)
I was using O3. Following on from this I found that if
I removed either O3 or
ffast-math the problem disappeared.
OK, so SSE/AVX/compiler-intrinsics are the issue at hand.
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 :(
You could also disable this for a specific function only
__attribute__((optimize("-fno-fast-math")))
void different_mysteries_here (float val)
{
...
}
ciao,
robin
4:18 p.m.
On Fri, Nov 23, 2018 at 02:09:02PM +0100, Robin Gareus wrote:
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Surely the answer is just to use 16-bit ints for everything, then...?
--
Gordonjcp
Thanks for going into this in such detail Robin.
I never realised fp stuff
could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely.. 
7:33 p.m.
In Fri, 23 Nov 2018 14:18:01 +0000
Gordonjcp <gordonjcp(a)gjcp.net> wrote:
On Fri, Nov 23, 2018 at 02:09:02PM +0100, Robin Gareus
wrote:
Let me note - it would be at least 32-bit int, since in 32bit float
significand is 23bit.
Of course, if clipping-resistance is not important (special precaucions
like keeping signal level in range of e.g. -6db or -12db could emulate
it).
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Surely the answer is just to use 16-bit ints for everything, then...?
Thanks for going into this in such detail Robin.
I never realised
fp stuff could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
7:56 p.m.
On Fri, Nov 23, 2018 at 10:33:24PM +0500, Nikita Zlobin wrote:
In Fri, 23 Nov 2018 14:18:01 +0000
Gordonjcp <gordonjcp(a)gjcp.net> wrote:
We don't need that much precision, we only have 10-bit ears.
--
Gordonjcp
On Fri, Nov 23, 2018 at 02:09:02PM +0100, Robin
Gareus wrote:
Let me note - it would be at least 32-bit int, since in 32bit float
significand is 23bit. On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Surely the answer is just to use 16-bit ints for everything, then...?
Thanks for going into this in such detail Robin.
I never realised
fp stuff could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
8:31 p.m.
In Fri, 23 Nov 2018 17:56:38 +0000
Gordonjcp <gordonjcp(a)gjcp.net> wrote:
On Fri, Nov 23, 2018 at 10:33:24PM +0500, Nikita
Zlobin wrote:
Oh... when you wrote "for everything" - i thought, it includes
processing too, not just listening.
In Fri, 23 Nov 2018 14:18:01 +0000
Gordonjcp <gordonjcp(a)gjcp.net> wrote:
We don't need that much precision, we only have 10-bit ears.
On Fri, Nov 23, 2018 at 02:09:02PM +0100, Robin
Gareus wrote:
Let me note - it would be at least 32-bit int, since in 32bit float
significand is 23bit. On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
> Thanks for going into this in such detail Robin. I never
> realised fp stuff could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not
matter at all, or may be off completely..
Surely the answer is just to use 16-bit ints for everything,
then...? 
8:37 p.m.
On 11/23/2018 09:18 AM, Gordonjcp wrote:
On Fri, Nov 23, 2018 at 02:09:02PM +0100, Robin Gareus
wrote:
I've said before: Bankers forbid fp and use BCD math,
otherwise pennies go missing which accumulate into dollars.
Unless Will's requirement is scientific or astronomical,
yet still requires a fairly wide range, I suggest 128-bit int math.
It may not be useful with plugin or jack audio paths, which are
always fp, but for other things like time it really helps.
Until several months ago MusE used fp to represent time, such as
wall-clock time and jack time etc. It caused many timing problems.
So I ripped up all of that and used 128-bit int math instead.
You just have to carefully determine your required range and
change the units that it represents. In MusE's case, the time
value now represents 128-bit int nanoseconds instead of fp seconds.
128-bit math is not supported fully yet on all platforms
especially 32-bit platforms. So I used a third-party library
'libdivide' which handles the cross-platform stuff.
I made a convenient routine which does A * B / C in 128-bit math
since that is by far the most common general usage in our app.
Tim.
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Surely the answer is just to use 16-bit ints for everything, then...?
Thanks for going into this in such detail Robin.
I never realised fp stuff
could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
7:06 p.m.
It could be lame question, but still... is it possible, that some
implementations (compiler/hardware) will print 1? (even 64bit doesn't
hold 20 decimal digits).?
In Fri, 23 Nov 2018 14:09:02 +0100
Robin Gareus <robin(a)gareus.org> wrote:
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Thanks for going into this in such detail Robin.
I never realised
fp stuff could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
float a = (1 + 1e20) - 1e20;
float b = 1 + 1e20; b -= 1e20;
printf ("%f %f\n", a, b);
Any guesses what this code prints without compiling it? :)
I was using O3. Following on from this I found
that if I removed
either O3 or ffast-math the problem disappeared.
OK, so SSE/AVX/compiler-intrinsics are the issue at hand.
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 :(
You could also disable this for a specific function only
__attribute__((optimize("-fno-fast-math")))
void different_mysteries_here (float val)
{
...
}
ciao,
robin
_______________________________________________
Linux-audio-dev mailing list
Linux-audio-dev(a)lists.linuxaudio.org
https://lists.linuxaudio.org/listinfo/linux-audio-dev
7:37 p.m.
It's probably gonna be silly question, but after short analysis i don't
see, what could be broken in this demo snippet, when float is standard
single-precision 32bit type.
I omit first case, as optimizing compiler could just optimize it to
just =1; though it could do it in second case too... (as assumed, i did
not try to compile it yet).
But even without underhood tricks - 32bit float has 23/24bit for
significand, thus
In Fri, 23 Nov 2018 14:09:02 +0100
Robin Gareus <robin(a)gareus.org> wrote:
On 11/23/2018 01:00 PM, Will Godfrey wrote:
[...]
Thanks for going into this in such detail Robin.
I never realised
fp stuff could be *quite* so, umm, approximate!
Depending on context and the maths, the difference may not matter at
all, or may be off completely..
float a = (1 + 1e20) - 1e20;
float b = 1 + 1e20; b -= 1e20;
printf ("%f %f\n", a, b);
Any guesses what this code prints without compiling it? :)
I was using O3. Following on from this I found
that if I removed
either O3 or ffast-math the problem disappeared.
OK, so SSE/AVX/compiler-intrinsics are the issue at hand.
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 :(
You could also disable this for a specific function only
__attribute__((optimize("-fno-fast-math")))
void different_mysteries_here (float val)
{
...
}
ciao,
robin
_______________________________________________
Linux-audio-dev mailing list
Linux-audio-dev(a)lists.linuxaudio.org
https://lists.linuxaudio.org/listinfo/linux-audio-dev
24 Nov
24 Nov
12:14 p.m.
On Thu, Nov 22, 2018 at 11:29:11PM +0100, Robin Gareus wrote:
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)
Actually 0.999999940 and 0.000000596.
The 1.000000 would imply that the fmodf (a, 1.0f) would be
plain wrong, but it isn't.
The examples shown in this thread all involve converting
floats to ints. Typical application of this in audio is
to convert a float index into a wavetable to an int index
and a float < 1 interpolation argument.
The dangerous thing to do is:
// given float p
int i = (int) floorf (p);
float f = fmodf (p, 1.0f);
as you could end up with i + f != p.
The safe way is of course:
int i = (int) floorf (p);
float f = p - i;
Ciao,
--
FA
12:50 p.m.
On Sat, 24 Nov 2018 11:14:06 +0100
Fons Adriaensen <fons(a)linuxaudio.org> wrote:
<snip>
The dangerous thing to do is:
// given float p
int i = (int) floorf (p);
float f = fmodf (p, 1.0f);
as you could end up with i + f != p.
The safe way is of course:
int i = (int) floorf (p);
float f = p - i;
Ciao,
Thanks very much for that Fons :)
I'd been mulling over *exactly* that point for some time. My reasoning being
that in the latter case, if the integer was slightly wrong then using it for
the subtraction should give a remainder that was also slightly wrong, but in a
direction that tends to compensate for the error.
The other thing, is why do we need the floorf? Or in the original example
(which was taken from working code) truncf?
A straight cast to int seems to give exactly the same result, and is at least
twice as fast with -O3 and about 4 times as fast unoptimised.
--
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.
3:19 p.m.
On Sat, Nov 24, 2018 at 10:49:46AM +0000, Will Godfrey wrote:
>The safe way is of course:
>
>int i = (int) floorf (p);
>float f = p - i;
I'd been mulling over *exactly* that point for
some time. My reasoning being
that in the latter case, if the integer was slightly wrong then using it for
the subtraction should give a remainder that was also slightly wrong, but in a
direction that tends to compensate for the error.
How can an int be 'slightly wrong' ? :-)
The advantage of the 'safe' way is that you always have p == i + f.
The other thing, is why do we need the floorf? Or in
the original example
(which was taken from working code) truncf?
A straight cast to int seems to give exactly the same result, and is at least
twice as fast with -O3 and about 4 times as fast unoptimised.
We want f >= 0, so rounding must be towards -inf. Casting will truncate
(i.e. round towards zero). This gives the correct result only if p >= 0.
That may be all you need, but I wouldn't like to depend on it.
There is a way to avoid all float to int conversions, at least outside
the per-sample loops.
Suppose you have a wavetable of size L, current position is float p,
and increment is float f. To generate N samples you'd have something
like:
for (i = 0; i < N; i++)
{
k = floorf (p);
u = p - k;
// use k, u to interpolate in wavetable
p += f;
if (p >= L) p -= L;
}
To avoid floorf() inside the loop, instead of maintaining p and f
as floats, split both of them from the start into an integer and
float part:
k = floorf (p);
u = p - k;
kf = floorf (f);
uf = f - kf;
for (i = 0; i < N; i++)
{
// use k, u to interpolate in wavetable
k += kf;
u += uf;
if (u >= 1.0f)
{
k += 1;
u -= 1;
}
if (k >= L) k -= L;
// or k &= L-1 if L is a power of 2.
}
Ciao,
--
FA
25 Nov
25 Nov
12:18 p.m.
On Sat, 24 Nov 2018 14:19:01 +0100
Fons Adriaensen <fons(a)linuxaudio.org> wrote:
On Sat, Nov 24, 2018 at 10:49:46AM +0000, Will Godfrey
wrote:
Ummm, bad wording on my part. I meant out by 1. The difference would be made up
next period.
>The safe way is of course:
>
>int i = (int) floorf (p);
>float f = p - i;
I'd been mulling over *exactly* that point
for some time. My reasoning being
that in the latter case, if the integer was slightly wrong then using it for
the subtraction should give a remainder that was also slightly wrong, but in a
direction that tends to compensate for the error.
How can an int be 'slightly wrong' ? :-) The advantage of the 'safe' way is that you
always have p == i + f.
Well in this particular case, it can't be less than zero - notwithstanding time
travel :)
The other thing, is why do we need the floorf? Or
in the original example
(which was taken from working code) truncf?
A straight cast to int seems to give exactly the same result, and is at least
twice as fast with -O3 and about 4 times as fast unoptimised.
We want f >= 0, so rounding must be towards -inf. Casting will truncate
(i.e. round towards zero). This gives the correct result only if p >= 0.
That may be all you need, but I wouldn't like to depend on it. There is a way to avoid all float to int conversions,
at least outside
the per-sample loops.
Suppose you have a wavetable of size L, current position is float p,
and increment is float f. To generate N samples you'd have something
like:
for (i = 0; i < N; i++)
{
k = floorf (p);
u = p - k;
// use k, u to interpolate in wavetable
p += f;
if (p >= L) p -= L;
}
To avoid floorf() inside the loop, instead of maintaining p and f
as floats, split both of them from the start into an integer and
float part:
k = floorf (p);
u = p - k;
kf = floorf (f);
uf = f - kf;
for (i = 0; i < N; i++)
{
// use k, u to interpolate in wavetable
k += kf;
u += uf;
if (u >= 1.0f)
{
k += 1;
u -= 1;
}
if (k >= L) k -= L;
// or k &= L-1 if L is a power of 2.
}
Ciao,
Interesting. I'll look into this.
--
Will J Godfrey
http://www.musically.me.uk
Just because you know what I'm talking about, it doesn't mean *I* do.
2218
days inactive
2221
days old
linux-audio-dev@lists.linuxaudio.org
25 comments
9 participants
participants (9)
-
David Runge -
Fons Adriaensen -
Gordonjcp -
Hermann Meyer -
Nikita Zlobin -
Robin Gareus -
Tim -
Will Godfrey -
Will J Godfrey