@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; }
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(); }
@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. 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 ) {} }
/** * 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); } }
/** * 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); } }
l = m; } else if( found == 1 ) { bisectionRoot(l,m,startIndex + root++); l = m; u = lastUpper = allUpper; lastFound = 0;
throw new RuntimeException("Too many iterations when searching center region"); boundEachRoot(-r,r,0,totalFound); int N = sturm.countRealRoots(upper-width,upper); if( N != 0 ) { boundEachRoot(upper-width,upper,totalFound,N); target -= N; totalFound += N; int N = sturm.countRealRoots(lower,lower+width); if( N != 0 ) { boundEachRoot(lower,lower+width,totalFound,N); target -= N; totalFound += N;
@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; }
/** * 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"); }
/** * 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); } }
l = m; } else if( found == 1 ) { bisectionRoot(l,m,startIndex + root++); l = m; u = lastUpper = allUpper; lastFound = 0;
throw new RuntimeException("Too many iterations when searching center region"); boundEachRoot(-r,r,0,totalFound); int N = sturm.countRealRoots(upper-width,upper); if( N != 0 ) { boundEachRoot(upper-width,upper,totalFound,N); target -= N; totalFound += N; int N = sturm.countRealRoots(lower,lower+width); if( N != 0 ) { boundEachRoot(lower,lower+width,totalFound,N); target -= N; totalFound += N;
/** * 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 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 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"); }