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 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; }
/** * Configures the algorithm. * * @param maxPolySize The maximum number of coefficients on the polynomial being processed. */ public SturmSequence( int maxPolySize ) { next = new Polynomial(maxPolySize); previous = new Polynomial(maxPolySize); result = new Polynomial(maxPolySize); sequence = new Polynomial[maxPolySize+1]; for( int i = 0; i < sequence.length; i++ ) sequence[i] = new Polynomial(maxPolySize); f = new double[ maxPolySize + 1]; }
/** * Configures the algorithm. * * @param maxPolySize The maximum number of coefficients on the polynomial being processed. */ public SturmSequence( int maxPolySize ) { next = new Polynomial(maxPolySize); previous = new Polynomial(maxPolySize); result = new Polynomial(maxPolySize); sequence = new Polynomial[maxPolySize+1]; for( int i = 0; i < sequence.length; i++ ) sequence[i] = new Polynomial(maxPolySize); f = new double[ maxPolySize + 1]; }
/** * 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; }
/** * 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; }
private Polynomial createRandom( int N ) { Polynomial ret = new Polynomial(N); for( int i = 0; i < N; i++ ) { ret.c[i] = 5*(rand.nextDouble()-0.5); } return ret; }
/** * 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; }
/** * 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; }
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)); } }
/** * Compare the computed sequence against randomly generated polynomials of different length */ @Test public void checkSequence_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; } double value = (rand.nextDouble()-0.5)*4; compareSequences( p, value); compareSequences( p, Double.POSITIVE_INFINITY); compareSequences( p, Double.NEGATIVE_INFINITY); } } }
@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); } } }
/** * 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)); } } }
@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)); }
@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); }
@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); } } }
@Test public void resize() { Polynomial a = new Polynomial(10); assertEquals(10,a.size()); a.resize(5); assertEquals(5,a.size()); a.resize(15); assertEquals(15,a.size()); assertEquals(15,a.c.length); }
@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 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); }