/*	File:	Romberg.cpp



			Christer Karlsson

This code implements Rombergs Algorithm.



    Allowed inputs are real numbers (long double).

    Define f in the function procedure



	NOTE!  All arrays have an upper limit of n=10.		*/



#include <iostream>

#include <cmath>

#include <iomanip>		// Needed too manipulate the display



using namespace std;





// the function

long double f(long double x)

{

    long double f;

    f=8*(sqrt(1-x*x)-x);

    return f;

}



// The Integration function

void Romberg(long double a, long double b, int n, long double R[10][10])

{

    int i, j, k;

    long double h, sum;



    h=b-a;

    R[0][0]=(h/2)*(f(a)+f(b));

    for (i=1; i<=n; i++)

    {

        h=h/2;

        sum=0;

        for (k=1; k<=int(pow(2.0,i)-1); k +=2)

            sum +=f(a+k*h);



        R[i][0]=R[i-1][0]/2+sum*h;

        for(j=1; j<=i; j++)

            R[i][j]=R[i][j-1]+(R[i][j-1]-R[i-1][j-1])/(pow(4.0,j)-1);

    }

}

int main(int argc, char *argv[])

{

    long double R[10][10];

    int n=5, digits=13;

    long double a=0, b=1/sqrt(2.0);



    // Call the integrating function

    Romberg(a, b, n, R);



    // Set a fixed amount of digits to show

    cout.setf(ios::fixed );

    cout.setf(ios::showpoint);

    cout <<  setprecision(digits) << endl;

    // Printout the diagram

    for(int i=0; i<n; i++)

    {

        for(int j=0; j<=i; j++)

        {

            cout << setw(digits+2) << R[i][j] << ' ';

        }



        cout << endl;

    }

    // Unset the fixed point

    cout.unsetf(ios::fixed);

    cout.unsetf(ios::showpoint);



    system("PAUSE");

    return 0;

}