Source code

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