#include <stdio.h>
#include <stdlib.h>
#include <math.h>  /* fmod() */
#include <time.h>  /* clock() */
#define KK 1597             /* can be changed               */
#define PP 2147483647       /* 2^31-1                       */
#define HH 1/(2.0*PP)
#define B_LCG 16807         /* for LCG, can be changed      */
#define SHIFT(b,x) ( ((x)>>(31-b)) + ((x)<<(b))&PP )

unsigned long MODP(unsigned long z) {return (((z)&PP)+((z)>>31));}

/* internal buffer and status, initialized in su_dx_k() */
static unsigned long  XX[KK];  /* buffer for u_dx_k */
static int II;                  /* index */
static int TT;                  /* index needed for dx-k-s */
static int K12;                 /* used by u_dx_k3 */
static int K13, K23;            /* used by u_dx_k4 */

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

    II = KK-1;                              /* running index */
    TT = tt-1;                              /* changes needed for dx*/
    K12 = KK/2-1;                           /* used by u_dx_k3 */
    K13 = 2*KK/3-1;    K23 = KK/3-1;        /* used by u_dx_k4 */
}