Fantastic!

Thanks,
Mike

----- Original Message ----- 
From: "Istvan Varga" 
To: 
Sent: Tuesday, September 27, 2005 7:11 PM
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC 
assembly code


> Here is a comparison of the time it takes to generate 100,000,000 random
> numbers with three different generators, compiled with various compilers
> on Linux.
>
> rand31: linear congruential generator (m = 0x7FFFFFFF, a = 742938285)
> rand:   ANSI C rand()
>         (note that it is much slower if not compiled with -static)
> MT:     Mersenne Twister (MT19937)
>
> ------------------------------------------------------------------------
>
> GCC 4.1 beta (-O2 -march=pentium3 -static)
>
>   rand31: 1.703, rand: 4.559, MT: 2.501
>
> GCC 4.1 beta (-O2 -march=pentium3 -fomit-frame-pointer -static)
>
>   rand31: 1.731, rand: 4.561, MT: 2.321
>
> GCC 4.1 beta (-O2 -march=pentium3)
>
>   rand31: 1.702, rand: 7.004, MT: 2.500
>
> GCC 4.0.0 (-O2 -march=pentium3 -static)
>
>   rand31: 1.702, rand: 4.561, MT: 2.405
>
> GCC 4.0.0 (-O2 -march=pentium3 -fomit-frame-pointer -static)
>
>   rand31: 1.732, rand: 4.562, MT: 2.339
>
> GCC 3.3.4 (-O2 -march=pentium3 -static)
>
>   rand31: 1.863, rand: 4.172, MT: 2.579
>
> GCC 3.3.4 (-O2 -march=pentium3 -fomit-frame-pointer -static)
>
>   rand31: 1.730, rand: 4.172, MT: 2.385
>
> ICC 9.0 (-O2 -xK -static)
>
>   rand31: 1.841, rand: 4.410, MT: 2.474
>
> ------------------------------------------------------------------------
>
> #include 
> #include 
> #include 
> #include 
>
> #define CS_PURE __attribute__ ((__pure__))
> #define CS_NOINLINE __attribute__ ((__noinline__))
> #define PUBLIC
>
> typedef struct CsoundRandMTState_ {
>     int         mti;
>     uint32_t    mt[624];
> } CsoundRandMTState;
>
> PUBLIC CS_PURE int csoundRand31(int seedVal)
> {
>     uint64_t  tmp1;
>     uint32_t  tmp2;
>
>     /* x = (742938285 * x) % 0x7FFFFFFF */
>     tmp1 = (uint64_t) ((int32_t) seedVal * (int64_t) 742938285);
>     tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
>     tmp2 += (uint32_t) (tmp1 >> 31);
>     if ((int32_t) tmp2 < (int32_t) 0)
>       tmp2 = (tmp2 + (uint32_t) 1) & (uint32_t) 0x7FFFFFFF;
>     return (int) tmp2;
> }
>
> /* Period parameters */
>
> #define N           624
> #define M           397
> #define MATRIX_A    0x9908B0DFU     /* constant vector a */
> #define UPPER_MASK  0x80000000U     /* most significant w-r bits */
> #define LOWER_MASK  0x7FFFFFFFU     /* least significant r bits */
>
> /* initializes mt[N] with a seed */
>
> PUBLIC void csoundSeedRandMT(CsoundRandMTState *p, uint32_t seedVal)
> {
>     int       i;
>     uint32_t  x;
>
>     p->mt[0] = x = seedVal;
>     for (i = 1; i < N; i++) {
>       /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
>       /* In the previous versions, MSBs of the seed affect   */
>       /* only MSBs of the array mt[].                        */
>       /* 2002/01/09 modified by Makoto Matsumoto             */
>       x = ((uint32_t) 1812433253 * (x ^ (x >> 30)) + (uint32_t) i);
>       p->mt[i] = x;
>     }
>     p->mti = N;
> }
>
> static CS_NOINLINE void MT_update_state(CsoundRandMTState *p)
> {
>     /* mag01[x] = x * MATRIX_A  for x=0,1 */
>     const uint32_t  mag01[2] = { (uint32_t) 0, (uint32_t) MATRIX_A };
>     int       i;
>     uint32_t  y;
>     for (i = 0; i < (N - M); i++) {
>       y = (p->mt[i] & UPPER_MASK) | (p->mt[i + 1] & LOWER_MASK);
>       p->mt[i] = p->mt[i + M] ^ (y >> 1) ^ mag01[y & (uint32_t) 1];
>     }
>     for ( ; i < (N - 1); i++) {
>       y = (p->mt[i] & UPPER_MASK) | (p->mt[i + 1] & LOWER_MASK);
>       p->mt[i] = p->mt[i + (M - N)] ^ (y >> 1) ^ mag01[y & (uint32_t) 1];
>     }
>     y = (p->mt[N - 1] & UPPER_MASK) | (p->mt[0] & LOWER_MASK);
>     p->mt[N - 1] = p->mt[M - 1] ^ (y >> 1) ^ mag01[y & (uint32_t) 1];
> }
>
> /* generates a random number on [0,0xffffffff]-interval */
>
> PUBLIC uint32_t csoundRandMT(CsoundRandMTState *p)
> {
>     int       i = p->mti;
>     uint32_t  y;
>
>     if (i >= N) {                   /* generate N words at one time */
>       MT_update_state(p);
>       i = 0;
>     }
>     y = p->mt[i];
>     p->mti = i + 1;
>     /* Tempering */
>     y ^= (y >> 11);
>     y ^= (y << 7) & (uint32_t) 0x9D2C5680U;
>     y ^= (y << 15) & (uint32_t) 0xEFC60000U;
>     y ^= (y >> 18);
>
>     return y;
> }
>
> volatile uint32_t seed_;
>
> int main(void)
> {
>     CsoundRandMTState st;
>     double t0, t1, t2, t3;
>     struct timeval tv;
>     int i;
>
>     gettimeofday(&tv, NULL);
>     t0 = (double) tv.tv_usec * 0.000001 + (double) tv.tv_sec;
>     seed_ = 12345;
>     for (i = 0; i < 100000000; i++)
>       seed_ = csoundRand31(seed_);
>     gettimeofday(&tv, NULL);
>     t1 = (double) tv.tv_usec * 0.000001 + (double) tv.tv_sec;
>     srand(314159);
>     for (i = 0; i < 100000000; i++)
>       seed_ = rand();
>     gettimeofday(&tv, NULL);
>     t2 = (double) tv.tv_usec * 0.000001 + (double) tv.tv_sec;
>     csoundSeedRandMT(&st, 5489);
>     for (i = 0; i < 100000000; i++) {
>       seed_ = csoundRandMT(&st);
> /*    if (i < 500) {
>         printf("%10u ", seed_);
>         if ((i % 5) == 4) printf("\n");
>       } */
>     }
>     gettimeofday(&tv, NULL);
>     t3 = (double) tv.tv_usec * 0.000001 + (double) tv.tv_sec;
>     printf("rand31: %.3f, rand: %.3f, MT: %.3f\n",
>            (t1 - t0), (t2 - t1), (t3 - t2));
>     return 0;
> }
> -- 
> Send bugs reports to this list.
> To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk
> 


-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk