// L-13 MCS 572 Wed 8 Feb 2012 : romberg4piqd.cpp
// Applies Romberg integration with quad double arithmetic
// to the integral of the rational function
// (16 x - 16)/(x^4 - 2*x^3 + 4*x - 4) over [0,1].
// The value of this integral equals pi.
// The code below uses OpenMP to distribute the function evaluations.

#include <iostream>
#include <iomanip>
#include <vector>
#include <qd/qd_real.h>
#include <cstdlib>
#include <ctime>
#include <omp.h>
using namespace std;

vector<qd_real> comptrap
  ( qd_real f ( qd_real x ),
    qd_real a, qd_real b, int n );
// The vector on return contains at position i
// the result of the composite trapezoidal rule 
// using i+1 subintervals,
// for i running from 0 to n-1.

void write_vector ( vector<qd_real> v );
// Writes the vector in scientific notation
// with 64 decimal places and the difference with Pi.

void romberg_extrapolation ( vector<qd_real>& t );
// Applies Romberg extrapolation to t,
// obtained with the composite trapezoidal rule.

qd_real f ( qd_real x );
// Returns (16 x - 16)/(x^4 - 2*x^3 + 4*x - 4) */

int main ( int argc, char *argv[] )
{
   int n;
   if(argc > 1)
      n = atoi(argv[1]);
   else
   {
      cout << "Give n : "; 
      cin >> n;
   }
   omp_set_num_threads(2);
   clock_t tstart,tstop;
   tstart = clock();
   vector<qd_real> t;
   t = comptrap(f,qd_real(0.0),qd_real(1.0),n);
   cout << "Initial approximations :" << endl;
   write_vector(t);
   cout << "Romberg integration :" << endl;
   romberg_extrapolation(t);
   write_vector(t);
   tstop = clock();
   cout << fixed << "elapsed time : " 
        << (tstop - tstart)/double(CLOCKS_PER_SEC)
        << " seconds" << endl;
   return 0;
}

vector<qd_real> comptrap
  ( qd_real f ( qd_real x ),
    qd_real a, qd_real b, int n )
{
   vector<qd_real> t(n);
   int id,jstart,jstop;
   qd_real val;

   qd_real h = b - a;

   t[0] = (f(a) + f(b))*h/2;
   for(int i=1, m=1; i<n; i++,m=m*2)
   {
      h = h/2;
      t[i] = 0.0;
      #pragma omp parallel private(id,jstart,jstop,val)
      {
         id = omp_get_thread_num();
         jstart = id*m/2;
         jstop = (id+1)*m/2;
         for(int j=jstart; j<jstop; j++)
            val += f(a+h+j*2*h);
         #pragma omp critical
            t[i] += val;
      }
      t[i] = t[i-1]/2 + h*t[i];
   }
   return t;
}

void write_vector ( vector<qd_real> v )
{
   qd_real pi = qd_real::_pi;
   cout << scientific;
   for(int i=0; i<v.size(); i++)
      cout << setprecision(64) << v[i]
           << "  " << setprecision(3)
           << v[i] - pi << endl;
 
}

void romberg_extrapolation ( vector<qd_real>& t )
{
   int n = t.size();
   qd_real e[n][n];
   int m = 0;
   for(int i=0; i<n; i++)
   {
      e[i][0] = t[i];
      for(int j=1; j<n; j++) e[i][j] = 0.0;
   }
   for(int i=1; i<n; i++)
   {
      for(int j=1,m=2; j<=i; j++,m=m+2)
      {
         qd_real r = pow(2.0,m);
         e[i][j] = (r*e[i][j-1] - e[i-1][j-1])/(r-1);
      }
   }
   for(int i=1; i<n; i++) t[i] = e[i][i];
}

qd_real f ( qd_real x )
{
   qd_real y = 16*x - 16;
   qd_real d = (((x-2)*x)*x + 4)*x - 4;
   return y/d;
}