#include <stdio.h>
#include <malloc.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <time.h>

// code output at http://www.rankyouragent.com/primes/patn.c.output.txt

const char* INFINITE_ZEROS = "";

long AbsoluteValue(const long in)
{
    if(in<0)
        return in*(-1);
    return in;
}

long FindFirstZero(const char* bitString )
{
    long firstZero = 1;
    long bitStringLength = strlen(bitString);
    long i=0;
    for (i = 0; i < bitStringLength; i++ )
      if ( bitString[i]=='0' ) return i + 1;
    return firstZero;
}

char* Elemental(long firstZero )
{
    long elementalSize = (2*firstZero) + 1;
    char* elemental = (char*) malloc(elementalSize+1);
    long i=0;
    for ( i = 0; i < elementalSize; i++ )
      if ( i == 0 )
        elemental[i] = '1';
      else
        elemental[i] = '0';
     elemental[elementalSize]='\0';
    return elemental;
}

const char* Merge(const char* lhsPattern, const char* rhsPattern ) {
    // Note 0... MERGE anyPattern = anyPattern
    long firstZeroLhsPattern = FindFirstZero( lhsPattern );
    long mergedPatternLength = strlen(lhsPattern) * strlen(rhsPattern);
    long lhsPatternLength = strlen(lhsPattern);
    long rhsPatternLength = strlen(rhsPattern);
    char* mergedPattern = (char*) malloc(mergedPatternLength+1);
    long i=0;
    time_t start=0,finish=0;
    start = time(NULL);
    if ( lhsPattern == INFINITE_ZEROS ) return rhsPattern;
    for ( i = 0; i < mergedPatternLength; i++ ) {
      char currLhsBit = '0';
      char currRhsBit = '0';
      currLhsBit = lhsPattern[i % ( lhsPatternLength )];
      // rhs pattern is shifted to align with first zero in LHS pattern
      currRhsBit = rhsPattern[AbsoluteValue( ( i + 1 - firstZeroLhsPattern ) ) % rhsPatternLength];
      // or the result to produce the merged pattern
      mergedPattern[i] = currLhsBit=='1' || currRhsBit=='1' ? '1':'0';
    }
    mergedPattern[mergedPatternLength]='\0';
    finish = time(NULL);
    //printf("Merge took %li seconds\n",finish-start);
    return mergedPattern;
  }
const char* currElemental = NULL;
const char* Pat(const long n)
{
    if(n==0)
        return INFINITE_ZEROS;
    else
    {
        const char* PatNMinusOne = Pat(n-1);
        long firstZeroLhsPattern = FindFirstZero(PatNMinusOne);
        char* elemental = Elemental(firstZeroLhsPattern);
        currElemental = elemental;
        const char * toReturn = Merge(PatNMinusOne,elemental);
        return toReturn;
    }

}
// a cache to speed things along
// initialize to -1
#define PRIME_CACHE_SIZE 200000
char primeCache[PRIME_CACHE_SIZE];

 // a basic implementation of isPrime
char isPrime( long n )
{
    //check cache first
    if(n<PRIME_CACHE_SIZE)
    {

      char cacheValue = primeCache[n];
      if(cacheValue>0)
        return cacheValue;
    }
    char prime = 1;
    long i=3;
    for (i = 3; i <= lround(sqrt((double)n)); i += 2 )
      if ( n % i == 0 ) {
        prime = 0;
        break;
      }
    if ( ( n % 2 != 0 && prime && n > 2 ) || n == 2 ) {
      primeCache[n]=1;
      return true;
    }
    else {
      primeCache[n]=0;
      return false;
    }
}

// return curr prime note I consider P(1) = 3
long P( int n )
{
     long i=0;
     long currPrime = 0;
     for ( i = 3; i < INT_MAX ; i += 2 ) {
         if ( isPrime( i ) )
            currPrime++;
         if ( currPrime == n ) return i;
     }
     return 1;
}


// returns the number of zeros in Pat(f(n),n)
long zeros( int n )
{
       if ( n == 1 ) return 2;
       return zeros( n - 1 ) * ( P( n ) - 1 );
}
long twins( int n )
{
     if ( n == 1 ) return 1;
     return twins( n - 1 ) * ( P( n ) - 2 );
}
long triplets( int n )
{
     if ( n == 2 ) return 2;
     return triplets( n - 1 ) * ( P( n ) - 3 );
}
long quadruplets( int n )
{
     if ( n == 2 ) return 1;
     return quadruplets( n - 1 ) * ( P( n ) - 4 );
}

// returns the number of ones in Pat(f(n),n)
long ones( int n )
{
     if ( n == 1 ) return 1;
     return ( ones( n - 1 ) * P( n ) ) + zeros( n - 1 );
}

// #P[P(n)->P(n)^2] = [zeros(n-1) * (P(n)^2 - P(n))] / P(n-1)#
  // #P_Corrected[P(n)->P(n)^2] = {[zeros(n-1) * (P(n)^2 - P(n))] / P(n-1)#} + 2 ^ (#P[P(n)->P(n)^2] / 100)
  // this is the number of primes between P(n) and P(n)^2
long numberPrimesBetweenPnAndPnSquared(int n, char addCorrection)
{
if(n<2)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
int primorialIndex = n-1;
double result = ((double)(zeros(n-1) * (PnSquared-Pn)));
// divide through by the primorial(n-1)
while(primorialIndex>0)
{
  result /=P(primorialIndex);
  primorialIndex--;
}
// my P(n) function omits 2 so we'll add it back into the primorial here
result /= 2;

if(addCorrection)
{
  double nDouble = (double)n;
  double correction = log(nDouble)*0.145348;
  nDouble *= correction;
  double nDoubleSquared = nDouble*nDouble;
  long nDoubleSquaredLong = lround(nDoubleSquared);
  result -= nDoubleSquaredLong;
}
return lround(result);
}

long numPrimePnMinusOneSquaredToPnSquared(int n)
{
if(n<2)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
long PnMinusOne = P(n-1);
long PnMinusOneSquared = PnMinusOne * PnMinusOne;
long numPrime = 0;
for (long i=PnMinusOneSquared+2; i<PnSquared; i+=2)
  if(isPrime(i))
    numPrime++;
return numPrime;
}

long numberPrimesBetweenPnMinusOneSquaredAndPnSquared(int n, char correct)
{
if(n<2)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
long PnMinusOne = P(n-1);
long PnMinusOneSquared = PnMinusOne * PnMinusOne;
int primorialIndex = n;
double result = ((double)(zeros(n) * (PnSquared-PnMinusOneSquared)));
// divide through by the primorial(n-1)
while(primorialIndex>0)
{
  result /=P(primorialIndex);
  primorialIndex--;
}
// my P(n) function omits 2 so we'll add it back into the primorial here
result /= 2;

if(correct)
{
  double nDouble = (double)n;
  double correction = log(nDouble)*0.145348;
  nDouble *= correction;
  double nDoubleSquared = nDouble*nDouble;
  long nDoubleSquaredLong = lround(nDoubleSquared);
  result -= nDoubleSquaredLong;
}

return lround(result);
}

long numberTwinsBetweenPnAndPnSquared(int n, char addCorrection)
{
if(n<2)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
int primorialIndex = n-1;
double result = ((double)(twins(n-1) * (PnSquared-Pn)));
// divide through by the primorial(n-1)
while(primorialIndex>0)
{
  result /=P(primorialIndex);
  primorialIndex--;
}
// my P(n) function omits 2 so we'll add it back into the primorial here
result /= 2;

if(addCorrection)
{
  double nDouble = (double)n;
  double correction = log(nDouble)*0.058652;
  nDouble *= correction;
  double nDoubleSquared = nDouble*nDouble;
  long nDoubleSquaredLong = lround(nDoubleSquared);
  result -= nDoubleSquaredLong;
}
return lround(result);
}

long numberTripletsBetweenPnAndPnSquared(int n, char addCorrection)
{
if(n<3)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
int primorialIndex = n-1;
double result = ((double)(triplets(n-1) * (PnSquared-Pn)));
// divide through by the primorial(n-1)
while(primorialIndex>0)
{
  result /=P(primorialIndex);
  primorialIndex--;
}
// my P(n) function omits 2 so we'll add it back into the primorial here
result /= 2;

if(addCorrection)
{
  double nDouble = (double)n;
  double correction = log(nDouble)*0.02881;
  nDouble *= correction;
  double nDoubleSquared = nDouble*nDouble;
  long nDoubleSquaredLong = lround(nDoubleSquared);
  result -= nDoubleSquaredLong;
}
return lround(result);
}

long numberQuadrupletsBetweenPnAndPnSquared(int n, char addCorrection)
{
if(n<3)
  return -1; //error
long Pn = P(n);
long PnSquared = Pn*Pn;
int primorialIndex = n-1;
double result = ((double)(quadruplets(n-1) * (PnSquared-Pn)));
// divide through by the primorial(n-1)
while(primorialIndex>0)
{
  result /=P(primorialIndex);
  primorialIndex--;
}
// my P(n) function omits 2 so we'll add it back into the primorial here
result /= 2;

if(addCorrection)
{
  double nDouble = (double)n;
  double correction = log(nDouble)*0.00718;
  nDouble *= correction;
  double nDoubleSquared = nDouble*nDouble;
  long nDoubleSquaredLong = lround(nDoubleSquared);
  result -= nDoubleSquaredLong;
}
return lround(result);
}

int main()
{
    time_t start=0,finish=0;
    long i=0,x=0;
    const char* patnDisplay = NULL;
    const char* patn = NULL;
    long calculatedOnes=0,calculatedZeros=0;
    start = time(NULL);

    // initialize primeCache to all -1;
    for(i=0;i<PRIME_CACHE_SIZE;i++)
      primeCache[i]=-1;

    // demonstrate Pat(n)
    for(i=0; i<8; i++)
    {
        long  n = i+1;
        long lenPatn=0;
        long numOnes=0,numZeros=0,singletonCandidates=0,twinCandidates=0,totalCandidates=0;
        int x=0;
        numOnes=0; numZeros=0;singletonCandidates=0;twinCandidates=0,totalCandidates=0;
        patn = Pat(n);
        lenPatn = strlen(patn);
        for(x=0;x<lenPatn;x++)
        {
            if(patn[x]=='1')
                numOnes++;
            else
                numZeros++;

            //singleton candidates
            if ( ( patn[x]=='0' ) && ( patn[x - 1]=='1' ) && ( patn[(x + 1) % lenPatn]=='1' ) ) singletonCandidates++;
            //twin candidates
            if ( ( patn[x]=='0' ) && ( patn[x - 1]=='0' ) ) twinCandidates++;

        }
        patnDisplay = patn;
        if(i>4)
          patnDisplay = "TOO BIG!";
        printf( "---------------------------------------------------------\n" );
        int elementalShiftedLength = strlen(currElemental)+1;
        char* elementalShifted = (char*)malloc(elementalShiftedLength);
        int f0 =(int)(P((int)n)-1)/2;
        for(int s=0;s<elementalShiftedLength-1;s++)
        {
          if(s+1==f0)
            elementalShifted[s]='1';
          else
            elementalShifted[s]='0';
        }
        elementalShifted[elementalShiftedLength-1]='\0'; // null terminate
        printf("Elemental(2f0+1,f0=%i) = %s -- shifted to align with f0 --> %s\n",f0,currElemental,elementalShifted);
        printf("Pat(2n+1,n=%i) = %s\n",n,patnDisplay);
        printf("Size = %li\n",lenPatn);
        printf("Number of 1's = %li; %li/%li = %Lf \n",numOnes,numOnes,lenPatn,(double)numOnes/(double)lenPatn);
        printf("Number of 0's = %li; %li/%li = %Lf \n",numZeros,numZeros,lenPatn,(double)numZeros/(double)lenPatn);
        totalCandidates = singletonCandidates+twinCandidates;
        printf("Total Candidates= %li\n",totalCandidates);
        printf("Singleton Candidates = %li; %li/%li = %Lf \n",singletonCandidates,singletonCandidates, singletonCandidates+twinCandidates,(double)singletonCandidates/(double)totalCandidates);
        printf("Twin Candidates = %li; %li/%li = %Lf \n",twinCandidates,twinCandidates, singletonCandidates+twinCandidates,(double)twinCandidates/(double)totalCandidates);
        printf( "---------------------------------------------------------\n" );
        free((void*)patn);
        free((void*)currElemental);
        free((void*)elementalShifted);

    }// for

    // examine the ones(n) and zeros(n) functions
    printf("%s","IT IS ALSO POSSIBLE TO ALGORITHMICALLY DETERMINE THE # OF ONES AND ZEROS AS OUTLINED ON MY WEBSITE.\n" );
    for ( x = 2; x < 9 ; x++ ) {
    calculatedOnes = ones( x );
    calculatedZeros = zeros( x );
    printf( "n=%li; Ones->%li Zeros->%li zeros/(ones+zeros) = %Lf\n",x,calculatedOnes,calculatedZeros,((double)calculatedZeros )/(double)( calculatedOnes + calculatedZeros ) );

    }
    printf( "%s", "---------------------------------------------------------\n" );
    // compare with actual values of primes:oddNumbers between P(n) and P(n)^2
    // System.out.println("Long.MAX_VALUE="+Long.MAX_VALUE);
    printf( "COMPARE THE ABOVE VALUES TO THE NUMBER OF PRIMES/ODD BETWEEN P(n)->P(n)^2\n" );
    printf( "NOTE: P(n)->P(n)^2 compares to n-1 above because (Pat(n-1)) is in effect between P(n)->P(n)^2)\n" );
    // an array to store n->numPrime for comparison later with #P[P(n)->P(n)^2]
    #define ITERATIONS 80
    long nToNumPrime[ITERATIONS+1];
    long nToNumTwin[ITERATIONS+1];
    long nToNumTrip[ITERATIONS+1];
    long nToNumQuad[ITERATIONS+1];

    for ( x = 1; x < ITERATIONS; x++ ) {
      long thisPrime = P( x );
      long thisPrimeSquared = thisPrime * thisPrime;
      long numPrime, numOdd, numTwin, numTrip,numQuad;
      numPrime = 0;
      numTwin = 0;
      numOdd = 0;
      numTrip=0;
      numQuad=0;
      for ( long j = thisPrime + 2; j < thisPrimeSquared; j += 2 ) {
        if ( isPrime( j ) ) numPrime++;
        if ( isPrime( j ) && isPrime(j-2) && (j>3) ) numTwin++;
        if ( isPrime( j ) && isPrime(j-4) && isPrime(j-6) && (j>5) ) numTrip++;
        if ( isPrime( j ) && isPrime(j-2) && isPrime(j-6) && isPrime(j-8) && (j>9) ) numQuad++;
        numOdd++;
      }
      nToNumPrime[x]=numPrime;
      nToNumTwin[x]=numTwin;
      nToNumTrip[x]=numTrip;
      nToNumQuad[x]=numQuad;
      printf( "P(%li)->P(%li)^2 : %li->%li numPrime/numOdd = %li/%li = %Lf\n",x,x,thisPrime,thisPrimeSquared,numPrime,numOdd,((double)numPrime/(double)numOdd) );
    }

    // examining the #P[P(n)->P(n)^2] formula
    printf( "---------------------------------------------------------\n" );
    printf( "EXAMINING THE #P[P(n)->P(n)^2] AND #P_Corrected[P(n)->P(n)^2] FORMULAE\n" );
    // anything more than 9 overflows double
    for ( x = 2; x < 9; x++ ) {
      printf( "#P[P(%li)->P(%li)^2] = %li; #P_Corrected[P(%li)->P(%li)^2] = %li; actual number primes in this interval = %li\n",x,x,numberPrimesBetweenPnAndPnSquared(x,0),x,x,numberPrimesBetweenPnAndPnSquared(x,1),nToNumPrime[x]);
    }
    // examining the #P[P(n-1)^2->P(n)^2] formula
    printf( "---------------------------------------------------------\n" );
    printf( "EXAMINING THE #P[P(n-1)^2->P(n)^2] AND #P_Corrected[P(n-1)^2->P(n)^2] FORMULAE\n" );
    // anything more than 8 overflows double
    for ( x = 2; x < 8; x++ ) {
      printf( "#P[P(%li)^2->P(%li)^2] = %li; #P_Corrected[P(%li)^2->P(%li)^2] = %li; actual number primes in this interval = %li\n",x-1,x,numberPrimesBetweenPnMinusOneSquaredAndPnSquared(x,0),x-1,x,numberPrimesBetweenPnMinusOneSquaredAndPnSquared(x,1),numPrimePnMinusOneSquaredToPnSquared(x));
    }
    // examining the #Twins[P(n)->P(n)^2] formula
    printf( "---------------------------------------------------------\n" );
    printf( "EXAMINING THE #Twins[P(n)->P(n)^2] AND #Twins_Corrected[P(n)->P(n)^2] FORMULAE\n" );
    // anything more than 9 overflows double
    for ( x = 2; x < 9; x++ ) {
      printf( "#Twins[P(%li)->P(%li)^2] = %li; #Twins_Corrected[P(%li)->P(%li)^2] = %li; actual number twins in this interval = %li\n",x,x,numberTwinsBetweenPnAndPnSquared(x,0),x,x,numberTwinsBetweenPnAndPnSquared(x,1),nToNumTwin[x]);
    }

    // examining the #Triplets[P(n)->P(n)^2] formula
    printf( "---------------------------------------------------------\n" );
    printf( "EXAMINING THE #Triplets[P(n)->P(n)^2] AND #Triplets_Corrected[P(n)->P(n)^2] FORMULAE\n" );
    // anything more than 9 overflows double
    for ( x = 3; x < 9; x++ ) {
      printf( "#Triplets[P(%li)->P(%li)^2] = %li; #Triplets_Corrected[P(%li)->P(%li)^2] = %li; actual number triplets in this interval = %li\n",x,x,numberTripletsBetweenPnAndPnSquared(x,0),x,x,numberTripletsBetweenPnAndPnSquared(x,1),nToNumTrip[x]);
    }
    // examining the #Quadruplets[P(n)->P(n)^2] formula
    printf( "---------------------------------------------------------\n" );
    printf( "EXAMINING THE #Quadruplets[P(n)->P(n)^2] AND #Quadruplets_Corrected[P(n)->P(n)^2] FORMULAE\n" );
    // anything more than 9 overflows double
    for ( x = 3; x < 9; x++ ) {
      printf( "#Quadruplets[P(%li)->P(%li)^2] = %li; #Quadruplets_Corrected[P(%li)->P(%li)^2] = %li; actual number quadruplets in this interval = %li\n",x,x,numberQuadrupletsBetweenPnAndPnSquared(x,0),x,x,numberQuadrupletsBetweenPnAndPnSquared(x,1),nToNumQuad[x]);
    }
    finish = time(NULL);
    printf( "---------------------------------------------------------\n" );
    printf("Execution time = %li seconds.",finish-start);
    return 0;
}