/** * Computes the polynomials output given the variable value Can handle infinite numbers * * @return Output */ public double evaluate( double variable ) { if( size == 0 ) { return 0; } else if( size == 1 ) { return c[0]; } else if( Double.isInfinite(variable)) { // Only the largest power with a non-zero coefficient needs to be evaluated int degree = computeDegree(); if( degree%2 == 0 ) variable = Double.POSITIVE_INFINITY; if( c[degree] < 0 ) variable = variable == Double.POSITIVE_INFINITY ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; return variable; } // Evaluate using Horner's Method. // This formulation produces more accurate results than the straight forward one and has fewer // multiplications double total = c[size-1]; for( int i = size-2; i >= 0; i-- ) { total = c[i] + total*variable; } return total; }
Polynomial deriv = new Polynomial(poly.size()); derivative(poly,deriv); double v = poly.evaluate(root); double d = deriv.evaluate(root);
/** * Wraps the polynomial around the array:<br> * f(x) = c[0] + c[1]*x + ... + c[n]*x<sup>n-1</sup> * @param coefficients Polynomial coefficients * @return new instance of a polyonimial which is identical to the input array */ public static Polynomial wrap( double ...coefficients ) { Polynomial p = new Polynomial(coefficients.length); p.setTo(coefficients,coefficients.length); return p; }
/** * Multiplies the two polynomials together. * @param a Polynomial * @param b Polynomial * @param result Optional storage parameter for the results. Must be have enough coefficients to store the results. * If null a new instance is declared. * @return Results of the multiplication */ public static Polynomial multiply( Polynomial a , Polynomial b , Polynomial result ) { int N = Math.max(0,a.size() + b.size() - 1); if( result == null ) { result = new Polynomial(N); } else { if( result.size < N ) throw new IllegalArgumentException("Unexpected length of 'result'"); result.zero(); } for( int i = 0; i < a.size; i++ ) { double coef = a.c[i]; int index = i; for( int j = 0; j < b.size; j++ ) { result.c[index++] += coef*b.c[j]; } } return result; }
/** * Compute the sturm sequence using a straight forward method */ private List<Double> computeSturm( Polynomial poly , double x ) { Polynomial d = new Polynomial(poly.size); PolynomialOps.derivative(poly, d); List<Double> found = new ArrayList<Double>(); found.add( poly.evaluate(x)); found.add( d.evaluate(x)); Polynomial q = new Polynomial(poly.size); Polynomial r = new Polynomial(poly.size); Polynomial p1 = new Polynomial(poly.size); Polynomial p2 = new Polynomial(poly.size); p1.setTo(poly); p2.setTo(d); do { PolynomialOps.divide(p1,p2,q,r); for( int i = 0; i < r.size; i++ ) { r.c[i] = -r.c[i]; } found.add(r.evaluate(x)); p1.setTo(p2); p2.setTo(r); } while( r.computeDegree() > 0 ); return found; }
@Test public void evaluate() { assertEquals(5,Polynomial.wrap(1,2).evaluate(2),1e-8); assertEquals(5,Polynomial.wrap(1,2,0).evaluate(2),1e-8); assertEquals(10,Polynomial.wrap(0,1,2).evaluate(2),1e-8); assertEquals(0,Polynomial.wrap(0).evaluate(2),1e-8); assertEquals(0, Polynomial.wrap().evaluate(Double.POSITIVE_INFINITY), 1e-8); assertEquals(2, Polynomial.wrap(2).evaluate(Double.POSITIVE_INFINITY), 1e-8); assertEquals(Double.POSITIVE_INFINITY,Polynomial.wrap(2,3).evaluate(Double.POSITIVE_INFINITY),1e-8); assertEquals(Double.POSITIVE_INFINITY,Polynomial.wrap(2,3,5).evaluate(Double.POSITIVE_INFINITY),1e-8); assertEquals(Double.NEGATIVE_INFINITY,Polynomial.wrap(-2,-3).evaluate(Double.POSITIVE_INFINITY),1e-8); assertEquals(Double.NEGATIVE_INFINITY,Polynomial.wrap(-2,-3,-5).evaluate(Double.POSITIVE_INFINITY),1e-8); assertEquals(Double.NEGATIVE_INFINITY,Polynomial.wrap(2,3).evaluate(Double.NEGATIVE_INFINITY),1e-8); assertEquals(Double.POSITIVE_INFINITY,Polynomial.wrap(2,3,5).evaluate(Double.NEGATIVE_INFINITY),1e-8); assertEquals(Double.POSITIVE_INFINITY,Polynomial.wrap(-2,-3).evaluate(Double.NEGATIVE_INFINITY),1e-8); assertEquals(Double.NEGATIVE_INFINITY,Polynomial.wrap(-2,-3,-5).evaluate(Double.NEGATIVE_INFINITY),1e-8); assertEquals(Double.POSITIVE_INFINITY,Polynomial.wrap(2,3,0,0).evaluate(Double.POSITIVE_INFINITY),1e-8); assertEquals(Double.NEGATIVE_INFINITY,Polynomial.wrap(2,3,0,0).evaluate(Double.NEGATIVE_INFINITY),1e-8); // Make sure it's using size and not length Polynomial p = new Polynomial(10); p.size = 2; p.c[0] = 2; p.c[1] = 3; p.c[4] = 7; assertEquals(8,p.evaluate(2),1e-8); }
sequence[0].setTo(poly); sequence[1].setTo(previous); sequenceLength = 2; return; int degree = result.computeDegree(); if( degree <= 0 ) { if( degree < 0 ) { sequence[i].setTo(next); sequenceLength = i+1; } else { sequence[i+1].setTo(result); PolynomialOps.divide(next, result, sequence[i], previous); sequenceLength = i+2;
@Test public void computeDegree() { assertEquals(-1,Polynomial.wrap().computeDegree()); assertEquals(-1,Polynomial.wrap(0).computeDegree()); assertEquals(-1,Polynomial.wrap(0,0).computeDegree()); assertEquals(0,Polynomial.wrap(1).computeDegree()); assertEquals(0,Polynomial.wrap(1,0).computeDegree()); assertEquals(0,Polynomial.wrap(1,0,0).computeDegree()); assertEquals(2,Polynomial.wrap(0,1,2).computeDegree()); assertEquals(3,Polynomial.wrap(0,1,2,1e-15).computeDegree()); assertEquals(2,Polynomial.wrap(-0,-1,-2,-0).computeDegree()); }
@Test public void rootsSmallReal() { for( int numCoef = 2; numCoef < 6; numCoef++ ) { Polynomial poly = new Polynomial(numCoef); for( int trial = 0; trial < 20; trial++ ) { for( int i = 0; i < numCoef; i++ ) { poly.c[i] = 10*(rand.nextDouble()-0.5); } List<Double> roots = computeRealRoots(poly); // the root of the polynomial should be zero. higher degree polynomials have // more stability problems for( double d : roots ) { assertEquals(0,poly.evaluate(d),1e-8); } int expectedRoots = PolynomialOps.countRealRoots(poly); assertTrue(roots.size() == expectedRoots); } } }
@Test public void derivative() { Polynomial p = Polynomial.wrap(2,3,4,5); p.size = 2; Polynomial d = new Polynomial(10); PolynomialOps.derivative(p,d); assertTrue(d.isIdentical(Polynomial.wrap(3), 1e-8)); p.size = 3; PolynomialOps.derivative(p,d); assertTrue(d.isIdentical(Polynomial.wrap(3, 8), 1e-8)); }
public static Polynomial derivative( Polynomial poly , Polynomial deriv ) { if( deriv == null ) { deriv = new Polynomial(poly.size-1); } else { deriv.size = poly.size - 1; } for( int i = 1; i < poly.size; i++ ) { deriv.c[i-1] = poly.c[i]*i; } return deriv; }
public static void main( String args[] ) { Polynomial a = Polynomial.wrap(5,-2,3,0.5,30.4); Polynomial b = Polynomial.wrap(-0.4,8.4,-2.3); System.out.println("a = "+a); System.out.println("b = "+b); System.out.println("a + b = "+PolynomialOps.add(a,b,null)); System.out.println("a * b = "+PolynomialOps.multiply(a, b, null)); // Declare storage for the quotient and the remainder Polynomial q = new Polynomial(10); Polynomial r = new Polynomial(10); PolynomialOps.divide(a, b, q,r); System.out.println("a / b = ( "+q+" , "+r+" )"); System.out.println("Derivative a = "+PolynomialOps.derivative(a,null)); } }
/** * Finds real and imaginary roots in a polynomial using the companion matrix and * Eigenvalue decomposition. The coefficients order is specified from smallest to largest. * Example, 5 + 6*x + 7*x^2 + 8*x^3 = [5,6,7,8] * * @param coefficients Polynomial coefficients from smallest to largest. * @return The found roots. */ @SuppressWarnings("ToArrayCallWithZeroLengthArrayArgument") public static Complex_F64[] polynomialRootsEVD(double... coefficients) { PolynomialRoots alg = new RootFinderCompanion(); if( !alg.process( Polynomial.wrap(coefficients)) ) throw new IllegalArgumentException("Algorithm failed, was the input bad?"); List<Complex_F64> coefs = alg.getRoots(); return coefs.toArray(new Complex_F64[0]); }
@Test public void truncateZeros() { Polynomial a = Polynomial.wrap(0,1,2,0); Polynomial b = Polynomial.wrap(0,1,2,1e-15); Polynomial c = Polynomial.wrap(0,1,2); Polynomial d = Polynomial.wrap(0,1,2,-1); a.truncateZeros(1e-15); b.truncateZeros(1e-15); c.truncateZeros(1e-15); d.truncateZeros(1e-15); assertTrue(a.isIdentical(Polynomial.wrap(0, 1, 2), 1e-8)); assertTrue(b.isIdentical(Polynomial.wrap(0, 1, 2), 1e-8)); assertTrue(c.isIdentical(Polynomial.wrap(0, 1, 2), 1e-8)); assertTrue(d.isIdentical(Polynomial.wrap(0, 1, 2, -1), 1e-8)); } }
@Test public void identical() { Polynomial a = Polynomial.wrap(0,1,2); assertTrue(a.isIdentical(Polynomial.wrap(0, 1, 2), 1e-8)); assertTrue(a.isIdentical(Polynomial.wrap(0, 1, 2, 0), 1e-8)); assertTrue(a.isIdentical(Polynomial.wrap(0, 1, 2, 1e-10), 1e-8)); assertTrue(Polynomial.wrap(0, 1, 1e-20).isIdentical(Polynomial.wrap(0, 1, 1e-14), 1e-8)); assertTrue(Polynomial.wrap(0, 1, 1e-20, 0).isIdentical(Polynomial.wrap(0, 1, 0, 1e-14), 1e-8)); assertFalse(a.isIdentical(Polynomial.wrap(0, 1, 3), 1e-8)); assertFalse(a.isIdentical(Polynomial.wrap(0, 1, 2 + 1e-5), 1e-8)); assertFalse(a.isIdentical(Polynomial.wrap(0, 1, 2, 3), 1e-8)); }
/** * Several consistency checks on random polynomials to the number of roots within a random interval */ @Test public void rootCountConsistency_Random() { for( int i = 3; i < 20; i++ ) { Polynomial p = new Polynomial(i); for( int trial = 0; trial < 20; trial++ ) { for( int j = 0; j < p.size; j++ ) { p.c[j] = (rand.nextDouble()-0.5)*2; } SturmSequence alg = new SturmSequence(p.size()); double low = (rand.nextDouble() - 0.5)*200; double high = low + rand.nextDouble()*200; double middle = (low+high)/2.0; alg.initialize(p); int every = alg.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); int all = alg.countRealRoots(low,high); int lowN = alg.countRealRoots(low,middle); int highN = alg.countRealRoots(middle,high); assertTrue(all >= 0); assertTrue(lowN >= 0); assertTrue(highN >= 0); assertEquals(all, lowN + highN); assertTrue(all <= every); } } }
/** * Randomly generate polynomials. First divide then multiply them back together to reconstruct * the original. */ @Test public void divide_then_multiply() { for( int i = 2; i < 10; i++ ) { Polynomial n = createRandom(i); for( int j = 1; j < i; j++ ) { Polynomial d = createRandom(j); Polynomial q = new Polynomial(10); Polynomial r = new Polynomial(10); PolynomialOps.divide(n,d,q,r); Polynomial a = PolynomialOps.multiply(d,q,null); Polynomial b = PolynomialOps.add(a,r,null); assertTrue(b.isIdentical(n, 1e-8)); } } }