001/*
002 * $Id$
003 */
004
005package edu.jas.root;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.LogManager;
012import org.apache.logging.log4j.Logger;
013
014import edu.jas.arith.BigDecimal;
015import edu.jas.arith.BigInteger;
016import edu.jas.arith.BigRational;
017import edu.jas.arith.Rational;
018import edu.jas.poly.GenPolynomial;
019import edu.jas.poly.GenPolynomialRing;
020import edu.jas.poly.PolyUtil;
021import edu.jas.structure.RingElem;
022import edu.jas.structure.RingFactory;
023import edu.jas.structure.UnaryFunctor;
024
025
026/**
027 * Real roots abstract class.
028 * @param <C> coefficient type.
029 * @author Heinz Kredel
030 */
031public abstract class RealRootsAbstract<C extends RingElem<C> & Rational> implements RealRoots<C> {
032
033
034    private static final Logger logger = LogManager.getLogger(RealRootsAbstract.class);
035
036
037    //private static final boolean debug = logger.isDebugEnabled();
038
039
040    /**
041     * Real root bound. With f(-M) * f(M) != 0.
042     * @param f univariate polynomial.
043     * @return M such that -M &lt; root(f) &lt; M.
044     */
045    public C realRootBound(GenPolynomial<C> f) {
046        if (f == null) {
047            return null;
048        }
049        RingFactory<C> cfac = f.ring.coFac;
050        C M = cfac.getONE();
051        if (f.isZERO()) {
052            return M;
053        }
054        if (f.isConstant()) {
055            M = f.leadingBaseCoefficient().abs().sum(cfac.getONE());
056            return M;
057        }
058        C a = f.leadingBaseCoefficient().abs();
059        for (C c : f.getMap().values()) {
060            C d = c.abs().divide(a);
061            if (M.compareTo(d) < 0) {
062                M = d;
063            }
064        }
065        BigRational r = M.getRational();
066        logger.info("rational root bound: {}", r);
067        BigInteger i = new BigInteger(r.numerator().divide(r.denominator()));
068        i = i.sum(BigInteger.ONE); // ceiling
069        M = cfac.fromInteger(i.getVal());
070        M = M.sum(f.ring.coFac.getONE());
071        logger.info("integer root bound: {}", M);
072        return M;
073    }
074
075
076    /**
077     * Magnitude bound.
078     * @param iv interval.
079     * @param f univariate polynomial.
080     * @return B such that |f(c)| &lt; B for c in iv.
081     */
082    @SuppressWarnings("unchecked")
083    public C magnitudeBound(Interval<C> iv, GenPolynomial<C> f) {
084        if (f == null) {
085            return null;
086        }
087        if (f.isZERO()) {
088            return f.ring.coFac.getONE();
089        }
090        if (f.isConstant()) {
091            return f.leadingBaseCoefficient().abs();
092        }
093        GenPolynomial<C> fa = f.map(new UnaryFunctor<C, C>() {
094
095
096            public C eval(C a) {
097                return a.abs();
098            }
099        });
100        //System.out.println("fa = " + fa);
101        C M = iv.left.abs();
102        if (M.compareTo(iv.right.abs()) < 0) {
103            M = iv.right.abs();
104        }
105        //System.out.println("M = " + M);
106        RingFactory<C> cfac = f.ring.coFac;
107        C B = PolyUtil.<C> evaluateMain(cfac, fa, M);
108        // works also without this case, only for optimization 
109        // to use rational number interval end points
110        // can fail if real root is in interval [r,r+1] 
111        // for too low precision or too big r, since r is approximation
112        //if ((Object) B instanceof RealAlgebraicNumber) {
113        //    RealAlgebraicNumber Br = (RealAlgebraicNumber) B;
114        //    BigRational r = Br.magnitude();
115        //    B = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator()));
116        //}
117        BigRational r = B.getRational();
118        B = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator()));
119        //System.out.println("B = " + B);
120        return B;
121    }
122
123
124    /**
125     * Real minimal root bound.
126     * @param f univariate polynomial.
127     * @return M such that abs(xi) &gt; M for f(xi) == 0.
128     */
129    public C realMinimalRootBound(GenPolynomial<C> f) {
130        if (f == null) {
131            return null;
132        }
133        RingFactory<C> cfac = f.ring.coFac;
134        // maxNorm root bound
135        BigRational mr = f.maxNorm().getRational().sum(BigRational.ONE);
136        BigRational di = mr.sum(BigRational.ONE).inverse();
137        C B = cfac.fromInteger(di.numerator()).divide(cfac.fromInteger(di.denominator()));
138        //System.out.println("B = " + B + ", sign(B) = " + B.signum());
139        return B;
140    }
141
142
143    /**
144     * Real minimal root separation.
145     * @param f univariate polynomial.
146     * @return M such that abs(xi-xj) &gt; M for roots xi, xj of f.
147     */
148    public C realMinimalRootSeparation(GenPolynomial<C> f) {
149        if (f == null) {
150            return null;
151        }
152        RingFactory<C> cfac = f.ring.coFac;
153        // sumNorm root bound
154        BigRational pr = f.sumNorm().getRational();
155        pr = pr.sum(BigRational.ONE);
156        BigRational sep = BigRational.ZERO;
157        long n = f.degree();
158        if (n > 0) {
159            sep = pr.power(2 * n).multiply(pr.fromInteger(n).power(n + 1)).inverse();
160        }
161        //System.out.println("sep = " + sep + ", sign(sep) = " + sep.signum());
162        C M = cfac.fromInteger(sep.numerator()).divide(cfac.fromInteger(sep.denominator()));
163        return M;
164    }
165
166
167    /**
168     * Bi-section point.
169     * @param iv interval with f(left) * f(right) != 0.
170     * @param f univariate polynomial, non-zero.
171     * @return a point c in the interval iv such that f(c) != 0.
172     */
173    public C bisectionPoint(Interval<C> iv, GenPolynomial<C> f) {
174        if (f == null) {
175            return null;
176        }
177        RingFactory<C> cfac = f.ring.coFac;
178        C two = cfac.fromInteger(2);
179        C c = iv.left.sum(iv.right);
180        c = c.divide(two);
181        if (f.isZERO() || f.isConstant()) {
182            return c;
183        }
184        C m = PolyUtil.<C> evaluateMain(cfac, f, c);
185        while (m.isZERO()) {
186            C d = iv.left.sum(c);
187            d = d.divide(two);
188            if (d.equals(c)) {
189                d = iv.right.sum(c);
190                d = d.divide(two);
191                if (d.equals(c)) {
192                    throw new RuntimeException("should not happen " + iv);
193                }
194            }
195            c = d;
196            m = PolyUtil.<C> evaluateMain(cfac, f, c);
197            //System.out.println("c = " + c);
198        }
199        //System.out.println("c = " + c);
200        return c;
201    }
202
203
204    /**
205     * Isolating intervals for the real roots.
206     * @param f univariate polynomial.
207     * @return a list of isolating intervals for the real roots of f.
208     */
209    public abstract List<Interval<C>> realRoots(GenPolynomial<C> f);
210
211
212    /**
213     * Isolating intervals for the real roots.
214     * @param f univariate polynomial.
215     * @param eps requested intervals length.
216     * @return a list of isolating intervals v such that |v| &lt; eps.
217     */
218    public List<Interval<C>> realRoots(GenPolynomial<C> f, C eps) {
219        return realRoots(f, eps.getRational());
220    }
221
222
223    /**
224     * Isolating intervals for the real roots.
225     * @param f univariate polynomial.
226     * @param eps requested intervals length.
227     * @return a list of isolating intervals v such that |v| &lt; eps.
228     */
229    public List<Interval<C>> realRoots(GenPolynomial<C> f, BigRational eps) {
230        List<Interval<C>> iv = realRoots(f);
231        return refineIntervals(iv, f, eps);
232    }
233
234
235    /**
236     * Sign changes on interval bounds.
237     * @param iv root isolating interval with f(left) * f(right) != 0.
238     * @param f univariate polynomial.
239     * @return true if f(left) * f(right) &lt; 0, else false
240     */
241    public boolean signChange(Interval<C> iv, GenPolynomial<C> f) {
242        if (f == null) {
243            return false;
244        }
245        RingFactory<C> cfac = f.ring.coFac;
246        C l = PolyUtil.<C> evaluateMain(cfac, f, iv.left);
247        C r = PolyUtil.<C> evaluateMain(cfac, f, iv.right);
248        return l.signum() * r.signum() < 0;
249    }
250
251
252    /**
253     * Number of real roots in interval.
254     * @param iv interval with f(left) * f(right) != 0.
255     * @param f univariate polynomial.
256     * @return number of real roots of f in I.
257     */
258    public abstract long realRootCount(Interval<C> iv, GenPolynomial<C> f);
259
260
261    /**
262     * Half interval.
263     * @param iv root isolating interval with f(left) * f(right) &lt; 0.
264     * @param f univariate polynomial, non-zero.
265     * @return a new interval v such that |v| &lt; |iv|/2.
266     */
267    public Interval<C> halfInterval(Interval<C> iv, GenPolynomial<C> f) {
268        if (f == null || f.isZERO()) {
269            return iv;
270        }
271        BigRational len = iv.rationalLength();
272        BigRational two = len.factory().fromInteger(2);
273        BigRational eps = len.divide(two);
274        return refineInterval(iv, f, eps);
275    }
276
277
278    /**
279     * Refine interval.
280     * @param iv root isolating interval with f(left) * f(right) &lt; 0.
281     * @param f univariate polynomial, non-zero.
282     * @param eps requested interval length.
283     * @return a new interval v such that |v| &lt; eps.
284     */
285    public Interval<C> refineInterval(Interval<C> iv, GenPolynomial<C> f, BigRational eps) {
286        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
287            return iv;
288        }
289        if (iv.rationalLength().compareTo(eps) < 0) {
290            return iv;
291        }
292        RingFactory<C> cfac = f.ring.coFac;
293        C two = cfac.fromInteger(2);
294        Interval<C> v = iv;
295        while (v.rationalLength().compareTo(eps) >= 0) {
296            C c = v.left.sum(v.right);
297            c = c.divide(two);
298            //System.out.println("c = " + c);
299            //c = RootUtil.<C>bisectionPoint(v,f);
300            if (PolyUtil.<C> evaluateMain(cfac, f, c).isZERO()) {
301                v = new Interval<C>(c, c);
302                break;
303            }
304            Interval<C> iv1 = new Interval<C>(v.left, c);
305            if (signChange(iv1, f)) {
306                v = iv1;
307            } else {
308                v = new Interval<C>(c, v.right);
309            }
310        }
311        return v;
312    }
313
314
315    /**
316     * Refine intervals.
317     * @param V list of isolating intervals with f(left) * f(right) &lt; 0.
318     * @param f univariate polynomial, non-zero.
319     * @param eps requested intervals length.
320     * @return a list of new intervals v such that |v| &lt; eps.
321     */
322    public List<Interval<C>> refineIntervals(List<Interval<C>> V, GenPolynomial<C> f, BigRational eps) {
323        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
324            return V;
325        }
326        List<Interval<C>> IV = new ArrayList<Interval<C>>();
327        for (Interval<C> v : V) {
328            Interval<C> iv = refineInterval(v, f, eps);
329            IV.add(iv);
330        }
331        return IV;
332    }
333
334
335    /**
336     * Invariant interval for algebraic number sign.
337     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
338     * @param f univariate polynomial, non-zero.
339     * @param g univariate polynomial, gcd(f,g) == 1.
340     * @return v with v a new interval contained in iv such that g(v) != 0.
341     */
342    public abstract Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g);
343
344
345    /**
346     * Real algebraic number sign.
347     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0,
348     *            with iv such that g(iv) != 0.
349     * @param f univariate polynomial, non-zero.
350     * @param g univariate polynomial, gcd(f,g) == 1.
351     * @return sign(g(iv)) .
352     */
353    public int realIntervalSign(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
354        if (g == null || g.isZERO()) {
355            return 0;
356        }
357        if (f == null || f.isZERO() || f.isConstant()) {
358            return g.signum();
359        }
360        if (g.isConstant()) {
361            return g.signum();
362        }
363        RingFactory<C> cfac = f.ring.coFac;
364        C c = iv.left.sum(iv.right);
365        c = c.divide(cfac.fromInteger(2));
366        C ev = PolyUtil.<C> evaluateMain(cfac, g, c);
367        //System.out.println("ev = " + ev);
368        return ev.signum();
369    }
370
371
372    /**
373     * Real algebraic number sign.
374     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
375     * @param f univariate polynomial, non-zero.
376     * @param g univariate polynomial, gcd(f,g) == 1.
377     * @return sign(g(v)), with v a new interval contained in iv such that g(v)
378     *         != 0.
379     */
380    public int realSign(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
381        if (g == null || g.isZERO()) {
382            return 0;
383        }
384        if (f == null || f.isZERO() || f.isConstant()) {
385            return g.signum();
386        }
387        if (g.isConstant()) {
388            return g.signum();
389        }
390        Interval<C> v = invariantSignInterval(iv, f, g);
391        return realIntervalSign(v, f, g);
392    }
393
394
395    /**
396     * Invariant interval for algebraic number magnitude.
397     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
398     * @param f univariate polynomial, non-zero.
399     * @param g univariate polynomial, gcd(f,g) == 1.
400     * @param eps length limit for interval length.
401     * @return v with v a new interval contained in iv such that |g(a) - g(b)|
402     *         &lt; eps for a, b in v in iv.
403     */
404    public Interval<C> invariantMagnitudeInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g,
405                    BigRational eps) {
406        Interval<C> v = iv;
407        if (g == null || g.isZERO()) {
408            return v;
409        }
410        if (g.isConstant()) {
411            return v;
412        }
413        if (f == null || f.isZERO() || f.isConstant()) { // ?
414            return v;
415        }
416        GenPolynomial<C> gp = PolyUtil.<C> baseDerivative(g);
417        //System.out.println("g  = " + g);
418        //System.out.println("gp = " + gp);
419        C B = magnitudeBound(iv, gp);
420        //System.out.println("B = " + B);
421        RingFactory<C> cfac = f.ring.coFac;
422        C two = cfac.fromInteger(2);
423        while (B.multiply(v.length()).getRational().compareTo(eps) >= 0) {
424            C c = v.left.sum(v.right);
425            c = c.divide(two);
426            Interval<C> im = new Interval<C>(c, v.right);
427            if (signChange(im, f)) {
428                v = im;
429            } else {
430                v = new Interval<C>(v.left, c);
431            }
432            //System.out.println("v = " + v.toDecimal());
433        }
434        return v;
435    }
436
437
438    /**
439     * Real algebraic number magnitude.
440     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0,
441     *            with iv such that |g(a) - g(b)| &lt; eps for a, b in iv.
442     * @param f univariate polynomial, non-zero.
443     * @param g univariate polynomial, gcd(f,g) == 1.
444     * @return g(iv) .
445     */
446    public C realIntervalMagnitude(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
447        if (g.isZERO() || g.isConstant()) {
448            return g.leadingBaseCoefficient();
449        }
450        RingFactory<C> cfac = f.ring.coFac;
451        C evl = PolyUtil.<C> evaluateMain(cfac, g, iv.left);
452        C evr = PolyUtil.<C> evaluateMain(cfac, g, iv.right);
453        C ev = evl;
454        if (evl.compareTo(evr) <= 0) {
455            ev = evr;
456        }
457        //System.out.println("ev = " + ev + ", evl = " + evl + ", evr = " + evr + ", iv = " + iv);
458        return ev;
459    }
460
461
462    /**
463     * Real algebraic number magnitude.
464     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
465     * @param f univariate polynomial, non-zero.
466     * @param g univariate polynomial, gcd(f,g) == 1.
467     * @param eps length limit for interval length.
468     * @return g(iv) .
469     */
470    public C realMagnitude(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g, BigRational eps) {
471        if (g.isZERO() || g.isConstant()) {
472            return g.leadingBaseCoefficient();
473        }
474        Interval<C> v = invariantMagnitudeInterval(iv, f, g, eps);
475        return realIntervalMagnitude(v, f, g);
476    }
477
478
479    /**
480     * Real algebraic number magnitude.
481     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0,
482     *            with iv such that |g(a) - g(b)| &lt; eps for a, b in iv.
483     * @param f univariate polynomial, non-zero.
484     * @param g univariate polynomial, gcd(f,g) == 1.
485     * @return Interval( g(iv.left), g(iv.right) ) .
486     */
487    public Interval<C> realIntervalMagnitudeInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
488        if (g.isZERO() || g.isConstant()) {
489            return iv;
490        }
491        RingFactory<C> cfac = f.ring.coFac;
492        C evl = PolyUtil.<C> evaluateMain(cfac, g, iv.left);
493        C evr = PolyUtil.<C> evaluateMain(cfac, g, iv.right);
494        Interval<C> ev = new Interval<C>(evr, evl);
495        if (evl.compareTo(evr) <= 0) {
496            ev = new Interval<C>(evl, evr);
497        }
498        System.out.println("ev = " + ev + ", iv = " + iv);
499        return ev;
500    }
501
502
503    /**
504     * Approximate real root.
505     * @param iv real root isolating interval with f(left) * f(right) &lt; 0.
506     * @param f univariate polynomial, non-zero.
507     * @param eps requested interval length.
508     * @return a decimal approximation d such that |d-v| &lt; eps, for f(v) = 0,
509     *         v real.
510     */
511    public BigDecimal approximateRoot(Interval<C> iv, GenPolynomial<C> f, BigRational eps)
512                    throws NoConvergenceException {
513        if (iv == null) {
514            throw new IllegalArgumentException("null interval not allowed");
515        }
516        BigDecimal d = iv.toDecimal();
517        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
518            return d;
519        }
520        if (iv.rationalLength().compareTo(eps) < 0) {
521            return d;
522        }
523        BigDecimal left = new BigDecimal(iv.left.getRational());
524        BigDecimal right = new BigDecimal(iv.right.getRational());
525        BigRational reps = eps.getRational();
526        BigDecimal e = new BigDecimal(reps);
527        BigDecimal q = new BigDecimal("0.25");
528        //System.out.println("left  = " + left);
529        //System.out.println("right = " + right);
530        e = e.multiply(d); // relative error
531        //System.out.println("e     = " + e);
532        BigDecimal dc = BigDecimal.ONE;
533        // polynomials with decimal coefficients
534        GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring);
535        GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f);
536        GenPolynomial<C> fp = PolyUtil.<C> baseDerivative(f);
537        GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp);
538        // Newton Raphson iteration: x_{n+1} = x_n - f(x_n)/f'(x_n)
539        int i = 0;
540        final int MITER = 50;
541        int dir = 0;
542        while (i++ < MITER) {
543            BigDecimal fx = PolyUtil.<BigDecimal> evaluateMain(dc, df, d); // f(d)
544            if (fx.isZERO()) {
545                return d;
546            }
547            BigDecimal fpx = PolyUtil.<BigDecimal> evaluateMain(dc, dfp, d); // f'(d)
548            if (fpx.isZERO()) {
549                throw new NoConvergenceException("zero derivative should not happen");
550            }
551            BigDecimal x = fx.divide(fpx);
552            BigDecimal dx = d.subtract(x);
553            //System.out.println("dx = " + dx + ", d = " + d);
554            if (d.subtract(dx).abs().compareTo(e) <= 0) {
555                return dx;
556            }
557            while (dx.compareTo(left) < 0 || dx.compareTo(right) > 0) {
558                // dx < left: dx - left < 0
559                // dx > right: dx - right > 0
560                //System.out.println("trying to leave interval");
561                if (i++ > MITER) { // dx > right: dx - right > 0
562                    throw new NoConvergenceException("no convergence after " + i + " steps");
563                }
564                if (i > MITER / 2 && dir == 0) {
565                    BigDecimal sd = new BigDecimal(iv.randomPoint().getRational());
566                    d = sd;
567                    x = sd.getZERO();
568                    logger.info("trying new random starting point {}", d);
569                    i = 0;
570                    dir = 1;
571                }
572                if (i > MITER / 2 && dir == 1) {
573                    BigDecimal sd = new BigDecimal(iv.randomPoint().getRational());
574                    d = sd;
575                    x = sd.getZERO();
576                    logger.info("trying new random starting point {}", d);
577                    //i = 0;
578                    dir = 2; // end
579                }
580                x = x.multiply(q); // x * 1/4
581                dx = d.subtract(x);
582                //System.out.println(" x = " + x);
583                //System.out.println("dx = " + dx);
584            }
585            d = dx;
586        }
587        throw new NoConvergenceException("no convergence after " + i + " steps");
588    }
589
590
591    /**
592     * Approximate real roots.
593     * @param f univariate polynomial, non-zero.
594     * @param eps requested interval length.
595     * @return a list of decimal approximations d such that |d-v| &lt; eps for
596     *         all real v with f(v) = 0.
597     */
598    public List<BigDecimal> approximateRoots(GenPolynomial<C> f, BigRational eps) {
599        List<Interval<C>> iv = realRoots(f);
600        List<BigDecimal> roots = new ArrayList<BigDecimal>(iv.size());
601        for (Interval<C> i : iv) {
602            BigDecimal r = null; //approximateRoot(i, f, eps); roots.add(r);
603            while (r == null) {
604                try {
605                    r = approximateRoot(i, f, eps);
606                    roots.add(r);
607                } catch (NoConvergenceException e) {
608                    // fall back to exact algorithm
609                    //System.out.println("" + e);
610                    BigRational len = i.rationalLength();
611                    len = len.divide(len.factory().fromInteger(1000));
612                    i = refineInterval(i, f, len);
613                    logger.info("fall back rootRefinement = {}", i);
614                }
615            }
616        }
617        return roots;
618    }
619
620
621    /**
622     * Test if x is an approximate real root.
623     * @param x approximate real root.
624     * @param f univariate polynomial, non-zero.
625     * @param eps requested interval length.
626     * @return true if x is a decimal approximation of a real v with f(v) = 0
627     *         with |d-v| &lt; eps, else false.
628     */
629    public boolean isApproximateRoot(BigDecimal x, GenPolynomial<C> f, C eps) {
630        if (x == null) {
631            throw new IllegalArgumentException("null root not allowed");
632        }
633        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
634            return true;
635        }
636        BigDecimal e = new BigDecimal(eps.getRational());
637        e = e.multiply(new BigDecimal("1000")); // relax
638        BigDecimal dc = BigDecimal.ONE;
639        // polynomials with decimal coefficients
640        GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring);
641        GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f);
642        GenPolynomial<C> fp = PolyUtil.<C> baseDerivative(f);
643        GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp);
644        //
645        return isApproximateRoot(x, df, dfp, e);
646    }
647
648
649    /**
650     * Test if x is an approximate real root.
651     * @param x approximate real root.
652     * @param f univariate polynomial, non-zero.
653     * @param fp univariate polynomial, non-zero, derivative of f.
654     * @param eps requested interval length.
655     * @return true if x is a decimal approximation of a real v with f(v) = 0
656     *         with |d-v| &lt; eps, else false.
657     */
658    public boolean isApproximateRoot(BigDecimal x, GenPolynomial<BigDecimal> f, GenPolynomial<BigDecimal> fp,
659                    BigDecimal eps) {
660        if (x == null) {
661            throw new IllegalArgumentException("null root not allowed");
662        }
663        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
664            return true;
665        }
666        BigDecimal dc = BigDecimal.ONE; // only for clarity
667        // f(x)
668        BigDecimal fx = PolyUtil.<BigDecimal> evaluateMain(dc, f, x);
669        //System.out.println("fx    = " + fx);
670        if (fx.isZERO()) {
671            return true;
672        }
673        // f'(x)
674        BigDecimal fpx = PolyUtil.<BigDecimal> evaluateMain(dc, fp, x); // f'(d)
675        //System.out.println("fpx   = " + fpx);
676        if (fpx.isZERO()) {
677            return false;
678        }
679        BigDecimal d = fx.divide(fpx);
680        if (d.isZERO()) {
681            return true;
682        }
683        if (d.abs().compareTo(eps) <= 0) {
684            return true;
685        }
686        System.out.println("x     = " + x);
687        System.out.println("d     = " + d);
688        return false;
689    }
690
691
692    /**
693     * Test if each x in R is an approximate real root.
694     * @param R ist of approximate real roots.
695     * @param f univariate polynomial, non-zero.
696     * @param eps requested interval length.
697     * @return true if each x in R is a decimal approximation of a real v with
698     *         f(v) = 0 with |d-v| &lt; eps, else false.
699     */
700    public boolean isApproximateRoot(List<BigDecimal> R, GenPolynomial<C> f, BigRational eps) {
701        if (R == null) {
702            throw new IllegalArgumentException("null root not allowed");
703        }
704        if (f == null || f.isZERO() || f.isConstant() || eps == null) {
705            return true;
706        }
707        BigDecimal e = new BigDecimal(eps.getRational());
708        e = e.multiply(new BigDecimal("1000")); // relax
709        BigDecimal dc = BigDecimal.ONE;
710        // polynomials with decimal coefficients
711        GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring);
712        GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f);
713        GenPolynomial<C> fp = PolyUtil.<C> baseDerivative(f);
714        GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp);
715        for (BigDecimal x : R) {
716            if (!isApproximateRoot(x, df, dfp, e)) {
717                return false;
718            }
719        }
720        return true;
721    }
722
723
724    /**
725     * Fourier sequence.
726     * @param f univariate polynomial.
727     * @return (f, f', ..., f(n)) a Fourier sequence for f.
728     */
729    public List<GenPolynomial<C>> fourierSequence(GenPolynomial<C> f) {
730        List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>();
731        if (f == null || f.isZERO()) {
732            return S;
733        }
734        long d = f.degree();
735        GenPolynomial<C> F = f;
736        S.add(F);
737        while (d-- > 0) {
738            GenPolynomial<C> G = PolyUtil.<C> baseDerivative(F);
739            F = G;
740            S.add(F);
741        }
742        //System.out.println("F = " + F);
743        return S;
744    }
745
746
747    /**
748     * Thom sign sequence.
749     * @param f univariate polynomial.
750     * @param v interval for a real root, f(v.left) * f(v.right) &lt; 0.
751     * @return (s1, s2, ..., sn) = (sign(f'(v)), .... sign(f(n)(v))) a Thom sign
752     *         sequence for the real root in v of f.
753     */
754    public List<Integer> signSequence(GenPolynomial<C> f, Interval<C> v) {
755        List<Integer> S = new ArrayList<Integer>();
756        GenPolynomial<C> fp = PolyUtil.<C> baseDerivative(f);
757        List<GenPolynomial<C>> fs = fourierSequence(fp);
758        for (GenPolynomial<C> p : fs) {
759            int s = realSign(v, f, p);
760            S.add(s);
761        }
762        return S;
763    }
764
765
766    /**
767     * Root number.
768     * @param f univariate polynomial.
769     * @param v interval for a real root, f(v.left) * f(v.right) &lt; 0.
770     * @return r the number of this root in the sequence a1 &lt; a2 &lt; ...,
771     *         &lt; am of all real roots of f
772     */
773    public Long realRootNumber(GenPolynomial<C> f, Interval<C> v) {
774        C M = realRootBound(f);
775        Interval<C> iv = new Interval<C>(M.negate(), v.right);
776        long r = realRootCount(iv, f);
777        if (r <= 0) {
778            logger.warn("no real root in interval {}", v);
779        }
780        return r;
781    }
782
783}