//
//  main.cpp
//  randomph
//
//  Created by Paul Beale on 12/4/14.
//  Copyright (c) 2014 Paul Beale. All rights reserved.
//
// Test randomph() using a single process

#include "randomph.h"
#include <iostream>
#include <iomanip>
#include <math.h>
#include <time.h>
#include <sys/time.h>


using namespace std;


//determine the one-sided p-value from the chi square test for a uniform distribution with a large number of degrees of freedom
double ChiSquareLargeNPValue(int* histogram,int n)
{
    double degreesoffreedom=n-1;
    double average=0.0,chisqr=0.0,pvalue,sigma;
    for (int i=0;i<n;i++) average += histogram[i];
    average /= n;
    for (int i=0;i<n;i++) chisqr += pow(histogram[i]-average,2)/average;
    sigma = (chisqr-degreesoffreedom)/sqrt(2.*degreesoffreedom);
    pvalue=0.5*erfc(fabs(sigma)/sqrt(2.0));
    return pvalue;
}

int main(int argc, char ** argv)
{
    

    long ntests = 100000000L; //number of pseudorandom numbers used in 3D, 4D, and 5D chi squared tests, twice as many used in 2D test, and five times as many used in 1D test
    long nrounds=10; //number of rounds of tests

    double time=0.0,timelast=0.0;
    time =((double) clock())/((double) CLOCKS_PER_SEC); timelast=time;

        cout << "Pohlig-Hellman pseudorandom number generator based on "<<endl;
        cout << "'A new class of scalable parallel pseudorandom number generators"<<endl;
        cout<< "based on Pohlig-Hellman exponentiation ciphers,'"<<endl;
        cout << "Paul D. Beale, (2014)."<<endl;
        cout << "http://arxiv.org/abs/1411.2484"<<endl;
        cout << endl;


        cout << "Intrastream tests of sequences of "<< 5*ntests*nrounds <<" pseudorandom numbers"<<endl;
        cout << "using 1D, 2D, 3D, 4D, and 5D histogram tests"<<endl;
        cout << "With these parameters, the run time should be roughly " << 5*ntests*nrounds/8e6/60. << " minutes."<<endl<<endl;

    
    //histograms for 1D, 2D, 3D, 4D, and 5D correlations tests
    int * hist1d = new int[1048576];
    int * hist2d = new int[1048576];
    int * hist3d = new int[1000000];
    int * hist4d = new int[1048576];
    int * hist5d = new int[1048576];
    
    //countpvalues[0]=total number of number of pvalue tests
    //countpvalues[1]=number with pvalue < 1e-1
    //countpvalues[2]=number with pvalue < 1e-2
    //countpvalues[3]=number with pvalue < 1e-3
    //countpvalues[4]=number with pvalue < 1e-4
    //countpvalues[5]=number with pvalue < 1e-5
    //countpvalues[6]=number with pvalue < 1e-6
    //countpvalues[7]=number with pvalue < 1e-7
    //countpvalues[8]=number with pvalue < 1e-8

    long countpvalues[9]={0};
    unsigned long integerparameters[7];
    unsigned long m0, s0, c0, n0, e0,p0, a0;
    double pvalue,log10pvalue;
    
    randomph();
    getparameters_randomph(integerparameters);
    cout << endl<< "initial state of the generator: m, s, c, n, e, p, a" << endl;
    for (int j=0;j<7;j++) cout << integerparameters[j] << " "; cout << endl;
    m0=integerparameters[0];
    s0=integerparameters[1];
    c0=integerparameters[2];
    n0=integerparameters[3];
    e0=integerparameters[4];
    p0=integerparameters[5];
    a0=integerparameters[6];

    for (int irounds=0;irounds<nrounds;irounds++)
    {
    
    
    for (int i=0;i<1048576;i++) hist1d[i]=hist2d[i]=hist4d[i]=hist5d[i]=0;
    for (int i=0;i<1000000;i++) hist3d[i]=0;
    
     
    double r1,r2,r3,r4,r5;
    
    for (long itest=0;itest<ntests;itest++)
    {
        r1=randomph();
        r2=randomph();
        r3=randomph();
        r4=randomph();
        r5=randomph();
        hist1d[(int)(1048576.0*r1)]++;
        hist1d[(int)(1048576.0*r2)]++;
        hist1d[(int)(1048576.0*r3)]++;
        hist1d[(int)(1048576.0*r4)]++;
        hist1d[(int)(1048576.0*r5)]++;
        hist2d[1024*(int) (1024.0*r2) +(int)(1024*r1)]++;
        hist2d[1024*(int) (1024.0*r4) +(int)(1024*r3)]++;
        hist3d[10000*(int) (100.0*r3)+100*(int) (100.0*r2) +(int)(100.0*r1)]++;
        hist4d[32768*(int) (32.0*r4)+1024*(int) (32.0*r3)+32*(int) (32.0*r2) +(int)(32.0*r1)]++;
        hist5d[65536*(int)(16.0*r5)+4096*(int) (16.0*r4)+256*(int) (16.0*r3)+16*(int) (16.0*r2) +(int)(16.0*r1)]++;
    }
    
    //analyze the chi square values from the histograms
    pvalue=ChiSquareLargeNPValue(hist1d,1048576);
    log10pvalue = log10(pvalue);
    countpvalues[0]++;
    for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;
    

    pvalue=ChiSquareLargeNPValue(hist2d,1048576);
    log10pvalue = log10(pvalue);
    countpvalues[0]++;
    for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;
 
    pvalue = ChiSquareLargeNPValue(hist3d,1000000);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;

    pvalue = ChiSquareLargeNPValue(hist4d,1048576);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;

    pvalue = ChiSquareLargeNPValue(hist5d,1048576);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;


    }
    
    getparameters_randomph(integerparameters);
    unsigned long mf, sf, cf, q, r;
    q=(5*ntests*nrounds)/(p0-1);
    r=(5*ntests*nrounds)%(p0-1);
    sf=s0;
    mf=(p0*(p0-1)/2)%n0;
    mf= (q*mf)%n0;
    mf=(mf+m0)%n0;  // m_{q (p-1)} = (m_0 + q p(p-1)/2 ) mod n
    // step forward r steps
    for (int j=1;j<=r;j++)
    {
        sf=(a0*sf)%p0;
        mf=(mf+sf)%n0;
    }
    cf=powermod(mf,e0,n0);
    cout << "final state of the generator: m, s, c, n, e, p, a" << endl;
    for (int j=0;j<7;j++) cout << integerparameters[j] << " "; cout << endl;
    
    cout << "expected final state of the generator: m, s, c, n, e, p, a" <<endl; // see equation (7) of paper above.
    cout << mf << " " << sf << " " << cf << " " << n0 << " " << e0 << " " << p0 << " " << a0 << endl;

    
 
    cout << "total number of chi square tests = "<<countpvalues[0]<<endl;
    for (int ip=1;ip<=8;ip+=1) cout << "number of p<1e-"<<ip<<" events = "<<countpvalues[ip]<<",   expected number = " << 2.0*countpvalues[0]*pow(0.1,ip) << " +/- " << sqrt(2.0*countpvalues[0]*pow(0.1,ip))<<endl;

    
    cout <<endl;
    time =((double) clock())/((double) CLOCKS_PER_SEC);
    cout << "time = "<<time-timelast << " seconds"<<endl;
    cout << "total number of pseudorandom numbers tested = "<< 5L*ntests*nrounds <<endl;
    //cout << "number generated and analyzed per second in process 0 = " << (double) ntests * (double) nrounds * 5.0 / time << endl;
    timelast=time;
    
    
    
    return 0;

    }