/** * Computes the values for each function in the sequence iterative starting at the end and working its way * towards the beginning.. */ protected void computeFunctions( double value ) { f[sequenceLength-1] = sequence[sequenceLength-1].c[0]; if( Double.isInfinite(value)) { f[sequenceLength-2] = multiplyInfinity(sequence[sequenceLength-2].evaluate(value),f[sequenceLength-1]); for( int i = sequenceLength-3; i > 0; i-- ) { // no need to consider the remainder since this will have a higher degree f[i] = multiplyInfinity(sequence[i].evaluate(value),f[i+1]); } } else { f[sequenceLength-2] = sequence[sequenceLength-2].evaluate(value)*f[sequenceLength-1]; for( int i = sequenceLength-3; i > 0; i-- ) { f[i] = sequence[i].evaluate(value)*f[i+1] - f[i+2]; } } f[0] = sequence[0].evaluate(value); }
/** * Computes the values for each function in the sequence iterative starting at the end and working its way * towards the beginning.. */ protected void computeFunctions( double value ) { f[sequenceLength-1] = sequence[sequenceLength-1].c[0]; if( Double.isInfinite(value)) { f[sequenceLength-2] = multiplyInfinity(sequence[sequenceLength-2].evaluate(value),f[sequenceLength-1]); for( int i = sequenceLength-3; i > 0; i-- ) { // no need to consider the remainder since this will have a higher degree f[i] = multiplyInfinity(sequence[i].evaluate(value),f[i+1]); } } else { f[sequenceLength-2] = sequence[sequenceLength-2].evaluate(value)*f[sequenceLength-1]; for( int i = sequenceLength-3; i > 0; i-- ) { f[i] = sequence[i].evaluate(value)*f[i+1] - f[i+2]; } } f[0] = sequence[0].evaluate(value); }
@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); }
public static void main( String args[] ) { // Select which algorithm to use PolynomialRoots finder = PolynomialOps.createRootFinder(6, RootFinderType.EVD); // Create an arbitrary 3rd order polynomial // f(x) = 2 + 0.2*x + 5*x^2 + 3*x^3 Polynomial poly = Polynomial.wrap( 2 , 0.2 , 5 , 3 ); // Find the roots if( !finder.process(poly) ) throw new RuntimeException("Failed to find solution!"); // Print the solution List<Complex_F64> roots = finder.getRoots(); System.out.println("Total roots found: "+roots.size()); for( Complex_F64 c : roots ) { if( !c.isReal() ) { System.out.println("root is imaginary: "+c); continue; } double value = poly.evaluate(c.real); System.out.println("Polynomial value at "+c.real+" is "+value); } } }
@Test public void cubicRealRoot() { Polynomial p = Polynomial.wrap(4,-5,2,1.3); double root = PolynomialOps.cubicRealRoot(p.c[0],p.c[1],p.c[2],p.c[3]); double eval = p.evaluate(root); assertEquals(0, eval, 1e-8); } }
/** * 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 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 rootsLargeReal() { for( int numCoef = 10; numCoef < 15; numCoef++ ) { Polynomial poly = new Polynomial(numCoef); for( int trial = 0; trial < 20; trial++ ) { for( int i = 0; i < numCoef; i++ ) { poly.c[i] = 2000*(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 basicTest() { for( int numCoef = 2; numCoef < 6; numCoef++ ) { Polynomial poly = new Polynomial(numCoef); for( int i = 0; i < numCoef; i++ ) { poly.c[i] = 10*(rand.nextDouble()-0.5); } PolynomialRoots alg = new RootFinderCompanion(); assertTrue(alg.process(poly)); List<Complex_F64> roots = alg.getRoots(); int numReal = 0; for( Complex_F64 c : roots ) { if( c.isReal() ) { assertEquals(0,poly.evaluate(c.real),1e-8); numReal++; } } int expectedRoots = PolynomialOps.countRealRoots(poly); assertTrue(numReal == expectedRoots); } } }