/*  dx-102.c : DX random number generators 
            produces uniform random numbers (0,1) using DX-102-s
       
Note: generic generators using double or 64-bit integers
    K =102 and P=2^31-1    
    see Table 1(a)(b)(c) 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 <time.h>  /* clock() */

/* Section 5.1: C directives for various implementations */ 

/*_int64/int64_t : machine dependent, see sys/types.h   */
/*_int64 used in MS-C, int64_t used in GNU-C            */
#if defined(_WIN32)
typedef _int64 XXTYPE;
#define DMOD(n, p) ((n) % (p))
#elif defined(__GNUC__)
typedef int64_t XXTYPE;
#define DMOD(n, p) ((n) % (p))
#else
#include <math.h>
typedef double XXTYPE;
#define DMOD(n, p) fmod((n), (p))
#endif

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

/* internal buffer and status, initialized in su_dx() */
static XXTYPE 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] = DMOD( 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.3: Generation */
#define BB1  1048554  /* from Table 1(c), (b) or (a) */
#define BB2  1047849
#define BB3   523905
#define BB4   524076
double u_dx1(void)
{
    int II0 = II;
    if(++II >= KK)  II = 0;     /*wrap around running index */
    XX[II] = DMOD(BB1 * XX[II] + XX[II0],  PP);
    return ((double) XX[II] /PP) + HH;
}
double u_dx2(void)
{
    int II0 = II;
    if(++II >= KK)  II = 0;     /*wrap around running index */
    XX[II] = DMOD(BB2 * (XX[II] + XX[II0]),  PP);
    return ((double) XX[II] /PP) + HH;
}
double u_dx3(void)
{
    int II0 = II;
    if(++II >= KK)  II = 0;     /*wrap around running index */
    if(++K12 >= KK) K12 = 0;    /*wrap around K12*/
    XX[II] = DMOD(BB3 * (XX[II] + XX[K12] + XX[II0]), PP);
    return ((double) XX[II] /PP) + HH;
}
double u_dx4(void)
{
    int II0 = II;
    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] = DMOD(BB4 * (XX[II]+XX[K13]+XX[K23]+XX[II0]), PP);
    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_dx1() with B=%d\n", runs, BB1 );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx1() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx2() with B=%d\n", runs, BB2 );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx2() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx3() with B=%d\n", runs, BB3 );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx3() );  if(i%5 == 4) printf("\n"); }
    
    printf("=====> The first %d random numbers of u_dx4() with B=%d\n", runs, BB4 );
    su_dx(12345);		/* initialization */
    for(i=0; i<runs; i++){ printf("%16.12lf", u_dx4() );  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-102-%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, BB1, u_dx1);
    verify(seed, nTry, 2, BB2, u_dx2);
    verify(seed, nTry, 3, BB3, u_dx3);
    verify(seed, nTry, 4, BB4, u_dx4);

#if defined(_WIN32)
    printf("NOTE: _int64 is used\n");
#elif defined(__GNUC__)
    printf("NOTE: int64_t is used\n");
#else
    printf("NOTE: double is used\n");
#endif
    return 0;
}
#endif