#include "newH2.h"
#include <math.h>
//int mynode, totalnodes; //only need these if using parallel computing

/* Algorithm for modular expoentiation.*/
long long ml_exp(long long b, int e, long long m)
{
   long long result = 1;
   while (e > 0)
    {
        /* Multiply in this bits $-1òù contribution while using
modulus to keep result small*/
        /* Note '&' is the bitwise 'and' operator*/
        if ((e & 1) == 1) result = (result * b) % m;
        e >>= 1;
        b = (b * b) % m;
   }
   return result;
}

/* Another method for modular expoentiation called on integers.*/
int m_expn(int b, int r, int num)
{
   return (int) ml_exp((long long) b, r, (long long) num);
}

/* Yet another simplified method for expoentiation that may be used
when the modulus (n) has already been defined. */
/* The inputs are the base b and the exponent r. */
int m_exp(int b, int r)
{
   return m_expn(b, r, n);
}

/* Compact method to find total 'Distance to Cycle' and 'Cycle Size'
based on previously computed arrays containing results from the
individual nodes. */
void computeResults(const unsigned short* distToCycle, const unsigned
short* cycleSize,
int* allToCycleSum, long long* allCycleLengthSum)
{
   int sumDistToCycle = 0;
   long long sumCycleSize = 0;
   int i = 0;
   for(i = 0; i < n; i++)
   {
       sumDistToCycle += distToCycle[i];
       sumCycleSize += cycleSize[i];
   }
   *allToCycleSum = sumDistToCycle;
   *allCycleLengthSum = sumCycleSize;
}

/* Implementation of the Rabin-Miller test for primality.*/
bool MillerRabin(int num, int k, int q, int a)
{
   int n1 = num-1;
   if(m_expn(a,q, num) == 1) return true;
   int i = 0;
   for(i = 0; i < k; i++)
       if(m_expn(a,(int)pow(2,i)*q, num) == n1) return true;
   return false;
}

/*Test for primiality using an attached array of primes in the header
file bn_prime.h.*/
bool isPrime(int num)
{
   int i = 0;
   for(i = 0; i < 50;i++)
   {
       if(primes[i] > (unsigned)num) return true;
       if((num % primes[i] == 0) && (num!=primes[i])) return false;
   }
   int k = 0;
   int q = num-1;
   while(q % 2 == 0)
   {
       k++;
       q >>= 1;
   }
   srand(time(0));
   int a;
   for(i = 0; i < 10; i++)
   {
       a = (rand() % (num-2)) + 1;
       if(!MillerRabin(num, k, q, a)) return false;
   }
   return true;
}

/* Method to test whether or not an element is a Primitive Root of an
inputted modulus. */
bool isPrimRoot(int base)
{
   if(!isPrime(n)) return false;
   int n_1 = n-1;
   if ((unsigned)n_1 > (primes[NUMPRIMES-1]*primes[NUMPRIMES-1]))
       printf("Error in Primitive Root Testing, n could have prime
factor too large for testing\n");
   int n1 = n_1;
   int index = 0;
   int p;
   while(n1 > 1 && index < NUMPRIMES)
   {
       /*find the primes that divide phi(n)*/
       if((n1 % primes[index]) == 0)
       {
           p = primes[index];
           /* divide out that prime all the way so it isn $-1òùt tested again*/
           while((n1 % primes[index] == 0)) n1/=primes[index];
           /*if base^phi(n)/p is 1, not a prim root*/
           if(m_exp(base,n_1/p) == 1) return false;
           /*if(isPrime(n1)) return m_exp(base,n_1/n1) == 1;*/
           if(n1 == 50021) return (m_exp(base,n_1/50021) != 1);
       }
       index++;
   }
   return true;
}

/*Method to find the greatest common divisors of two elements.*/
int gcd(int a, int b)
{
   if(a== 0) return b;
   if(b==0) return a;
   int r = a % b;
   int d = b;
   int c;
   while (r > 0)
   {
       c = d;
       d = r;
       r = c % d;
   }
   return d;
}

/*Method to determine whether an element is relatively prime to the modulus.*/
bool isRelPrime(int base)
{
   return gcd(base, n-1) == 1;
}

/*Method to determine the number of integers relatively prime, that
is, Euler Phi Function. This is used to find
the number of m-ary graphs that will be computers thereby giving the
number of trials.*/
double euler(int numb)
{
 double value;
 int i=0;
 int result = numb;
 for(i=2;i*i <= numb;i++)
 {
     if (numb % i == 0) result -= result / i;
     while (numb % i == 0) numb /= i;
 }
 if (numb > 1) result -= result / numb;
 value=result;
 return value;
}



/*Method that initializes the arrays based on the size of the modulus.*/
/*come back here to take a look at this*/
void setArrays(unsigned short * cycleSize, bool* visit, unsigned
short* distToCycle, bool* image)
{
   int i = 0;
   for(i = 0; i < n; i++)
   {
       visit[i] = false;
       cycleSize[i] = 0;
       distToCycle[i] = 0;
       image[i] = false;
   }
}

//zeroing the list...yeah
void zeroList(int * listArray)
{
   int i = 0;
   for(i = 0; i < n; i++)
   listArray[i] = 0;
}

//this is where the data is put into files
void writeTotalResults(
   int maxTAll,
   int maxCAll,
   int terminalAll,
   int allComponents,
   int allCyclicNodes,
   int allToCycleSum,
   long long allCycleLengthSum,
   int oneCycles,
   int twoCycles,
   int threeCycles,
   int fourCycles,
   int fiveCycles,
   int sixCycles,
   int sevenCycles,
   int eightCycles,
   int nineCycles,
   int tenCycles)
{
   char fileStr[20];
   sprintf(fileStr, "%d_selfpower.dat", n);
   FILE * out = fopen(fileStr, "w");
   /* cycles base i
   sum of cycle size seen from nodes in base i
   sum of distance to cycle from nodes in base i
   terminal nodes for base i
   max cycle for base i
   max tail for base i
   cyclic nodes for base i */
   fprintf(out, "%d %ld %d %d %d %d %d\n", allComponents, allCycleLengthSum,
           allToCycleSum, terminalAll, maxCAll, maxTAll, allCyclicNodes);
   fclose(out);

   char fileStr2[20];
   sprintf(fileStr2, "distrib_%d_selfpower.dat", n);
   FILE * dis = fopen(fileStr2, "w");
   //this file contains counts on the numbers of cycles with sizes 1-10
   fprintf(dis, "%d %d %d %d %d %d %d %d %d %d\n", oneCycles,
twoCycles, threeCycles, fourCycles, fiveCycles, sixCycles,
sevenCycles, eightCycles, nineCycles, tenCycles);
   fclose(dis);

   //unused variable from previous code that are used to calculate
statistics over multiple graphs

   // double cComponents = 0;
   // double cComponentsSquared = 0;
   // double cCyclicNodes = 0;
   // double cCyclicNodesSquared = 0;
   // double cImageNodes = 0;
   /*variance = 0*/
   // double cMaxCycle = 0;
   // double cMaxCycleSquared = 0;
   // double cMaxTail = 0;
   // double cMaxTailSquared = 0;
   // double cWeightedCycle = 0;
   // double cWeightedCycleSquared = 0;
   // double cWeightedTail = 0;
   //    double cWeightedTailSquared = 0;


   // cComponents += ((double)allComponents);
   // cComponentsSquared += (double)allComponents * (double)allComponents;
   // cCyclicNodes += (double)allCyclicNodes;
   // cCyclicNodesSquared += (double)allCyclicNodes * (double)allCyclicNodes;
   // long double cycle = (long double)allCycleLengthSum / (double)(n-1);
   // cWeightedCycle += (long double)(cycle);
   // cWeightedCycleSquared += (long double)(cycle*cycle);
   // double tail = (double)allToCycleSum / (double)(n-1);
   // cWeightedTail += tail;
   // cWeightedTailSquared += (double)(tail*tail);
   // if (terminalAll > 0)
     // cImageNodes += ((double)(n-1) - terminalAll);
   // cMaxCycle += (double)maxCAll;
   // cMaxCycleSquared += (double)maxCAll*(double)maxCAll;
   // cMaxTail += (double)maxTAll;
   // cMaxTailSquared += (double)maxTAll*(double)maxTAll;

   //calls for parallel computing procedures that aren't needed for
the self-power graph

   //MPI_Reduce( &cComponents, &cComponentsTot, 1, MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cComponentsSquared, &cComponentsSquaredTot, 1,
MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   //MPI_Reduce( &cCyclicNodes, &cCyclicNodesTot, 1, MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
   //MPI_Reduce( &cCyclicNodesSquared, &cCyclicNodesSquaredTot, 1,
MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   //MPI_Reduce( &cImageNodes, &cImageNodesTot, 1, MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
   //MPI_Reduce( &cMaxCycle, &cMaxCycleTot, 1, MPI_DOUBLE, MPI_SUM,
0, MPI_COMM_WORLD);
   //    MPI_Reduce( &cMaxCycleSquared, &cMaxCycleSquaredTot, 1,
MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cMaxTail, &cMaxTailTot, 1, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD);
   // MPI_Reduce( &cMaxTailSquared, &cMaxTailSquaredTot, 1,
MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cWeightedCycle, &cWeightedCycleTot, 1, MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cWeightedCycleSquared, &cWeightedCycleSquaredTot,
1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cWeightedTail, &cWeightedTailTot, 1, MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
   // MPI_Reduce( &cWeightedTailSquared, &cWeightedTailSquaredTot, 1,
MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

   //third file that is not needed for the self-power map. Contains
some things that are already in the first file and some averages if
working over multiple maps, look at Hoffman's version for a better
example of the way it works

   // if (mynode == 0)
   // {
   //  char res[20];
   //  sprintf(res, "results_%d.dat", n);
   //  FILE * r = fopen(res, "w");
   //  fprintf(r, "%lf ", cycle);
   //  fprintf(r, "%lf\n", tail);
   //  fclose(r);
       // }
}
void run()
{
 //MPI_Comm_size(MPI_COMM_WORLD, &totalnodes);
 //MPI_Comm_rank(MPI_COMM_WORLD, &mynode);
 FILE * s;
   // if (mynode == 0)
   // {
   //  s = fopen(STATUS, "w");
   //  fprintf(s, "Allocating...\n");
   //  fclose(s);
   //}
   /*I think these indicate if a node has been visited*/
 bool visit[n];
 bool image[n];
   /*might be able to change all of these from arrays to just one variable*/
   /*Maximum tail length for base [i]*/
 int maxTAll;
   /*Maximum cycle lenghth for base [i]*/
 int maxCAll;
   /*Terminal nodes for base [i]*/
 int terminalAll;
   /*Size of cycle for the component this node is a part of*/
 unsigned short cycleSize[n];
   /*Distance to cycle from node n (0 if node n is cyclic)*/
 unsigned short distToCycle[n];
   /*Number of components for base i*/
 int allComponents;
   /*Number of image nodes for base i*/
 int allCyclicNodes;
   /*Sum of all nodes $-1òù distance to cycle for base i*/
 int allToCycleSum;
   /*Sum of each node $-1òùs cycle length for base i*/
 long long allCycleLengthSum;
   /*Record the number of cycles of arbitrary lengths.*/
 int oneCycles;
 int twoCycles;
 int threeCycles;
 int fourCycles;
 int fiveCycles;
 int sixCycles;
 int sevenCycles;
 int eightCycles;
 int nineCycles;
 int tenCycles;
   /*probably don't need marker for my code*/
   /*Initialize variables that will store max cycle length (mC) and
max tail length (mT).*/
 int mC, mT;
 int next, loc, baseTail, cycleLength, terminal;
 int exp;
 int listArray[n];
 int listSize = 0;
   //  if (mynode == 0)
   // {
   //  s = fopen(STATUS, "w");
   //  fprintf(s, "zeroing...\n");
   //  fclose(s);
   // }
   /*initialize arrays to 0 */
 maxTAll = 0;
 maxCAll = 0;
 terminalAll = 0;
 allComponents = 0;
 allCyclicNodes = 0;
 allToCycleSum = 0;
 allCycleLengthSum = 0;
 oneCycles = 0;
 twoCycles = 0;
 threeCycles = 0;
 fourCycles = 0;
 fiveCycles = 0;
 sixCycles = 0;
 sevenCycles=0;
 eightCycles=0;
 nineCycles=0;
 tenCycles = 0;
 //double t;
   // if (mynode == 0)
   //   t = MPI_Wtime();
   // double tt;
   //double expTime = 0;
   //double tailTime = 0;
   //double intoCycleTime = 0;
   //double cycleTime = 0;
   //double resultsTime = 0;
   /*Recall mC is max cycle, mT is max tail.*/
 mC = 0;
 mT = 0;
   /*0 out everything*/
 setArrays(cycleSize, visit, distToCycle, image);
   /*begin making graph, using gamma(i) = base^i mod n*/
   /*need to change this because I want one graph, x^x mod n */
 int i=0;
 //here's where the meat of the work is done
 for(i = 1; i < n; i++)
   {
     //this first section finds a node that hasn't been visited yet
     if(visit[i])
       continue;
     next = i;
     listArray[0] = next;
     listSize = 1;
     //        tt = MPI_Wtime();
     while(!visit[next])
       {
         visit[next] = true;
         /*found a place to change*/
         next = m_exp(next,next);
         image[next] = true;
         listArray[listSize] = next;
         listSize++;
       }
     //        expTime += MPI_Wtime() - tt;
     int j = 0;
     if(cycleSize[next] != 0)
       {
         if(distToCycle[next] == 0)
           {
             /*all tail into cycle*/
             //        tt = MPI_Wtime();
             cycleLength = cycleSize[listArray[listSize-1]];
             if(listSize - 1 > mT) mT = listSize - 1;
             for(j = 0; j < listSize-1; j++)
               {
                 distToCycle[listArray[j]] = listSize - 1 - j;
                 cycleSize[listArray[j]] = cycleLength;
               }
             //        intoCycleTime += MPI_Wtime() - tt;
           }
         else
           {
             /*extension of tail*/
             //        tt = MPI_Wtime();
             baseTail = distToCycle[listArray[listSize-1]];
             cycleLength = cycleSize[listArray[listSize-1]];
             if(listSize-1 + baseTail > mT) mT = listSize-1 + baseTail;
             for(j = 0; j < listSize-1; j++)
               {
                 distToCycle[listArray[j]] = baseTail + listSize - 1 - j;
                 cycleSize[listArray[j]] = cycleLength;
               }
             //        tailTime += MPI_Wtime() - tt;
           }
       }
     else
       {
         /*new cycle found*/
         // tt = MPI_Wtime();
         /*loc will be the first node in the cycle we ran in to*/
         int repeat = listArray[listSize-1];
         for(j = 0; listArray[j] != repeat; j++);
         int firstCycle = j;
         cycleLength = listSize - (j+1);
         if(cycleLength==1) oneCycles++;
         if(cycleLength==2) twoCycles++;
         if(cycleLength==3) threeCycles++;
         if(cycleLength==4) fourCycles++;
           if(cycleLength==5) fiveCycles++;
           if(cycleLength==6) sixCycles++;
           if(cycleLength==7) sevenCycles++;
           if(cycleLength==8) eightCycles++;
           if(cycleLength==9) nineCycles++;
           if(cycleLength==10) tenCycles++;
           if(cycleLength > mC) mC = cycleLength;
           if(firstCycle > mT) mT = firstCycle;
           /*mark each tail node along the way with how far it is to
             the cycle (marked as a negative number)*/
           for(j = 0; j < firstCycle; j++)
             {
               distToCycle[listArray[j]] = firstCycle - j;
               cycleSize[listArray[j]] = cycleLength;
             }
           /*mark each cycle node with how big the cycle is*/
           for(j = firstCycle; j < listSize - 1; j++)
             cycleSize[listArray[j]] = cycleLength;
           allComponents++;
           allCyclicNodes += cycleLength;
           //  cycleTime += MPI_Wtime() - tt;
         }
     }
   // tt = MPI_Wtime();
   terminal=0;
   for(i = 1; i < n; i++)
     if(!image[i]) terminal++;
   maxTAll = mT;
   maxCAll = mC;
   terminalAll = terminal;
   computeResults(distToCycle, cycleSize, &allToCycleSum, &allCycleLengthSum);
   // resultsTime += MPI_Wtime() - tt;

   // if(mynode == 0)
   // {
   //  s = fopen(STATUS, "w");
   //  fprintf(s, "Writing Results...\n");
   //  fclose(s);
   // }
   writeTotalResults(
                     maxTAll,
                     maxCAll,
                     terminalAll,
                     allComponents,
                     allCyclicNodes,
                     allToCycleSum,
                     allCycleLengthSum,
                     oneCycles,
                     twoCycles,
                     threeCycles,
                     fourCycles,
                     fiveCycles,
                     sixCycles,
                     sevenCycles,
                     eightCycles,
                     nineCycles,
                     tenCycles);
   // if(mynode == 0)
     // {
       //      s = fopen(STATUS, "w");
       //      fprintf(s, "%lf minutes...\n Exiting...\n%lf %lf %lf %lf %lf\n",
(MPI_Wtime() - t)/60, expTime, tailTime);
       //      fclose(s);
       // }
   //  printf("%lf\n", expTime);
}