// ************************************************************************** // // // // eses eses // // eses eses // // eses eseses esesese eses Embedded Systems Group // // ese ese ese ese ese // // ese eseseses eseseses ese Department of Computer Science // // eses eses ese eses // // eses eseses eseseses eses University of Kaiserslautern // // eses eses // // // // ************************************************************************** // // The module below computes all roots of a given polynomial. The algorithm // // applies Newton-Raphson iteration to find an approximation root[i] of a root// // and divides then the polynomial by the monomial (x-root[i]). To this end, // // the algorithm works with complex numbers to that we know that there are N // // complex roots. As soon as all roots have been computed with precision eps, // // we polish them by another round of Newton-Raphson iterations starting with // // the computed roots, but this time with the original polynomial. // // The algorithm below is not recommended for use in scientific computing; // // to that purpose, it is better to use methods that reduce the problem to the// // computation of the eigenvalues of the companion matrix of the polynomial. // // ************************************************************************** // macro N = 5; // macros on real numbers macro rroot(x) = exp(x,0.5); macro rsign(x) = (x<0 ? -1 : +1); // arithmetic on complex numbers macro cadd(x,y) = (x.0+y.0,x.1+y.1); macro csub(x,y) = (x.0-y.0,x.1-y.1); macro cmul(x,y) = (x.0*y.0-x.1*y.1,x.1*y.0+x.0*y.1); macro cdiv(x,y) = ((x.0*y.0+x.1*y.1)/(y.0*y.0+y.1*y.1), (x.1*y.0-x.0*y.1)/(y.0*y.0+y.1*y.1)); macro sqabs(x) = (x.1*x.1 + x.0*x.0); macro cabs(x) = rroot(x.1*x.1 + x.0*x.0); macro croot(x) = (rsign(x.1) * rroot(0.5 * (cabs(x) + x.0)), rroot(0.5 * (cabs(x) - x.0))); macro Re(x) = x.0; macro Im(x) = x.1; // numeric precision macro eps = 1.0e-10; module PolynomialRootsNewton([N+1]real ?c,[N](real*real) root,event !rdy) { [N+1](real*real) cc,dd; // coefficients as complex numbers (real*real) pp,dp; // value of p and its derivative (real*real) x; // value used in Newton iteration // initialization: convert given coefficients to complex numbers for(i=0..N) { cc[i].0 = c[i]; cc[i].1 = 0.0; } for(i=0..N-1) { // apply Newton-Raphson iteration to find a complex root do { // evaluate p(x) and its derivative at x in O(N) time, which can be // improved to log(N) time but does not matter for small polynomials next(pp) = cc[N]; next(dp) = (0.0,0.0); pause; for(j=0..N-1) { next(dp) = cadd(cmul(dp,root[i]),pp); next(pp) = cadd(cmul(pp,root[i]),cc[N-1-j]); pause; } // to avoid division by zero, we cheat a bit let(dp1 = (dp==(0.0,0.0) ? (eps,eps) : dp)) next(root[i]) = csub(root[i],cdiv(pp,dp1)); next(x) = root[i]; pause; } while(sqabs(csub(x,root[i])) > eps); // deflation step: divide polynomial cc by (x-root[i]) to continue with // the computation of the remaining roots using a polynomial of lower degree dd[N] = (0.0,0.0); for(j=1..N) dd[j-1] = (j==N-i ? cc[N-i] : cadd(cc[j],cmul(root[i],dd[j]))); pause; // now assign cc := dd to proceed with the computation of remaining roots for(j=0..N) cc[j] = dd[j]; } pause; // At this point, we have computed all the N roots. However, due to the deflation // the precision of the computed roots might not be that good. For this reason, // we polish them by another round of Newton-Raphson iterations starting with the // so far obtained roots. for(i=0..N) { cc[i].0 = c[i]; cc[i].1 = 0.0; } for(i=0..N-1) { // apply Newton-Raphson iteration starting with the current approximations do { // evaluate p(x) and its derivative at x in O(N) time, which can be // improved to log(N) time but does not matter for small polynomials next(pp) = cc[N]; next(dp) = (0.0,0.0); pause; for(j=0..N-1) { next(dp) = cadd(cmul(dp,root[i]),pp); next(pp) = cadd(cmul(pp,root[i]),cc[N-1-j]); pause; } // to avoid division by zero, we cheat a bit let(dp1 = (dp==(0.0,0.0) ? (eps*eps,eps*eps) : dp)) next(root[i]) = csub(root[i],cdiv(pp,dp1)); next(x) = root[i]; pause; } while(sqabs(csub(x,root[i])) > eps*eps); } emit(rdy); } drivenby p1 { // note: x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120 = (x-1)(x-2)(x-3)(x-4)(x-5) c = [-120.0,+274.0,-225.0,+85.0,-15.0,1.0]; await(rdy); } drivenby p2 { // note: x^5-2x^4-3x^3+6x^2+2x-4 = (x-2)(x-1)(x+1)(x-sqrt(2))(x+sqrt(2)) c = [-4.0,+2.0,+6.0,-3.0,-2.0,1.0]; await(rdy); } drivenby p3 { // note: x^5+6x^4+14x^3+24x^2+33x+18 = (x+3)(x+2)(x+1)(x^2+3) c = [+18.0,33.0,24.0,14.0,6.0,1.0]; await(rdy); } drivenby p4 { // note: x^5-2x^4+2x^3-4x^2-3x+6 = (x-2)(x-1)(x+1)(x^2+3) c = [+6.0,-3.0,-4.0,+2.0,-2.0,1.0]; await(rdy); }