 andy0130tw
7/9/2016 - 3:38 PM

## Newton method to numerically solve all roots of a polynomial of degree of 5.

Newton method to numerically solve all roots of a polynomial of degree of 5.

``````#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), fabs(poly->coeff));
}

// 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;
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);
}
``````