public static int countRealRoots( Polynomial poly ) { SturmSequence sturm = new SturmSequence(poly.size); sturm.initialize(poly); return sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); }
/** * Determines the number of real roots there are in the polynomial within the specified bounds. Must call * {@link #initialize(Polynomial)} first. * * @param lower lower limit on the bound. * @param upper Upper limit on the bound * @return Number of real roots */ public int countRealRoots( double lower , double upper ) { // There are no roots for constant equations if( sequenceLength <= 1 ) return 0; computeFunctions(lower); int numLow = countSignChanges(); computeFunctions(upper); int numHigh = countSignChanges(); return numLow-numHigh; }
/** * Creates a generic polynomial root finding class which will return all real or all real and complex roots * depending on the algorithm selected. * * @param type Which algorithm is to be returned. * @param maxDegree Maximum degree of the polynomial being considered. * @return Root finding algorithm. */ public static PolynomialRoots createRootFinder( RootFinderType type , int maxDegree ) { switch ( type ) { case EVD: return new RootFinderCompanion(); case STURM: FindRealRootsSturm sturm = new FindRealRootsSturm(maxDegree,-1,1e-10,30,20); return new WrapRealRootsSturm(sturm); } throw new IllegalArgumentException("Unknown type"); }
@Override public List<Double> computeRealRoots(Polynomial poly) { FindRealRootsSturm alg = new FindRealRootsSturm(poly.size,-1,1e-16,500,500); alg.process(poly); int N = alg.getNumberOfRoots(); double[] roots = alg.getRoots(); List<Double> ret = new ArrayList<Double>(); for( int i = 0; i < N; i++ ) ret.add(roots[i]); return ret; }
/** * Find real roots for the specified polynomial. * * @param poly Polynomial which has less than or equal to the number of coefficients specified in this class's * constructor. */ public void process( Polynomial poly ) { sturm.initialize(poly); if( searchRadius <= 0 ) numRoots = sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); else numRoots = sturm.countRealRoots(-searchRadius,searchRadius); if( numRoots == 0 ) return; if( searchRadius <= 0 ) handleAllRoots(); else { boundEachRoot(-searchRadius,searchRadius,0,numRoots); } // improve the accuracy for( int i = 0; i < numRoots; i++ ) { roots[i] = PolynomialOps.refineRoot(poly, roots[i], maxRefineIterations); } }
@Override public boolean process(Polynomial poly) { try { alg.process(poly); } catch( RuntimeException e ) { return false; } roots.clear(); double found[] = alg.getRoots(); int N = alg.getNumberOfRoots(); for( int i = 0; i < N; i++ ) { Complex_F64 c = storage[i]; c.real = found[i]; c.imaginary = 0; roots.add(c); } return true; }
/** * 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); }
private void compareSequences( Polynomial p, double value) { List<Double> expected = computeSturm(p,value); SturmSequence alg = new SturmSequence(p.size); alg.initialize(p); alg.computeFunctions(value); assertEquals(expected.size(),alg.sequenceLength); for( int j = 0; j < expected.size(); j++ ) { if( Double.isInfinite(expected.get(j)) ) { assertTrue(expected.get(j) == alg.f[j]); } else { assertEquals(expected.get(j),alg.f[j],Math.abs(alg.f[j])*1e-6); } } }
@Test public void checkCount() { SturmSequence alg = new SturmSequence(10); alg.f = new double[]{0,1,1,-1,0,-1,-1,0,1,1}; alg.sequenceLength = alg.f.length; assertEquals(2,alg.countSignChanges()); alg.f = new double[]{1,1,1,-1,0,-1,-1,0,1,1,0,-1,-1,1}; alg.sequenceLength = alg.f.length; assertEquals(4,alg.countSignChanges()); }
public WrapRealRootsSturm(FindRealRootsSturm alg) { this.alg = alg; storage = new Complex_F64[alg.getMaxRoots()]; for( int i = 0; i < storage.length; i++ ) storage[i] = new Complex_F64(); }
/** * Creates different polynomial root finders. * * @param maxCoefficients The maximum number of coefficients that will be processed. This is the order + 1 * @param which 0 = Sturm and 1 = companion matrix. * @return PolynomialRoots */ public static PolynomialRoots createRootFinder( int maxCoefficients , RootFinderType which ) { switch( which ) { case STURM: FindRealRootsSturm sturm = new FindRealRootsSturm(maxCoefficients,-1,1e-10,200,200); return new WrapRealRootsSturm(sturm); case EVD: return new RootFinderCompanion(); default: throw new IllegalArgumentException("Unknown algorithm: "+which); } }
public static int countRealRoots( Polynomial poly ) { SturmSequence sturm = new SturmSequence(poly.size); sturm.initialize(poly); return sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); }
/** * Find real roots for the specified polynomial. * * @param poly Polynomial which has less than or equal to the number of coefficients specified in this class's * constructor. */ public void process( Polynomial poly ) { sturm.initialize(poly); if( searchRadius <= 0 ) numRoots = sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY); else numRoots = sturm.countRealRoots(-searchRadius,searchRadius); if( numRoots == 0 ) return; if( searchRadius <= 0 ) handleAllRoots(); else { boundEachRoot(-searchRadius,searchRadius,0,numRoots); } // improve the accuracy for( int i = 0; i < numRoots; i++ ) { roots[i] = PolynomialOps.refineRoot(poly, roots[i], maxRefineIterations); } }
/** * Determines the number of real roots there are in the polynomial within the specified bounds. Must call * {@link #initialize(Polynomial)} first. * * @param lower lower limit on the bound. * @param upper Upper limit on the bound * @return Number of real roots */ public int countRealRoots( double lower , double upper ) { // There are no roots for constant equations if( sequenceLength <= 1 ) return 0; computeFunctions(lower); int numLow = countSignChanges(); computeFunctions(upper); int numHigh = countSignChanges(); return numLow-numHigh; }
public WrapRealRootsSturm(FindRealRootsSturm alg) { this.alg = alg; storage = new Complex_F64[alg.getMaxRoots()]; for( int i = 0; i < storage.length; i++ ) storage[i] = new Complex_F64(); }
/** * Creates different polynomial root finders. * * @param maxCoefficients The maximum number of coefficients that will be processed. This is the order + 1 * @param which 0 = Sturm and 1 = companion matrix. * @return PolynomialRoots */ public static PolynomialRoots createRootFinder( int maxCoefficients , RootFinderType which ) { switch( which ) { case STURM: FindRealRootsSturm sturm = new FindRealRootsSturm(maxCoefficients,-1,1e-10,200,200); return new WrapRealRootsSturm(sturm); case EVD: return new RootFinderCompanion(); default: throw new IllegalArgumentException("Unknown algorithm: "+which); } }
public int countRealRoots( Polynomial p , double low , double high ) { SturmSequence alg = new SturmSequence(p.size); alg.initialize(p); return alg.countRealRoots(low,high); }
/** * Creates a generic polynomial root finding class which will return all real or all real and complex roots * depending on the algorithm selected. * * @param type Which algorithm is to be returned. * @param maxDegree Maximum degree of the polynomial being considered. * @return Root finding algorithm. */ public static PolynomialRoots createRootFinder( RootFinderType type , int maxDegree ) { switch ( type ) { case EVD: return new RootFinderCompanion(); case STURM: FindRealRootsSturm sturm = new FindRealRootsSturm(maxDegree,-1,1e-10,30,20); return new WrapRealRootsSturm(sturm); } throw new IllegalArgumentException("Unknown type"); }