[SciPy-dev] fft code for testing

Chuck Harris Chuck.Harris at sdl.usu.edu
Sat May 31 14:06:22 CDT 2003


Hi,


> -----Original Message-----
> From:	Pearu Peterson [mailto:pearu at scipy.org]
> Sent:	Sat 5/31/2003 11:49 AM
> To:	scipy-dev at 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 at scipy.net
http://www.scipy.net/mailman/listinfo/scipy-dev



-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 3897 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/scipy-dev/attachments/20030531/efb46607/attachment.bin 


More information about the Scipy-dev mailing list