View previous topic :: View next topic 
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(). 

Back to top 


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.");
// BitReverse operation
S = 32cast(int)wtemp;
for( ii=0; ii<data.length; ii++ ){
jj = (bitReverse(ii)>>S);
if( jj>ii ){
temp = data[ii];
data[ii] = data[jj];
data[jj] = temp;
}
}
while( nb>0 ){
...
Tiago 

Back to top 


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.");
// BitReverse operation
S = 32cast(int)wtemp;
for( ii=0; ii<data.length; ii++ ){
jj = (bitReverse(ii)>>S);
if( jj>ii ){
temp = data[ii];
data[ii] = data[jj];
data[jj] = temp;
}
}
while( nb>0 ){
...
Tiago 

Back to top 


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 twoscomplement 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? 

Back to top 


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 twoscomplement 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 clockcycles.
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 

Back to top 


tgasiba
Joined: 25 Nov 2005 Posts: 5 Location: Germany

Posted: Tue Nov 29, 2005 3:24 am Post subject: Fast WalshHadamard Transform 


By the way, I forgot to mention in my previous message that yesterday night I have finished coding a Fast WalshHadamard 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: 
/**
WalshHadamard 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<<L) )
throw new Error("Input data length must be a power of 2.");
S = (32L);
for( ii=0; ii<data.length/2; ii++ ){
jj = (bitReverse(ii)>>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<k3; i3++ ){
ii = i3 + L1 1;
jj = ii + k3;
temp1 = data[ii];
temp2 = data[jj];
if( (i2?2) == 0 ){
data[ii] = temp1  temp2;
data[jj] = temp1 + temp2;
} else {
data[ii] = temp1 + temp2;
data[jj] = temp1  temp2;
}
}
L1 += k1;
}
k1 >>= 1;
k2 <<= 1;
k3 >>= 1;
}
if( inverse==true )
for( ii=0; ii<data.length; ii++ )
data[ii] /= data.length;
}
}
alias _FastHadamardTransform!(double[]) FastHadamardTransform;
alias _FastHadamardTransform!(cdouble[]) FastHadamardTransform;
alias _FastHadamardTransform!(float[]) FastHadamardTransform;
alias _FastHadamardTransform!(cfloat[]) FastHadamardTransform;
alias _FastHadamardTransform!(int[]) FastHadamardTransform;
alias _FastHadamardTransform!(uint[]) FastHadamardTransform;

( I have just discovered how to include code hehehe)
Best,
Tiago 

Back to top 


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 opensource 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! 

Back to top 




You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum

Powered by phpBB © 2001, 2005 phpBB Group
