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

void main()
{
    //program d9r4
    //driver for routine svdfit
    int i,mp,np,npt = 100;
    double chisq;
	double spread = 0.02;
    int npol = 5;
    double x[101], y[101], sig[101], a[6], cvm[26];
    double u[501], v[26], w[6];
    //polynomial fit
    long idum = -911;
    mp = npt;
    np = npol;
    for (i = 1; i<=npt; i++)
	{
        x[i] = 0.02 * i;
        y[i] = 1.0 + x[i] * (2.0 + x[i] * (3.0 + x[i] * (4.0 + x[i] * 5.0)));
        y[i] = y[i] + spread * gasdev(idum);
		sig[i] = spread;
    }
    svdfit(x, y, sig, npt, a, npol, u, v, w, mp, np, chisq, "fpoly");
    svdvar(v, npol, np, w, cvm, npol);
    cout<<endl;
    cout<<"Polynomial fit"<<endl;
    cout<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
    for (i = 1; i<=npol; i++)
	{
		cout.width(12);
		cout<<a[i]<<"      +-";
		cout.width(12);
        cout<<sqrt(cvm[(i-1)*npol+i])<<endl;
    }
    cout<<endl;
    cout<<"Chi-squared  ";
	cout<<chisq<<endl;
    svdfit(x, y, sig, npt, a, npol, u, v, w, mp, np, chisq, "fleg");
    svdvar(v, npol, np, w, cvm, npol);
    cout<<endl;
    cout<<"Legendre polynomial fit"<<endl;
    for (i = 1; i<=npol; i++)
	{
		cout.width(12);
		cout<<a[i]<<"      +-";
		cout.width(12);
        cout<<sqrt(cvm[(i-1)*npol+i])<<endl;
    }
    cout<<endl;
    cout<<"Chi-squared  "<<chisq<<endl;
}