Author Message
Don Clugston

Joined: 05 Oct 2005
Posts: 91
Location: Germany (expat Australian)

 Posted: Fri Nov 25, 2005 9:52 am    Post subject: Bit operations They've got no business being in MathExtra, but I had to put them somewhere... bitoperations.d They supplement the ones in std.intrinsic. The most useful being bitcount().
tgasiba

Joined: 25 Nov 2005
Posts: 5
Location: Germany

 Posted: Mon Nov 28, 2005 11:10 am    Post subject: Re: Bit operations Ok. Here goes a new function. uint bitReverse( uint x ){ asm { mov EDX, EAX; shr EAX, 1; and EDX, 0x5555_5555; and EAX, 0x5555_5555; shl EDX, 1; or EAX, EDX; mov EDX, EAX; shr EAX, 2; and EDX, 0x3333_3333; and EAX, 0x3333_3333; shl EDX, 2; or EAX, EDX; mov EDX, EAX; shr EAX, 4; and EDX, 0x0f0f_0f0f; and EAX, 0x0f0f_0f0f; shl EDX, 4; or EAX, EDX; bswap EAX; } } which rewrites an integer with the bits reversed. Using this routine, the FastFourierTransform now becomes slightly more optimized: ... wtemp = floor(log2(cast(double)data.length)); if( data.length != cast(int)pow(2.0,wtemp) ) throw new Error("Input data length must be a power of 2."); // Bit-Reverse operation S = 32-cast(int)wtemp; for( ii=0; ii>S); if( jj>ii ){ temp = data[ii]; data[ii] = data[jj]; data[jj] = temp; } } while( nb>0 ){ ... Tiago
tgasiba

Joined: 25 Nov 2005
Posts: 5
Location: Germany

 Posted: Mon Nov 28, 2005 11:11 am    Post subject: Re: Bit operations Ok. Here goes a new function. uint bitReverse( uint x ){ asm { mov EDX, EAX; shr EAX, 1; and EDX, 0x5555_5555; and EAX, 0x5555_5555; shl EDX, 1; or EAX, EDX; mov EDX, EAX; shr EAX, 2; and EDX, 0x3333_3333; and EAX, 0x3333_3333; shl EDX, 2; or EAX, EDX; mov EDX, EAX; shr EAX, 4; and EDX, 0x0f0f_0f0f; and EAX, 0x0f0f_0f0f; shl EDX, 4; or EAX, EDX; bswap EAX; } } which rewrites an integer with the bits reversed. Using this routine, the FastFourierTransform now becomes slightly more optimized: ... int S; ... wtemp = floor(log2(cast(double)data.length)); if( data.length != cast(int)pow(2.0,wtemp) ) throw new Error("Input data length must be a power of 2."); // Bit-Reverse operation S = 32-cast(int)wtemp; for( ii=0; ii>S); if( jj>ii ){ temp = data[ii]; data[ii] = data[jj]; data[jj] = temp; } } while( nb>0 ){ ... Tiago
Don Clugston

Joined: 05 Oct 2005
Posts: 91
Location: Germany (expat Australian)

 Posted: Tue Nov 29, 2005 2:45 am    Post subject: Nice. Did you know that this: wtemp = floor(log2(cast(double)data.length)); if( data.length != cast(int)pow(2.0,wtemp) ) throw new Error("Input data length must be a power of 2."); can be done with: // returns 2^b such that 2^b <= x < 2^b+1, where b is an integer // Relies on twos-complement arithmetic. int nextlowerpowerof2(int x) { return x & -x; } if (data.length !=nextlowerpowerof2(data.length)) { throw new Error("Input data length must be a power of 2."); } Can you think of a better name for that function?
tgasiba

Joined: 25 Nov 2005
Posts: 5
Location: Germany

Posted: Tue Nov 29, 2005 3:15 am    Post subject: nextlowerpowerof2

 Don Clugston wrote: can be done with: // returns 2^b such that 2^b <= x < 2^b+1, where b is an integer // Relies on twos-complement arithmetic. int nextlowerpowerof2(int x) { return x & -x; } if (data.length !=nextlowerpowerof2(data.length)) { throw new Error("Input data length must be a power of 2."); } Can you think of a better name for that function?

That is very interesting. Never thought of that. I have also fixed it like this:

<snip>
L = cast(int)floor(log2(cast(double)data.length));
if( data.length != (1<<L) )
throw new Error("Input data length must be a power of 2.");
<snip>

But I think that your routine is faster, because the shift takes a few more clock-cycles.
The name "nextlowerpowerof2" does not seem very good because of the word "next" and "lower". I would suggest the name "nextpow2" as in Matlab.

Best,
Tiago
tgasiba

Joined: 25 Nov 2005
Posts: 5
Location: Germany

Posted: Tue Nov 29, 2005 3:24 am    Post subject: Fast Walsh-Hadamard Transform

By the way, I forgot to mention in my previous message that yesterday night I have finished coding a Fast Walsh-Hadamard Transform. It might also be nice to include into the standard library.
Another comment I have is the following:
The MathExtra project might get contributions that are away advanced for a standard library. I have never saw a programming language that comes already packed with FFT, FWHT or Galois Field arithmetic etc... The question would be the following. Are all the routines intended to be part of the standard library, will some routines be considered to be there and others (as the name indicated) be extra? If so, then which ones?
In any case, I think that all the routines should belong to the standard language. Is this a good idea?

 Code: /** Walsh-Hadamard Transform */ template _FastHadamardTransform( T : T[] ){   void _FastHadamardTransform( inout T[] data, bool inverse=false ) {     int    L, ii, jj;     int    k1, k2, k3;     int    i1, i2, i3;     int    L1, S;     T      temp1, temp2;     L = cast(int)floor(log2(cast(double)data.length));     if( data.length != (1<>S);       if( jj>ii ){         temp1    = data[ii];         data[ii] = data[jj];         data[jj] = temp1;       }     }     k1 = data.length;     k2 = 1;     k3 = (k1>>1);     for( i1=1; i1<=L; i1++ ){       for( L1=1, i2=1; i2<=k2; i2++ ){         for( i3=0; i3>= 1;       k2 <<= 1;       k3 >>= 1;     }     if( inverse==true )       for( ii=0; ii

( I have just discovered how to include code hehehe)

Best,
Tiago
Don Clugston

Joined: 05 Oct 2005
Posts: 91
Location: Germany (expat Australian)

 Posted: Tue Nov 29, 2005 7:49 am    Post subject: What should be included [quote]The MathExtra project might get contributions that are away advanced for a standard library. I have never saw a programming language that comes already packed with FFT, FWHT or Galois Field arithmetic etc... The question would be the following. Are all the routines intended to be part of the standard library, will some routines be considered to be there and others (as the name indicated) be extra? If so, then which ones? [/quote] Personally, I use the "Excel Test". If Microsoft Excel has it, it is reasonable for it to be in a standard library. This includes: * the statistical functions * FFT, but only for powers of 2. * basic matrix operations. The C standard library is quite impoverished mathematically. Since D has exceptional support for mathematics, and also because the standard library is open-source and redistributable, I think it's easy to argue that it should at least match Excel. There's actually some wierd things in the C standard library. Why on earth are Bessel functions in there?? I can think of a lot of things that are far more useful. I'm not picky. I'll accept anything into MathExtra, provided it's under some sort of unrestricted license. When parts of it stabilise, I'll try to get them into Phobos. But so far, all I've succeeded into getting into Phobos is the gamma function (which actually went in the C++ compiler, my D code disappeared!), and my feqrel() function. BTW, if you like the bit manipulation tricks, you should take a look at std.math.feqrel. I'm rather proud of it. It was also my very first D function. I must confess that I'm still a D newbie. I've had a big impact on the last couple of compiler releases, and seem to have invented some revolutionary template metaprogramming techniques, but ... (shock) I've never actually created a class in D!
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 6 Hours Page 1 of 1