c_cpp 牛顿法数值求解5次多项式的所有根。

Posted 2021-05-14

#include<math.h>
#include<string.h>
#include<stdio.h>
#include<stdlib.h>

typedef struct {int n_order; double *coeff; } ebao;

#define FP_OUTPUT_FORMAT "%.5lf"
#define SIGN_CHAR(x) (x >= 0 ? '+' : '-')

// allocate the poly and specify its degree
ebao* newpoly(int n) {
    ebao* poly = (ebao*) malloc(sizeof(ebao));
    poly->n_order = n;
    poly->coeff = (double*) malloc(sizeof(double) * (n + 1));
    return poly;
}

// deallocate the poly
void deletepoly(ebao* poly) {
    free(poly->coeff);
    poly->coeff = NULL;
    free(poly);
}

// make a copy the poly
ebao* copypoly(ebao* poly) {
    int n = poly->n_order;
    ebao* poly2 = newpoly(n);
    memcpy(poly2->coeff, poly->coeff, sizeof(double) * (n + 1));
    return poly2;
}

// read the coefficients (in ascending powers) from stdin
ebao* inputpoly(int deg) {
    ebao* poly = newpoly(deg);
    for (int i = 0; i <= deg; i++)
        scanf("%lf", poly->coeff + i);
    return poly;
}

// pretty-print the poly
void printpoly(ebao* poly) {
    int n = poly->n_order;
    if (n > 0)
        printf(FP_OUTPUT_FORMAT " x^%d", poly->coeff[n], n);
    for (int i = n - 1; i > 0; i--) {
        printf(" %c " FP_OUTPUT_FORMAT " x^%d",
            SIGN_CHAR(poly->coeff[i]), fabs(poly->coeff[i]), i);
    }
    printf(" %c " FP_OUTPUT_FORMAT, SIGN_CHAR(poly->coeff[0]), fabs(poly->coeff[0]));
}

// differentiate the poly
ebao* diffpoly(ebao* poly) {
    int n = poly->n_order;
    ebao* diff;
    if (poly->n_order == 0) {
        diff = newpoly(0);
        diff-> coeff[0] = 0;
        return diff;
    } else {
        diff = newpoly(n - 1);
        for (int i = 1; i <= n; i++) {
            diff->coeff[i - 1] = i * poly->coeff[i];
        }
    }
    return diff;
}

// substitute the x and yield the result
double evalpoly(double x, ebao* poly) {
    int n = poly->n_order;
    double val = poly->coeff[n];
    for (int i = n - 1; i >= 0; i--) {
        val = val * x + poly->coeff[i];
    }
    return val;
}

// divide the polynomial by (x - k); the remainder is dropped
ebao* divpoly(double k, ebao* poly) {
    int n = poly->n_order;
    ebao* quot = newpoly(n - 1);
    double val = 0;
    quot->coeff[n - 1] = poly->coeff[n];
    for (int i = n - 1; i > 0; i--) {
        quot->coeff[i - 1] = poly->coeff[i] + quot->coeff[i] * k;
    }
    return quot;
}

// get the zero of poly by Newton's method
double getzeropoly(double x, ebao* poly, double eps, int max_iter) {
    ebao* diff = diffpoly(poly);
    double xbest = x, mbest = INFINITY;
    while (max_iter) {
        double evalf = evalpoly(x, poly);
        double evald = evalpoly(x, diff);
        double evalm = evalf / evald;
        // printf("    x = %lf\tf = %lf\tf' = %lf\tf/f' = %lf\n", x, evalf, evald, evalm);

        x -= evalm;

        if (fabs(evalm) <= mbest) {
            mbest = fabs(evalm);
            xbest = x;
        }

        if (isnan(x) || mbest <= eps)
            break;
        --max_iter;
    }
    deletepoly(diff);
    if (isnan(x)) return NAN;
    return xbest;
}

void test() {
    // we manually pass predefined coefficients
    ebao* p5 = (ebao*) malloc(sizeof(ebao));
    p5->n_order = 5;
    double coeff[] = {1, 5, 10, 10, 5, 1};
    p5->coeff = coeff;
    ebao* dp5 = diffpoly(p5);
    ebao* p5dm1 = divpoly(-1, p5);

    free(p5);
    deletepoly(dp5);
    deletepoly(p5dm1);
}

int main() {
    printf("Enter coeffs: ");
    ebao* p5 = inputpoly(5);

    // test();

    printf("Input: ");
    printpoly(p5);
    puts("");

    printf("Now solving the roots\n");

    ebao* p5dup = copypoly(p5);

    for (int i = 0; i < 5; i++) {
        double sol;
        double try_init[] = {0, 1, -1, evalpoly(1, p5dup), evalpoly(-1, p5dup)};

        for (int i = 0; i < sizeof(try_init) / sizeof(double); i++) {
            printf("Trying %lf\n", try_init[i]);
            sol = getzeropoly(try_init[i], p5dup, 1e-7, 1000);
            if (!isnan(sol)) break;
        }
        if (isnan(sol)) break;

        printf("  ROOT #%d: " FP_OUTPUT_FORMAT "\n", i + 1, sol);

        ebao* fact = divpoly(sol, p5dup);
        if (fact->n_order == 0) {
            deletepoly(fact);
            break;
        }

        printf("  Factored form: (x %c " FP_OUTPUT_FORMAT ")", SIGN_CHAR(-sol), fabs(sol));
        printf("(");
        printpoly(fact);
        printf(")\n");

        deletepoly(p5dup);
        p5dup = fact;
    }

    if (p5dup->n_order == 1) {
        printf("ALL REAL ROOTS ARE SOLVED\n");
    } else {
        printf("CANNOT SOLVE ALL; MAY HAVE %d IMAGARY ROOT(S)\n", p5dup->n_order);
    }

    deletepoly(p5dup);
    deletepoly(p5);
}