d13r3.cpp

使用VC++编写的大量数学算法的源代码

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

void sort(int n, double ra[])
{
	int i,j,l,ir;
	double rra;
	l = 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 = j + 1;
				}
			}
			if (rra < ra[j])
			{
				ra[i] = ra[j];
				i = j;
				j = j + j;
			}
			else
			{
				j = ir + 1;
			}
		}
		ra[i] = rra;
	}while(1);
}

void mdian1(double x[], int n, double& xmed)
{
	int n2;
    sort(n, x);
    n2 = n / 2;
    if (2 * n2 == n)
	{
        xmed = 0.5 * (x[n2] + x[n2 + 1]);
	}
    else
	{
        xmed = x[n2 + 1];
    }
}

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 mdian2(double x[], int n, double& xmed)
{
    double ap,am,big = 1e+30;
    double dum,bbb,aa,afac = 1.5;
    double amp = 1.5;
    double a = 0.5 * (x[1] + x[n]);
    double eps = fabs(x[n] - x[1]);
    ap = big;
    am = -big;
r1: double xx,sum = 0.0;
    double sumx = 0.0;
    int j,np = 0;
    int nm = 0;
    double xp = big;
    double xm = -big;
    for (j = 1; j<=n; j++)
	{
        xx = x[j];
        if (xx != a)
		{
            if (xx > a)
			{
                np = np + 1;
                if (xx < xp)
				{
					xp = xx;
				}
			}
            else
			{
				if (xx < a)
				{
					nm = nm + 1;
					if (xx > xm)
					{
						xm = xx;
					}
				}
            }
            dum = 1.0 / (eps + fabs(xx - a));
            sum = sum + dum;
            sumx = sumx + xx * dum;
        }
    }
    if (np - nm >= 2)
	{
        am = a;
        if (sumx / sum - a < 0)
		{
            bbb = 0;
		}
        else
		{
            bbb = sumx / sum - a;
        }
        aa = xp + bbb * amp;
        if (aa > ap)
		{
			aa = 0.5 * (a + ap);
		}
        eps = afac * fabs(aa - a);
        a = aa;
        goto r1;
	}
    else
	{
		if (nm - np >= 2)
		{
			ap = a;
			if (sumx / sum - a < 0)
			{
				bbb = sumx / sum - a;
			}
			else
			{
				bbb = 0;
			}
			aa = xm + bbb * amp;
			if (aa < am)
			{
				aa = 0.5 * (a + am);
			}
			eps = afac * fabs(aa - a);
			a = aa;
			goto r1;
		}
		else
		{
			if ((n % 2) == 0)
			{
				if (np == nm)
				{
					xmed = 0.5 * (xp + xm);
				}
				else
				{
					if (np > nm)
					{
						xmed = 0.5 * (a + xp);
					}
					else
					{
						xmed = 0.5 * (xm + a);
					}
				}
			}
			else
			{
				if (np == nm)
				{
					xmed = a;
				}
				else
				{
					if (np > nm)
					{
						xmed = xp;
					}
					else
					{
						xmed = xm;
					}
				}
			}
		}
	}
}

void main()
{
    //program d13r3
    //driver for routine mdian2
    int i,npts = 50;
	double xmed;
    double data[51];
    long idum = -5;
    for (i = 1; i<=npts; i++)
	{
        data[i] = gasdev(idum);
    }
    mdian2(data, npts, xmed);
	cout.setf(ios::fixed|ios::left);
	cout.precision(6);
    cout<<endl;
    cout<<"Gaussian distrib., zero mean, unit variance"<<endl;
    cout<<"Median according to mdian2 is   ";
	cout.width(12);
    cout<<xmed<<endl;
    mdian1(data, npts, xmed);
    cout<<"Median according to mdian1 is   ";
	cout.width(12);
    cout<<xmed<<endl;
}