mdian2.cpp

此程序为VC++常用数值算法这本书附赠的光盘中包含了本书中全部的源代码

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;
					}
				}
			}
		}
	}
}