int main() {
array double x1[3], p1[4] = {1, -6, -72, -27}; /* x^3-6x^2-72x-27 =0*/
array double x2[4], p2[5] = {1, -12, 0, 25, 116}; /* x^4-12x^3-25x+116 =0 */
array double complex z2[4];
array double complex z3[4], p3[5] = {complex(3,4), complex(4,2),
complex(5,3), complex(2,4), complex(1,5)};
int status, i;
roots(x1, p1);
printf("x1 from roots(x1, p1) = %.3f\n", x1);
for(i=0; i<3; i++)
printf("polyeval(p1, x1[%d]) = %.3f\n", i, polyeval(p1, x1[i]));
/* x^4-12x^3-25x+116 =0 has two complex roots and two real roots */
roots(x2, p2);
printf("x2 from roots(x2, p2) = %.3f\n", x2);
roots(z2, p2);
printf("z2 from roots(z2, p2) = \n%.3f\n", z2);
status = roots(z3, p3);
if(status == 0)
printf("z3 from roots(z3, p3) = \n%.3f\n", z3);
else
printf("roots(z3, p3) failed\n");
for(i=0; i<4; i++)
printf("cpolyeval(p3, z3[%d]) = %.3f\n", i, cpolyeval(p3, z3[i]));
}
x1 from roots(x1, p1) = 12.123 -5.735 -0.388
polyeval(p1, x1[0]) = 0.000
polyeval(p1, x1[1]) = -0.000
polyeval(p1, x1[2]) = 0.000
x2 from roots(x2, p2) = 11.747 2.703 NaN NaN
z2 from roots(z2, p2) =
complex(11.747,0.000) complex(2.703,0.000) complex(-1.225,1.467) complex(-1.225,-1.467)
z3 from roots(z3, p3) =
complex(0.226,1.281) complex(0.431,-0.728) complex(-0.754,-0.708) complex(-0.703,0.554)
cpolyeval(p3, z3[0]) = complex(0.000,-0.000)
cpolyeval(p3, z3[1]) = complex(0.000,-0.000)
cpolyeval(p3, z3[2]) = complex(0.000,-0.000)
cpolyeval(p3, z3[3]) = complex(0.000,-0.000)