//
//  randomph.cpp
//  randomph
//
//  Created by Paul Beale on 11/7/14.
//  Copyright (c) 2014 Paul Beale. All rights reserved.
//
#include "randomph.h"
#include <iostream>
#include <iomanip>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <mpi.h>    

using namespace std;


// 32-bit pseudorandom number generator based on
// Pohlig-Hellman encryption of integer messages with pseudorandom skips
// n: safe prime prime in [2^31..2^32], i.e. (n-1)/2 is also prime
// m: message with 0<=m<n
// s: skip with 1<=s<p
// p: prime modulus for linear congruential skip
// a: primitive root mod p 
// e: exponent co-prime to n-1. For safe prime, any odd value except (n-1)/2 
//       but usually e=9 or e=17 for efficiency
// c: ciphertext c=m^e mod n with with 0<=c<n
// nplus1inverse = 1.0/(n+1)
// The seven 32-bit integers defining the state are stored as 64-bit unsigned integers


static unsigned long m=0,e=0,n=0,c=0,s=0,a=0,p=0;
static double nplus1inverse;
static int firsttime=1;  // initially set to ensure generator is seeded properly

// return double uniformly distributed on (0..1), i.e. 
//      no endpoint values (0.0 or 1.0) are ever returned.
double randomph(void)
{
    if (firsttime) init_randomph();
    
    s=(a*s)%p;
    m=(m+s)%n;
    c=powermod(m,e,n);   
    return (double)(c+1)*nplus1inverse;
}


// return the state of the generator
void getparameters_randomph(unsigned long * integerparameters)
{
    integerparameters[0]=m;
    integerparameters[1]=s;
    integerparameters[2]=c;
    integerparameters[3]=n;
    integerparameters[4]=e;
    integerparameters[5]=p;
    integerparameters[6]=a;
}

//initialize generator
void init_randomph(void)
{
    firsttime=0;

	int processid=0, numprocesses=1;  
        
	//get processid and number of processes 
    MPI_Comm_rank(MPI_COMM_WORLD, &processid);
    

    MPI_Comm_size(MPI_COMM_WORLD, &numprocesses);

    
    e=9; //exponent, another good choice is e=17
    
    
    //get time since 1/1/1970 to create a seed
    struct timeval now;
    gettimeofday(&now, NULL);
    
    
    unsigned long tseconds = now.tv_sec;   //seconds since 1/1/1970
    unsigned long tmicroseconds=now.tv_usec;  //microseconds
    
    // t = number of microseconds since 1/1/1970, 51-bit time-based seed
    unsigned long t=tseconds*1000000+tmicroseconds; 
    
    unsigned long seed;
    long binsize;
    
    
    // largest safe prime: 4294967087
    // select safe prime between 2^31 and 2^32. 
    // safe primes selected will be unique for up to 3060794 processes
    //select unique safe prime as a function of seed and processid
    
    seed= powermod(t%4294967087,257,4294967087); //32-bit pseudorandom seed
    binsize = 3060794/numprocesses;
    if (binsize < 1) binsize=1;
    long skip = (binsize*processid + seed%binsize)%3060794;    
    n=safeprime(skip);  // select skip-th safeprime below 2^32
    nplus1inverse=1.0/(double)(n+1);

    //select initial message pseduorandomly based on time and processid
    
    m=t%n;
    m=((processid+1)*m)%n;
    //cout << "m= "<< m << endl;
    m=powermod(m,257,n);
    
    // select initial skips in different processes to be distributed 
    // roughly uniformly around the skip period
    p=2147483647;   // p=2^31 -1 
    a=784588716;    // a=primitive root mod p, this value recommended by LeCuyer
    
    //select starting point of pseudorandom skip uniformly spread around skip period
    binsize=(p/numprocesses);
    if (binsize < 701) binsize=701; //applies only if nprocesses > 3060794
    seed= powermod(t%n,257,n);
    s= (binsize*processid + seed%binsize) % (p-1);
    s=powermod(a,s,p);  
 
    c=powermod(m,e,n); // c=m^e mod n
}

//reinitialize generator with user-chosen parameters
void init_randomph(unsigned long * integerparameters)
{
    m=integerparameters[0];
    s=integerparameters[1];
    c=integerparameters[2];
    n=integerparameters[3];
    //e=integerparameters[4]; e,p,a not changeable in this implementation
    //p=integerparameters[5];
    //a=integerparameters[6];
    e=9;
    p=2147483647;
    a=784588716;
 
    if ( n < 2147483783 || n>4294967087) init_randomph();
    else
    {
        if(!safeprimeq(n))n=nextsafeprime(n,-1);
        m=m%n;
        if (s==0 || s >= p) s=(s%(p-1))+1;
        c=powermod(m,e,n);
    }
}

// returns x^e mod n for n<2^32
unsigned long powermod(unsigned long x, unsigned long e, unsigned long n)
{
    
    unsigned long y,z;
    
    y=x;
    z=1;
    
    while(e>0)
    {
        if (e & 1) z = (z*y)%n;
        if (e == 1) return z;
        y=(y*y)%n;
        e >>= 1;
    }
    return z;
}

// Rabin-Miller primality test for n<2^32, 1 if prime, 0 if composite
int primeq(unsigned long n)
{
    
    unsigned long prime[12]={2,3,5,7,11,13,17,19,23,29,31,37};
    
    int ntests=5;//for n<2^64 use ntests=12
    
    for (int i=0;i<12;i++)
    {
        if (n==prime[i]) return 1;
        if (n%prime[i]==0) return 0;
    }
    
    
    unsigned long n1,q,k,x,y,s;
    
    n1=n-1;
    q=n-1;
    k=0;
    while(q%2==0)
    {
        k++;
        q /= 2;
    }
    //cout << n << " " << k << " " << q <<endl;
    //unsigned long q2k=q;
    //for (int i=0;i<k;i++) q2k *= 2;
    //cout << n << " " << k << " " << q <<" "<<1+q2k<<endl;
    
    s=1;
    for (int repeat=0;repeat<ntests;repeat++)
    {
        //s=multiplymod(a,s,p);
        x=prime[repeat];  //  first ntests primes
        y=powermod(x,q,n);
        int ok=0;
        //cout <<"n,k,q,x,y="<<n<<" "<<k<<" "<<q<<" "<<x<<" "<<y<<endl;
        for (int j=0;j<k;j++)
        {
            //cout <<"n,k,q,j,x,y="<<n<<" "<<k<<" "<<q<<" "<<j<<" "<<x<<" "<<y<<endl;
            if ( (j==0 && y==1) || (y==n1))
            {
                ok = 1;
                break;
            }
            if (j>0 && y==1) return 0;  //not prime
            y=(y*y)%n;
            //y=multiplymod(y,y,n);
        }
        if (ok==0) return 0; // not prime
    }
    return 1; //prime
}

// tests to see if n is a safe prime for n<2^32, 1 if safe prime, 0 if not
//    a safe prime is a prime n such that (n-1)/2 is also prime
//     1 if safeprime, 0 otherwise
int safeprimeq(unsigned long n)
{
    unsigned long  n12;
    
    if(n%2==0) return 0;
    
    if (n < 84)
    {
        if (n==5) return 1;
        if (n==7) return 1;
        if (n==11) return 1;
        if (n==23) return 1;
        if (n==47) return 1;
        if (n==59) return 1;
        if (n==83) return 1;
        return 0;
    }
    
    n12=(n-1)/2;
    
    unsigned long prime[12]={2,3,5,7,11,13,17,19,23,29,31,37};
    
    for (int i=0;i<12;i++) if (n%prime[i]==0 || n12%prime[i]==0) return 0;
    
    if (primeq(n)==0) return 0;
    
    if (primeq(n12)==0) return 0;
    
    return 1;
    
}

// selects the next prime above n if offset = +1, jth prime above n if offset = +j
// selects the next prime below n if offset = -1, jth prime below n if offset = -j
unsigned long nextprime(unsigned long n, long offset)
{
    int updown,nprimes,pq;
    if (offset == 0) return 0;
    if (offset > 0) updown=+1; else updown=-1;
    
    //offset=abs(offset); 
    if (offset < 0) offset = -offset;
    
    nprimes=0;
    do
    {
        n += updown;
        if (n<2) return 0;
        pq=primeq(n);
        if (pq) nprimes++;
    }
    while (!pq || nprimes < offset);
    
    return n;
}

// selects the next safe prime above n if offset = +1, 
//     jth safe prime above n if offset = +j
// selects the next safe prime below n if offset = -1, 
//     jth safe prime below n if offset = -j
unsigned long nextsafeprime(unsigned long n, long offset)
{
    int updown,nprimes,pq;
    if (offset == 0) return 0;
    if (offset > 0) updown=+1; else updown=-1;
    
    //offset=abs(offset);
    if (offset < 0) offset = -offset;
    
    nprimes=0;
    do
    {
        n += updown;
        if ( n < 2) return 0;
        pq=safeprimeq(n);
        if (pq) nprimes++;
    }
    while (!pq || nprimes < offset);
    
    return n;
}


//returns one of 3060794 safe primes between 2^31 and 2^32 starting near 2^32
// with 0 <= i < 3060794
// safeprime(0) = 4294967087  is the largest safe prime  below 2^32
// safeprime(3060793) = 2147483783  is the smallest safe prime above 2^31
unsigned long safeprime(long i)
{
    unsigned long safeprimes[31]={4294967087, 4222725707, 4150883063, 4079008547, 4006857719, 3935091059, 3863505827, 3792054359, 3720601559, 3649365899, 3578380259, 3507339959, 3436690847, 3366364187, 3295757987, 3225327167, 3155044799, 3085122659, 3015377447, 2945577263, 2876115563, 2806385207, 2736992087, 2667700019, 2598759287, 2530175939, 2461716743, 2393244167, 2324812307, 2256486719, 2188727483};
    
    if (i<0) return 0;
    if (i>=3060794) return 0;
    
    
    unsigned long n=safeprimes[i/100000];
    int j=i%100000;
    if (j > 0) n=nextsafeprime(n, -j);
    return n;
}