source: branches/lib/GDE/CLUSTALW/random.c

/*
*       
*       Rand.c
*       
*       -       linear and additive congruential random number generators
*               (see R. Sedgewick, Algorithms, Chapter 35)
*
*       Implementation: R. Fuchs, EMBL Data Library, 1991
*       
*/
#include <stdio.h>

unsigned long linrand(unsigned long r);
unsigned long addrand(unsigned long r);
void addrandinit(unsigned long s);

static unsigned long mult(unsigned long p,unsigned long q);


#define m1      10000
#define m       100000000

static unsigned long mult(unsigned long p, unsigned long q);

/* linear congruential method
*       
*       linrand() returns an unsigned long random number in the range 0 to r-1
*/


unsigned long linrand(unsigned long r)
{
        static unsigned long a=1234567;
        
        a = (mult(a,31415821)+1) % m;
        return( ( (a / m1) * r) / m1 );
}

static unsigned long mult(unsigned long p, unsigned long q)
{
        unsigned long p1,p0,q1,q0;
        
        p1 = p/m1; p0 = p % m1;
        q1 = q/m1; q0 = q % m1;
        return((((p0*q1 + p1*q0) % m1) * m1 + p0*q0) % m);
}


/* additive congruential method
*       
*       addrand() returns an unsigned long random number in the range 0 to r-1
*       The random number generator is initialized by addrandinit()
*/

static unsigned long j;
static unsigned long a[55];

unsigned long addrand(unsigned long r)
{
int x,y;
/*        fprintf(stdout,"\n j = %d",j);  */
        j = (j + 1) % 55;
/*        fprintf(stdout,"\n j = %d",j);  */
        x = (j+23)%55;
        y = (j+54)%55;
        a[j] = (a[x] + a[y]) % m;
/*      a[j] = (a[(j+23)%55] + a[(j+54)%55]) % m;  */
/*        fprintf(stdout,"\n a[j] = %d",a[j]);     */
        return( ((a[j] / m1) * r) / m1 );
}

void addrandinit(unsigned long s)
{
        a[0] = s;
        j = 0;
        do {
                ++j;
                a[j] = (mult(31,a[j-1]) + 1) % m;
        } while (j<54);
}