#include "./include/includes.h"
#include "./include/random.h"


/* 64 bit MT PRNG from mt19937-64.c */
#define NN 312
#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

extern int mti;
extern unsigned long long mt[NN];
void init_genrand64(unsigned long long seed);
double genrand64_real1(void);

double myrand()
{
    return genrand64_real1();
}

/* generate new rng.dat/ from rng.dat.init */
void make_rngdat()
{
    int idum=4711;
    unsigned long long seed;

    /* start MT */
    mti=NN+1;
    if (idum<0) idum=-idum;
    seed=(unsigned long long)idum;
    init_genrand64( seed );

    g_sync_KE();
    ranend();
}

void ranstart()
{
    make_rngdat();
}

void ranend()
{
}

void random_su2( su2_matr_comp * r_su2, double eps )
{
    double arg_cs, cs, sn, su2_select;

    arg_cs = 2 * M_PI * eps * ( myrand( 0 ) - 0.5 );
    cs = cos( arg_cs );
    sn = sin( arg_cs );
    r_su2->a[0] = cs;
    su2_select = myrand( 0 );
    if( su2_select < 0.3333333 )
    {
	r_su2->a[1] = sn;
	r_su2->a[2] = 0;
	r_su2->a[3] = 0;
    }
    else if( su2_select < 0.6666667 )
    {
	r_su2->a[1] = 0;
	r_su2->a[2] = sn;
	r_su2->a[3] = 0;
    }
    else
    {
	r_su2->a[1] = 0;
	r_su2->a[2] = 0;
	r_su2->a[3] = sn;
    }
}

void random_su3_KE( su3_matrix * r_su3, double eps )
{
    su2_matr_comp t;
    int a, b, index, i, j;

    /*  pick out an SU(2) subgroup */
    index = 3.0 * myrand( 0 );
    a = ( index + 1 ) % 3;
    b = ( index + 2 ) % 3;
    if( a > b )
    {
	i = a;
	a = b;
	b = i;
    }
    for ( i = 0; i < 3; i++ )
	for ( j = 0; j < 3; j++ )
	{
	    if( i == j )
	    {
		r_su3->ROWCOL( i, j ) = cmplx_KE( 1, 0 );
	    }
	    else
	    {
		r_su3->ROWCOL( i, j ) = cmplx_KE( 0, 0 );
	    }
	}
    random_su2( &t, eps );
    r_su3->ROWCOL( a, a ) = cmplx_KE( t.a[0], t.a[3] );
    r_su3->ROWCOL( a, b ) = cmplx_KE( t.a[2], t.a[1] );
    r_su3->ROWCOL( b, a ) = cmplx_KE( -t.a[2], t.a[1] );
    r_su3->ROWCOL( b, b ) = cmplx_KE( t.a[0], -t.a[3] );
}