#include <iostream>
#include "PhotosRandom.h"
#include "Photos.h"
#include "Log.h"

namespace Photospp
{

bool         PhotosRandom::init    = false;
int          PhotosRandom::iseed[2]= { 1802, 9373 };
int          PhotosRandom::i97     = 96;
int          PhotosRandom::j97     = 32;
double       PhotosRandom::uran[97]= { 0.0 };
double       PhotosRandom::cran    = 362436.0  /16777216.0;
const double PhotosRandom::cdran   = 7654321.0 /16777216.0;
const double PhotosRandom::cmran   = 16777213.0/16777216.0;

void PhotosRandom::setSeed(int s1,int s2)
{
	if(s1<0 || s1>31327) Log::Fatal("PhotosRandom::setSeed(): Seed(1) out of range [0,31327]",8);
	if(s2<0 || s2>30080) Log::Fatal("PhotosRandom::setSeed(): Seed(2) out of range [0,30080]",9);
	iseed[0]=s1;
	iseed[1]=s2;
}

/*******************************************************************************
  PHORIN:   PHOton radiation  in decays RANdom number generator init

  Purpose:  Initialse PHORAN  with  the user  specified seeds in the
            array iseed.  For details  see also:  F. James  CERN DD-
            Report November 1988.

  Author(s):  B. van Eijk and F. James        Created at:  27/09/89
                                              Last Update: 22/02/90
                                         Rewritten to C++: 18/10/10
                                         by T. Przedzinski  (tprzedzi@cern.ch)
*******************************************************************************/
void PhotosRandom::initialize()
{
	long IS1,IS2,IS3,IS4,IS5;
	double S,T;

// Calculate Marsaglia and Zaman seeds (by F. James)
	IS1=(iseed[0]/177)%177+2;
	IS2= iseed[0]%177+2;
	IS3=(iseed[1]/169)%178+1;
	IS4= iseed[1]%169;
	for(int i=0;i<97;i++)
	{
		S=0.0;
		T=0.5;
		for(int j=0;j<24;j++)
		{
			IS5=( ((IS1*IS2)%179)*IS3 )%179;
			IS1=IS2;
			IS2=IS3;
			IS3=IS5;
			IS4=(53*IS4+1)%169;
			if( (IS4*IS5)%64>=32) S=S+T;
			T=0.5*T;
		}
		uran[i]=S;
	}
	init=true;
	Log::Debug(0)<<"PhotosRandom::inititalize(): seed: "<<iseed[0]<<", "<<iseed[1]<<std::endl;
}

/*******************************************************************************
  PHORAN:   PHOton radiation in decays ret number generator based
            on Marsaglia Algorithm

  Purpose:  Generate  uniformly  distributed  random numbers between
            0 and 1.  Super long period:  2**144.  See also:
            G. Marsaglia and A. Zaman,  FSU-SCR-87-50,  for seed mo-
            difications  to  this version  see:  F. James DD-Report,
            November 1988.  The generator  has  to be initialized by
            a call to PHORIN ( C++ version: initialize() ).

  Author(s):  B. van Eijk, G. Marsaglia and   Created at:  27/09/89
              A. Zaman                        Last Update: 27/09/89
                                         Rewritten to C++: 18/10/10
                                         by T. Przedzinski  (tprzedzi@cern.ch)
*******************************************************************************/
double PhotosRandom::randomReal()
{
	if(!init) Log::Fatal("PhotosRandom::randomReal(): generator not initialized",1);
	double ret=0.0;
	while(true)
	{
		ret = uran[i97]-uran[j97];
		if(ret<0.0) ret+=1.;
		uran[i97]=ret;
		i97--;
		if(i97<0) i97=96;
		j97--;
		if(j97<0) j97=96;
		cran-=cdran;
		if(cran<0.0) cran+=cmran;
		ret-=cran;
		if(ret<0.0) ret+=1.0;
		if(ret>0.0) break;
	}
	return ret;
}

} // namespace Photospp