Source code

001/*
002 * $Id: RealRootsSturm.java 4111 2012-08-19 12:30:30Z kredel $
003 */
004
005package edu.jas.root;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.log4j.Logger;
012
013import edu.jas.arith.BigRational;
014import edu.jas.arith.Rational;
015import edu.jas.poly.GenPolynomial;
016import edu.jas.poly.PolyUtil;
017import edu.jas.structure.RingElem;
018import edu.jas.structure.RingFactory;
019
020
021/**
022 * Real root isolation using Sturm sequences.
023 * @param <C> coefficient type.
024 * @author Heinz Kredel
025 */
026public class RealRootsSturm<C extends RingElem<C> & Rational> extends RealRootsAbstract<C> {
027
028
029    private static final Logger logger = Logger.getLogger(RealRootsSturm.class);
030
031
032    private static boolean debug = logger.isDebugEnabled();
033
034
035    /**
036     * Sturm sequence.
037     * @param f univariate polynomial.
038     * @return a Sturm sequence for f.
039     */
040    public List<GenPolynomial<C>> sturmSequence(GenPolynomial<C> f) {
041        List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>();
042        if (f == null || f.isZERO()) {
043            return S;
044        }
045        if (f.isConstant()) {
046            S.add(f.monic());
047            return S;
048        }
049        GenPolynomial<C> F = f;
050        S.add(F);
051        GenPolynomial<C> G = PolyUtil.<C> baseDeriviative(f);
052        while (!G.isZERO()) {
053            GenPolynomial<C> r = F.remainder(G);
054            F = G;
055            G = r.negate();
056            S.add(F/*.monic()*/);
057        }
058        //System.out.println("F = " + F);
059        if (F.isConstant()) {
060            return S;
061        }
062        // make squarefree
063        List<GenPolynomial<C>> Sp = new ArrayList<GenPolynomial<C>>(S.size());
064        for (GenPolynomial<C> p : S) {
065            p = p.divide(F);
066            Sp.add(p);
067        }
068        return Sp;
069    }
070
071
072    /**
073     * Isolating intervals for the real roots.
074     * @param f univariate polynomial.
075     * @return a list of isolating intervals for the real roots of f.
076     */
077    @Override
078    public List<Interval<C>> realRoots(GenPolynomial<C> f) {
079        List<Interval<C>> R = new ArrayList<Interval<C>>();
080        if (f == null || f.isConstant()) {
081            return R;
082        }
083        if (f.isZERO()) {
084            C z = f.ring.coFac.getZERO();
085            R.add(new Interval<C>(z));
086            return R;
087        }
088        if (f.degree(0) == 1L) {
089            C z = f.monic().trailingBaseCoefficient().negate();
090            R.add(new Interval<C>(z));
091            return R;
092        }
093        GenPolynomial<C> F = f;
094        C M = realRootBound(F); // M != 0, since >= 2
095        Interval<C> iv = new Interval<C>(M.negate(), M);
096        //System.out.println("iv = " + iv);
097        List<GenPolynomial<C>> S = sturmSequence(F);
098        //System.out.println("S = " + S);
099        //System.out.println("f_S = " + S.get(0));
100        List<Interval<C>> Rp = realRoots(iv, S);
101        if (logger.isInfoEnabled() && !(((Object) f.ring.coFac) instanceof BigRational)) {
102            //logger.info("realRoots bound: " + iv);
103            logger.info("realRoots: " + Rp);
104        }
105        R.addAll(Rp);
106        return R;
107    }
108
109
110    /**
111     * Isolating intervals for the real roots.
112     * @param iv interval with f(left) * f(right) != 0.
113     * @param S sturm sequence for f and I.
114     * @return a list of isolating intervals for the real roots of f in I.
115     */
116    public List<Interval<C>> realRoots(Interval<C> iv, List<GenPolynomial<C>> S) {
117        List<Interval<C>> R = new ArrayList<Interval<C>>();
118        GenPolynomial<C> f = S.get(0); // squarefree part
119        if (f.isZERO()) {
120            C z = f.leadingBaseCoefficient();
121            if (!iv.contains(z)) {
122                throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv
123                                + ", z = " + z);
124            }
125            Interval<C> iv1 = new Interval<C>(z);
126            R.add(iv1);
127            return R;
128        }
129        if (f.isConstant()) {
130            return R;
131            //throw new IllegalArgumentException("f has no root: f = " + f + ", iv = " + iv);
132        }
133        if (f.degree(0) == 1L) {
134            C z = f.monic().trailingBaseCoefficient().negate();
135            if (!iv.contains(z)) {
136                return R;
137                //throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv + ", z = " + z);
138            }
139            Interval<C> iv1 = new Interval<C>(z);
140            R.add(iv1);
141            return R;
142        }
143        //System.out.println("iv = " + iv);
144        // check sign variations at interval bounds
145        long v = realRootCount(iv, S);
146        //System.out.println("v = " + v);
147        if (v == 0) {
148            return R;
149        }
150        if (v == 1) {
151            R.add(iv);
152            return R;
153        }
154        // now v &gt; 1
155        // bi-sect interval, such that f(c) != 0
156        C c = bisectionPoint(iv, f);
157        //System.out.println("c = " + c);
158        // recursion on both sub-intervals
159        Interval<C> iv1 = new Interval<C>(iv.left, c);
160        Interval<C> iv2 = new Interval<C>(c, iv.right);
161        List<Interval<C>> R1 = realRoots(iv1, S);
162        //System.out.println("R1 = " + R1);
163        if (debug) {
164            logger.info("R1 = " + R1);
165        }
166        List<Interval<C>> R2 = realRoots(iv2, S);
167        //System.out.println("R2 = " + R2);
168        if (debug) {
169            logger.info("R2 = " + R2);
170        }
171
172        // refine isolating intervals if adjacent 
173        if (R1.isEmpty()) {
174            R.addAll(R2);
175            return R;
176        }
177        if (R2.isEmpty()) {
178            R.addAll(R1);
179            return R;
180        }
181        iv1 = R1.get(R1.size() - 1); // last
182        iv2 = R2.get(0); // first
183        if (iv1.right.compareTo(iv2.left) < 0) {
184            R.addAll(R1);
185            R.addAll(R2);
186            return R;
187        }
188        // now iv1.right == iv2.left
189        //System.out.println("iv1 = " + iv1);
190        //System.out.println("iv2 = " + iv2);
191        R1.remove(iv1);
192        R2.remove(iv2);
193        while (iv1.right.equals(iv2.left)) {
194            C d1 = bisectionPoint(iv1, f);
195            C d2 = bisectionPoint(iv2, f);
196            Interval<C> iv11 = new Interval<C>(iv1.left, d1);
197            Interval<C> iv12 = new Interval<C>(d1, iv1.right);
198            Interval<C> iv21 = new Interval<C>(iv2.left, d2);
199            Interval<C> iv22 = new Interval<C>(d2, iv2.right);
200
201            boolean b11 = signChange(iv11, f);
202            boolean b12 = signChange(iv12, f); // TODO check unnecessary
203            //boolean b21 = signChange(iv21, f); // TODO check unused or unnecessary
204            boolean b22 = signChange(iv22, f);
205            if (b11) {
206                iv1 = iv11;
207                if (b22) {
208                    iv2 = iv22;
209                } else {
210                    iv2 = iv21;
211                }
212                break; // done, refine
213            }
214            if (b22) {
215                iv2 = iv22;
216                if (b12) {
217                    iv1 = iv12;
218                } else {
219                    iv1 = iv11;
220                }
221                break; // done, refine
222            }
223            iv1 = iv12;
224            iv2 = iv21;
225            //System.out.println("iv1 = " + iv1);
226            //System.out.println("iv2 = " + iv2);
227        }
228        R.addAll(R1);
229        R.add(iv1);
230        R.add(iv2);
231        R.addAll(R2);
232        return R;
233    }
234
235
236    /**
237     * Number of real roots in interval.
238     * @param iv interval with f(left) * f(right) != 0.
239     * @param S sturm sequence for f and I.
240     * @return number of real roots of f in I.
241     */
242    public long realRootCount(Interval<C> iv, List<GenPolynomial<C>> S) {
243        // check sign variations at interval bounds
244        GenPolynomial<C> f = S.get(0); // squarefree part
245        //System.out.println("iv = " + iv);
246        RingFactory<C> cfac = f.ring.coFac;
247        List<C> l = PolyUtil.<C> evaluateMain(cfac, S, iv.left);
248        List<C> r = PolyUtil.<C> evaluateMain(cfac, S, iv.right);
249        long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r);
250        //System.out.println("v = " + v);
251        if (v < 0L) {
252            v = -v;
253        }
254        return v;
255    }
256
257
258    /**
259     * Number of real roots in interval.
260     * @param iv interval with f(left) * f(right) != 0.
261     * @param f univariate polynomial.
262     * @return number of real roots of f in I.
263     */
264    @Override
265    public long realRootCount(Interval<C> iv, GenPolynomial<C> f) {
266        if (f == null || f.isConstant()) { // ? 
267            return 0L;
268        }
269        if (f.isZERO()) {
270            C z = f.leadingBaseCoefficient();
271            if (!iv.contains(z)) {
272                return 0L;
273            }
274            return 1L;
275        }
276        List<GenPolynomial<C>> S = sturmSequence(f);
277        return realRootCount(iv, S);
278    }
279
280
281    /**
282     * Invariant interval for algebraic number sign.
283     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
284     * @param f univariate polynomial, non-zero.
285     * @param g univariate polynomial, gcd(f,g) == 1.
286     * @return v with v a new interval contained in iv such that g(w) != 0 for w
287     *         in v.
288     */
289    @Override
290    public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
291        Interval<C> v = iv;
292        if (g == null || g.isZERO()) {
293            //throw new IllegalArgumentException("g == 0");
294            return v;
295        }
296        if (g.isConstant()) {
297            return v;
298        }
299        if (f == null || f.isZERO()) { // ? || f.isConstant()
300            throw new IllegalArgumentException("f == 0");
301            //return v;
302        }
303        List<GenPolynomial<C>> Sg = sturmSequence(g.monic());
304        Interval<C> ivp = invariantSignInterval(iv, f, Sg);
305        return ivp;
306    }
307
308
309    /**
310     * Invariant interval for algebraic number sign.
311     * @param iv root isolating interval for f, with f(left) * f(right) &lt; 0.
312     * @param f univariate polynomial, non-zero.
313     * @param Sg Sturm sequence for g, a univariate polynomial with gcd(f,g) ==
314     *            1.
315     * @return v with v a new interval contained in iv such that g(w) != 0 for w
316     *         in v.
317     */
318    public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, List<GenPolynomial<C>> Sg) {
319        Interval<C> v = iv;
320        GenPolynomial<C> g = Sg.get(0);
321        if (g == null || g.isZERO()) {
322            return v;
323        }
324        if (g.isConstant()) {
325            return v;
326        }
327        if (f == null || f.isZERO()) { // ? || f.isConstant()
328            return v;
329        }
330        RingFactory<C> cfac = f.ring.coFac;
331        C two = cfac.fromInteger(2);
332
333        while (true) {
334            long n = realRootCount(v, Sg);
335            logger.debug("n = " + n);
336            if (n == 0) {
337                return v;
338            }
339            C c = v.left.sum(v.right);
340            c = c.divide(two);
341            Interval<C> im = new Interval<C>(c, v.right);
342            if (signChange(im, f)) {
343                v = im;
344            } else {
345                v = new Interval<C>(v.left, c);
346            }
347        }
348        // return v;
349    }
350
351}