/** * 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]); }
/** * 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 add() { Polynomial a = PolynomialOps.add(Polynomial.wrap(0.5,1,1),Polynomial.wrap(0.5,1,2),null); Polynomial b = PolynomialOps.add(Polynomial.wrap(0.5,1),Polynomial.wrap(0.5,1,2),null); Polynomial c = PolynomialOps.add(Polynomial.wrap(0.5,1,1),Polynomial.wrap(0.5,1),null); assertTrue(Polynomial.wrap(1, 2, 3).isIdentical(a, 1e-8)); assertTrue(Polynomial.wrap(1, 2, 2).isIdentical(b, 1e-8)); assertTrue(Polynomial.wrap(1, 2, 1).isIdentical(c, 1e-8)); }
@Test public void multiply() { Polynomial a = PolynomialOps.multiply(Polynomial.wrap(2), Polynomial.wrap(3), null); Polynomial b = PolynomialOps.multiply(Polynomial.wrap(),Polynomial.wrap(3),null); Polynomial c = PolynomialOps.multiply(Polynomial.wrap(),Polynomial.wrap(),null); Polynomial d = PolynomialOps.multiply(Polynomial.wrap(4),Polynomial.wrap(2,3),null); Polynomial e = PolynomialOps.multiply(Polynomial.wrap(2,3),Polynomial.wrap(4),null); Polynomial f = PolynomialOps.multiply(Polynomial.wrap(2,3),Polynomial.wrap(4,8),null); assertTrue(Polynomial.wrap(6).isIdentical(a, 1e-8)); assertTrue(b.size == 0); assertTrue(c.size == 0); assertTrue(Polynomial.wrap(8, 12).isIdentical(d, 1e-8)); assertTrue(Polynomial.wrap(8, 12).isIdentical(e, 1e-8)); assertTrue(Polynomial.wrap(8, 28, 24).isIdentical(f, 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)); }
@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)); } }
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)); } }
@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 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)); }
@Test public void divide() { // numerator and denominator, intentionally add trailing zeros to skew things up Polynomial n = Polynomial.wrap(1,2,3,0); Polynomial d = Polynomial.wrap(-3,1,0); // quotient and remainder Polynomial q = new Polynomial(10); Polynomial r = new Polynomial(10); // expected solutions Polynomial expectedQ = Polynomial.wrap(11,3,0); Polynomial expectedR = Polynomial.wrap(34); PolynomialOps.divide(n,d,q,r); assertEquals(3,q.size); assertEquals(1,r.size); assertTrue(expectedQ.isIdentical(q, 1e-8)); assertTrue(expectedR.isIdentical(r, 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); } } }
/** * Examine a case which was found to cause problems. It would get stuck in a loop indefinitely */ @Test public void checkSpecificPoly02() { Polynomial poly = Polynomial.wrap(-0.3497655753671151,2.9621784756210587,-12.9131723419964320, 8.0038345403612960,31.0841473414231300,-90.1765840283034800,62.2630323196965500, -17.4634213565573200,1.4562328842432635,0.1227493996590678,-0.0133645583045475); FindRealRootsSturm sturm = new FindRealRootsSturm(11,-1,1e-10,20,20); try { sturm.process(poly); } catch( RuntimeException ignore ) {} }
@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); }
@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); } }
@Test public void countRealRoots() { Polynomial poly = Polynomial.wrap(-1,0,3,1); assertEquals(1, countRealRoots(poly, -3, -2)); assertEquals(0, countRealRoots(poly, 2, 3)); assertEquals(3, countRealRoots(poly, -3, 3)); assertEquals(3, countRealRoots(poly, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); // see if it handles small polynomials correctly assertEquals(0, countRealRoots(Polynomial.wrap(2), Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); assertEquals(1, countRealRoots(Polynomial.wrap(2, 3), Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); assertEquals(0, countRealRoots(Polynomial.wrap(2, 3), 0, Double.POSITIVE_INFINITY)); assertEquals(0, countRealRoots(Polynomial.wrap(2, 3), Double.NEGATIVE_INFINITY, -10)); assertEquals(1, countRealRoots(Polynomial.wrap(2, 3), -2, -0.5)); assertEquals(0, countRealRoots(Polynomial.wrap(2, 3, 4), Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); assertEquals(2, countRealRoots(Polynomial.wrap(2, -1, -4), Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); // more difficult case where special code is required for infinite bounds poly = Polynomial.wrap(-1.322309e+02 , 3.713984e+02 , -5.007874e+02 , 3.744386e+02 ,-1.714667e+02 , 4.865014e+01 ,-1.059870e+01 , 1.642273e+00 ,-2.304341e-01,2.112391e-03,-2.273737e-13); assertEquals(2, countRealRoots(poly, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); }
/** * Check sequence against a hand selected sequence */ @Test public void checkSequence() { Polynomial poly = Polynomial.wrap(-1,0,3,1,5,-3); SturmSequence alg = new SturmSequence(10); alg.initialize(poly); alg.computeFunctions(-3); List<Double> expected = computeSturm(poly,-3); assertEquals(expected.size(),alg.sequenceLength); for( int i = 0; i < expected.size(); i++ ) { assertEquals(expected.get(i),alg.f[i],1e-8); } }
@Test public void setTo() { Polynomial a = new Polynomial(10); Polynomial b = Polynomial.wrap(1,2,3,4); a.setTo(b); assertEquals(b.size(),a.size()); for( int i = 0; i < a.size(); i++ ) assertEquals(a.c[i],b.c[i],1e-8); }
/** * Examine a case which was found to cause problems */ @Test public void checkSpecificPoly01() { Polynomial poly = Polynomial.wrap(-41.118263303597175,-120.95384505825373,-417.8477600492497,-634.5308297409192, -347.7885168491812,6.771313016808563,79.70258790927392,31.68212813610444,5.0248961592587875, 0.2879701466217739,0.0); // Compare computed sequence to the standard compareSequences( poly, -500); compareSequences( poly, Double.NEGATIVE_INFINITY); SturmSequence alg = new SturmSequence(poly.size); alg.initialize(poly); int N = alg.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); int M = alg.countRealRoots(-500,500); assertTrue(M <= N); }