28 Jan
2011
28 Jan
'11
1:35 a.m.
Hi All,
I have been working on the Neil tracker program recently, and hit a
weird bug that seems to affect only my computer! I get a segfault when
trying to use fft.h (which is used by the oomek 303 emulator, which has
incidentally now been released as open source - a really good sound
IMHO, but I digress). The following test code results in a segfault on
my computer - an AMD XP2500+ (but no-one elses as far as I can tell).
bt:
Program received signal SIGSEGV, Segmentation fault.
0x080486d7 in IFFT (fftBuffer=0x8049b00, fftFrameSize=2048, sign=1) at
fft.h:51
51 tr = *p2r * ur - *p2i * ui;
Any thoughts as to what might be happening?
James
test.cpp:
#include <cstdio>
#include <cmath>
#include "fft.h"
float a[4096];
int main ()
{
for(int i=0; i<4096; i++)
{
a[i]=0;
}
IFFT(a,2048,1);
printf("%f",a[1]);
}
fft.h:
#define M_PI 3.14159265358979323846
void IFFT(float *fftBuffer, long fftFrameSize, long sign)
/*
FFT routine, (C)1996 S.M.Sprenger. Sign = -1 is FFT, 1 is iFFT (inverse)
Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the
time domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes
and returns the cosine and sine parts in an interleaved manner, ie.
fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize
must be a power of 2. It expects a complex input signal (see
footnote 2),
ie. when working with 'common' audio signals our input signal has to be
passed as {in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the
transform
of the frequencies of interest is in fftBuffer[0...fftFrameSize].
*/
{
float wr, wi, arg, *p1, *p2, temp;
float tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k;
for (i = 2; i < 2 * fftFrameSize - 2; i += 2) {
for (bitm = 2, j = 0; bitm < 2 * fftFrameSize; bitm <<= 1) {
if (i & bitm) j++;
j <<= 1;
}
if (i < j)
{
p1 = fftBuffer + i;
p2 = fftBuffer + j;
temp = *p1;
*(p1++) = *p2;
*(p2++) = temp;
temp = *p1;
*p1 = *p2;
*p2 = temp;
}
}
for (k = 0, le = 2; k < log(fftFrameSize)/log(2.); k++) {
le <<= 1;
le2 = le>>1;
ur = 1.0;
ui = 0.0;
arg = (float)M_PI / (le2>>1);
wr = (float)cos(arg);
wi = sign*(float)sin(arg);
for (j = 0; j < le2; j += 2) {
p1r = fftBuffer+j;
p1i = p1r+1;
p2r = p1r+le2;
p2i = p2r+1;
for (i = j; i < 2 * fftFrameSize; i += le) {
tr = *p2r * ur - *p2i * ui;
ti = *p2r * ui + *p2i * ur;
*p2r = *p1r - tr;
*p2i = *p1i - ti;
*p1r += tr;
*p1i += ti;
p1r += le;
p1i += le;
p2r += le;
p2i += le;
}
tr = ur*wr - ui*wi;
ui = ur*wi + ui*wr;
ur = tr;
}
}
}
28 Jan
28 Jan
1:40 a.m.
On Thu, Jan 27, 2011 at 7:35 PM, James Stone <jamesmstone(a)gmail.com> wrote:
Hi All,
I have been working on the Neil tracker program recently, and hit a weird
bug that seems to affect only my computer! I get a segfault when trying to
use fft.h (which is used by the oomek 303 emulator, which has incidentally
now been released as open source - a really good sound IMHO, but I digress).
The following test code results in a segfault on my computer - an AMD
XP2500+ (but no-one elses as far as I can tell).
use valgrind. its the right tool for the job.
4:01 a.m.
Hi James,
You always know how to nerd snipe me!
On Thursday, January 27, 2011 06:35:08 pm James Stone wrote:
I have been working on the Neil tracker program
recently,
and hit a weird bug that seems to affect only my
computer! I get a segfault when trying to use fft.h
Not just you. It segfaults for me, too. Ubuntu 10.04 with Core Duo
processor.
With the patch below (to fft.h) I detect that p2r and p2i overrun the
buffer when k==11, j==0, i==0. However, k was supposed to STOP when
it reached 11 = log(2048)/log(2).
My guess is that is that the integer k is promoted to a float and the
comparison is performed. When I replace:
for (k = 0, le = 2; k < log(fftFrameSize)/log(2.); k++) {
with:
long ITERS = log(fftFrameSize)/log(2.0) + 0.5;
for (k = 0, le = 2; k < ITERS; k++) {
...the code doesn't crash on me. Here's a sample program that
illustrates what's happening.
/* BEGIN */
#include <cmath>
#include <iostream>
using namespace std;
#define ITERS (log(2048)/log(2.0))
int main(void)
{
long k;
for( k=0 ; k < ITERS ; ++k ) {
cout << "k = " << k
<< " ITERS = " << ITERS
<< " diff = " << (ITERS - k)
<< endl;
}
return 0;
}
/* END */
Output on my machine:
k = 0 ITERS = 11 diff = 11
k = 1 ITERS = 11 diff = 10
k = 2 ITERS = 11 diff = 9
k = 3 ITERS = 11 diff = 8
k = 4 ITERS = 11 diff = 7
k = 5 ITERS = 11 diff = 6
k = 6 ITERS = 11 diff = 5
k = 7 ITERS = 11 diff = 4
k = 8 ITERS = 11 diff = 3
k = 9 ITERS = 11 diff = 2
k = 10 ITERS = 11 diff = 1
k = 11 ITERS = 11 diff = 4.80518e-16
Thanks,
Gabriel
Here's the testing patch....
--- a/fft.h
+++ b/fft.h
@@ -1,5 +1,7 @@
#define M_PI 3.14159265358979323846
+#include <assert.h>
+
void IFFT(float *fftBuffer, long fftFrameSize, long sign)
/*
FFT routine, (C)1996 S.M.Sprenger. Sign = -1 is FFT, 1 is iFFT
(inverse)
@@ -17,6 +19,7 @@ void IFFT(float *fftBuffer, long fftFrameSize, long
sign)
float tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k;
+
for (i = 2; i < 2 * fftFrameSize - 2; i += 2) {
for (bitm = 2, j = 0; bitm < 2 * fftFrameSize; bitm <<= 1) {
if (i & bitm) j++;
@@ -48,6 +51,10 @@ void IFFT(float *fftBuffer, long fftFrameSize, long
sign)
p2r = p1r+le2;
p2i = p2r+1;
for (i = j; i < 2 * fftFrameSize; i += le) {
+ assert(p1r < (fftBuffer + 2 * fftFrameSize));
+ assert(p1i < (fftBuffer + 2 * fftFrameSize));
+ assert(p2r < (fftBuffer + 2 * fftFrameSize));
+ assert(p2i < (fftBuffer + 2 * fftFrameSize));
tr = *p2r * ur - *p2i * ui;
ti = *p2r * ui + *p2i * ur;
*p2r = *p1r - tr;
12:02 p.m.
On Fri, Jan 28, 2011 at 3:00 AM, Gabriel M. Beddingfield
<gabrbedd(a)gmail.com> wrote:
Hi James,
You always know how to nerd snipe me!
On Thursday, January 27, 2011 06:35:08 pm James Stone wrote:
Wow! Gabriel, thanks for this! :) Excellent and very speedy work!!
James
I have been working on the Neil tracker program
recently,
and hit a weird bug that seems to affect only my
computer! I get a segfault when trying to use fft.h
Not just you. It segfaults for me, too. Ubuntu 10.04 with Core Duo
processor.
With the patch below (to fft.h) I detect that p2r and p2i overrun the
buffer when k==11, j==0, i==0. However, k was supposed to STOP when
it reached 11 = log(2048)/log(2).
My guess is that is that the integer k is promoted to a float and the
comparison is performed. When I replace:
for (k = 0, le = 2; k < log(fftFrameSize)/log(2.); k++) {
with:
long ITERS = log(fftFrameSize)/log(2.0) + 0.5;
for (k = 0, le = 2; k < ITERS; k++) {
...the code doesn't crash on me. Here's a sample program that
illustrates what's happening.
12:28 p.m.
Excerpts from Gabriel M. Beddingfield's message of 2011-01-28 04:00:51 +0100:
Hi James,
You always know how to nerd snipe me!
On Thursday, January 27, 2011 06:35:08 pm James Stone wrote:
Interesting.. would you mind explaining how this can be?
How can 11-11 yield 4.80518e-16?
Do I understand it correctly: log(2048)/log(2.0) yields a float and for comparison k is
implicitly cast to float as well. Why does the output of your
test-Program still show "k = 11"?
And why does it only happen on some machines? I'd play with it and try
to figure out what's going on and how this kind of problem can be
avoided, but the error doesn't appear here.
Here's the output on my machine:
k = 0 ITERS = 11 diff = 11
k = 1 ITERS = 11 diff = 10
k = 2 ITERS = 11 diff = 9
k = 3 ITERS = 11 diff = 8
k = 4 ITERS = 11 diff = 7
k = 5 ITERS = 11 diff = 6
k = 6 ITERS = 11 diff = 5
k = 7 ITERS = 11 diff = 4
k = 8 ITERS = 11 diff = 3
k = 9 ITERS = 11 diff = 2
k = 10 ITERS = 11 diff = 1
I have been working on the Neil tracker program
recently,
and hit a weird bug that seems to affect only my
computer! I get a segfault when trying to use fft.h
Not just you. It segfaults for me, too. Ubuntu 10.04 with Core Duo
processor.
With the patch below (to fft.h) I detect that p2r and p2i overrun the
buffer when k==11, j==0, i==0. However, k was supposed to STOP when
it reached 11 = log(2048)/log(2).
My guess is that is that the integer k is promoted to a float and the
comparison is performed. When I replace:
for (k = 0, le = 2; k < log(fftFrameSize)/log(2.); k++) {
with:
long ITERS = log(fftFrameSize)/log(2.0) + 0.5;
for (k = 0, le = 2; k < ITERS; k++) {
...the code doesn't crash on me. Here's a sample program that
illustrates what's happening.
/* BEGIN */
#include <cmath>
#include <iostream>
using namespace std;
#define ITERS (log(2048)/log(2.0))
int main(void)
{
long k;
for( k=0 ; k < ITERS ; ++k ) {
cout << "k = " << k
<< " ITERS = " << ITERS
<< " diff = " << (ITERS - k)
<< endl;
}
return 0;
}
/* END */
Output on my machine:
k = 0 ITERS = 11 diff = 11
k = 1 ITERS = 11 diff = 10
k = 2 ITERS = 11 diff = 9
k = 3 ITERS = 11 diff = 8
k = 4 ITERS = 11 diff = 7
k = 5 ITERS = 11 diff = 6
k = 6 ITERS = 11 diff = 5
k = 7 ITERS = 11 diff = 4
k = 8 ITERS = 11 diff = 3
k = 9 ITERS = 11 diff = 2
k = 10 ITERS = 11 diff = 1
k = 11 ITERS = 11 diff = 4.80518e-16
Thanks,
Gabriel 
1:04 p.m.
On Fri, 2011-01-28 at 12:28 +0100, Philipp Überbacher wrote:
Interesting.. would you mind explaining how this can
be?
How can 11-11 yield 4.80518e-16?
Because, rather than 11, log(2048)/log(2.0) is actually (to 40 d.p.)
7.6246189861593984035895533360399422488305 /
0.6931471805599453094172321214581765680755
But doubles don't contain 40 d.p.; it is approximately 16. Thus some
rounding errors may occur...
Peter.
1:09 p.m.
On Fri, Jan 28, 2011 at 12:03 PM, Peter Nelson <peter(a)fuzzle.org> wrote:
On Fri, 2011-01-28 at 12:28 +0100, Philipp Überbacher
wrote:
.. I wonder if this bugfix should be added to some other projects? I
think the same fft.h is used in one of the functions in audacity?
Perhaps someone who is on LAD could cross-post this?
James
Interesting.. would you mind explaining how this
can be?
How can 11-11 yield 4.80518e-16?
Because, rather than 11, log(2048)/log(2.0) is actually (to 40 d.p.)
7.6246189861593984035895533360399422488305 /
0.6931471805599453094172321214581765680755
But doubles don't contain 40 d.p.; it is approximately 16. Thus some
rounding errors may occur...
1:28 p.m.
On Fri, Jan 28, 2011 at 6:09 AM, James Stone <jamesmstone(a)gmail.com> wrote:
.. I wonder if this bugfix should be added to some
other projects? I
think the same fft.h is used in one of the functions in audacity?
Perhaps someone who is on LAD could cross-post this?
First, should probably check if the value should be floor()'d or
round()'ed. The "fix" I have above is rounded... but I haven't dug in
to the algo see if it should actually be floored().
-gabriel
1:41 p.m.
On Fri, 2011-01-28 at 06:28 -0600, Gabriel Beddingfield wrote:
On Fri, Jan 28, 2011 at 6:09 AM, James Stone
<jamesmstone(a)gmail.com> wrote:
I wouldn't use the log operation at all:
for (k = 0, le = 2; (1 << k) < fftFrameSize; k++)
Peter.
.. I wonder if this bugfix should be added to
some other projects? I
think the same fft.h is used in one of the functions in audacity?
Perhaps someone who is on LAD could cross-post this?
First, should probably check if the value should be floor()'d or
round()'ed. The "fix" I have above is rounded... but I haven't dug in
to the algo see if it should actually be floored().
3:31 p.m.
On Fri, 28 Jan 2011, Peter Nelson wrote:
FTW!
Not only is this a cool trick, but I think it's the Right
Answer.
Peace,
Gabriel
First, should
probably check if the value should be floor()'d or
round()'ed. The "fix" I have above is rounded... but I haven't dug in
to the algo see if it should actually be floored().
I wouldn't use the log operation at all:
for (k = 0, le = 2; (1 << k) < fftFrameSize; k++)
2:02 p.m.
Excerpts from Peter Nelson's message of 2011-01-28 13:03:56 +0100:
On Fri, 2011-01-28 at 12:28 +0100, Philipp Überbacher
wrote:
Thanks Peter, this section was supposed to refer to the output of the
test program. cout printed 11 for both values but 4.80518e-16 for one
value minus the other, so this at least a very nasty inconsistency.
I still don't understand why this results in 11 on some machines and
something else on others. 32/64 bit? Compiler differences? Luck?
rant_begin
Why can't log mean the same thing everywhere? Why does it need to be
base e here and base 10 there? Why is there no consistency?
And why is there no proper logarithmus dualis function? Because you
can simply do log(n)/log(2)? We've just seen how well this works.
How about:
log() - base 10
ln() - base e - logarithmus naturalis
ld() - base 2 - logarithmus dualis
rant_end
The next obvious question is: Does the inaccuracy reliably result in
values bigger than 11?
Interesting.. would you mind explaining how this
can be?
How can 11-11 yield 4.80518e-16?
Because, rather than 11, log(2048)/log(2.0) is actually (to 40 d.p.)
7.6246189861593984035895533360399422488305 /
0.6931471805599453094172321214581765680755
But doubles don't contain 40 d.p.; it is approximately 16. Thus some
rounding errors may occur...
Peter.
3:39 p.m.
On Fri, 28 Jan 2011, Philipp Überbacher wrote:
Thanks Peter, this section was supposed to refer to
the output of the
test program. cout printed 11 for both values but 4.80518e-16 for one
value minus the other, so this at least a very nasty inconsistency.
I still don't understand why this results in 11 on some machines and
something else on others. 32/64 bit? Compiler differences? Luck?
Truncation error related to (a) how the compiler ordered
the operations and (b) the log() implementation. I also
noticed if you assign `double iters = log(N)/log(2.0)` that
`k < iters` works as expected on my machine.
But the thing is, he violated one of the priciples taught in
Programming 101: floating point comparisons are dangerous!
Why can't log mean the same thing everywhere?
Why does it need to be
AFAIK, nearly every programming language is consistent that
log() is the natural log and log10() is the common log.
This only conflicts with normal mathematical notation.
However, in the case of log(N)/log(2.0)... it doesn't matter
if it's log-e, log-10, log-2, log-666.... the answer is the
same. This is a well-known formula to calculate the number
of bits required to represent N.
The next obvious question is: Does the inaccuracy
reliably result in
values bigger than 11?
I like Peter Nelson's solution to do (1<<k) and skip log()
alltogether.
-gabriel
4:11 p.m.
On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp Überbacher wrote:
rant_begin
Why can't log mean the same thing everywhere? Why does it need to be
base e here and base 10 there? Why is there no consistency?
And why is there no proper logarithmus dualis function? Because you
can simply do log(n)/log(2)? We've just seen how well this works.
How about:
log() - base 10
ln() - base e - logarithmus naturalis
ld() - base 2 - logarithmus dualis
rant_end
Libm has log(), log10, and log2().
The next obvious question is: Does the inaccuracy
reliably result in
values bigger than 11?
No.
If the input is a power of two, and you expect an integer as
a result, just do
k = (int)(log2(x) + 1e-6)
or
k = (int)(log(x)/log(2) + 1e-6)
or
int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Ciao,
--
FA
There are three of them, and Alleline.
29 Jan
29 Jan
12:07 p.m.
Excerpts from fons's message of 2011-01-28 16:11:52 +0100:
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant.
On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp
Überbacher wrote:
Took me a while to figure out that libm is part of glibc :)
Good to know that those functions are available on probably pretty much
all linux systems.
rant_begin
Why can't log mean the same thing everywhere? Why does it need to be
base e here and base 10 there? Why is there no consistency?
And why is there no proper logarithmus dualis function? Because you
can simply do log(n)/log(2)? We've just seen how well this works.
How about:
log() - base 10
ln() - base e - logarithmus naturalis
ld() - base 2 - logarithmus dualis
rant_end
Libm has log(), log10, and log2(). The next
obvious question is: Does the inaccuracy reliably result in
values bigger than 11?
No.
If the input is a power of two, and you expect an integer as
a result, just do
k = (int)(log2(x) + 1e-6) or
k = (int)(log(x)/log(2) + 1e-6)
or
int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
Anyway, it's a neat way to avoid the problem and the rounding properties
of mult/div in case of not power of 2 could be useful as well.
Regards,
Philipp
12:46 p.m.
On Sat, Jan 29, 2011 at 12:07:12PM +0100, Philipp Überbacher wrote:
Everything floating point 'suffers from the same problem'.
Even addition. (1/3.0 + 1/3.0 + 1/3.0 == 1.0) doesn't have
to result in true.
It's a single CPU cycle on most processors, and on some (e.g. ARM)
it can even be combined with other instructions.
Ciao,
--
FA
There are three of them, and Alleline.
Took me a while to figure out that libm is part of
glibc :)
It isn't. But it lives in /lib, not /usr/lib
k =
(int)(log2(x) + 1e-6)
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant. int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
1:22 p.m.
Excerpts from fons's message of 2011-01-29 12:46:49 +0100:
Everything floating point 'suffers from the same problem'.
Even addition. (1/3.0 + 1/3.0 + 1/3.0 == 1.0) doesn't have
to result in true.
I thought that log2() might use a different routine if the argument is a
power of 2. I guess that even if it does it's nothing to rely on because
another libm implementation might not do the same thing.
Why did you choose 1e-6 specifically?
It's a single CPU cycle on most processors, and on some (e.g. ARM)
it can even be combined with other instructions.
Really nice, but I wonder whether it really matters at the end of the
day, simply because of optimisation. I tried to compare your variant
with the one from Peter Nelson using a nice little line I got from
Torben Hohn at some point. Without optimisation the output is quite
different, but beginning with -O2 the machine code appears all the same.
# simple variant using *= 2
echo "void a() { int m, k, j = 2048; for (k = 0, m = 1; m < j; k++, m *=
2); }" | gcc -O2 -xc -S - -o -
# Fons's variant
echo "void a() { int m, k, j = 2048; for (k = 0, m = 1; m < j; k++, m
<<= 1); }" | gcc -O2 -xc -S - -o -
# Peter Nelsons variant
echo "void a() { int m, k, j = 2048; for (k = 0; (1 << k) < j; k++);
}"
| gcc -O2 -xc -S - -o -
To be honest I don't know what the line does exactly and hence whether
this is a valid comparison. Also, my knowledge of assembly is so minimal
that I have no idea where the computation takes place at all in the
optimised versions.
If the optimised code really is the same is it worth it to use funky bit
shifting operations? Especially Peter Nelsons version is non-obvious to
me (although it's probably the most efficient out of the three if no
optimisation is used).
On Sat, Jan 29, 2011 at 12:07:12PM +0100, Philipp
Überbacher wrote:
The glibc package seems to own it on my system. Anyway, I was searching
the web for libm, then looking for a header and so on, but it's pretty
clear now, it's installed by default and the header for C is math.h.
Took me a while to figure out that libm is part
of glibc :)
It isn't. But it lives in /lib, not /usr/lib k = (int)(log2(x) + 1e-6)
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant. int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
1:46 p.m.
On Sat, Jan 29, 2011 at 01:22:10PM +0100, Philipp Überbacher wrote:
I thought that log2() might use a different routine if
the argument is a
power of 2. I guess that even if it does it's nothing to rely on because
another libm implementation might not do the same thing.
There are very good reasons why math routines don't do such things.
Why did you choose 1e-6 specifically?
It does the job. 1e-6f will also work for single precision.
... Without optimisation the output is quite
different, but beginning with -O2 the machine code appears all the same.
That could well be completely different on an another CPU and compiler,
and even more so if tested ***in context***. Such tests are really
pointless, in particular if you consider the actual use of this code.
It isn't executed 100 times for each audio sample.
If the optimised code really is the same is it worth
it to use funky bit
shifting operations?
There's nothing funky about them, they are part of C and C++.
Ciao,
--
FA
There are three of them, and Alleline.
30 Jan
30 Jan
2:46 p.m.
Excerpts from fons's message of 2011-01-29 13:46:19 +0100:
Thanks, I see this now after Peter did another test with quite different
results. I looked at it out of context because I simply don't understand
the context. I only know what a ft does but not how it works, nor have I
implemented a fft, so the code with its two-letter variables is
meaningless to me.
They were funky to me mainly because I had to look them up but also
because operating on the bit representation requires a different mode of
thinking. I bet it also has its own set of pitfalls.
It's quite interesting that it still seems to make a difference
performance-wise. Thanks for your insights.
On Sat, Jan 29, 2011 at 01:22:10PM +0100, Philipp
Überbacher wrote:
What are those?
I thought that log2() might use a different
routine if the argument is a
power of 2. I guess that even if it does it's nothing to rely on because
another libm implementation might not do the same thing.
There are very good reasons why math routines don't do such things. Why did you
choose 1e-6 specifically?
It does the job. 1e-6f will also work for single precision.
... Without optimisation the output is quite
different, but beginning with -O2 the machine code appears all the same.
That could well be completely different on an another CPU and compiler,
and even more so if tested ***in context***. Such tests are really
pointless, in particular if you consider the actual use of this code.
It isn't executed 100 times for each audio sample. If the
optimised code really is the same is it worth it to use funky bit
shifting operations?
There's nothing funky about them, they are part of C and C++.
3:29 p.m.
On Sun, Jan 30, 2011 at 02:46:21PM +0100, Philipp Überbacher wrote:
If your machine uses BCD integers you could have some surprises...
If you don't like m <<= 1, just write m *= 2. The compiler will
very probably use a bit shift to do this even at low optimisation
levels.
Ciao,
--
FA
There are three of them, and Alleline.
Thanks, I see this now after Peter did another test
with quite different
results. I looked at it out of context because I simply don't understand
the context. I only know what a ft does but not how it works, nor have I
implemented a fft, so the code with its two-letter variables is
meaningless to me.
It's not clear if Peter's test included the '+ 1e-6'. I'd expect the
result to be correct then even with -O3. If not, increasing the added
constant should fix it.
There's
nothing funky about them, they are part of C and C++.
They were funky to me mainly because I had to look them up but also
because operating on the bit representation requires a different mode of
thinking. I bet it also has its own set of pitfalls.
29 Jan
29 Jan
1:54 p.m.
On Sat, 2011-01-29 at 13:22 +0100, Philipp Überbacher wrote:
Without optimisation the output is quite
different, but beginning with -O2 the machine code appears all the same.
That's because your 2048 is static, and nothing is happening in the
loop, therefore gcc can optimise way more stuff away than in the real
code, because it "knows" that buffer size is always 2048 and that k is
not used inside the loop.
Peter.
1:47 p.m.
On Sat, 2011-01-29 at 12:07 +0100, Philipp Überbacher wrote:
Excerpts from fons's message of 2011-01-28
16:11:52 +0100:
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant.
Well, now I'm just being pedantic :-), but as a quick test using rdtsc
(i.e. profiling to be taken with a grain of salt):
1: log(x)/log(2)
2: (1 << k) < x
3: m < x & m <<= 1
This is 1000 iterations; cycling through x of 512, 1024, 2048, 4096 (to
prevent the compiler optimizing log(x)/log(2) to a single call
throughout the whole test). The fourth line is the sum of the iterator
in each loop.
./a.out (unoptimized)
1 3132000 cycles
2 261468 cycles
3 285273 cycles
1 50500, 2 50500, 3 50500
./a.out (-O3)
1 2598345 cycles
2 170055 cycles
3 173673 cycles
1 50500, 2 50500, 3 50500
I must admit, I'm surprised that fons way shows slightly slower in my
test!
As an aside, here's the result with (-O3 -ffast-math)
./a.out
1 150894 cycles,
2 169983 cycles,
3 158850 cycles,
1 47750, 2 50500, 3 50500
Yeah, faster, but wrong :-)
Peter.
On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp
Überbacher wrote:
Took me a while to figure out that libm is part of glibc :)
Good to know that those functions are available on probably pretty much
all linux systems.
rant_begin
Why can't log mean the same thing everywhere? Why does it need to be
base e here and base 10 there? Why is there no consistency?
And why is there no proper logarithmus dualis function? Because you
can simply do log(n)/log(2)? We've just seen how well this works.
How about:
log() - base 10
ln() - base e - logarithmus naturalis
ld() - base 2 - logarithmus dualis
rant_end
Libm has log(), log10, and log2(). The next
obvious question is: Does the inaccuracy reliably result in
values bigger than 11?
No.
If the input is a power of two, and you expect an integer as
a result, just do
k = (int)(log2(x) + 1e-6) or
k = (int)(log(x)/log(2) + 1e-6)
or
int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
Anyway, it's a neat way to avoid the problem and the rounding properties
of mult/div in case of not power of 2 could be useful as well.
30 Jan
30 Jan
2:29 p.m.
Excerpts from Peter Nelson's message of 2011-01-29 13:47:07 +0100:
On Sat, 2011-01-29 at 12:07 +0100, Philipp Überbacher
wrote:
Is it really surprising (for the non-optimised variant)? Your variant
uses one less variable and assignment. That this isn't optimised away
completely with -O3 does surprise me too.
Excerpts from fons's message of 2011-01-28
16:11:52 +0100:
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant.
Well, now I'm just being pedantic :-), but as a quick test using rdtsc
(i.e. profiling to be taken with a grain of salt):
1: log(x)/log(2)
2: (1 << k) < x
3: m < x & m <<= 1
This is 1000 iterations; cycling through x of 512, 1024, 2048, 4096 (to
prevent the compiler optimizing log(x)/log(2) to a single call
throughout the whole test). The fourth line is the sum of the iterator
in each loop.
./a.out (unoptimized)
1 3132000 cycles
2 261468 cycles
3 285273 cycles
1 50500, 2 50500, 3 50500
./a.out (-O3)
1 2598345 cycles
2 170055 cycles
3 173673 cycles
1 50500, 2 50500, 3 50500
I must admit, I'm surprised that fons way shows slightly slower in my
test! On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp
Überbacher wrote:
Took me a while to figure out that libm is part of glibc :)
Good to know that those functions are available on probably pretty much
all linux systems.
rant_begin
Why can't log mean the same thing everywhere? Why does it need to be
base e here and base 10 there? Why is there no consistency?
And why is there no proper logarithmus dualis function? Because you
can simply do log(n)/log(2)? We've just seen how well this works.
How about:
log() - base 10
ln() - base e - logarithmus naturalis
ld() - base 2 - logarithmus dualis
rant_end
Libm has log(), log10, and log2(). The next
obvious question is: Does the inaccuracy reliably result in
values bigger than 11?
No.
If the input is a power of two, and you expect an integer as
a result, just do
k = (int)(log2(x) + 1e-6) or
k = (int)(log(x)/log(2) + 1e-6)
or
int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
Anyway, it's a neat way to avoid the problem and the rounding properties
of mult/div in case of not power of 2 could be useful as well. As an aside, here's the result with (-O3
-ffast-math)
./a.out
1 150894 cycles,
2 169983 cycles,
3 158850 cycles,
1 47750, 2 50500, 3 50500
Yeah, faster, but wrong :-)
Peter.
After a brief look at man gcc -ffast-math doesn't seem like a very good
idea to use it in general, so no surprise it's not part of -O3 or
something.
Quite interesting stuff, thanks Peter.
1 Feb
1 Feb
2 p.m.
On Sun, Jan 30, 2011 at 02:29:13PM +0100, Philipp Überbacher wrote:
Excerpts from Peter Nelson's message of 2011-01-29
13:47:07 +0100:
i like -ffast-math... this whole thread is about float not being exact.
sure its worse with -ffast-math ... but if you dont assume correctness,
and do appropriate workarounds in case where you need correct results,
its the way to go.
On Sat, 2011-01-29 at 12:07 +0100, Philipp
Überbacher wrote:
Is it really surprising (for the non-optimised variant)? Your variant
uses one less variable and assignment. That this isn't optimised away
completely with -O3 does surprise me too.
Excerpts from fons's message of 2011-01-28
16:11:52 +0100:
Well, now I'm just being pedantic :-), but as a quick test using rdtsc
(i.e. profiling to be taken with a grain of salt):
1: log(x)/log(2)
2: (1 << k) < x
3: m < x & m <<= 1
This is 1000 iterations; cycling through x of 512, 1024, 2048, 4096 (to
prevent the compiler optimizing log(x)/log(2) to a single call
throughout the whole test). The fourth line is the sum of the iterator
in each loop.
./a.out (unoptimized)
1 3132000 cycles
2 261468 cycles
3 285273 cycles
1 50500, 2 50500, 3 50500
./a.out (-O3)
1 2598345 cycles
2 170055 cycles
3 173673 cycles
1 50500, 2 50500, 3 50500
I must admit, I'm surprised that fons way shows slightly slower in my
test! On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp
Überbacher wrote:
> rant_begin
> Why can't log mean the same thing everywhere? Why does it need to be
> base e here and base 10 there? Why is there no consistency?
> And why is there no proper logarithmus dualis function? Because you
> can simply do log(n)/log(2)? We've just seen how well this works.
> How about:
> log() - base 10
> ln() - base e - logarithmus naturalis
> ld() - base 2 - logarithmus dualis
> rant_end
Libm has log(), log10, and log2().
Took me a while to figure out that libm is part of glibc :)
Good to know that those functions are available on probably pretty much
all linux systems.
> The next obvious question is: Does the
inaccuracy reliably result in
> values bigger than 11?
No.
If the input is a power of two, and you expect an integer as
a result, just do
k = (int)(log2(x) + 1e-6)
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant.
or
k = (int)(log(x)/log(2) + 1e-6)
or
int m, k;
for (k = 0, m = 1; m < x; k++, m <<= 1);
which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
Anyway, it's a neat way to avoid the problem and the rounding properties
of mult/div in case of not power of 2 could be useful as well. As an aside, here's the result with (-O3
-ffast-math)
./a.out
1 150894 cycles,
2 169983 cycles,
3 158850 cycles,
1 47750, 2 50500, 3 50500
Yeah, faster, but wrong :-)
Peter.
After a brief look at man gcc -ffast-math doesn't seem like a very good
idea to use it in general, so no surprise it's not part of -O3 or
something.
Quite interesting stuff, thanks Peter.
_______________________________________________
Linux-audio-user mailing list
Linux-audio-user(a)lists.linuxaudio.org
http://lists.linuxaudio.org/listinfo/linux-audio-user
--
torben Hohn
13 Feb
13 Feb
4:30 p.m.
Excerpts from torbenh's message of 2011-02-01 14:00:17 +0100:
On Sun, Jan 30, 2011 at 02:29:13PM +0100, Philipp
Überbacher wrote:
Since --ffast-math is a collection of other options I would have thought
that you'd pick the ones that make sense for your use case instead of
--ffast-math.
Excerpts from Peter Nelson's message of
2011-01-29 13:47:07 +0100:
i like -ffast-math... this whole thread is about float not being exact.
sure its worse with -ffast-math ... but if you dont assume correctness,
and do appropriate workarounds in case where you need correct results,
its the way to go. On Sat, 2011-01-29 at 12:07 +0100, Philipp
Überbacher wrote:
Is it really surprising (for the non-optimised variant)? Your variant
uses one less variable and assignment. That this isn't optimised away
completely with -O3 does surprise me too.
Excerpts from fons's message of 2011-01-28
16:11:52 +0100:
> On Fri, Jan 28, 2011 at 02:02:36PM +0100, Philipp Überbacher wrote:
>
> > rant_begin
> > Why can't log mean the same thing everywhere? Why does it need to be
> > base e here and base 10 there? Why is there no consistency?
> > And why is there no proper logarithmus dualis function? Because you
> > can simply do log(n)/log(2)? We've just seen how well this works.
> > How about:
> > log() - base 10
> > ln() - base e - logarithmus naturalis
> > ld() - base 2 - logarithmus dualis
> > rant_end
>
> Libm has log(), log10, and log2().
Took me a while to figure out that libm is part of glibc :)
Good to know that those functions are available on probably pretty much
all linux systems.
> > The next obvious question is: Does the inaccuracy reliably result in
> > values bigger than 11?
>
> No.
>
> If the input is a power of two, and you expect an integer as
> a result, just do
>
> k = (int)(log2(x) + 1e-6)
log2() suffers from the same problem? I somewhat dislike the idea of
adding a constant.
> or
>
> k = (int)(log(x)/log(2) + 1e-6)
>
> or
>
> int m, k;
> for (k = 0, m = 1; m < x; k++, m <<= 1);
>
> which will round up if x is not a power of 2.
Neat. I thought about it myself yesterday but my ideas weren't exactly
brilliant. One idea was to divide by 2, the other to use a small
lookup table for powers of 2. I don't really know about efficiency, but
I guess bit shifting is as efficient as it gets?
Anyway, it's a neat way to avoid the problem and the rounding properties
of mult/div in case of not power of 2 could be useful as well.
Well, now I'm just being pedantic :-), but as a quick test using rdtsc
(i.e. profiling to be taken with a grain of salt):
1: log(x)/log(2)
2: (1 << k) < x
3: m < x & m <<= 1
This is 1000 iterations; cycling through x of 512, 1024, 2048, 4096 (to
prevent the compiler optimizing log(x)/log(2) to a single call
throughout the whole test). The fourth line is the sum of the iterator
in each loop.
./a.out (unoptimized)
1 3132000 cycles
2 261468 cycles
3 285273 cycles
1 50500, 2 50500, 3 50500
./a.out (-O3)
1 2598345 cycles
2 170055 cycles
3 173673 cycles
1 50500, 2 50500, 3 50500
I must admit, I'm surprised that fons way shows slightly slower in my
test! As an aside, here's the result with (-O3
-ffast-math)
./a.out
1 150894 cycles,
2 169983 cycles,
3 158850 cycles,
1 47750, 2 50500, 3 50500
Yeah, faster, but wrong :-)
Peter.
After a brief look at man gcc -ffast-math doesn't seem like a very good
idea to use it in general, so no surprise it's not part of -O3 or
something.
4821
days inactive
4837
days old
linux-audio-user@lists.linuxaudio.org
23 comments
8 participants
participants (8)
-
fons＠kokkinizita.net -
Gabriel Beddingfield -
Gabriel M. Beddingfield
-
James Stone -
Paul Davis -
Peter Nelson -
Philipp Überbacher -
torbenh