qroot.cpp

Visual C++ 常用数值算法集 源代码

void value1(double qq[],double rem1[],double d[],double q[])
{
	int i;
	for(i=1; i<=20; i++)
	{
	   qq[i] = 0.0;
	   rem1[i] = 0.0;
	   q[i] = 0.0;
	}
	for(i=1; i<=3; i++)
	{
		d[i] = 0.0;
	}
}

void qroot(double p[8],int n,double& b,double& c,double eps)
{
	int itmax,iter,i;
	double tiny,s,r,sc,rc,sb,rb,div,delb,delc,db,dc;
    itmax = 20;
    tiny = 0.000001;
    double q[21], d[4], rem1[21], qq[21];
    d[3] = 1.0;
    for (iter = 1; iter<=itmax; iter++)
	{
        d[2] = b;
        d[1] = c;
        poldiv(p, n, d, 3, q, rem1);
        s = rem1[1];
        r = rem1[2];
        poldiv(q, n - 1, d, 3, qq, rem1);
        sc = -rem1[1];
        rc = -rem1[2];
        for (i = n - 1; i>=1; i--)
            q[i + 1] = q[i];
        q[1] = 0.0;
        poldiv(q, n, d, 3, qq, rem1);
        sb = -rem1[1];
        rb = -rem1[2];
        div = 1.0/ (sb * rc - sc * rb);
        delb = (r * sc - s * rc) * div;
        delc = (-r * sb + s * rb) * div;
        b = b + delb;
        c = c + delc;
        db = fabs(delb) - eps * fabs(b);
        dc = fabs(delc) - eps * fabs(c);
        if (((db <= 0.0)|| (fabs(b) < tiny)) && ((dc <= 0) || (fabs(c) < tiny)))
		{
            value1( qq, rem1, d, q);
            return;
		}
    }
    cout<< "too many iterations in qroot"<<endl;
}