#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 polint(double xa[], double ya[], double n, double x, double& y, double& dy)
{
  double c[11], d[11];
  int i,m,ns = 1;
  double dif,dift,ho,hp,den,w;
  dif= fabs(x - xa[1]);
  for (i = 1; i<=n; i++)
  {
      dift = fabs(x - xa[i]);
      if (dift < dif)
	  {
          ns = i;
          dif = dift;
      }
      c[i] = ya[i];
      d[i] = ya[i];
  }
  y = ya[ns];
  ns = ns - 1;
  for (m = 1; m<=n-1; m++)
  {
      for (i = 1; i<=n - m; i++)
	  {
          ho = xa[i] - x;
          hp = xa[i + m] - x;
          w = c[i + 1] - d[i];
          den = ho - hp;
          if (den == 0.0)
		  {
             cout<< "pause"<<endl;
             return;
          }
          den = w / den;
          d[i] = hp * den;
          c[i] = ho * den;
      }
      if (2 * ns < n - m)
	  {
          dy = c[ns + 1];
	  }
      else
	  {
          dy = d[ns];
          ns = ns - 1;
      }
      y = y + dy;
  }
}

void qromb(double a, double b, double& ss)
{
    double dss,eps = 0.000001;
    const int jmax = 20;
    int jmaxp = jmax + 1;
    int j,k = 5;
    int km = k - 1;
    double s[21], h[21];
    h[1] = 1.0;
    for (j = 1; j<=jmax; j++)
	{
        trapzd(a, b, s[j], j);
        if (j >= k)
		{
            polint(h, s, k, 0.0, ss, dss);
            if (fabs(dss) < eps * fabs(ss))
			{
				return;
			}
        }
        s[j + 1] = s[j];
        h[j + 1] = 0.25 * h[j];
    }
    cout<<"too many steps."<<endl;
}

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