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

double func(double x)
{
    return x * exp(-x);
}

void gauleg(double x1, double x2, double x[], double w[], int n)
{
	int i,j,m;
    double p1,p2,p3,z,z1,pp,xm,xl,eps = 0.00000000000003;
    m = (n + 1) / 2;
    xm = 0.5 * (x2 + x1);
    xl = 0.5 * (x2 - x1);
    for (i = 1; i<=m; i++)
	{
        z = cos(3.141592654 * (i - 0.25) / (n + 0.5));
		do
		{
	        p1 = 1.0;
			p2 = 0.0;
			for (j = 1; j<=n; j++)
			{
				p3 = p2;
				p2 = p1;
				p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
			}
			pp = n * (z * p1 - p2) / (z * z - 1.0);
			z1 = z;
			z = z1 - p1 / pp;
		}while (fabs(z - z1) > eps);
		x[i] = xm - xl * z;
        x[n + 1 - i] = xm + xl * z;
        w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
        w[n + 1 - i] = w[i];
    }
}

void main()
{
    //program d3r8
    //driver for routine gauleg
    int i,npoint = 10;
    double x1 = 0.0;
    double x2 = 1.0;
    double x3 = 10.0;
    double x[11], w[11];
    gauleg(x1, x2, x, w, npoint);
    cout<<endl;
    cout<<"      #       x(i)        w(i)"<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
    for (i = 1; i<=npoint; i++)
	{
		cout.width(2);
        cout<<"      "<<i;
		cout.width(12);
        cout<<x[i];
		cout.width(12);
        cout<<w[i]<<endl;
    }
    //demonstrate the use of gauleg for an integral
    gauleg(x1, x3, x, w, npoint);
    double xx = 0.0;
    for (i = 1; i<=npoint; i++)
	{
        xx = xx + w[i] * func(x[i]);
    }
    cout<<endl;
    cout<<"integral from gauleg:  "<<xx<<endl;
    cout<<"actual value:  ";
    cout<<(1.0 - (1.0 + x3) * exp(-x3))<<endl;
}