#include "newL2.h"

int mynode, totalnodes;

/* 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. 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 int* distToCycle, const int* cycleSize,
long long* allToCycleSum, long long* allCycleLengthSum) 
{
    long long 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.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.*/
void setArrays(int * cycleSize, bool* visit, int* distToCycle, bool* image)
{
    int i = 0;
    for(i = 0; i < n; i++)
    {
        visit[i] = false;
        cycleSize[i] = 0;
        distToCycle[i] = 0;
        image[i] = false;
    }
}

void zeroList(int * listArray) 
{
    int i = 0;
    for(i = 0; i < n; i++)
    listArray[i] = 0;
}

void writeTotalResults(
    int* maxTAll,
    int* maxCAll,
    int* terminalAll,
    int* imageAll,
    int* allComponents,
    int* allCyclicNodes,
    long long* allToCycleSum,
    long long* allCycleLengthSum,
    int* oneCycles,
    int* twoCycles,
    int* threeCycles,
    int* fourCycles,
    int* fiveCycles,
    int* sevenCycles,
    int* tenCycles,
    int* twentyfiveCycles,
    int* marker,
    int* orderMarker)
    /*int ccMinTotal,
    double expectedComponents,
    double expectedCyclicNodes,
    double expectedTerminalNodes,
    double expectedImageNodes,
    double expectedTailLength,
    double expectedCycleLength,
    double expectedFixedPoints,
    double expectedMaxCycle,
    double expectedMaxTail)*/
{
    double trials=n-3;//euler((n-1)/M_ARY);
    char fileStr[20];
    sprintf(fileStr, "%d_%d_%d.dat", n, M_ARY, mynode);
    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 */
    double ccMin, numGraphs;
    int i = 2;
    for(i = 2; i < n-1; i++) 
    {
        fprintf(out, "%d %lld %lld %d %d %d %d\n", allComponents[i], allCycleLengthSum[i],
        allToCycleSum[i], terminalAll[i], maxCAll[i], maxTAll[i], allCyclicNodes[i]);
    }
    fclose(out);

    char fileStr2[20];
    sprintf(fileStr2, "distrib_%d_%d.dat", n, M_ARY);
    FILE * dis = fopen(fileStr2, "w");
    for(i=2; i < n-1; i++)
    {
       fprintf(dis, "%d %d %d %d %d %d %d %d %d\n", marker[i], oneCycles[i], twoCycles[i], threeCycles[i], fourCycles[i], fiveCycles[i], sevenCycles[i], tenCycles[i], twentyfiveCycles[i]); 
    }
    fclose(dis);
	
    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;
    double cFixedPoints = 0;
    for (i = 2; i < n-1; i++)
    {   
        cComponents += ((double)allComponents[i]) / trials;
        cComponentsSquared += (double)allComponents[i] * (double)allComponents[i] / trials;
        cCyclicNodes += (double)allCyclicNodes[i] / trials;
        cCyclicNodesSquared += (double)allCyclicNodes[i] * (double)allCyclicNodes[i] / trials;
        long double cycle = (long double)allCycleLengthSum[i] / (double)(n-1);
        cWeightedCycle += (long double)(cycle) / trials;
        cWeightedCycleSquared += (long double)(cycle*cycle) / trials;
        double tail = (double)allToCycleSum[i] / (double)(n-1);
        cWeightedTail += tail / trials;
        cWeightedTailSquared += (double)(tail*tail) / trials;
        cImageNodes += ((double)(n-1) - terminalAll[i]) / trials;
        cMaxCycle += (double)maxCAll[i] / trials;
        cMaxCycleSquared += (double)maxCAll[i]*(double)maxCAll[i] / trials;
        cMaxTail += (double)maxTAll[i] / trials;
        cMaxTailSquared += (double)maxTAll[i]*(double)maxTAll[i] / trials;
        cFixedPoints += (double) oneCycles[i] / trials;
    }
    

    double cComponentsTot = 0;
    double cComponentsSquaredTot = 0;
    double cCyclicNodesTot = 0;
    double cCyclicNodesSquaredTot = 0;
    double cImageNodesTot = 0;
    double cMaxCycleTot = 0;
    double cMaxCycleSquaredTot = 0;
    double cMaxTailTot = 0;
    double cMaxTailSquaredTot = 0;
    double cWeightedCycleTot = 0;
    double cWeightedCycleSquaredTot = 0;
    double cWeightedTailTot = 0;
    double cWeightedTailSquaredTot = 0;
    
    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);

    if (mynode == 0)
    {
        char res[20];
        char resAvg[25];
        char resExp[25];
        char resAll[25];
        sprintf(res, "results_%d.dat", n);
        sprintf(resAvg, "avg_results_%d.csv", n);
        sprintf(resExp, "exp_results_%d.csv", n);
        sprintf(resAll, "all_results_%d.csv", n);
        FILE* r = fopen(res, "w");
        FILE* r1 = fopen(resAvg, "w");
        FILE* r2 = fopen(resExp, "w");
        FILE* r3 = fopen(resAll, "w");
         
        fprintf(r1, "Prime\nOrder, numGraphs, avgCC, avgCyclicN, avgTermN, avgImgN, avgFP, avgTotCycLgth, avgDistToCyc, avgMCyc, avgMTail\n");
        fprintf(r2, "Prime\nOrder, numGraphs, expCC, expCyclicN, expTermN, expImgN, expFP, expTotCycLgth, expDistToCyc\n");
        fprintf(r3, "Prime\nOrder, numGraphs, avgCC, expCC, avgCyclicN, expCyclicN, avgTermN, expTermN, avgImgN, expImgN, avgFP, expFP, avgTotCycLgth, expTotCycLgth, avgDistToCyc, expDistToCyc, avgMCyc, avgMTail\n");
        
        fprintf(r1, "%d\n", n);
        fprintf(r2, "%d\n", n);
        fprintf(r3, "%d\n", n);
        double ComponentsVariance = cComponentsSquaredTot - cComponentsTot*cComponentsTot;
        double CyclicNodesVariance = cCyclicNodesSquaredTot - cCyclicNodesTot*cCyclicNodesTot;
        double WeightedCycleVariance = cWeightedCycleSquaredTot - cWeightedCycleTot*cWeightedCycleTot;
        double WeightedTailVariance = cWeightedTailSquaredTot - cWeightedTailTot*cWeightedTailTot;
        double MaxCycleVariance = cMaxCycleSquaredTot - cMaxCycleTot*cMaxCycleTot;
        double MaxTailVariance = cMaxTailSquaredTot - cMaxTailTot*cMaxTailTot;

        double avgCC, expCC;
        double avgCyclicN, expCyclicN;
        double avgTermN, expTermN;
        double avgImgN, expImgN;
        double avgFP, expFP;
        double avgTotCycLgth, expTotCycLgth;
        double avgDistToCyc, expDistToCyc;
        double avgMCyc, expMCyc;
        double avgMTail, expMTail;

        for (i=2; i<=(n-1)/2;++i)
        {
	   if (orderMarker[i] != 0)
           {
              ccMin = (n-1)/((double) i);
              fprintf(r, "Order: %d\n", i);
              fprintf(r, "Number of graphs: %d\n", orderMarker[i]);
              fprintf(r, "ccMin: %f\n", ccMin);
              numGraphs = ((double) orderMarker[i]);
              
              fprintf(r1, "%d,%d,", i, orderMarker[i]);
              fprintf(r2, "%d,%d,", i, orderMarker[i]);
              fprintf(r3, "%d,%d,", i, orderMarker[i]);

              avgCC = allComponents[i]/numGraphs;
              expCC = ccMin*(0.5)*log(i);
              fprintf(r, "Average connected components: %f\n", avgCC);
              fprintf(r, "Expected connected components: %f\n", expCC);
              fprintf(r1, "%f,", avgCC);
              fprintf(r2, "%f,", expCC);
              fprintf(r3, "%f,%f,", avgCC, expCC); 
              
              avgCyclicN = allCyclicNodes[i]/numGraphs;
              expCyclicN = ccMin*(sqrt((M_PI*i)/2));
              fprintf(r, "Average cyclic nodes: %f\n", avgCyclicN);
              fprintf(r, "Expected cyclic nodes: %f\n", expCyclicN); 
              fprintf(r1, "%f,", avgCyclicN);
              fprintf(r2, "%f,", expCyclicN);
              fprintf(r3, "%f,%f,", avgCyclicN, expCyclicN); 

              avgTermN = terminalAll[i]/numGraphs;
              expTermN = ccMin*((1.0/M_E)*i);
              fprintf(r, "Average terminal nodes: %f\n", avgTermN);
              fprintf(r, "Expected terminal nodes: %f\n", expTermN); 
              fprintf(r1, "%f,", avgTermN);
              fprintf(r2, "%f,", expTermN);
              fprintf(r3, "%f,%f,", avgTermN, expTermN); 

              avgImgN = imageAll[i]/numGraphs;
              expImgN = ccMin*((1-(1.0/M_E))*i);
              fprintf(r, "Average image nodes: %f\n", avgImgN);
              fprintf(r, "Expected image nodes: %f\n", expImgN); 
              fprintf(r1, "%f,", avgImgN);
              fprintf(r2, "%f,", expImgN);
              fprintf(r3, "%f,%f,", avgImgN, expImgN); 
              
              avgFP = oneCycles[i]/numGraphs;
              expFP = ccMin*1;
              fprintf(r, "Average fixed points: %f\n", oneCycles[i]/numGraphs);
              fprintf(r, "Expected fixed points: %f\n", expFP); 
              fprintf(r1, "%f,", avgFP);
              fprintf(r2, "%f,", expFP);
              fprintf(r3, "%f,%f,", avgFP, expFP); 

              avgTotCycLgth = allCycleLengthSum[i]/numGraphs;
              expTotCycLgth = ccMin*i*(sqrt((M_PI*i)/8));
              fprintf(r, "Average total cycle length: %f\n", avgTotCycLgth);
              fprintf(r, "Expected total cycle length: %f\n", expTotCycLgth); 
              fprintf(r1, "%f,", avgTotCycLgth);
              fprintf(r2, "%f,", expTotCycLgth);
              fprintf(r3, "%f,%f,", avgTotCycLgth, expTotCycLgth); 

              avgDistToCyc = allToCycleSum[i]/numGraphs;
              expDistToCyc = ccMin*i*(sqrt((M_PI*i)/8));
              fprintf(r, "Average distance to cycle: %f\n", avgDistToCyc);
              fprintf(r, "Expected distance to cycle: %f\n", expDistToCyc); 
              fprintf(r1, "%f,", avgDistToCyc);
              fprintf(r2, "%f\n", expDistToCyc);
              fprintf(r3, "%f,%f,", avgDistToCyc, expDistToCyc); 

              avgMCyc = maxCAll[i]/numGraphs;
              avgMTail = maxTAll[i]/numGraphs;
              fprintf(r, "Average max cycle: %f\n", avgMCyc);
              fprintf(r, "Average max tail: %f\n\n", avgMTail);
              fprintf(r1, "%f,%f\n", avgMCyc, avgMTail);
              fprintf(r3, "%f,%f\n", avgMCyc, avgMTail); 
           }
        }

        i = n-1;
	if (orderMarker[n-1] != 0)
        {
              ccMin = (n-1)/((double) i);
              fprintf(r, "Order: %d\n", i);
              fprintf(r, "Number of graphs: %d\n", orderMarker[i]);
              fprintf(r, "ccMin: %f\n", ccMin);
              numGraphs = ((double) orderMarker[i]);
              
              fprintf(r1, "%d,%d,", i, orderMarker[i]);
              fprintf(r2, "%d,%d,", i, orderMarker[i]);
              fprintf(r3, "%d,%d,", i, orderMarker[i]);

              avgCC = allComponents[i]/numGraphs;
              expCC = ccMin*(0.5)*log(i);
              fprintf(r, "Average connected components: %f\n", avgCC);
              fprintf(r, "Expected connected components: %f\n", expCC);
              fprintf(r1, "%f,", avgCC);
              fprintf(r2, "%f,", expCC);
              fprintf(r3, "%f,%f,", avgCC, expCC); 
              
              avgCyclicN = allCyclicNodes[i]/numGraphs;
              expCyclicN = ccMin*(sqrt((M_PI*i)/2));
              fprintf(r, "Average cyclic nodes: %f\n", avgCyclicN);
              fprintf(r, "Expected cyclic nodes: %f\n", expCyclicN); 
              fprintf(r1, "%f,", avgCyclicN);
              fprintf(r2, "%f,", expCyclicN);
              fprintf(r3, "%f,%f,", avgCyclicN, expCyclicN); 

              avgTermN = terminalAll[i]/numGraphs;
              expTermN = ccMin*((1.0/M_E)*i);
              fprintf(r, "Average terminal nodes: %f\n", avgTermN);
              fprintf(r, "Expected terminal nodes: %f\n", expTermN); 
              fprintf(r1, "%f,", avgTermN);
              fprintf(r2, "%f,", expTermN);
              fprintf(r3, "%f,%f,", avgTermN, expTermN); 

              avgImgN = imageAll[i]/numGraphs;
              expImgN = ccMin*((1-(1.0/M_E))*i);
              fprintf(r, "Average image nodes: %f\n", avgImgN);
              fprintf(r, "Expected image nodes: %f\n", expImgN); 
              fprintf(r1, "%f,", avgImgN);
              fprintf(r2, "%f,", expImgN);
              fprintf(r3, "%f,%f,", avgImgN, expImgN); 
              
              avgFP = oneCycles[i]/numGraphs;
              expFP = ccMin*1;
              fprintf(r, "Average fixed points: %f\n", avgFP);
              fprintf(r, "Expected fixed points: %f\n", expFP); 
              fprintf(r1, "%f,", avgFP);
              fprintf(r2, "%f,", expFP);
              fprintf(r3, "%f,%f,", avgFP, expFP); 

              avgTotCycLgth = allCycleLengthSum[i]/numGraphs;
              expTotCycLgth = ccMin*i*(sqrt((M_PI*i)/8));
              fprintf(r, "Average total cycle length: %f\n", avgTotCycLgth);
              fprintf(r, "Expected total cycle length: %f\n", expTotCycLgth); 
              fprintf(r1, "%f,", avgTotCycLgth);
              fprintf(r2, "%f,", expTotCycLgth);
              fprintf(r3, "%f,%f,", avgTotCycLgth, expTotCycLgth); 

              avgDistToCyc = allToCycleSum[i]/numGraphs;
              expDistToCyc = ccMin*i*(sqrt((M_PI*i)/8));
              fprintf(r, "Average distance to cycle: %f\n", avgDistToCyc);
              fprintf(r, "Expected distance to cycle: %f\n", expDistToCyc); 
              fprintf(r1, "%f,", avgDistToCyc);
              fprintf(r2, "%f\n", expDistToCyc);
              fprintf(r3, "%f,%f,", avgDistToCyc, expDistToCyc); 

              avgMCyc = maxCAll[i]/numGraphs;
              avgMTail = maxTAll[i]/numGraphs;
              fprintf(r, "Average max cycle: %f\n", avgMCyc);
              fprintf(r, "Average max tail: %f\n\n", avgMTail);
              fprintf(r1, "%f,%f\n", avgMCyc, avgMTail);
              fprintf(r3, "%f,%f\n", avgMCyc, avgMTail);       
        } 

        fclose(r);
        fclose(r1);
        fclose(r2);
        fclose(r3);
          
        /*
        //fprintf(r, "components average: %lf \n", cComponentsTot);
        fprintf(r, "components total: %lf \n", cComponentsTot*trials);
        fprintf(r, "expected components: %lf \n\n", expectedComponents);
        //fprintf(r, "components variance: %lf \n", ComponentsVariance);
        
        //fprintf(r, "cyclic nodes average: %lf \n", cCyclicNodesTot);
        fprintf(r, "cyclic nodes total: %lf \n", cCyclicNodesTot*trials);
        fprintf(r, "expected cyclic nodes: %lf \n\n", expectedCyclicNodes);
        //fprintf(r, "cyclic nodes variance: %lf\n", CyclicNodesVariance);

        //fprintf(r, "terminal nodes average: %lf\n",(n-1)-cImageNodesTot);
        fprintf(r, "terminal nodes total: %lf\n", ((n-1)-cImageNodesTot)*trials);
        fprintf(r, "expected terminal nodes: %lf\n\n", expectedTerminalNodes);

        //fprintf(r, "image nodes average: %lf\n", cImageNodesTot);
        fprintf(r, "image nodes total: %lf\n", cImageNodesTot*trials);
        fprintf(r, "expected image nodes: %lf\n\n", expectedImageNodes);
         
        fprintf(r, "total cycle length (as seen from a node): %lf\n", cWeightedCycle*(trials*(n-1)));
        fprintf(r, "expected total cycle length: %lf\n\n", expectedCycleLength);
       
        fprintf(r, "total distance to a cycle: %lf\n", cWeightedTail*(trials*(n-1)));
        fprintf(r, "expected total distance to a cycle: %lf\n\n", expectedTailLength);
        
        fprintf(r, "avg cycle (as seen from a node): %lf\n", cWeightedCycleTot);
        fprintf(r, "expected avg cycle: %lf\n\n", expectedCycleLength/ccMinTotal);
        //fprintf(r, "avg cycle variance: %lf\n", WeightedCycleVariance);

        fprintf(r, "avg tail (distance from cycle as seen from a node): %lf\n", cWeightedTailTot);
        fprintf(r, "expected avg tail: %lf\n\n", expectedTailLength/ccMinTotal);
        //fprintf(r, "avg tail variance: %lf\n", WeightedTailVariance);
        
        fprintf(r, "fixed points total: %lf \n",cFixedPoints*trials);
        fprintf(r, "expected fixed points: %lf \n\n",expectedFixedPoints);
        
        fprintf(r, "max cycle average: %lf\n", cMaxCycleTot);
        fprintf(r, "expected max cycle average: %lf\n\n", expectedMaxCycle);
        //fprintf(r, "max cycle variance: %lf\n", MaxCycleVariance);

        fprintf(r, "max tail average: %lf\n", cMaxTailTot);
        fprintf(r, "expected max tail average: %lf\n\n", expectedMaxTail);
        //fprintf(r, "max tail variance: %lf\n", MaxTailVariance);
        */
        //fclose(r);
        //fclose(r1);
        //fclose(r2);
        //fclose(r3);
    }
}
void run() 
{
    MPI_Comm_size(MPI_COMM_WORLD, &totalnodes);
    MPI_Comm_rank(MPI_COMM_WORLD, &mynode);

    FILE * s;
    if (mynode == 0)
    {
    	s = fopen(STATUS, "a");
        fprintf(s, "Allocating...\n");
        fclose(s);                                      
    }	
    bool visit[n];
    bool image[n];
    /*Maximum tail length for base [i]*/
    int maxTAll[n];
    /*Maximum cycle lenghth for base [i]*/
    int maxCAll[n];
    /*Terminal nodes for base [i]*/
    int terminalAll[n];
    int imageAll[n];
    /*Size of cycle for the component this node is a part of*/
    int cycleSize[n];
    /*Distance to cycle from node n (0 if node n is cyclic)*/
    int distToCycle[n];
    /*Number of components for base i*/
    int allComponents[n+1];
    /*Number of image nodes for base i*/
    int allCyclicNodes[n+1];
    /*Sum of all nodes. distance to cycle for base i*/
    long long allToCycleSum[n+1];
    /*Sum of each node.s cycle length for base i*/
    long long allCycleLengthSum[n+1];
    /*Record the number of cycles of arbitrary lengths.*/
    int oneCycles[n];
    int twoCycles[n];
    int threeCycles[n];
    int fourCycles[n];
    int fiveCycles[n];
    int sevenCycles[n];
    int tenCycles[n];
    int twentyfiveCycles[n];
    int marker[n];
    int orderMarker[n];
    /*Initialize variables that will store max cycle length (mC) and max tail length (mT).*/
    int mC, mT;
    int next, loc, baseTail, cycleLength, terminal;
    int root, exp, base;
    int listArray[n];
    int listSize = 0;
    int rootOrder, ccMin, ccMinTotal;
    double expectedComponents=0, expectedCyclicNodes=0, expectedTerminalNodes=0, expectedImageNodes=0, expectedTailLength=0, expectedCycleLength=0;
    double expectedFixedPoints=0, expectedMaxCycle=0, expectedMaxTail=0;

    if (mynode == 0)
    {
        s = fopen(STATUS, "a");
        fprintf(s, "zeroing...\n");
        fclose(s);
    }
    /*initialize arrays to 0 */
    int i = 0;
    for(i = 0; i < n; i++)
    {
        maxTAll[i] = 0;
        maxCAll[i] = 0;
        terminalAll[i] = 0;
        imageAll[i] = 0;
        allComponents[i] = 0;
        allCyclicNodes[i] = 0;
        allToCycleSum[i] = 0;
        allCycleLengthSum[i] = 0;
        oneCycles[i] = 0;
        twoCycles[i] = 0;
        threeCycles[i] = 0;
        fourCycles[i] = 0;
        fiveCycles[i] = 0;
        sevenCycles[i]=0;
        tenCycles[i] = 0;
        twentyfiveCycles[i] = 0;
        marker[i] = 0;
        orderMarker[i] = 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;
    
    if (mynode == 0)
    {
        s = fopen(STATUS, "w");
        fprintf(s, "Finding a PR...\n");
        fclose(s);
    }

    //find the smallest primitive root
    for(root = 1; !isPrimRoot(root); root++);
    if (mynode == 0)
    {
        s = fopen(STATUS, "a");
        fprintf(s, "Prim root is %d...\n", root);
        fclose(s);
    }

    //int count = -1;

    for(exp = 1; exp < n-1; exp++) 
    {
        if (exp == ((n-1)/2)) continue;

	base = m_exp(root,exp);
        rootOrder = (n-1)/gcd(exp,n-1);
        orderMarker[rootOrder] += 1;

       /* ccMin = (n-1)/rootOrder;
        ccMinTotal += ccMin;

        expectedComponents += ccMin*(0.5*log(rootOrder));
        expectedCyclicNodes += ccMin*(sqrt((M_PI*rootOrder)/2));
        expectedTerminalNodes += ccMin*((1/M_E)*rootOrder);
        expectedImageNodes += ccMin*((1-(1/M_E))*rootOrder);
        expectedTailLength += ccMin*(sqrt((M_PI*rootOrder)/8));
        expectedCycleLength += ccMin*(sqrt((M_PI*rootOrder)/8));
        expectedFixedPoints += ccMin*(1);
        expectedMaxCycle += ccMin*(0.78248*sqrt(rootOrder));
        expectedMaxTail += ccMin*(sqrt(2*M_PI)*log(2)*sqrt(rootOrder));*/
        

        if(exp % 100 == 0 && mynode == 0) 
        {
            s = fopen(STATUS, "a");
            fprintf(s, "Exp is %d\n", exp);
            fclose(s);
        }
        //printf("g=%d\n",base);
        /*discard all but the bases which will make the target M-ARY graphs*/
        /*if(gcd(exp,n-1) != M_ARY) continue;
        count++;
        if(count % totalnodes != mynode) continue;
        base = m_exp(root,exp);*/
        /*Mark each base that will produce target graph.*/
       // marker[base]=1;
        /*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*/
        for(i = 1; i < n; i++) 
        {
            if(visit[i])
                continue;
            next = i;
            listArray[0] = next;
            listSize = 1;
            tt = MPI_Wtime();
            while(!visit[next])
            {
                visit[next] = true;
                next = (next * m_exp(base,next)) % n;
                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;
                            //printf("All. Node %d, Cycle Length %d\n",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;
                        //printf("Ext. Node %d, Cycle Length %d\n",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[rootOrder]++;
                if(cycleLength==2) twoCycles[rootOrder]++;
                if(cycleLength==3) threeCycles[rootOrder]++;
                if(cycleLength==3) fourCycles[rootOrder]++;
                if(cycleLength==5) fiveCycles[rootOrder]++;
                if(cycleLength==7) sevenCycles[rootOrder]++;
                if(cycleLength==10) tenCycles[rootOrder]++;
                if(cycleLength==20) twentyfiveCycles[rootOrder]++;
                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;
                    //printf("Node %d, Cycle Length %d\n",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[rootOrder]++;
                allCyclicNodes[rootOrder] += cycleLength;
                cycleTime += MPI_Wtime() - tt;
            }
        }
        tt = MPI_Wtime();
        terminal=0;
        for(i = 1; i < n; i++)
            if(!image[i]) terminal++;
        maxTAll[rootOrder] += mT;
        maxCAll[rootOrder] += mC;
        terminalAll[rootOrder] += terminal;
        imageAll[rootOrder] += (n-1)-terminal;
        computeResults(distToCycle, cycleSize, &allToCycleSum[rootOrder], &allCycleLengthSum[rootOrder]);
        //printf("all cyle length sum %lld\n",allCycleLengthSum[base]);
        resultsTime += MPI_Wtime() - tt;
    }
    if(mynode == 0)
    {
        s = fopen(STATUS, "a");
        fprintf(s, "Writing Results...\n");
        fclose(s);
    }
    writeTotalResults(
    maxTAll,
    maxCAll,
    terminalAll,
    imageAll,
    allComponents,
    allCyclicNodes,
    allToCycleSum,
    allCycleLengthSum,
    oneCycles,
    twoCycles,
    threeCycles,
    fourCycles, 
    fiveCycles,
    sevenCycles,
    tenCycles,
    twentyfiveCycles,
    marker,
    orderMarker);
    /*ccMinTotal,
    expectedComponents,
    expectedCyclicNodes,
    expectedTerminalNodes,
    expectedImageNodes,
    expectedTailLength,
    expectedCycleLength,
    expectedFixedPoints,
    expectedMaxCycle,
    expectedMaxTail);*/
    if(mynode == 0)
    {
        s = fopen(STATUS, "a");
        fprintf(s, "%lf minutes...\n Exiting...\n%lf %lf %lf %lf %lf\n", (MPI_Wtime() - t)/60, expTime, tailTime, intoCycleTime, cycleTime, resultsTime);
        fclose(s);
    }
    printf("%lf\n", expTime);
}