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

char choose[10];

double funcl(double x)
{
    return sqrt(x) / sin(x);
}

double funcu(double x)
{
    double pi = 3.1415926;
    return sqrt(pi - x) / sin(x);
}

double funcinf(double x)
{
    return sin(x) / (pow(x, 2));
}

double funcend(double x)
{
    return exp(-x) / sqrt(x);
}

double func(double x)
{
    //dim choose as string
    if (strcmp(choose,"funcl")==0)
	{
		return funcl(x);
	}
    if (strcmp(choose,"funcu")==0)
	{
		return funcu(x);
	}
    if (strcmp(choose,"funcinf")==0)
	{
		return funcinf(x);
	}
    if (strcmp(choose,"funcend")==0)
	{
		return funcend(x);
	}
	return -1;
}

double inf(double x)
{
    return func(1 / x) / pow(x , 2);
}


void midinf(double aa, double bb, double& s, int n)
{
    int ii,tnm,j;
	double del,ddel,x,sum;
	double b = 1.0 / aa;
    double a = 1.0 / bb;
    if (n == 1)
	{
        s = (b - a) * inf(0.5 * (a + b));
        ii = 1;
	}
    else
	{
        ii = (int)pow(3 , n - 2);
        tnm = ii;
        del = (b - a) / (3.0 * tnm);
        ddel = del + del;
        x = a + 0.5 * del;
        sum = 0.0;
        for (j = 1; j<=ii; j++)
		{
            sum = sum + inf(x);
            x = x + ddel;
            sum = sum + inf(x);
            x = x + del;
        }
        s = (s + (b - a) * sum / tnm) / 3.0;
    }
}

double sql(double x, double aa)
{
    return 2 * x * func(aa + pow(x , 2));
}


void midsql(double aa, double bb, double& s, int n)
{
	int it,tnm;
	double del,ddel,x,sum;
	double b,a;
    b = sqrt(bb - aa);
    a = 0.0;
    if (n == 1)
	{
        s = (b - a) * sql(0.5 * (a + b), aa);
        it = 1;
	}
    else
	{
        it = (int)pow(3, n - 2);
        tnm = it;
        del = (b - a) / (3.0 * tnm);
        ddel = del + del;
        x = a + 0.5 * del;
        sum = 0.0;
        for (int j = 1; j<=it; j++)
		{
            sum = sum + sql(x, aa);
            x = x + ddel;
            sum = sum + sql(x, aa);
            x = x + del;
        }
        s = (s + (b - a) * sum / tnm) / 3.0;
    }
}

double squ(double x, double bb)
{
    return 2 * x * func(bb - pow(x , 2));
}

void midsqu(double aa, double bb, double& s, int n)
{
	int it,tnm;
	double del,ddel,x,sum;
	double b,a;
    b = sqrt(bb - aa);
    a = 0.0;
    if (n == 1)
	{
        s = (b - a) * squ(0.5 * (a + b), bb);
        it = 1;
	}
    else
	{
        it = (int)pow(3, n - 2);
        tnm = it;
        del = (b - a) / (3.0 * tnm);
        ddel = del + del;
        x = a + 0.5 * del;
        sum = 0.0;
        for (int j = 1; j<=it; j++)
		{
            sum = sum + squ(x, bb);
            x = x + ddel;
            sum = sum + squ(x, bb);
            x = x + del;
        }
        s = (s + (b - a) * sum / tnm) / 3.0;
    }
}

void midpnt(double a, double b, double& s, int n)
{
	int ii,tnm,j;
	double del,ddel,x,sum;

    if (n == 1)
	{
        s = (b - a) * func(0.5 * (a + b));
        ii = 1;
	}
    else
	{
        ii = (int)pow(3 , n - 2);
        tnm = ii;
        del = (b - a) / (3.0 * tnm);
        ddel = del + del;
        x = a + 0.5 * del;
        sum = 0.0;
        for (j = 1; j<=ii; j++)
		{
            sum = sum + func(x);
            x = x + ddel;
            sum = sum + func(x);
            x = x + del;
        }
        s = (s + (b - a) * sum / tnm) / 3.0;
    }
}

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 qromo(double a, double b, double& ss, char pick[])
{
    double dss,eps = 0.00003;
    int jmaxp,jmax = 14;
    jmaxp = jmax + 1;
    int k = 7;
    double s[15], h[15];
    h[1] = 1.0;
    for (int j = 1; j<=jmax; j++)
	{
        if (strcmp(pick,"midpnt")==0)
		{
			midpnt(a, b, s[j], j);
		}
        if (strcmp(pick,"midinf")==0)
		{
			midinf(a, b, s[j], j);
		}
        if (strcmp(pick,"midsql")==0)
		{
			midsql(a, b, s[j], j);
		}
        if (strcmp(pick,"midsqu")==0)
		{
			midsqu(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] = h[j] / 9.0;
    }
    cout<< "too many steps."<<endl;
}


void main()
{
    //program d3r6
    //driver for routine qromo
    double res1,res2,result,x1 = 0.0;
    double x2 = 1.5707963;
    double x3 = 3.1415926;
    double ainf = 1e+20;
    cout<<endl;
    cout<<"improper integrals:"<<endl;
    cout<<endl;
    strcpy(choose,"funcl");
    qromo(x1, x2, result, "midsql");
    cout<<"function: sqr(x)/sin(x)      interval: (0,pi/2)"<<endl;
    cout<<"using: midsql                result: ";
    cout<<result<<endl;
    cout<<endl;
    strcpy(choose,"funcu");    
    qromo(x2, x3, result, "midsqu");
    cout<<"function: sqr(pi-x)/sin(x)   interval: (pi/2,pi)"<<endl;
    cout<<"using: midsqu                result: ";
    cout<<result<<endl;
    cout<<endl;
	strcpy(choose,"funcinf");
    qromo(x2, ainf, result, "midinf");
    cout<<"function: sin(x)/x^2        interval: (pi/2,infty)"<<endl;
    cout<<"using: midinf                result: ";
    cout<<result<<endl;
    cout<<endl;
	strcpy(choose,"funcinf");
    qromo(-ainf, -x2, result, "midinf");
    cout<<"function: sin(x)/x^2        interval: (-infty,-pi/2)"<<endl;
    cout<<"using: midinf                result: ";
    cout<<result<<endl;
    cout<<endl;
	strcpy(choose,"funcend");
    qromo(x1, x2, res1, "midsql");
    qromo(x2, ainf, res2, "midinf");
    cout<<"function: exp(-x)/sqr(x)     interval: (0,infty)"<<endl;
    cout<<"using: midsql,midinf         result: ";
    cout<<(res1 + res2)<<endl;
}