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

int ndatat;
double xt[1001], yt[1001], arr[1001], aa, abdevt;

int sgn(double x)
{
	if (x>0)
	{
		return 1;
	}
	else
	{
		if(x<0)
		{
			return -1;
		}
	}
	return 0;
}

double gammln(double xx)
{
	double cof[7];
	cof[1] = 76.18009173;
	cof[2] = -86.50532033;
	cof[3] = 24.01409822;
	cof[4] = -1.231739516;
	cof[5] = 0.00120858003;
	cof[6] = -0.00000536382;
	double ser,x,tmp,stp = 2.50662827465;
	double half = 0.5;
	double one = 1.0;
	double fpf = 5.5;
	int j;
	x = xx - one;
	tmp = x + fpf;
	tmp = (x + half) * log(tmp) - tmp;
	ser = one;
	for (j = 1; j<=6; j++)
	{
		x = x + one;
		ser = ser + cof[j] / x;
	}
	return tmp + log(stp * ser);
}

void gcf(double& gammcf, double a, double x, double gln)
{
	int an,n,itmax = 100;
	double anf,g,ana,eps = 0.0000003;
	gln = gammln(a);
	double gold = 0.0;
	double a1,a0 = 1.0;
	a1 = x;
	double b0 = 0.0;
	double b1 = 1.0;
	double fac = 1.0;
	for (n = 1; n<=itmax; n++)
	{
		an = n;
		ana = an - a;
		a0 = (a1 + a0 * ana) * fac;
		b0 = (b1 + b0 * ana) * fac;
		anf = an * fac;
		a1 = x * a0 + anf * a1;
		b1 = x * b0 + anf * b1;
		if (a1 != 0.0)
		{
			fac = 1.0 / a1;
			g = b1 * fac;
			if (fabs((g - gold) / g) < eps)
			{
				gammcf = exp(-x + a * log(x) - gln) * g;
				return;
			}
			gold = g;
		}
	}
	cout<< "a too large, itmax too small"<<endl;
	gammcf = exp(-x + a * log(x) - gln) * g;
}

void gser(double& gamser, double a, double x, double gln)
{
	int n,itmax = 100;
	double ap,sum,del,eps = 0.0000003;
	gln = gammln(a);
	if (x <= 0.0)
	{
		if (x < 0.0)
		{
			cout<<"pause in gser"<<endl;
			return;
		}
		gamser = 0.0;
		return;
	}
	ap = a;
	sum = 1.0 / a;
	del = sum;
	for (n = 1; n<=itmax; n++)
	{
		ap = ap + 1.0;
		del = del * x / ap;
		sum = sum + del;
		if (fabs(del) < fabs(sum) * eps)
		{
			gamser = sum * exp(-x + a * log(x) - gln);
			return;
		}
	}
	cout<< "a too large, itmax too small";
	gamser = sum * exp(-x + a * log(x) - gln);
}

double gammq(double a, double x)
{
	double gamser,gammcf,gln=0;
	if (x < 0.0 || a <= 0.0)
	{
		cout<< "pause"<<endl;
		exit(1);
	}
	if (x < a + 1.0)
	{
		gser(gamser, a, x, gln);
		return 1.0 - gamser;
	}
	else
	{
		gcf(gammcf, a, x, gln);
		return gammcf;
	}
}

double ran1(long& idum)
{
    int j,iff=-1;
	static long ix1,ix2,ix3;
	static double r[98];
    long m1 = 259200; long m2 = 134456; long m3 = 243000;
	long ia1 = 7141; long ia2 = 8121; long ia3 = 4561;
	long ic1 = 54773; long ic2 = 28411; long ic3 = 51349;
	double rm1 = 0.0000038580247; double rm2 = 0.0000074373773;    
    if (idum < 0 || iff == 0)
	{
        iff = 1;
        ix1 = (ic1 - idum) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix2 = (ia2 * ix2 + ic2) % m2;
            r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
        }
        idum = 1;
    }
    ix1 = (ia1 * ix1 + ic1) % m1;
    ix2 = (ia2 * ix2 + ic2) % m2;
    ix3 = (ia3 * ix3 + ic3) % m3;
    j = 1 + int((97 * ix3) / m3);
    if (j > 97 || j < 1)
	{
		cout<<"abnormal exit in ran1"<<endl;
		exit(1);
	}
	double temp=r[j];
    r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
	return temp;

}

double gasdev(long& idum)
{
    static int iset;
	static double gset;
	double v1,v2,r,fac;
    if (iset == 0)
	{
		do
		{
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			r = v1 * v1 + v2 * v2;
		}while (r >= 1.0 || r == 0);
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		iset = 1;
		return v2 * fac;
	}
    else
	{
		iset = 0;
		return gset;
    }
}

void fit(double x[], double y[], int ndata, double sig[], int mwt, double& a,
		double& b, double& siga, double& sigb, double& chi2, double& q)
{
	int i;
    double sigdat,t,sxoss,wt,ss,sx = 0.0;
    double sy = 0.0;
    double st2 = 0.0;
    b = 0.0;
    if (mwt != 0)
	{
        ss = 0.0;
        for (i = 1; i<=ndata; i++)
		{
            wt = 1.0 / (sig[i] * sig[i]);
            ss = ss + wt;
            sx = sx + x[i] * wt;
            sy = sy + y[i] * wt;
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            sx = sx + x[i];
            sy = sy + y[i];
        }
        ss = ndata;
    }
    sxoss = sx / ss;
    if (mwt != 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            t = (x[i] - sxoss) / sig[i];
            st2 = st2 + t * t;
            b = b + t * y[i] / sig[i];
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            t = x[i] - sxoss;
            st2 = st2 + t * t;
            b = b + t * y[i];
		}
    }
    b = b / st2;
    a = (sy - sx * b) / ss;
    siga = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
    sigb = sqrt(1.0 / st2);
    chi2 = 0.0;
    if (mwt == 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow((y[i] - a - b * x[i]) , 2);
        }
        q = 1.0;
        sigdat = sqrt(chi2 / (ndata - 2));
        siga = siga * sigdat;
        sigb = sigb * sigdat;
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
        }
        q = gammq(0.5 * (ndata - 2), 0.5 * chi2);
    }
}

void sort(int n, double ra[])
{
    int i,j,l,ir;
	double rra;
	l = int(n / 2) + 1;
    ir = n;
    do
	{
        if (l > 1)
		{
            l = l - 1;
            rra = ra[l];
		}
        else
		{
            rra = ra[ir];
            ra[ir] = ra[1];
            ir = ir - 1;
            if (ir == 1)
			{
              ra[1] = rra;
              return;
            }
        }
        i = l;
        j = l + l;
        while (j <= ir)
		{
            if (j < ir)
			{
                if (ra[j] < ra[j + 1])
				{
					j+=1;
				}
            }
            if (rra < ra[j])
			{
                ra[i] = ra[j];
                i = j;
                j += j;
			}
            else
			{
                j = ir + 1;
            }
        }
        ra[i] = rra;
    }while(1);
}

double rofunc(double b)
{
    int n1,j;
	double nml,nmh;
	double d;
	n1 = ndatat + 1;
    nml = (double)n1 / 2;
    nmh = n1 - nml;
    for (j = 1; j<=ndatat; j++)
	{
        arr[j] = yt[j] - b * xt[j];
    }
    sort(ndatat, arr);
    aa = 0.5 * (arr[(int)nml] + arr[(int)nmh]);
    double sum = 0.0;
    abdevt = 0.0;
    for (j = 1; j<=ndatat; j++)
	{
        d = yt[j] - (b * xt[j] + aa);
        abdevt = abdevt + fabs(d);
        sum = sum + xt[j] * sgn(d);
    }
    return sum;
}

void medfit(double x[], double y[], int ndata, double& a, double& b, double& abdev)
{
	int j;
    double bb,del,chisq,f,sx = 0.0;
    double sigb,b1,b2,f1,f2,sy = 0.0;
    double sxy = 0.0;
    double sxx = 0.0;
    for (j = 1; j<=ndata; j++)
	{
        xt[j] = x[j];
        yt[j] = y[j];
        sx = sx + x[j];
        sy = sy + y[j];
        sxy = sxy + x[j] * y[j];
        sxx = sxx + x[j] * x[j];
    }
    ndatat = ndata;
    del = ndata * sxx - sx * sx;
    aa = (sxx * sy - sx * sxy) / del;
    bb = (ndata * sxy - sx * sy) / del;
    chisq = 0.0;
    for (j = 1; j<=ndata; j++)
	{
        chisq = chisq + pow((y[j] - (aa + bb * x[j])) , 2);
    }
    sigb = sqrt(chisq / del);
    b1 = bb;
    f1 = rofunc(b1);
    b2 = bb + fabs(3 * sigb) * sgn(f1);
    f2 = rofunc(b2);
    while (f1 * f2 > 0.0)
	{
        bb = 2.0 * b2 - b1;
        b1 = b2;
        f1 = f2;
        b2 = bb;
        f2 = rofunc(b2);
	}
    sigb = 0.01 * sigb;
    while (fabs(b2 - b1) > sigb)
	{
        bb = 0.5 * (b1 + b2);
        if (bb == b1 || bb == b2)
		{
			break;
		}
        f = rofunc(bb);
        if (f * f1 >= 0.0)
		{
            f1 = f;
            b1 = bb;
		}
        else
		{
            f2 = f;
            b2 = bb;
        }
    }
    a = aa;
    b = bb;
    abdev = abdevt / ndata;
}

void main()
{
    //program d9r10
    //driver for routine medfit
    int i,npt = 100;
    double a,b,siga,sigb,chi2,q,abdev,spread = 0.1;
    double x[101], y[101], sig[101];
    long idum = -1984;
    for (i = 1; i<=npt; i++)
	{
        x[i] = 0.1 * i;
        y[i] = -2.0 * x[i] + 1.0 + spread * gasdev(idum);
        sig[i] = spread;
    }
    int mwt = 1;
    fit(x, y, npt, sig, mwt, a, b, siga, sigb, chi2, q);
	cout.setf(ios::fixed|ios::left);
	cout.precision(4);
    cout<<"According to routine fit the result is:"<<endl;
    cout<<endl;
    cout<<"a = ";
	cout.width(12);
	cout<<a;
	cout<<"Uncertainty: ";
	cout.width(12);
    cout<<siga<<endl;
	cout<<endl;
    cout<<"b = ";
	cout.width(12);
	cout<<b;
	cout<<"Uncertainty: ";
    cout<<sigb<<endl;
    cout<<endl;
    cout<<"Chi-squared: ";
	cout.width(12);
	cout<<chi2<<" for   ";
	cout.width(5);
	cout<<npt<<"points";
    cout<<endl;cout<<endl;
    cout<<"Goodness-of-fit: ";
	cout.width(12);
	cout<<q;
    cout<<endl;cout<<endl;
    cout<<"According to routine medfit the result is:";
    cout<<endl;cout<<endl;
    medfit(x, y, npt, a, b, abdev);
    cout<<"a = ";
	cout.width(12);
	cout<<a;
    cout<<"b = ";
	cout.width(12);
	cout<<b;
    cout<<endl;cout<<endl;
    cout<<"Absolute deviation (per data point): ";
	cout.width(12);
	cout<<abdev;
    cout<<endl;cout<<endl;
    cout<<"(note: Gaussian spread is  ";
	cout<<spread;
	cout<<")";
	cout<<endl;cout<<endl;
}