#include "iostream.h"
#include "math.h"

double func(double x)
{
    return pow(x, 2) * (pow(x,2) - 2.0) * sin(x);
}

double fint(double x)
{
    //integral of func
	double aaa;
    aaa = 4.0 * x * (pow(x,2) - 7.0) * sin(x);
    return aaa - (pow(x , 4) - 14.0 * pow(x , 2) + 28.0) * cos(x);
}

void trapzd(double a, double b, double& s, int n)
{
	double del,x,sum;
	int j,it, tnm;
    if (n ==1)
	{
        s = 0.5 * (b - a) * (func(a) + func(b));
        it = 1;
	}
    else
	{
        it = (int)pow(2 , n - 2);
        tnm = it;
        del = (b - a) / tnm;
        x = a + 0.5 * del;
        sum = 0.0;
        for (j = 1; j<=it; j++)
		{
            sum = sum + func(x);
            x = x + del;
        }
        s = 0.5 * (s + (b - a) * sum / tnm);
    }
}

void qsimp(double a, double b, double& s)
{
    double eps = 0.000001;
    const int jmax = 20;
    double ost = -1e+30;
    double st, os = -1e+30;
    for (int j = 1; j<=jmax; j++)
	{
        trapzd(a, b, st, j);
        s = (4.0 * st - ost) / 3.0;
        if (fabs(s - os) < eps * fabs(os))
		{
			return;
		}
        os = s;
        ost = st;
    }
    cout<<"too many steps."<<endl;
}


void main()
{
    //program d3r3
    //driver for routine qsimp
    double pio2 = 1.5707963;
    double s, b,a = 0.0;
    b = pio2;
    cout<<endl;
    cout<<"integral of func computed with qsimp"<<endl;
    cout<<endl;
    cout<<"actual value of integral is:";
    cout<<(fint(b) - fint(a))<<endl;
    cout<<endl;
    qsimp(a, b, s);
    cout<<"result from routine qtrap is:";
    cout<<s<<endl;
}