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; }
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); } } }
/** * 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); } }
@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()); }
while( u-l > boundTolerance*Math.abs(l) && iter++ < maxBoundIterations) { double m = (l+u)/2.0; int numRoots = sturm.countRealRoots(m,u); if( numRoots == 1 ) { l = m;
/** * Configures search parameters. * * @param maxCoefficients The maximum number of coefficients a polynomial will have. * @param searchRadius If {@code >} 0 then roots are searched for inside of -searchRadius ≤ x ≤ searchRadius. If * ≤ 0 then it will find all roots. * @param boundTolerance How accurately roots are found using bisection. Should be close enough that refinement can * accurately converge to the correct root. * @param maxBoundIterations Maximum number of iterations each step can perform when bounding. * @param maxRefineIterations Maximum number of iterations when refining. */ public FindRealRootsSturm(int maxCoefficients, double searchRadius , double boundTolerance, int maxBoundIterations , int maxRefineIterations ) { if( Double.isInfinite(searchRadius) ) searchRadius = -1; sturm = new SturmSequence(maxCoefficients); this.searchRadius = searchRadius; this.boundTolerance = boundTolerance; this.maxBoundIterations = maxBoundIterations; this.maxRefineIterations = maxRefineIterations; roots = new double[maxCoefficients]; }
/** * 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); }
/** * 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); } }
/** * 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); } }
while( u-l > boundTolerance*Math.abs(l) && iter++ < maxBoundIterations) { double m = (l+u)/2.0; int numRoots = sturm.countRealRoots(m,u); if( numRoots == 1 ) { l = m;
/** * Configures search parameters. * * @param maxCoefficients The maximum number of coefficients a polynomial will have. * @param searchRadius If {@code >} 0 then roots are searched for inside of -searchRadius ≤ x ≤ searchRadius. If * ≤ 0 then it will find all roots. * @param boundTolerance How accurately roots are found using bisection. Should be close enough that refinement can * accurately converge to the correct root. * @param maxBoundIterations Maximum number of iterations each step can perform when bounding. * @param maxRefineIterations Maximum number of iterations when refining. */ public FindRealRootsSturm(int maxCoefficients, double searchRadius , double boundTolerance, int maxBoundIterations , int maxRefineIterations ) { if( Double.isInfinite(searchRadius) ) searchRadius = -1; sturm = new SturmSequence(maxCoefficients); this.searchRadius = searchRadius; this.boundTolerance = boundTolerance; this.maxBoundIterations = maxBoundIterations; this.maxRefineIterations = maxRefineIterations; roots = new double[maxCoefficients]; }
/** * 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); }
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; }
while( iter++ < maxBoundIterations && (totalFound=sturm.countRealRoots(-r,r)) <= 0 ) { r = 2*r*r; if( totalFound < numRoots && (target = sturm.countRealRoots(Double.NEGATIVE_INFINITY,-r)) > 0 ) { double upper = -r; double width = r; int N = sturm.countRealRoots(upper-width,upper); if( N != 0 ) { boundEachRoot(upper-width,upper,totalFound,N); if( totalFound < numRoots && (target = sturm.countRealRoots(r,Double.POSITIVE_INFINITY)) > 0 ) { double lower = r; double width = r; int N = sturm.countRealRoots(lower,lower+width); if( N != 0 ) { boundEachRoot(lower,lower+width,totalFound,N);