#include <iostream>
#include <iomanip>

#include <math.h>
#include <time.h>


#define FIXED_X 0


int main()
{
   /* initialize random generator */
   int seed = time(NULL);
   std::cout << "Reseeding the random number generator (t="
             << seed << ")." << std::endl;
   srand (seed);

   int N = 100000;
   std::cout << "Running " << N
             << " random test cases with:" << std::endl;
   
   double h_t = 0.0162688; /* (3*10^(-6)/cos(0.8))^(1/3) evaluated in MatLab */
   double h_e = 0.0062500; /* Experimental value for h                       */

   std::cout << "     h_e = " << h_e << std::endl;   
   std::cout << "     h_t = " << h_t << std::endl;

   double m_t = 0, m_e = 0;
   double x = 0.9;

   for(int i=0; i<N; i++)
   {
      if (!FIXED_X)
         x = 0.8+(0.2)*((double)rand() / ((double)(RAND_MAX)));
      double r1 = (-1+2*((double)rand() / ((double)(RAND_MAX))))/pow(10,6);
      double r2 = (-1+2*((double)rand() / ((double)(RAND_MAX))))/pow(10,6);

      if (fabs(((sin(x+h_t)+r1)-(sin(x-h_t)+r2))/(2*h_t)-cos(x)) > m_t)
         m_t = fabs(((sin(x+h_t)+r1)-(sin(x-h_t)+r2))/(2*h_t)-cos(x));

      if (fabs(((sin(x+h_e)+r1)-(sin(x-h_e)+r2))/(2*h_e)-cos(x)) > m_e)
         m_e = fabs(((sin(x+h_e)+r1)-(sin(x-h_e)+r2))/(2*h_e)-cos(x));
   }

   std::cout << "Maximal error with h_e : m_e = " << m_e << std::endl;
   std::cout << "Maximal error with h_t : m_t = " << m_t << std::endl;   
}