/****************************************************************
Copyright (C) :    University of Michigan
Project:           NERS 544 - Class example
File:              tt.cpp      
Purpose:           Torture test of critical functions used in
                   Monte Carlo codes
Compiler:          g++ -ansi -pedantic -Wconversion
                   Use -O2 as well
Exacution          time a.out                   
Programmer:        Alex Bielajew
Revision Date:     09/15/2016
****************************************************************/

#include <iostream>
#include <cmath>

using namespace std;

int rand(int in) {return 69069*in;}
// The magic multiplier 69069 is cited by Marsaglia as the "best" for 32 bits 

int main(void)
{  //Starting seed = ";
   int iFirst = 987654321; // Starting seed: You can substitute any odd number 
                           // and the result should be the same. Try it!

   int Ncycles   = 1073741824;  // 2^(30), the sequence length of this generator
   int Imax      = 2147483647;  // maximum int
   int Imin      = -2147483648; // minimum int
   double halfPi = 1.570796326794897; // From Matlab
   double scale  = 1.000000000000003; // To put all algs on an == footing
 
   cout << "This code will torture test various functions used in MC\n";
   cout << "Input: 0 for mcrng32()\n"
        << "       1 for log()\n"
        << "       2 for sqrt()\n"
        << "       3 for sin()\n"
        << "       4 for cos()\n"
        << "       5 for exp()\n";
   int choice;
   cin >> choice;
   
   if (choice == 0) // rng < 2 s
   {  int ixx = iFirst;
      double r;
      for (int i = 0; i < Ncycles; i = i + 1)
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = scale*(0.5 + ixx*2.328306436538696e-10);
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final rn  : " << r   << endl; // Forces the optimizer to execute
   }

   else if (choice == 1) // the log function < 48 seconds
   {  int ixx = iFirst;
      double r = 1;
      for (int i = 0; i < Ncycles; i = i + 1) // A
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = scale*(0.5 + ixx*2.328306436538696e-10);
         r = r + log(r); // r + fools the optimizer
         
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final double: " << r << endl; // Forces the optimizer to execute
   }

   else if (choice == 2) // the sqrt function < 6 seconds for scale = 1
   {  int ixx = iFirst;
      double r = 1;
      for (int i = 0; i < Ncycles; i = i + 1) // 
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = scale*(0.5 + ixx*2.328306436538696e-10);
         r = r + sqrt(r); // r + fools the optimizer
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final double: " << r << endl; // Forces the optimizer to execute
   }

   else if (choice == 3) // the sin function < 26 seconds 
   {  int ixx = iFirst;
      double r = 1;
      for (int i = 0; i < Ncycles; i = i + 1) //
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = halfPi*(0.5 + ixx*2.328306436538696e-10);
         r = r + sin(r); // r+ fools the optimizer
         
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final double: " << r << endl; // Forces the optimizer to execute
   }

   else if (choice == 4) // the cos function < 31 seconds 
   {  int ixx = iFirst;
      double r = 1;
      for (int i = 0; i < Ncycles; i = i + 1)
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = halfPi*(0.5 + ixx*2.328306436538696e-10);
         r = r + cos(r); // r+ fools the optimizer
         
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final double: " << r << endl; // Forces the optimizer to execute
   }

   else if (choice == 5) // the exp function < 30 seconds for scale = 1
   {  int ixx = iFirst;
      double r = 1;
      for (int i = 0; i < Ncycles; i = i + 1) // < 30 seconds
      {  ixx = 69069*ixx; // Knuth's famous multiplier
         r = scale*(0.5 + ixx*2.328306436538696e-10);
         r = r + exp(r); // r+ fools the optimizer
         
      }
      cout << "Final seed: " << ixx << endl; // Should equal iFirst!
      cout << "Final double: " << r << endl; // Forces the optimizer to execute
   }

   return 0;
}