d10r12.cpp

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

#include <math.h>
#include <iomanip.h>
#include <iostream.h>
#include <process.h>

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 poldiv(double u[], int n, double v[], int nv, double q[], double r[])
{
	int j,k;
    for (j = 1;j<=n;j++)
	{
        r[j] = u[j];
        q[j] = 0.0;
    }
    for( k = n - nv ;k>=0;k--)
	{
        q[k + 1] = r[nv + k] / v[nv];
        for (j = nv + k - 1 ;j>= k + 1; j--)
		{
            r[j] = r[j] - q[k + 1] * v[j - k];
        }
    }
    r[nv] = 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;
}

void main()
{
    //program d10r12
    //driver for routine qroot
	int n,ntry,nroot,nflag,i,j;
	double eps,tiny,p[8],b[11],c[11],aaa,bbb;
    n = 7;
    eps = 0.000001;
    ntry = 10;
    tiny = 0.00001;
    p[1] = 10.0;    p[2] = -18.0;
    p[3] = 25.0;    p[4] = -24.0;
    p[5] = 16.0;    p[6] = -6.0;
    p[7] = 1.0;
    cout<<endl;
    cout<< setw(5)<<"P[x]=x^6-6x^5+16x^4-24x^3+25x^2-18x+10"<<endl;
    cout<<endl;
    cout<< setw(5)<<"Quadratic factors x^2+bx+c"<<endl;
    cout<<endl;
    cout<<  setw(5)<<" factor       b              c"<<endl;
    nroot = 0;
	cout<<setprecision(5)<<setiosflags(ios::fixed);
    for (i = 1; i<=ntry; i++)
	{
        c[i] = 0.5 * i;
        b[i] = -0.5 * i;
        qroot(p, n, b[i], c[i], eps);
        if (nroot == 0)
		{
            cout<<setw(3)<<nroot;
            cout<<setw(15)<<b[i];
            cout<<setw(15)<<c[i]<<endl;
            nroot = 1;
		}
        else
		{
            nflag = 0;
            for (j = 1; j<=nroot; j++)
			{
                aaa = fabs(b[i] - b[j]);
                bbb = fabs(c[i] - c[j]);
                if ((aaa < tiny) && (bbb < tiny))  nflag = 1;
            }
            if (nflag == 0)
			{
                cout<< setw(3)<<nroot;
                cout<< setw(15)<<b[i];
                cout<< setw(15)<<c[i]<<endl;
                nroot = nroot + 1;
            }
		}
    }
}