FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Bit operations

 
Post new topic   Reply to topic     Forum Index -> MathExtra
View previous topic :: View next topic  
Author Message
Don Clugston



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

PostPosted: Fri Nov 25, 2005 9:52 am    Post subject: Bit operations Reply with quote

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
View user's profile Send private message
tgasiba



Joined: 25 Nov 2005
Posts: 5
Location: Germany

PostPosted: Mon Nov 28, 2005 11:10 am    Post subject: Re: Bit operations Reply with quote

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<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
View user's profile Send private message Visit poster's website AIM Address MSN Messenger
tgasiba



Joined: 25 Nov 2005
Posts: 5
Location: Germany

PostPosted: Mon Nov 28, 2005 11:11 am    Post subject: Re: Bit operations Reply with quote

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<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
View user's profile Send private message Visit poster's website AIM Address MSN Messenger
Don Clugston



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

PostPosted: Tue Nov 29, 2005 2:45 am    Post subject: Reply with quote

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?
Back to top
View user's profile Send private message
tgasiba



Joined: 25 Nov 2005
Posts: 5
Location: Germany

PostPosted: Tue Nov 29, 2005 3:15 am    Post subject: nextlowerpowerof2 Reply with quote

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. Smile

Best,
Tiago
Back to top
View user's profile Send private message Visit poster's website AIM Address MSN Messenger
tgasiba



Joined: 25 Nov 2005
Posts: 5
Location: Germany

PostPosted: Tue Nov 29, 2005 3:24 am    Post subject: Fast Walsh-Hadamard Transform Reply with quote

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<<L) )
      throw new Error("Input data length must be a power of 2.");

    S = (32-L);
    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;

( Very Happy I have just discovered how to include code hehehe)

Best,
Tiago
Back to top
View user's profile Send private message Visit poster's website AIM Address MSN Messenger
Don Clugston



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

PostPosted: Tue Nov 29, 2005 7:49 am    Post subject: What should be included Reply with quote

[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!
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic     Forum Index -> MathExtra All times are GMT - 6 Hours
Page 1 of 1

 
Jump to:  
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