/*  dx-120d.c : DX random number generators 
            produces uniform random numbers (0,1) using DX-120-s
       
Note: specific efficient generators using bit operations on 32-bit unsigned long integers
    static unsigned long XX[KK];
    K =120 and P=2^31-1    
    see Table 1(d) of Deng and Xu (2002)
    
Authors: Lih-Yuan Deng (lihdeng@memphis.edu) and Hongquan Xu (hqxu@stat.ucla.edu)
Date: April 25, 2003

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  /* fmod() */
#include <time.h>  /* clock() */


/* Section 5.2: Initialization */ 
#define KK 120          
#define PP 2147483647   /* 2^31-1  */
#define HH 1/(2.0*PP)
#define B_LCG 16807     /* for LCG, can be changed      */

/* internal buffer and status, initialized in su_dx() */
static unsigned long XX[KK];           /* buffer */
static int II;                  /* index */
static int K12;                 /* used by u_dx3 */
static int K13, K23;            /* used by u_dx4 */

/* su_dx : Initialization of dx-k-s, using an LCG */
void su_dx(unsigned int seed)
{
    int i;
    if(seed==0) seed = 12345;           /* seed can't be zero */
    XX[0] = seed;
    for(i=1; i<KK; i++) XX[i] = fmod( (double)B_LCG *  XX[i-1],  PP);

    II = KK-1;                          /* running index */
    K12 = KK/2-1;                       /* used by u_dx3 */
    K13 = KK/3-1;    K23 = 2*KK/3-1;    /* used by u_dx4 */
}

/* Section 5.6: More efficient generators */
unsigned long MODP(unsigned long z) {return (((z)&PP)+((z)>>31));}

/* dx-120-1(d): (r,w)=(26, -22) */
#define MUL26(x) ( ((x)>>5) + ((x)<<26)&PP )
#define MUL22(x)  ( ((x)>>9) + ((x)<<22) &PP )
double u_dx1d()
{
    int II0 = II; unsigned long X ;
    if(++II >= KK)  II = 0;
    X = MODP( MUL26(XX[II]) +(PP-MUL22(XX[II])) );
    XX[II] = MODP(X + XX[II0]);
    return ((double) XX[II] /PP) + HH;
}

/* dx-120-1(e): (r,w)=(28, -16) */
#define MUL28(x) ( ((x)>>3) + ((x)<<28)&PP )
#define MUL16(x)  ( ((x)>>15) + ((x)<<16) &PP )
double u_dx1e()
{
    int II0 = II; unsigned long X ;
    if(++II >= KK)  II = 0;
    X = MODP( MUL28(XX[II]) +(PP-MUL16(XX[II])) );
    XX[II] = MODP(X + XX[II0]);
    return ((double) XX[II] /PP) + HH;
}

/* dx-120-2(d): (r,w)=(20, 9) */
#define MUL20(x) ( ((x)>>11) + ((x)<<20)&PP )
#define MUL9(x)  ( ((x)>>22) + ((x)<<9) &PP )
double u_dx2d()
{
    int II0 = II; unsigned long S ;
    if(++II >= KK)  II = 0;
    S = MODP(XX[II] + XX[II0]);
    XX[II] = MODP( MUL20(S) + MUL9(S));
    return ((double) XX[II] /PP) + HH;
}

/* dx-120-2(e): (r,w)=(28, -16) */ /* MUL28 and MUL16 have already been defined */
//#define MUL28(x) ( ((x)>>3) + ((x)<<28)&PP )  
//#define MUL16(x)  ( ((x)>>15) + ((x)<<16) &PP )
double u_dx2e()
{
    int II0 = II; unsigned long S ;
    if(++II >= KK)  II = 0;
    S = MODP(XX[II] + XX[II0]);
    XX[II] = MODP( MUL28(S) + (PP-MUL16(S)));
    return ((double) XX[II] /PP) + HH;
}

/* dx-120-3(d): (r,w)=(21, -8) */
#define MUL21(x) ( ((x)>>10) + ((x)<<21)&PP )  
#define MUL8(x)  ( ((x)>>23) + ((x)<<8) &PP )
double u_dx3d()
{
    int II0 = II; unsigned long S ;
    if(++II >= KK)  II = 0;     /*wrap around running index */
    if(++K12 >= KK) K12 = 0;    /*wrap around K12*/
    S = MODP(XX[II] + XX[K12]);   S = MODP(S + XX[II0]);
    XX[II] = MODP( MUL21(S) + (PP-MUL8(S)));
    return ((double) XX[II] /PP) + HH;
}

/* dx-120-4(d): (r,w)=(21, 17) */
//#define MUL21(x) ( ((x)>>10) + ((x)<<21)&PP )   /* MUL21 has already been defined */
#define MUL17(x)  ( ((x)>>14) + ((x)<<17) &PP )
double u_dx4d()
{
    int II0 = II; unsigned long S ;
    if(++II >= KK)   II = 0;    /*wrap around running index */
    if(++K13 >= KK) K13 = 0;    /*wrap around K13*/
    if(++K23 >= KK) K23 = 0;    /*wrap around K23*/
    XX[II] = MODP(XX[II]+XX[K13]); S = MODP(XX[K23]+XX[II0]);  S = MODP(S + XX[II]);
    XX[II] = MODP( MUL21(S) + MUL17(S));
    return ((double) XX[II] /PP) + HH;
}

#define _VERIFY_
#ifndef _VERIFY_
/* sample main() function */
int main()
{
    int i, runs=10;
    printf("=====> The first %d random numbers of u_dx1d() with B=2^26 - 2^22\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx1d() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx1e() with B=2^28 - 2^16\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx1e() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx2d() with B=2^20 + 2^9\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx2d() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx2e() with B=2^28 - 2^16\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx2e() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx3d() with B=2^21 - 2^8\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx3d() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx4d() with B=2^21 - 2^17\n", runs );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx4d() );  if(i%5 == 4) printf("\n"); }

    printf("\n");
    return 0;
}           
#else
void verify(int seed, int nTry, int IdxS, int BB, double (*u_dx)() )
{
    int i;
    clock_t t0, t1;
    su_dx(seed); t0 = clock();     for(i=0; i<nTry; i++) u_dx();     t1 = clock();
    printf("DX-120-%d (B=%d, %d numbers), time =%.2lfs\n", IdxS, BB, nTry, (double)(t1-t0)/CLOCKS_PER_SEC );
    for(i=0; i<5; i++) printf("%16.12lf", u_dx());
    printf("\n");
}

int main(int argc, char*argv[])
{
    int i, seed, nTry = 2000000;
    
    if(argc < 2){
        printf("[Usage] %s seed [nTry]\n", argv[0]);
        return 0;
        }
    seed = atoi(argv[1]);
    if(argc > 2) nTry = atoi(argv[2]);

    printf("=====> The next 5 numbers after %d trys\n", nTry );
    verify(seed, nTry, 1, (1<<26) - (1<<22), u_dx1d);
    verify(seed, nTry, 1, (1<<28) - (1<<16), u_dx1e);
    verify(seed, nTry, 2, (1<<20) + (1<<9), u_dx2d);
    verify(seed, nTry, 2, (1<<28) - (1<<16), u_dx2e);
    verify(seed, nTry, 3, (1<<21) - (1<<8), u_dx3d);
    verify(seed, nTry, 4, (1<<21) + (1<<17), u_dx4d);
    return 0;
}
#endif