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.BigRational;
015import edu.jas.arith.Rational;
016import edu.jas.poly.ExpVector;
017import edu.jas.poly.GenPolynomial;
018import edu.jas.poly.GenPolynomialRing;
019import edu.jas.poly.PolyUtil;
020import edu.jas.structure.RingElem;
021import edu.jas.structure.RingFactory;
022
023
024/**
025 * Real root isolation using Sturm sequences.
026 * @param <C> coefficient type.
027 * @author Heinz Kredel
028 */
029public class RealRootsSturm<C extends RingElem<C> & Rational> extends RealRootsAbstract<C> {
030
031
032    private static final Logger logger = LogManager.getLogger(RealRootsSturm.class);
033
034
035    private static final boolean debug = logger.isDebugEnabled();
036
037
038    /**
039     * Sturm sequence.
040     * @param f univariate polynomial.
041     * @return a Sturm sequence for f.
042     */
043    public List<GenPolynomial<C>> sturmSequence(GenPolynomial<C> f) {
044        List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>();
045        if (f == null || f.isZERO()) {
046            return S;
047        }
048        if (f.isConstant()) {
049            S.add(f.monic());
050            return S;
051        }
052        GenPolynomial<C> F = f;
053        S.add(F);
054        GenPolynomial<C> G = PolyUtil.<C> baseDerivative(f);
055        while (!G.isZERO()) {
056            GenPolynomial<C> r = F.remainder(G);
057            F = G;
058            G = r.negate();
059            S.add(F/*.monic()*/);
060        }
061        //System.out.println("F = " + F);
062        if (F.isConstant()) {
063            return S;
064        }
065        // make squarefree
066        List<GenPolynomial<C>> Sp = new ArrayList<GenPolynomial<C>>(S.size());
067        for (GenPolynomial<C> p : S) {
068            p = p.divide(F);
069            Sp.add(p);
070        }
071        return Sp;
072    }
073
074
075    /**
076     * Isolating intervals for the real roots.
077     * @param f univariate polynomial.
078     * @return a list of isolating intervals for the real roots of f.
079     */
080    @SuppressWarnings("cast")
081    @Override
082    public List<Interval<C>> realRoots(GenPolynomial<C> f) {
083        List<Interval<C>> R = new ArrayList<Interval<C>>();
084        if (f == null) {
085            return R;
086        }
087        GenPolynomialRing<C> pfac = f.ring;
088        if (f.isZERO()) {
089            C z = pfac.coFac.getZERO();
090            R.add(new Interval<C>(z));
091            return R;
092        }
093        // check trailing degree
094        ExpVector et = f.trailingExpVector();
095        if (!et.isZERO()) {
096            GenPolynomial<C> tr = pfac.valueOf(et);
097            logger.info("trailing term = {}", tr);
098            f = PolyUtil.<C> basePseudoDivide(f, tr);
099            R.add(new Interval<C>(pfac.coFac.getZERO()));
100        }
101        if (f.isConstant()) {
102            return R;
103        }
104        if (f.degree(0) == 1L) {
105            C z = f.monic().trailingBaseCoefficient().negate();
106            R.add(new Interval<C>(z));
107            return R;
108        }
109        //if (f.degree(0) == 2L) ... 
110        GenPolynomial<C> F = f;
111        C M = realRootBound(F); // M != 0, since >= 2
112        Interval<C> iv = new Interval<C>(M.negate(), M);
113        //System.out.println("iv = " + iv);
114        List<GenPolynomial<C>> S = sturmSequence(F);
115        //System.out.println("S = " + S);
116        //System.out.println("f_S = " + S.get(0));
117        List<Interval<C>> Rp = realRoots(iv, S);
118        if (!(((Object) f.ring.coFac) instanceof BigRational)) {
119            //logger.info("realRoots bound: {}", iv);
120            logger.info("realRoots: {}", Rp);
121        }
122        R.addAll(Rp);
123        return R;
124    }
125
126
127    /**
128     * Isolating intervals for the real roots.
129     * @param iv interval with f(left) * f(right) != 0.
130     * @param S sturm sequence for f and I.
131     * @return a list of isolating intervals for the real roots of f in I.
132     */
133    public List<Interval<C>> realRoots(Interval<C> iv, List<GenPolynomial<C>> S) {
134        List<Interval<C>> R = new ArrayList<Interval<C>>();
135        GenPolynomial<C> f = S.get(0); // squarefree part
136        if (f.isZERO()) {
137            C z = f.leadingBaseCoefficient();
138            if (!iv.contains(z)) {
139                throw new IllegalArgumentException(
140                                                   "root not in interval: f = " + f + ", iv = " + iv + ", z = " + z);
141            }
142            Interval<C> iv1 = new Interval<C>(z);
143            R.add(iv1);
144            return R;
145        }
146        if (f.isConstant()) {
147            return R;
148            //throw new IllegalArgumentException("f has no root: f = " + f + ", iv = " + iv);
149        }
150        if (f.degree(0) == 1L) {
151            C z = f.monic().trailingBaseCoefficient().negate();
152            if (!iv.contains(z)) {
153                return R;
154                //throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv + ", z = " + z);
155            }
156            Interval<C> iv1 = new Interval<C>(z);
157            R.add(iv1);
158            return R;
159        }
160        //System.out.println("iv = " + iv);
161        // check sign variations at interval bounds
162        long v = realRootCount(iv, S);
163        //System.out.println("v = " + v);
164        if (v == 0) {
165            return R;
166        }
167        if (v == 1) {
168            iv = excludeZero(iv, S);
169            R.add(iv);
170            return R;
171        }
172        // now v &gt; 1
173        // bi-sect interval, such that f(c) != 0
174        C c = bisectionPoint(iv, f);
175        //System.out.println("c = " + c);
176        // recursion on both sub-intervals
177        Interval<C> iv1 = new Interval<C>(iv.left, c);
178        Interval<C> iv2 = new Interval<C>(c, iv.right);
179        List<Interval<C>> R1 = realRoots(iv1, S);
180        //System.out.println("R1 = " + R1);
181        if (debug) {
182            logger.info("R1 = {}", R1);
183        }
184        List<Interval<C>> R2 = realRoots(iv2, S);
185        //System.out.println("R2 = " + R2);
186        if (debug) {
187            logger.info("R2 = {}", R2);
188        }
189
190        // refine isolating intervals if adjacent 
191        if (R1.isEmpty()) {
192            R.addAll(R2);
193            return R;
194        }
195        if (R2.isEmpty()) {
196            R.addAll(R1);
197            return R;
198        }
199        iv1 = R1.get(R1.size() - 1); // last
200        iv2 = R2.get(0); // first
201        if (iv1.right.compareTo(iv2.left) < 0) {
202            R.addAll(R1);
203            R.addAll(R2);
204            return R;
205        }
206        // now iv1.right == iv2.left
207        //System.out.println("iv1 = " + iv1);
208        //System.out.println("iv2 = " + iv2);
209        R1.remove(iv1);
210        R2.remove(iv2);
211        while (iv1.right.equals(iv2.left)) {
212            C d1 = bisectionPoint(iv1, f);
213            C d2 = bisectionPoint(iv2, f);
214            Interval<C> iv11 = new Interval<C>(iv1.left, d1);
215            Interval<C> iv12 = new Interval<C>(d1, iv1.right);
216            Interval<C> iv21 = new Interval<C>(iv2.left, d2);
217            Interval<C> iv22 = new Interval<C>(d2, iv2.right);
218
219            boolean b11 = signChange(iv11, f);
220            boolean b12 = signChange(iv12, f); // TODO check unnecessary
221            //boolean b21 = signChange(iv21, f); // TODO check unused or unnecessary
222            boolean b22 = signChange(iv22, f);
223            if (b11) {
224                iv1 = iv11;
225                if (b22) {
226                    iv2 = iv22;
227                } else {
228                    iv2 = iv21;
229                }
230                break; // done, refine
231            }
232            if (b22) {
233                iv2 = iv22;
234                if (b12) {
235                    iv1 = iv12;
236                } else {
237                    iv1 = iv11;
238                }
239                break; // done, refine
240            }
241            iv1 = iv12;
242            iv2 = iv21;
243            //System.out.println("iv1 = " + iv1);
244            //System.out.println("iv2 = " + iv2);
245        }
246        R.addAll(R1);
247        R.add(iv1);
248        R.add(iv2);
249        R.addAll(R2);
250        return R;
251    }
252
253
254    /**
255     * Number of real roots in interval.
256     * @param iv interval with f(left) * f(right) != 0.
257     * @param S sturm sequence for f and I.
258     * @return number of real roots of f in I.
259     */
260    public long realRootCount(Interval<C> iv, List<GenPolynomial<C>> S) {
261        // check sign variations at interval bounds
262        GenPolynomial<C> f = S.get(0); // squarefree part
263        //System.out.println("iv = " + iv);
264        RingFactory<C> cfac = f.ring.coFac;
265        List<C> l = PolyUtil.<C> evaluateMain(cfac, S, iv.left);
266        List<C> r = PolyUtil.<C> evaluateMain(cfac, S, iv.right);
267        long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r);
268        //System.out.println("v = " + v);
269        if (v < 0L) {
270            v = -v;
271        }
272        return v;
273    }
274
275
276    /**
277     * Number of real roots in interval.
278     * @param iv interval with f(left) * f(right) != 0.
279     * @param f univariate polynomial.
280     * @return number of real roots of f in I.
281     */
282    @Override
283    public long realRootCount(Interval<C> iv, GenPolynomial<C> f) {
284        if (f == null || f.isConstant()) { // ? 
285            return 0L;
286        }
287        if (f.isZERO()) {
288            C z = f.leadingBaseCoefficient();
289            if (!iv.contains(z)) {
290                return 0L;
291            }
292            return 1L;
293        }
294        List<GenPolynomial<C>> S = sturmSequence(f);
295        return realRootCount(iv, S);
296    }
297
298
299    /**
300     * Invariant interval for algebraic number sign.
301     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
302     * @param f univariate polynomial, non-zero.
303     * @param g univariate polynomial, gcd(f,g) == 1.
304     * @return v with v a new interval contained in iv such that g(w) != 0 for w
305     *         in v.
306     */
307    @Override
308    public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
309        Interval<C> v = iv;
310        if (g == null || g.isZERO()) {
311            //throw new IllegalArgumentException("g == 0");
312            return v;
313        }
314        if (g.isConstant()) {
315            return v;
316        }
317        if (f == null || f.isZERO()) { // ? || f.isConstant()
318            throw new IllegalArgumentException("f == 0");
319            //return v;
320        }
321        List<GenPolynomial<C>> Sg = sturmSequence(g.monic());
322        Interval<C> ivp = invariantSignInterval(iv, f, Sg);
323        return ivp;
324    }
325
326
327    /**
328     * Invariant interval for algebraic number sign.
329     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
330     * @param f univariate polynomial, non-zero.
331     * @param Sg Sturm sequence for (f,g), a univariate polynomial with gcd(f,g) ==
332     *            1.
333     * @return v with v a new interval contained in iv such that g(w) != 0 for w
334     *         in v.
335     */
336    public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, List<GenPolynomial<C>> Sg) {
337        Interval<C> v = iv;
338        GenPolynomial<C> g = Sg.get(0);
339        if (g == null || g.isZERO()) {
340            return v;
341        }
342        if (g.isConstant()) {
343            return v;
344        }
345        if (f == null || f.isZERO()) { // ? || f.isConstant()
346            return v;
347        }
348        RingFactory<C> cfac = f.ring.coFac;
349        C two = cfac.fromInteger(2);
350
351        while (true) {
352            long n = realRootCount(v, Sg);
353            logger.debug("n = {}", n);
354            if (n == 0) {
355                return v;
356            }
357            C c = v.left.sum(v.right);
358            c = c.divide(two);
359            Interval<C> im = new Interval<C>(c, v.right);
360            if (signChange(im, f)) {
361                v = im;
362            } else {
363                v = new Interval<C>(v.left, c);
364            }
365        }
366        // return v;
367    }
368
369
370    /**
371     * Exclude zero, old version.
372     * @param iv root isolating interval with f(left) * f(right) &lt; 0.
373     * @param S sturm sequence for f and I.
374     * @return a new interval v such that v &lt; 0 or v &gt; 0.
375     */
376    public Interval<C> excludeZeroOld(Interval<C> iv, List<GenPolynomial<C>> S) {
377        if (S == null || S.isEmpty()) {
378            return iv;
379        }
380        C zero = S.get(0).ring.coFac.getZERO();
381        if (!iv.contains(zero)) {
382            return iv;
383        }
384        Interval<C> vn = new Interval<C>(iv.left, zero);
385        if (realRootCount(vn,S) == 1) {
386            return vn;
387        }
388        vn = new Interval<C>(zero, iv.right);
389        return vn;
390    }
391
392
393    /**
394     * Exclude zero v2.
395     * @param iv root isolating interval with f(left) * f(right) &lt; 0.
396     * @param S sturm sequence for f and I.
397     * @return a new interval v such that v &lt; 0 or v &gt; 0 or v == 0.
398     */
399    public Interval<C> excludeZero(Interval<C> iv, List<GenPolynomial<C>> S) {
400        if (S == null || S.isEmpty()) {
401            return iv;
402        }
403        GenPolynomial<C> f = S.get(0);
404        C zero = f.ring.coFac.getZERO();
405        if (!iv.contains(zero)) { // left <= 0 <= right
406            return iv;
407        }
408        if (iv.left.isZERO() && iv.right.isZERO()) { // (0, 0)
409            return iv;
410        }
411        C m = realMinimalRootBound(f);
412        Interval<C> vi = iv;
413        Interval<C> vn;
414        if (vi.left.isZERO()) {
415            vi = new Interval<C>(m, vi.right);
416        } else if (vi.right.isZERO()) {
417            vi = new Interval<C>(vi.left, m.negate());
418        }
419        vn = new Interval<C>(vi.left, m.negate()); // l != 0
420        if (realRootCount(vn, S) == 1) {
421            return vn;
422        }
423        vn = new Interval<C>(m, vi.right); // r != 0
424        if (realRootCount(vn, S) == 1) {
425            return vn;
426        }
427        vn = new Interval<C>(zero);
428        logger.warn("interval is zero: iv = {}, trail = {}, vi = {}, vn = {}", iv, f.trailingExpVector().degree(),
429                     vi, vn);
430        return vn;
431    }
432
433}