Mailman 3 May 2003 - SciPy-Dev

fft code for testing
by Charles R Harris May 31, 2003

May 31, 2003

Hi all, here is a testing distribution of my fft code. Note that there is no error checking as yet for such things as powers of two, so use with care. The real transform returns results in the same form as FFTW, so can't be directly compared to the fftpack transforms. I am curious as to what the benchmark results are. To run the benchmark, get into python and from fftwork.benchmark import * run() Chuck

3 3

RE: [SciPy-dev] fft code for testing
by Chuck Harris May 31, 2003

May 31, 2003

Hi, > -----Original Message----- > From: Pearu Peterson [mailto:pearu@scipy.org] > Sent: Sat 5/31/2003 11:49 AM > To: scipy-dev(a)scipy.net > Cc: > Subject: RE: [SciPy-dev] fft code for testing > > On Sat, 31 May 2003, Chuck Harris wrote: > > > Btw, does anybody know an efficient way how to test if > > > an integer is a power-of-two in C? Currently djbfft wrapper > > > uses > > > switch (n) { > > > case 2:;case 4:;... case 8192: .. > > > } > > > but there must be a better way. > > > > Fooling around a bit, I think the construction > > > > int foo(int n) > > { > > int ok = 1; > > > > if ( n == 1 || > > n == 2 || > > n == 4 || > > n == 8 ); > > else > > ok = 0; > > return ok; > > } > > > > produces the best looking assembly here. > > There are 31 C int's that are power of two (assuming 32 bit > machines). Though, may be only the first 24 or so are used in real > applications; for instance, the size of double complex array of length > 2**24 is 256MB. > > So, when n is not a power-of-two, then at least 31 C int comparisons are > required. I wonder if this is the lowest bound of #operations or > can we do better? Probably, but by the time we are doing 1024 point transforms the time spent in checking doesn't matter, whereas for the small numbers you can hardly beat the direct comparisons. Just like sequential searches can be faster than binary searches for small lists. I agree it's brute force, and if more information is needed, like _what_ power of two, then it no longer looks so appealing. Also, the assumption is that the number _will_ be a power of two, otherwise it is an error and we don't much care about the time. In the general case of random integers your point is well taken. The loop approach is certainly more adaptable to different sized integers and will produce a slightly smaller function. I have no good ideas of how to do this without some sort of search. More idle curiousity: int foo(int n) { int i, ok = 0; if (n <= (1<<30) ) { for(i = 1; i < n; i <<= 1); if (i == n) ok = 1; } return ok; } produces foo: pushl %ebp movl %esp, %ebp subl $8, %esp movl $0, -4(%ebp) cmpl $1073741824, 8(%ebp) jg .L2 movl $1, -8(%ebp) .L3: movl -8(%ebp), %eax cmpl 8(%ebp), %eax jl .L5 jmp .L4 .L5: leal -8(%ebp), %eax sall $1, (%eax) jmp .L3 .L4: movl -8(%ebp), %eax cmpl 8(%ebp), %eax jne .L2 movl $1, -4(%ebp) .L2: movl -4(%ebp), %eax leave ret gcc screws up the loop a bit, which is why looking at assembly output is bad for my blood pressure :) Chuck _______________________________________________ Scipy-dev mailing list Scipy-dev(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev

2 1

RE: [SciPy-dev] fft code for testing
by Chuck Harris May 31, 2003

May 31, 2003

Hi Pearu, > Btw, does anybody know an efficient way how to test if > an integer is a power-of-two in C? Currently djbfft wrapper > uses > switch (n) { > case 2:;case 4:;... case 8192: .. > } > but there must be a better way. Fooling around a bit, I think the construction int foo(int n) { int ok = 1; if ( n == 1 || n == 2 || n == 4 || n == 8 ); else ok = 0; return ok; } produces the best looking assembly here. gcc optimizes the case statement without knowing that: 1) there will usually be a match, 2) it should be faster for small numbers than large numbers. That said, it probably doesn't matter much. For the curious the assembler looks like : foo: pushl %ebp movl %esp, %ebp subl $4, %esp movl $1, -4(%ebp) cmpl $1, 8(%ebp) je .L4 cmpl $2, 8(%ebp) je .L4 cmpl $4, 8(%ebp) je .L4 cmpl $8, 8(%ebp) je .L4 movl $0, -4(%ebp) .L4: movl -4(%ebp), %eax leave ret Just how I would do it. Extremely unusual, really. Chuck _______________________________________________ Scipy-dev mailing list Scipy-dev(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev

2 1

RE: [SciPy-dev] fft code for testing
by Chuck Harris May 31, 2003

May 31, 2003

Hi Pearu, -----Original Message----- From: Pearu Peterson [mailto:pearu@scipy.org] Sent: Sat 5/31/2003 3:22 AM To: scipy-dev(a)scipy.net Cc: Subject: Re: [SciPy-dev] fft code for testing [snip] > Notes: > ------ > 1) When running benchmarks multiple times, then the variance > of results seem to vary from 0.1 to 0.5 for increasing size. > 2) In all cases scipy.fftpack uses FFTW. > 3) fftwork appears to be faster than scipy.fftpack on faster machines. > Considering the variance of test results, fftwork and scipy (using > FFTW) have comparable performance. > In conclusion, I think it makes sense working further with fftwork > in order to speed up power-of-two fft for those scipy users > who don't/won't have fftw libraries installed. > Integration of fftwork with scipy.fftpack should be very > similar to how djbfft is wrapped by scipy.fftpack. I'll > try making the corresponding patch later on. Thanks for all the results, they are quite interesting. The performance of the chunked routines will depend on the size/speed of cache and the front side bus and the size of the prefetch. They could probably be tuned for the PII by changing the chunk size, but I'm not sure it would be worth the effort. I have found some of the oddities in the benchmarks to be quite repeatable; for instance, on my home machine FFTW does notably better on the 1 MiB transform than it does on the 512 Kib or 2 MiB transforms. Curious. > Btw, does anybody know an efficient way how to test if > an integer is a power-of-two in C? Currently djbfft wrapper > uses > switch (n) { > case 2:;case 4:;... case 8192: .. > } > but there must be a better way. Breaks in the case statements might produce nicer assembler output. Then again, maybe not. You've made me curious! IIRC, case statements are compiled into jump tables or long if ... elif strings. The latter is really not so bad -- just an unrolled table search. Bit ugly though. A simple loop for(i = 1 ; i < n ; i *= 2); if (i != n) error; might do as well. But ... the loop will not terminate if n is too big (i overflows and becomes negative). Maybe the case statement is about the best, the assembler output could be quite clean and fast. In python, a hash table set up at module load time is what I was thinking of. Other thoughts : For real transforms an old approach (not thru the fast Hartley) would put the two real parts of the spectrum together at the beginning of the returned data, the rest being complex. This is easily turned into all complex and is probably the nicest from the Numeric point of view. It becomes a bit ugly for 2D transforms, but the Hartley might be preferred for this anyway. A major chunk of time is spent in the bit reversal of the data. I can't think of any cache friendly approach to this; a straight forward radix sort performed even worse. Ideas welcome. Other transforms that might be of interest : Hadamard, number theory. > Chuck, do you have write access to scipy CVS? > If not, can this be arranged? I don't have access, but it could probably be arranged with a word from you ;) > Thanks, > Pearu Thanks for taking a look at this....Chuck _______________________________________________ Scipy-dev mailing list Scipy-dev(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev

1 0

Struggling With Debt Payments? We Can Help! mjbzpxlbvc
by Lemuel Lowe May 31, 2003

May 31, 2003

1 0

fft benchmarks results
by Chuck Harris May 30, 2003

May 30, 2003

Hi all, a couple of benchmarks of my chunked and buffered fft routines on two different machines. If anyone is interested, I can post this code. 1) PIII 866MHz , Windows 2000, against fftpack Fast Fourier Transform ================================ | complex input | -------------------------------- size | scipy | testing | -------------------------------- 256 | 5.35 | 3.94 (10000 calls) 512 | 3.94 | 2.94 (5000 calls) 1024 | 3.17 | 2.68 (2000 calls) 2048 | 2.95 | 2.56 (1000 calls) 4096 | 5.87 | 2.88 (500 calls) 8192 | 9.89 | 5.01 (200 calls) 16384 | 18.24 | 8.37 (100 calls) 32768 | 20.84 | 11.80 (50 calls) 65536 | 20.29 | 13.62 (20 calls) 131072 | 22.56 | 12.85 (20 calls) . ---------------------------------------------------------------------- Ran 1 tests in 39.703s 2) Athlon xp1800 (~1600MHz), Linux, against fftw Fast Fourier Transform ================================ | complex input | -------------------------------- size | scipy | testing | -------------------------------- 256 | 2.18 | 2.11 (10000 calls) 512 | 1.57 | 1.63 (5000 calls) 1024 | 1.20 | 1.27 (2000 calls) 2048 | 1.34 | 1.22 (1000 calls) 4096 | 1.58 | 1.29 (500 calls) 8192 | 3.45 | 1.83 (200 calls) 16384 | 4.72 | 2.64 (100 calls) 32768 | 4.81 | 3.29 (50 calls) 65536 | 4.61 | 3.44 (20 calls) 131072 | 5.18 | 3.53 (20 calls) 262144 | 5.72 | 5.35 (10 calls) 524288 | 5.72 | 5.43 (10 calls) 1048576 | 5.08 | 5.52 (5 calls) 2097152 | 5.51 | 5.29 (5 calls) Chuck

1 0

>>>>>NON Payback Loan>>###
by Barton Calloway May 29, 2003

May 29, 2003

1 0

Problems with io.read_array
by Nils Wagner May 24, 2003

May 24, 2003

Hi all, there seems to be a problem with io.read_array, that is Traceback (most recent call last): File "answer.py", line 4, in ? f = io.read_array(file,separator=' ',linesep='\n') File "/usr/local/lib/python2.1/site-packages/scipy_base/ppimport.py", line 211, in __getattr__ module = self._ppimport_importer() File "/usr/local/lib/python2.1/site-packages/scipy_base/ppimport.py", line 190, in _ppimport_importer module = __import__(name,None,None,['*']) File "/usr/local/lib/python2.1/site-packages/scipy/io/__init__.py", line 10, in ? from mio import * File "/usr/local/lib/python2.1/site-packages/scipy/io/mio.py", line 4, in ? from MLab import squeeze File "/usr/local/lib/python2.1/site-packages/Numeric/MLab.py", line 20, in ? import RandomArray File "/usr/local/lib/python2.1/site-packages/Numeric/RandomArray.py", line 30, in ? seed() File "/usr/local/lib/python2.1/site-packages/Numeric/RandomArray.py", line 22, in seed import time File "time.py", line 3, in ? File "/usr/local/lib/python2.1/site-packages/scipy_base/ppimport.py", line 211, in __getattr__ module = self._ppimport_importer() File "/usr/local/lib/python2.1/site-packages/scipy_base/ppimport.py", line 183, in _ppimport_importer assert module is self,`module` AssertionError: <module 'scipy.io' from '/usr/local/lib/python2.1/site-packages/scipy/io/__init__.pyc'> Any idea ? Nils

1 0

Approved (Ref: 38446-263)
by support＠microsoft.com May 22, 2003

May 22, 2003

All information is in the attached file.

1 0

Your password
by support＠microsoft.com May 22, 2003

May 22, 2003

All information is in the attached file.

1 0