### Polynomial Roots

#### Introduction

This Java program implements the Durand-Kerner-Weierstrass method to find the roots of a univariate polynomial with complex coefficients (listing):

Durand-Kerner method
The method of Weierstrass
Univariate Polynomial Root-Finding with Lower Computational Precision and Higher Convergence Rates

For better numerical stability, the implementation uses Horner's scheme to evaluate the polynomial. A similar implentation in Ada may be found here. Sample output is shown below. The implementation uses the JScience library for complex arithmetic.

#### Implementation

```import org.jscience.mathematics.function.Polynomial;
import org.jscience.mathematics.function.Term;
import org.jscience.mathematics.function.Variable;
import org.jscience.mathematics.number.Complex;

/**
* An implementation of the Durand-Kerner-Weierstrass method to
* determine the roots of a complex univariate polynomial, as described
* <a href="http://en.wikipedia.org/wiki/Durand-Kerner_method">here</a>.
*
* @author John B. Matthews; distribution per LGPL.
*/
final class ComplexPolynomial {

private static final int MAX_COUNT = 999;
private static final boolean DEBUG = false;
private static double epsilon = 1E-15;

private ComplexPolynomial() {}

/**
* Create a complex polynomial from a number of Complex values.
* The array should have the highest order coefficient first.
*
* @param  a a variable number of Complex coefficients.
* @return a Polynomial having those Complex coefficients.
*/
public static Polynomial<Complex> create(Complex... a) {
Variable<Complex> x = new Variable.Local<Complex>("x");
Polynomial<Complex> px = Polynomial.valueOf(Complex.ZERO, x);
for (int i = 0, e = a.length - 1; i < a.length; i++, e--) {
px = px.plus(Polynomial.valueOf(a[i], Term.valueOf(x, e)));
}
return px;
}

/**
* Create a complex array from an array of double.
* The array should have the highest order coefficient first.
*
* @param  a a variable number of double coefficients.
* @return an array of the corresponding Complex coefficients.
*/
public static Complex[] complexArray(double... a) {
Complex[] ca = new Complex[a.length];
for (int i = 0; i < a.length; i++) {
ca[i] = Complex.valueOf(a[i], 0);
}
return ca;
}

/**
* This implementation uses the Durand-Kerner-Weierstrass method
* to find the roots of a polynomial with complex coefficients.
* The method requires a monic polynomial; some error may occur
* when dividing by the leading coefficient.
* The array should have the highest order coefficient first.
*
* @param  ca a variable number of Complex polynomial coefficients.
* @return an array of the Complex roots of the polynomial.
*/
public static Complex[] roots(Complex[] ca) {
Complex[] a0 = new Complex[ca.length - 1];
Complex[] a1 = new Complex[ca.length - 1];

// Divide by leading coefficient if not monic
if (!ca[0].equals(Complex.ONE)) {
for (int i = 0; i < ca.length; i++) {
}
}

// Initialize a0
Complex result = Complex.valueOf(0.4, 0.9);
a0[0] = Complex.ONE;
for (int i = 1; i < a0.length; i++) {
a0[i] = a0[i - 1].times(result);
}

// Iterate
int count = 0;
while (true) {
for (int i = 0; i < a0.length; i++) {
result = Complex.ONE;
for (int j = 0; j < a0.length; j++) {
if (i != j) {
result = a0[i].minus(a0[j]).times(result);
}
}
a1[i] = a0[i].minus(ComplexPolynomial.
eval(ca, a0[i]).divide(result));
}
count++;
if (count > MAX_COUNT || done(a0, a1)) break;
System.arraycopy(a1, 0, a0, 0, a1.length); // a0 := a1
}
if (DEBUG) {
System.out.println("Iterations: " + count);
for (Complex c : a0) System.out.println(c);
}
return a1;
}

// Determine if the components of two arrays are unchanging
private static boolean done(Complex[] a, Complex[] b) {
boolean unchanged = true;
Complex delta;
for (int i = 0; i < a.length; i++) {
delta = a[i].minus(b[i]);
unchanged &= Math.abs(delta.getReal()) < epsilon
&& Math.abs(delta.getImaginary()) < epsilon;
}
return unchanged;
}

/**
* Evaluate the polynomial at x using
* <a href="http://en.wikipedia.org/wiki/Horner_scheme">Horner's scheme</a>.
* The array should have the highest order coefficient first.
*
* @param  ca an array of Complex polynomial coefficients.
* @param  x the function argument.
* @return the Complex value of the function at x.
*/
public static Complex eval(Complex[] ca, Complex x) {
Complex result = ca[0];
for (int i = 1; i < ca.length; i++) {
result = result.times(x).plus(ca[i]);
}
return result;
}

/** Return the nominal tolerance value. */
public static double getEpsilon() { return epsilon; }

/** Set the nominal tolerance value. */
public static void setEpsilon(double e) { epsilon = e; }
}
```

#### Sample test output

The test driver produces the following output using `double` complex arithmetic. On the test platform, `double` offers more than 15 decimal digits of precision, and a corresponding value of epsilon is used. For each polynomial P and root r, the maximal value of abs(P(r) – 0) is shown. The first few entries are easily factored polynomials with integer coefficients. These are followed by a series of all one polynomials. Finally, there are a few examples with complex coefficients. On PowerPC (G4), the largest useful order appears to be about 38; on Intel (Pentium 4) and AMD (Athlon), extended precision may give correct results up to order 48. Beyond that, round-off errors cause the roots to converge to zero.

```Poly: [1.0 + 0.0i]x^2 + [-3.0 + 0.0i]x + [2.0 + 0.0i]
Real:  1.00000000000000E0
Real:  2.00000000000000E0
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^2 + [-2.0 + 0.0i]x + [1.0 + 0.0i]
Real:  1.00000000000000E0
Real:  1.00000000000000E0
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^2 + [4.0 + 0.0i]
Comp:  0.00000000000000E0 +- 2.00000000000000E0i
Error: < 1.0E-15

Poly: [4.0 + 0.0i]x^2 + [1.0 + 0.0i]
Comp:  0.00000000000000E0 +- 5.00000000000000E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^3 + [-3.0 + 0.0i]x^2 + [2.0 + 0.0i]x
Real:  0.00000000000000E0
Real:  1.00000000000000E0
Real:  2.00000000000000E0
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^4 + [-13.0 + 0.0i]x^2 + [36.0 + 0.0i]
Real: -3.00000000000000E0
Real: -2.00000000000000E0
Real:  2.00000000000000E0
Real:  3.00000000000000E0
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^5 + [-13.0 + 0.0i]x^3 + [36.0 + 0.0i]x
Real: -3.00000000000000E0
Real: -2.00000000000000E0
Real:  0.00000000000000E0
Real:  2.00000000000000E0
Real:  3.00000000000000E0
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp:  0.00000000000000E0 +- 1.00000000000000E0i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp:  3.09016994374947E-1 +- 9.51056516295154E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp:  5.00000000000000E-1 +- 8.66025403784439E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.00968867902419E-1 +- 4.33883739117558E-1i
Comp: -2.22520933956314E-1 +- 9.74927912181824E-1i
Comp:  6.23489801858734E-1 +- 7.81831482468030E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -7.07106781186548E-1 +- 7.07106781186548E-1i
Comp:  0.00000000000000E0 +- 1.00000000000000E0i
Comp:  7.07106781186548E-1 +- 7.07106781186548E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.39692620785908E-1 +- 3.42020143325669E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp:  1.73648177666930E-1 +- 9.84807753012208E-1i
Comp:  7.66044443118978E-1 +- 6.42787609686539E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp: -3.09016994374947E-1 +- 9.51056516295154E-1i
Comp:  3.09016994374947E-1 +- 9.51056516295154E-1i
Comp:  8.09016994374948E-1 +- 5.87785252292473E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.59492973614497E-1 +- 2.81732556841430E-1i
Comp: -6.54860733945285E-1 +- 7.55749574354258E-1i
Comp: -1.42314838273285E-1 +- 9.89821441880933E-1i
Comp:  4.15415013001886E-1 +- 9.09631995354518E-1i
Comp:  8.41253532831181E-1 +- 5.40640817455598E-1i
Error:  1.01904927555717E-15

Poly: [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -8.66025403784439E-1 +- 5.00000000000000E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp:  0.00000000000000E0 +- 1.00000000000000E0i
Comp:  5.00000000000000E-1 +- 8.66025403784439E-1i
Comp:  8.66025403784439E-1 +- 5.00000000000000E-1i
Error:  1.14304456355485E-15

Poly: [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.70941817426052E-1 +- 2.39315664287558E-1i
Comp: -7.48510748171101E-1 +- 6.63122658240795E-1i
Comp: -3.54604887042536E-1 +- 9.35016242685415E-1i
Comp:  1.20536680255323E-1 +- 9.92708874098054E-1i
Comp:  5.68064746731156E-1 +- 8.22983865893656E-1i
Comp:  8.85456025653210E-1 +- 4.64723172043769E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.00968867902419E-1 +- 4.33883739117558E-1i
Comp: -6.23489801858734E-1 +- 7.81831482468030E-1i
Comp: -2.22520933956314E-1 +- 9.74927912181824E-1i
Comp:  2.22520933956314E-1 +- 9.74927912181824E-1i
Comp:  6.23489801858734E-1 +- 7.81831482468030E-1i
Comp:  9.00968867902419E-1 +- 4.33883739117558E-1i
Error: < 1.0E-15

Poly: [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.78147600733806E-1 +- 2.07911690817759E-1i
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: -1.04528463267653E-1 +- 9.94521895368273E-1i
Comp:  3.09016994374947E-1 +- 9.51056516295154E-1i
Comp:  6.69130606358858E-1 +- 7.43144825477394E-1i
Comp:  9.13545457642601E-1 +- 4.06736643075800E-1i
Error:  1.25607396694702E-15

Poly: [1.0 + 0.0i]x^15 + [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.23879532511287E-1 +- 3.82683432365090E-1i
Comp: -7.07106781186548E-1 +- 7.07106781186548E-1i
Comp: -3.82683432365090E-1 +- 9.23879532511287E-1i
Comp:  0.00000000000000E0 +- 1.00000000000000E0i
Comp:  3.82683432365090E-1 +- 9.23879532511287E-1i
Comp:  7.07106781186548E-1 +- 7.07106781186548E-1i
Comp:  9.23879532511287E-1 +- 3.82683432365090E-1i
Error:  2.71039981919101E-15

Poly: [1.0 + 0.0i]x^37 + [1.0 + 0.0i]x^36 + [1.0 + 0.0i]x^35 + [1.0 + 0.0i]x^34 + [1.0 + 0.0i]x^33 + [1.0 + 0.0i]x^32 + [1.0 + 0.0i]x^31 + [1.0 + 0.0i]x^30 + [1.0 + 0.0i]x^29 + [1.0 + 0.0i]x^28 + [1.0 + 0.0i]x^27 + [1.0 + 0.0i]x^26 + [1.0 + 0.0i]x^25 + [1.0 + 0.0i]x^24 + [1.0 + 0.0i]x^23 + [1.0 + 0.0i]x^22 + [1.0 + 0.0i]x^21 + [1.0 + 0.0i]x^20 + [1.0 + 0.0i]x^19 + [1.0 + 0.0i]x^18 + [1.0 + 0.0i]x^17 + [1.0 + 0.0i]x^16 + [1.0 + 0.0i]x^15 + [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.86361303402722E-1 +- 1.64594590280734E-1i
Comp: -9.45817241700635E-1 +- 3.24699469204684E-1i
Comp: -8.79473751206489E-1 +- 4.75947393037074E-1i
Comp: -7.89140509396394E-1 +- 6.14212712689668E-1i
Comp: -6.77281571625741E-1 +- 7.35723910673132E-1i
Comp: -5.46948158122427E-1 +- 8.37166478262529E-1i
Comp: -4.01695424652969E-1 +- 9.15773326655057E-1i
Comp: -2.45485487140799E-1 +- 9.69400265939330E-1i
Comp: -8.25793454723323E-2 +- 9.96584493006670E-1i
Comp:  8.25793454723323E-2 +- 9.96584493006670E-1i
Comp:  2.45485487140799E-1 +- 9.69400265939330E-1i
Comp:  4.01695424652970E-1 +- 9.15773326655058E-1i
Comp:  5.46948158122427E-1 +- 8.37166478262528E-1i
Comp:  6.77281571625741E-1 +- 7.35723910673132E-1i
Comp:  7.89140509396394E-1 +- 6.14212712689668E-1i
Comp:  8.79473751206489E-1 +- 4.75947393037074E-1i
Comp:  9.45817241700635E-1 +- 3.24699469204683E-1i
Comp:  9.86361303402722E-1 +- 1.64594590280734E-1i
Error:  1.36549533872468E-14

Poly: [1.0 + 3.0i]x^2 + [2.0 + 2.0i]x + [3.0 + 1.0i]
Comp: -8.00000000000000E-1 - 6.00000000000000E-1i
Comp:  0.00000000000000E0 + 1.00000000000000E0i
Error: < 1.0E-15

Poly: [1.0 + 1.0i]x^3 + [2.0 + 1.0i]x^2 + [3.0 + 1.0i]x + [4.0 + 1.0i]
Comp: -1.40135938330275E0 + 2.88269653138075E-1i
Comp: -2.84985631785343E-1 - 1.30378640290474E0i
Comp:  1.86345015088091E-1 + 1.51551674976666E0i
Error: < 1.0E-15

Poly: [1.0 + 1.0i]x^3 + [2.0 + 2.0i]x^2 + [3.0 + 3.0i]x + [4.0 + 4.0i]
Real: -1.65062919143939E0
Comp: -1.74685404280306E-1 +- 1.54686888723140E0i
Error: < 1.0E-15
```