//  RAN2.C
//  This contains a slightly modified version of
//  the ran2() function given in Numerical Recipes, p.282.
//  © 2021 Peter J. Meyer

//  ran2() has a period of > 2 x 10^18 whereas
//  ran1() has a period only of c. 2 x 10^9.

#define true  1
#define false 0

#define IM1  (2147483563)
#define IM2  (2147483399)
#define AM   (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1  (40014)
#define IA2  (40692)
#define IQ1  (53668)
#define IQ2  (52774)
#define IR1  (12211)
#define IR2  (3791)
#define NTAB (32)
#define NDIV (1+IMM1/NTAB)
#define EPS  (1e-16)            //  Changed from 1.2e-7       
#define RNMX (1.0-EPS)

static double num_calls = 0.0;
static double period_length = 2e18;
static int period_length_exceeded = false;
static long xn32 = -1;

long iy=0;

void reset_num_ran2_calls(void);

//  Returns a "uniform random deviate" between 0.0 and 1.0 (exclusive).
/*-------------*/
double ran2(void)
{
long j, k;
static long ix=123456789;
//  static long iy=0;
static long iv[NTAB];
double result;

num_calls += 1.0;

if ( num_calls > period_length )
    period_length_exceeded = true;    

if ( xn32 <= 0 )
    {
    //  Initialize
    if ( xn32 == 0 )
        xn32 = 1;
    else    // xn32 < 0 
        xn32 = -xn32;

    ix = xn32;

    //  Load the shuffle table (after 8 warm-ups).
    for ( j=NTAB+7; j>=0; j-- )
        {
        k = xn32/IQ1;
        xn32 = IA1*(xn32-k*IQ1) - IR1*k;
        if ( xn32 < 0 )
            xn32 += IM1;
        if ( j < NTAB )
            iv[j] = xn32;
        }

    iy = iv[0];
    }

k = xn32/IQ1;                       //  Start here when not initializing.
xn32 = IA1*(xn32 - k*IQ1) - IR1*k;  //  Compute xn32 = (IA1*xn32)%IM1 without overflow.
if ( xn32 < 0 )
    xn32 += IM1;

k = ix/IQ2;                         //  Compute ix%IM2 likewise.
ix = IA2*(ix - k*IQ2) - IR2*k;
if ( ix < 0 )
    ix += IM2;

j = iy/NDIV;                        //  Will be in the range 0 through NTAB-1.
iy = iv[j] - ix;                    //  Here xn32 is shuffled and xn32
                                    //  and ix are combined to generate output.
iv[j] = xn32;                       //  and refill the shuffle table.

if ( iy < 1 )
    iy += IMM1;

result = AM*iy;  

return ( result > RNMX ? RNMX : result );   //  So as not to return 1.0.
}

/*-----------------------------*/
void init_ran2(unsigned int seed)
{
int i;

xn32 = seed;

xn32 = xn32*(xn32+1) + 7;               //  Make sure xn32 is odd.
if ( xn32 > 0 )                         //  Make sure xn32 is negative
    xn32 = -xn32;

ran2();                                 //  Initialise

for ( i=0; i<64; i++ )                  //  Run ran2() a bit.
    ran2();

reset_num_ran2_calls();                 //  Reset to 0.
}

/*---------------------------*/
void reset_num_ran2_calls(void)
{
num_calls = 0.0;
}

/*-----------------------*/
double num_ran2_calls(void)
{
return ( num_calls );
}

/*--------------------------*/
int ran2_period_exceeded(void)
{
return ( period_length_exceeded );
}