多项式求根

ACM模版

多项式求根（牛顿法）

/* 
 *  牛顿法解多项式的根
 *  输入:多项式系数c[],多项式度数n,求在[a,b]间的根
 *  输出:根 要求保证[a,b]间有根
 */

double fabs(double x)
{
    return (x < 0) ? -x : x;
}

double f(int m, double c[], double x)
{
    int i;
    double p = c[m];
    for (i = m; i > 0; i--)
    {
        p = p * x + c[i - 1];
    }
    return p;
}

int newton(double x0, double *r, double c[], double cp[], int n, double a, double b, double eps)
{
    int MAX_ITERATION = 1000;
    int i = 1;
    double x1, x2, fp, eps2 = eps / 10.0;
    x1 = x0;
    while (i < MAX_ITERATION)
    {
        x2 = f(n, c, x1);
        fp = f(n - 1, cp, x1);
        if ((fabs(fp) < 0.000000001) && (fabs(x2) > 1.0))
        {
            return 0;
        }
        x2 = x1 - x2 / fp;
        if (fabs(x1 - x2) < eps2)
        {
            if (x2 < a || x2 > b)
            {
                return 0;
            }
            *r = x2;
            return 1;
        }
        x1 = x2;
        i++;
    }
    return 0;
}

double Polynomial_Root(double c[], int n, double a, double b, double eps)
{
    double *cp;
    int i;
    double root;
    cp = (double *)calloc(n, sizeof(double));
    for (i = n - 1; i >= 0; i--)
    {
        cp[i] = (i + 1) * c[i + 1];
    }
    if (a > b)
    {
        root = a;
        a = b;
        b = root;
    }
    if ((!newton(a, &root, c, cp, n, a, b, eps)) && (!newton(b, &root, c, cp, n, a, b, eps)))
    {
        newton((a + b) * 0.5, &root, c, cp, n, a, b, eps);
    }
    free(cp);
    if (fabs(root) < eps)
    {
        return fabs(root);
    }
    else
        return root;
}