#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <string.h>

#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/function.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>


using namespace std;

bool DEBUG = 0;

void TimeRoutine(string name,boost::function<int(int,int)>fun, 
            int arg1,int arg2, bool doTime);
int OldSchool(int a, int b);
int Euclid(int a, int b);
int Max(int, int);
int Min(int,int);
int Brute(int, int);
int School(int, int,const vector <int> &);
void SetAB(int, int & , int & );
bool Test(int a, int b, int d);

void args(string progname) {
   cout << progname << " [-l n -r -R n -t -d -a n -b n -h]" << endl;
   cout << endl;
   cout << "\t" << " -l n " << "\t" << " use a preset input" << endl;
   cout << "\t\t" << 0 << "\t  Book's Example " << endl;
   cout << "\t\t" << 1 << "\t  Two factorable numbers " << endl;
   cout << "\t\t" << 2 << "\t  Book's Problem " << endl;
   cout << "\t\t" << 3 << "\t  Medium Sized Numbers " << endl;
   cout << "\t\t" << 4 << "\t  Large Numbers " << endl;
   cout << "\t\t" << 5 << "\t  Very Large Numbers " << endl;
   cout << "\t" << " -r " << "\t" << " Generate random numbers" << endl;
   cout << "\t" << " -R n "<< "\t" << " Generate random numbers 2 to n" << endl;
   cout << "\t" << " -d " << "\t" << " Debug output" << endl;
   cout << "\t" << " -h " << "\t" << " This menu" << endl;
   cout << "\t" << " -t " << "\t" << " Show Timing" << endl;
   cout << "\t" << " -a n " << "\t" << " Set first number to be n" << endl;
   cout << "\t" << " -b n " << "\t" << " Set second number to be n" << endl;
   return;
}
 
int main(int argc, char * argv[]) {
    int a = 12;
    int b = 15;
    bool doTime = false;
    int arg;

    int i = 1; 
    while (i < argc) {
        if (strcmp(argv[i],"-l")==0) {
	   i++;
	   if (i < argc) {
	       arg = -1;
               arg = atoi(argv[i]); 
               SetAB(arg, a,b);
	       i++;
	   } else {
	       cout << "-l requires an integer level, using default " << endl;
	   }
	} else if (!strcmp(argv[i],"-a")) {
	   i++;
	   a = atoi(argv[i]);
	   if (a < 2) {
	      a = 12;
	   }
	   i++;
	} else if (!strcmp(argv[i],"-b")) {
	   i++;
	   b = atoi(argv[i]);
	   if (b < 2) {
	      b = 15;
	   }
	   i++;
	} else if (!strcmp(argv[i],"-t")) {
	   doTime = true;
	   i++;
	} else if (!strcmp(argv[i],"-r")) {

	   boost::mt19937 rng;
	   boost::uniform_int<> dist(2,1<<30);
	   boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
	         bigDie(rng,dist);

	   rng.seed(static_cast<unsigned int>(std::time(0)));

	   a = bigDie(); 
	   b = bigDie();
	   i++;
	} else if (!strcmp(argv[i],"-R")) {
	   int upperLimit;

	   i++;
	   upperLimit = atoi(argv[i]);
	   i++;

	   if (upperLimit < 2) {
	      upperLimit = 1 << 30;
	   }

	   boost::mt19937 rng;
	   boost::uniform_int<> dist(2,upperLimit);
	   boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
	         bigDie(rng,dist);

	   rng.seed(static_cast<unsigned int>(std::time(0)));

	   a = bigDie(); 
	   b = bigDie();

	} else if (!strcmp(argv[i],"-h")) {
	   args(argv[0]);
	   i++;
	} else if (!strcmp(argv[i],"-d")) {
	   DEBUG = true;
	   i++;
	} else {
	   cout << "unknown argument " << argv[i] << endl;
	   args(argv[0]);
	   i++;
	}
    }

    cout << endl;
    cout << "A = " << a << " and has " << int(log(a)/log(2)) << " bits" << endl;
    cout << "B = " << b << " and has " << int(log(b)/log(2)) << " bits" << endl;
    cout << endl;

    TimeRoutine("Euclid", Euclid, a,b,doTime);
    TimeRoutine("Brute ", Brute, a,b,doTime);
    TimeRoutine("School", OldSchool, a,b,doTime);

    return 0;
}

int Euclid(int a, int b) {
    if (b==0) {
       return a;
    } else  {
       return Euclid(b, a%b);
    }
}

int Max(int a, int b) {
   if (a > b) {
      return a;
   } else {
      return b;
   }
}

int Min(int a, int b) {
   if (a < b) {
      return a;
   } else {
      return b;
   }
}

int Brute(int a, int b) {
    int n;

    n = Min(a,b);
    while (a%n != 0 or b%n !=0) {
       n--;
    }

    return n;
}

void Sieve(vector<int> & primes, int n) {
    int i,j;
    vector<bool> A;

    primes.erase(primes.begin(), primes.end());

    for(i=0;i<=n;i++) {
       A.push_back(true);
    }

    for(i=2;i<=sqrt(n);i++) {
       if(A[i]) {
          primes.push_back(i);
	  if (DEBUG > 0)  {
	      cout << "inserting " << i << endl;
	  }
          for(j=i*i;j<=n;j+= i) {
	     if (DEBUG > 0) {
	          cout << "crossing off " << j << endl;
	     }
	     A[j] = false;
	  }
       }
    }

    for(j=i;j<=n;j++) {
        if(A[j]) {
	   primes.push_back(j);
	}
    }

    if( DEBUG > 0) {
        cout << "The primes between 2 and " << n << " are :" << endl;
        for(int i=0;i<primes.size();i++) {
           cout << primes[i] << endl;
        }
    }

    return;
}

void Factor(int n, vector<int> & factors, const vector<int> & primes) {
    int i;
    int num;

    factors.erase(factors.begin(), factors.end());

    if(DEBUG) {
        cout<< "Factoring : "  << n << endl;
    }

    i = 0;
    num = n;

    while(num > 1) {
       if (num%primes[i] == 0) {
           factors.push_back(primes[i]);
	   num = num/primes[i];
       } else {
          i++;
       }
    }

    if (DEBUG) {
        cout << "Factors of " << n << ": ";
        for(int i =0;i<factors.size();i++) {
           cout << factors[i] << " " ;
        }
        cout << endl;
    }


    return;
}

void CommonFactors(const vector<int> & A,
                   const vector<int> & B, 
		   vector<int>  & common){

     int i,j;

     common.erase(common.begin(),common.end());

     i = 0;
     j = 0;
     while (i < A.size() and j < B.size()) {
        if (A[i] == B[j]) {
	   // same, common, push it and move both
           common.push_back(A[i]);
	   i++;
	   j++;
	} else if (A[i] > B[j]) {
	   // A is too big, so bump B;
	   j++;
	} else {
	   i++;
	}
     }

     return;
}

int School(int a, int b, const vector<int> & primes) {
    vector<int> aFactors, bFactors, common;
    int gcd=1;
    int i;

    Factor(a, aFactors,primes);
    Factor(b, bFactors,primes);
    CommonFactors(aFactors, bFactors, common);

    for(i=0;i<common.size();i++) {
      if (DEBUG) {
         cout << common[i] << " is a common factor " << endl;
      }
      gcd *= common[i];
    }

    return gcd;

}

void SetAB(int arg, int & a, int & b) {
       switch(arg) {
          case 0:
             // the book's example
             a = 60;
	     b = 24;
	     break;
	  case 1:
             // 2^2 x 3^3 x 5^5 x 7   and 2^2 x 3^2 x 5^4 (to test factoring)
             a = 45000;
	     b = 787500;
	     break;
	  case 2:
             // an exercise in the book
             a=31415;
	     b= 14142;
	     break;
	  case 3:
             // not to big prime number, look out old school!
             a = 80147*2*3*3;
             b = 21011*2*3*3;
	     break;
	  case 4:
             // two bigger prime numbers
             a=  4095036*2*3*5;
             b = 4095043*2*3*5;
	     break;
	  default:
	     a = 40950391 * 2 * 3;
	     b = 40950401 * 2 * 3;
       }
       return;
}


int OldSchool(int a, int b) {
    int g;
    vector<int> primeList;
    Sieve(primeList, Max(a,b));
    return School(a,b,primeList);
}

void TimeRoutine(string name, boost::function < int (int, int)> fun, 
                   int arg1,int arg2, bool doTime){

    int gcd;
    boost::posix_time::ptime startTime, endTime;
    boost::posix_time::time_duration totalTime;

    startTime = boost::posix_time::microsec_clock::local_time();
    gcd = fun(arg1, arg2);
    endTime = boost::posix_time::microsec_clock::local_time();
    cout <<"\t" << name << ": " << setw(15) << gcd;
    if (doTime) {
       totalTime = endTime-startTime;
       cout << "\tTime: " << totalTime;
    }
    cout << endl;
    if (!Test(arg1, arg2, gcd)) {
        cout << endl;
        cout << "ERROR ERROR ERROR ERROR ERROR ERROR " << endl;
        cout << endl << endl << endl;
    }

    return;
}

bool Test(int a, int b, int d){
     int r1, r2, d1, d2, p1, p2;

     r1 = a%d;
     r2 = b%d;
     d1 = a/d;
     d2 = b/d;
     p1 = d1*d;
     p2 = d2*d;

     if (r1 != 0) {
         cout<< "Error " << a << "/" << d << " = " << d1 << " R " << r1 << endl;
	 cout << "The remainder should be 0  but is " << r1 << endl;

	 return false;
     }

     if (r2 != 0) {
         cout<< "Error " << b << "/" << d << " = " << d2 << " R " << r2 << endl;
	 cout << "The remainder should be 0  but is " << r2 << endl;

	 return false;
     }

     if (p1 != a) {
        cout << "Error, the proposed GCD is " << d << " but (" 
	     << a << "/" << d << ") * " << d << "  = " <<  p1 << endl;
	cout << "    and " << p1 << " !=  " << a << endl;
	return false;
     }
     if (p2 != b) {
        cout << "Error, the proposed GCD is " << d << " but (" 
	     << b << "/" << d << ") * " << d << "  = " <<  p2 << endl;
	cout << "    and " << p2 << " !=  " << b << endl;
	return false;
     }
     return true;

}