Source code

001/*
002 * $Id: RootFactory.java 5828 2018-05-13 15:52:01Z kredel $
003 */
004
005package edu.jas.root;
006
007
008import java.util.ArrayList;
009import java.util.LinkedList;
010import java.util.List;
011import java.util.Map;
012
013import edu.jas.arith.BigDecimal;
014import edu.jas.arith.BigRational;
015import edu.jas.arith.Rational;
016import edu.jas.poly.Complex;
017import edu.jas.poly.ComplexRing;
018import edu.jas.poly.GenPolynomial;
019import edu.jas.poly.GenPolynomialRing;
020import edu.jas.poly.PolyUtil;
021import edu.jas.structure.GcdRingElem;
022import edu.jas.ufd.FactorAbstract;
023import edu.jas.ufd.FactorFactory;
024import edu.jas.ufd.SquarefreeAbstract;
025import edu.jas.ufd.SquarefreeFactory;
026
027
028/**
029 * Roots factory. Generate real and complex algebraic numbers from root
030 * intervals or rectangles.
031 * @author Heinz Kredel
032 */
033public class RootFactory {
034
035
036    /**
037     * Is real algebraic number a root of a polynomial.
038     * @param f univariate polynomial.
039     * @param r real algebraic number.
040     * @return true, if f(r) == 0, else false;
041     */
042    public static <C extends GcdRingElem<C> & Rational> boolean isRoot(GenPolynomial<C> f,
043                    RealAlgebraicNumber<C> r) {
044        RealAlgebraicRing<C> rr = r.factory();
045        GenPolynomialRing<RealAlgebraicNumber<C>> rfac = new GenPolynomialRing<RealAlgebraicNumber<C>>(rr,
046                        f.factory());
047        GenPolynomial<RealAlgebraicNumber<C>> p;
048        p = PolyUtilRoot.<C> convertToRealCoefficients(rfac, f);
049        RealAlgebraicNumber<C> a = PolyUtil.<RealAlgebraicNumber<C>> evaluateMain(rr, p, r);
050        return a.isZERO();
051    }
052
053
054    /**
055     * Real algebraic numbers.
056     * @param f univariate polynomial.
057     * @return a list of different real algebraic numbers.
058     */
059    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbers(
060                    GenPolynomial<C> f) {
061        RealRoots<C> rr = new RealRootsSturm<C>();
062        SquarefreeAbstract<C> engine = SquarefreeFactory.<C> getImplementation(f.ring.coFac);
063        Map<GenPolynomial<C>, Long> SF = engine.squarefreeFactors(f);
064        //Set<GenPolynomial<C>> S = SF.keySet();
065        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
066        for (Map.Entry<GenPolynomial<C>, Long> me : SF.entrySet()) {
067            GenPolynomial<C> sp = me.getKey();
068            List<Interval<C>> iv = rr.realRoots(sp);
069            for (Interval<C> I : iv) {
070                RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(sp, I);
071                RealAlgebraicNumber<C> rn = rar.getGenerator();
072                long mult = me.getValue();
073                for (int i = 0; i < mult; i++) {
074                    list.add(rn);
075                }
076            }
077        }
078        return list;
079    }
080
081
082    /**
083     * Real algebraic numbers.
084     * @param f univariate polynomial.
085     * @param eps rational precision.
086     * @return a list of different real algebraic numbers.
087     */
088    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbers(
089                    GenPolynomial<C> f, BigRational eps) {
090        RealRoots<C> rr = new RealRootsSturm<C>();
091        SquarefreeAbstract<C> engine = SquarefreeFactory.<C> getImplementation(f.ring.coFac);
092        Map<GenPolynomial<C>, Long> SF = engine.squarefreeFactors(f);
093        //Set<GenPolynomial<C>> S = SF.keySet();
094        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
095        for (Map.Entry<GenPolynomial<C>, Long> me : SF.entrySet()) {
096            GenPolynomial<C> sp = me.getKey();
097            List<Interval<C>> iv = rr.realRoots(sp, eps);
098            for (Interval<C> I : iv) {
099                RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(sp, I);
100                rar.setEps(eps);
101                RealAlgebraicNumber<C> rn = rar.getGenerator();
102                long mult = me.getValue();
103                for (int i = 0; i < mult; i++) {
104                    list.add(rn);
105                }
106            }
107        }
108        return list;
109    }
110
111
112    /**
113     * Real algebraic numbers from a field.
114     * @param f univariate polynomial.
115     * @return a list of different real algebraic numbers from a field.
116     */
117    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbersField(
118                    GenPolynomial<C> f) {
119        RealRoots<C> rr = new RealRootsSturm<C>();
120        FactorAbstract<C> engine = FactorFactory.<C> getImplementation(f.ring.coFac);
121        Map<GenPolynomial<C>, Long> SF = engine.baseFactors(f);
122        //Set<GenPolynomial<C>> S = SF.keySet();
123        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
124        for (Map.Entry<GenPolynomial<C>, Long> me : SF.entrySet()) {
125            GenPolynomial<C> sp = me.getKey();
126            List<Interval<C>> iv = rr.realRoots(sp);
127            for (Interval<C> I : iv) {
128                RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(sp, I, true);//field
129                RealAlgebraicNumber<C> rn = rar.getGenerator();
130                long mult = me.getValue();
131                for (int i = 0; i < mult; i++) {
132                    list.add(rn);
133                }
134            }
135        }
136        return list;
137    }
138
139
140    /**
141     * Real algebraic numbers from a field.
142     * @param f univariate polynomial.
143     * @param eps rational precision.
144     * @return a list of different real algebraic numbers from a field.
145     */
146    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbersField(
147                    GenPolynomial<C> f, BigRational eps) {
148        RealRoots<C> rr = new RealRootsSturm<C>();
149        FactorAbstract<C> engine = FactorFactory.<C> getImplementation(f.ring.coFac);
150        Map<GenPolynomial<C>, Long> SF = engine.baseFactors(f);
151        //Set<GenPolynomial<C>> S = SF.keySet();
152        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
153        for (Map.Entry<GenPolynomial<C>, Long> me : SF.entrySet()) {
154            GenPolynomial<C> sp = me.getKey();
155            List<Interval<C>> iv = rr.realRoots(sp, eps);
156            for (Interval<C> I : iv) {
157                RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(sp, I, true);//field
158                rar.setEps(eps);
159                RealAlgebraicNumber<C> rn = rar.getGenerator();
160                long mult = me.getValue();
161                for (int i = 0; i < mult; i++) {
162                    list.add(rn);
163                }
164            }
165        }
166        return list;
167    }
168
169
170    /**
171     * Real algebraic numbers from a irreducible polynomial.
172     * @param f univariate irreducible polynomial.
173     * @return a list of different real algebraic numbers from a field.
174     */
175    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbersIrred(
176                    GenPolynomial<C> f) {
177        RealRoots<C> rr = new RealRootsSturm<C>();
178        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
179        List<Interval<C>> iv = rr.realRoots(f);
180        for (Interval<C> I : iv) {
181            RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(f, I, true);//field
182            RealAlgebraicNumber<C> rn = rar.getGenerator();
183            list.add(rn);
184        }
185        return list;
186    }
187
188
189    /**
190     * Real algebraic numbers from a irreducible polynomial.
191     * @param f univariate irreducible polynomial.
192     * @param eps rational precision.
193     * @return a list of different real algebraic numbers from a field.
194     */
195    public static <C extends GcdRingElem<C> & Rational> List<RealAlgebraicNumber<C>> realAlgebraicNumbersIrred(
196                    GenPolynomial<C> f, BigRational eps) {
197        RealRoots<C> rr = new RealRootsSturm<C>();
198        List<RealAlgebraicNumber<C>> list = new ArrayList<RealAlgebraicNumber<C>>();
199        List<Interval<C>> iv = rr.realRoots(f, eps);
200        for (Interval<C> I : iv) {
201            RealAlgebraicRing<C> rar = new RealAlgebraicRing<C>(f, I, true);//field
202            rar.setEps(eps);
203            RealAlgebraicNumber<C> rn = rar.getGenerator();
204            list.add(rn);
205        }
206        return list;
207    }
208
209
210    /**
211     * Is complex algebraic number a root of a polynomial.
212     * @param f univariate polynomial.
213     * @param r complex algebraic number.
214     * @return true, if f(r) == 0, else false;
215     */
216    public static <C extends GcdRingElem<C> & Rational> boolean isRoot(GenPolynomial<C> f,
217                    ComplexAlgebraicNumber<C> r) {
218        ComplexAlgebraicRing<C> cr = r.factory();
219        GenPolynomialRing<ComplexAlgebraicNumber<C>> cfac = new GenPolynomialRing<ComplexAlgebraicNumber<C>>(
220                        cr, f.factory());
221        GenPolynomial<ComplexAlgebraicNumber<C>> p;
222        p = PolyUtilRoot.<C> convertToComplexCoefficients(cfac, f);
223        ComplexAlgebraicNumber<C> a = PolyUtil.<ComplexAlgebraicNumber<C>> evaluateMain(cr, p, r);
224        return a.isZERO();
225    }
226
227
228    /**
229     * Is complex algebraic number a root of a complex polynomial.
230     * @param f univariate complex polynomial.
231     * @param r complex algebraic number.
232     * @return true, if f(r) == 0, else false;
233     */
234    public static <C extends GcdRingElem<C> & Rational> boolean isRootComplex(GenPolynomial<Complex<C>> f,
235                    ComplexAlgebraicNumber<C> r) {
236        ComplexAlgebraicRing<C> cr = r.factory();
237        GenPolynomialRing<ComplexAlgebraicNumber<C>> cfac = new GenPolynomialRing<ComplexAlgebraicNumber<C>>(
238                        cr, f.factory());
239        GenPolynomial<ComplexAlgebraicNumber<C>> p;
240        p = PolyUtilRoot.<C> convertToComplexCoefficientsFromComplex(cfac, f);
241        ComplexAlgebraicNumber<C> a = PolyUtil.<ComplexAlgebraicNumber<C>> evaluateMain(cr, p, r);
242        return a.isZERO();
243    }
244
245
246    /**
247     * Is complex algebraic number a real root of a polynomial.
248     * @param f univariate polynomial.
249     * @param c complex algebraic number.
250     * @param r real algebraic number.
251     * @return true, if f(c) == 0 and c == r, else false;
252     */
253    public static <C extends GcdRingElem<C> & Rational> boolean isRealRoot(GenPolynomial<C> f,
254                    ComplexAlgebraicNumber<C> c, RealAlgebraicNumber<C> r) {
255        boolean t = isRoot(f, c) && isRoot(f, r);
256        if (!t) {
257            return t;
258        }
259        Rectangle<C> rc = c.ring.root;
260        Interval<C> ivci = new Interval<C>(rc.getSW().getIm(), rc.getNE().getIm());
261        t = ivci.contains(f.ring.coFac.getZERO());
262        if (!t) {
263            return t;
264        }
265        //System.out.println("imag = 0");
266        Interval<C> ivc = new Interval<C>(rc.getSW().getRe(), rc.getNE().getRe());
267        Interval<C> ivr = r.ring.root;
268        // disjoint intervals
269        if (ivc.right.compareTo(ivr.left) < 0 || ivr.right.compareTo(ivc.left) < 0) {
270            //System.out.println("disjoint: ivc = " + ivc + ", ivr = " + ivr);
271            return false;
272        }
273        //System.out.println("not disjoint");
274        // full containement
275        t = ivc.contains(ivr) || ivr.contains(ivc);
276        if (t) {
277            return t;
278        }
279        //System.out.println("with overlap");
280        // overlap, refine to smaller interval
281        C left = ivc.left;
282        if (left.compareTo(ivr.left) > 0) {
283            left = ivr.left;
284        }
285        C right = ivc.right;
286        if (right.compareTo(ivr.right) < 0) {
287            right = ivr.right;
288        }
289        Interval<C> ref = new Interval<C>(left, right);
290        //System.out.println("refined interval " + ref);
291        RealRoots<C> reng = r.ring.engine; //new RealRootsSturm<C>(); 
292        long z = reng.realRootCount(ref, f);
293        if (z != 1) {
294            return false;
295        }
296        ComplexRing<C> cr = rc.getSW().ring;
297        Complex<C> sw = new Complex<C>(cr, left, rc.getSW().getIm());
298        Complex<C> ne = new Complex<C>(cr, right, rc.getNE().getIm());
299        //Complex<C> sw = new Complex<C>(cr, left, cr.ring.getZERO());
300        //Complex<C> ne = new Complex<C>(cr, right, cr.ring.getZERO());
301        Rectangle<C> rec = new Rectangle<C>(sw, ne);
302        //System.out.println("refined rectangle " + rec);
303        ComplexRoots<C> ceng = c.ring.engine; //new ComplexRootsSturm<C>(); 
304        GenPolynomial<Complex<C>> p = PolyUtilRoot.<C> complexFromAny(f);
305        try {
306            z = ceng.complexRootCount(rec, p);
307        } catch (InvalidBoundaryException e) {
308            System.out.println("should not happen, rec = " + rec + ", p = " + p);
309            e.printStackTrace();
310            z = 0;
311        }
312        //System.out.println("complexRootCount: " + z);
313        if (z != 1) {
314            return false;
315        }
316        return true;
317    }
318
319
320    /**
321     * Is complex decimal number a real root of a polynomial.
322     * @param f univariate polynomial.
323     * @param c complex decimal number.
324     * @param r real decimal number.
325     * @param eps desired precision.
326     * @return true, if f(c) == 0 and c == r, else false;
327     */
328    public static <C extends GcdRingElem<C> & Rational> boolean isRealRoot(GenPolynomial<C> f,
329                    Complex<BigDecimal> c, BigDecimal r, BigRational eps) {
330        BigDecimal e = new BigDecimal(eps);
331        if (c.getIm().abs().compareTo(e) > 0) {
332            return false;
333        }
334        if (c.getRe().subtract(r).abs().compareTo(e) > 0) {
335            return false;
336        }
337        GenPolynomialRing<Complex<C>> cfac = new GenPolynomialRing<Complex<C>>(
338                        new ComplexRing<C>(f.ring.coFac), f.ring);
339        GenPolynomial<Complex<C>> cf = PolyUtil.<C> complexFromAny(cfac, f);
340        ComplexRing<BigDecimal> cd = new ComplexRing<BigDecimal>(e);
341        GenPolynomialRing<Complex<BigDecimal>> rfac = new GenPolynomialRing<Complex<BigDecimal>>(cd, cf.ring);
342        GenPolynomial<Complex<BigDecimal>> cdf = PolyUtil.<C> complexDecimalFromRational(rfac, cf);
343        Complex<BigDecimal> z = PolyUtil.evaluateMain(cd, cdf, c);
344        if (z.isZERO()) {
345            return true;
346        }
347        System.out.println("z != 0: " + z);
348        return false;
349    }
350
351
352    /**
353     * Complex algebraic numbers.
354     * @param f univariate polynomial.
355     * @return a list of different complex algebraic numbers.
356     */
357    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(
358                    GenPolynomial<Complex<C>> f) {
359        ComplexRoots<C> cr = new ComplexRootsSturm<C>(f.ring.coFac);
360        SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory
361                        .<Complex<C>> getImplementation(f.ring.coFac);
362        Map<GenPolynomial<Complex<C>>, Long> SF = engine.squarefreeFactors(f);
363        //Set<GenPolynomial<Complex<C>>> S = SF.keySet();
364        List<ComplexAlgebraicNumber<C>> list = new ArrayList<ComplexAlgebraicNumber<C>>();
365        for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : SF.entrySet()) {
366            GenPolynomial<Complex<C>> sp = me.getKey();
367            List<Rectangle<C>> iv = cr.complexRoots(sp);
368            for (Rectangle<C> I : iv) {
369                ComplexAlgebraicRing<C> car = new ComplexAlgebraicRing<C>(sp, I);
370                ComplexAlgebraicNumber<C> cn = car.getGenerator();
371                long mult = me.getValue();
372                for (int i = 0; i < mult; i++) {
373                    list.add(cn);
374                }
375            }
376        }
377        return list;
378    }
379
380
381    /**
382     * Complex algebraic numbers.
383     * @param f univariate polynomial.
384     * @param eps rational precision.
385     * @return a list of different complex algebraic numbers.
386     */
387    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(
388                    GenPolynomial<Complex<C>> f, BigRational eps) {
389        ComplexRoots<C> cr = new ComplexRootsSturm<C>(f.ring.coFac);
390        SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory
391                        .<Complex<C>> getImplementation(f.ring.coFac);
392        Map<GenPolynomial<Complex<C>>, Long> SF = engine.squarefreeFactors(f);
393        //Set<GenPolynomial<Complex<C>>> S = SF.keySet();
394        List<ComplexAlgebraicNumber<C>> list = new ArrayList<ComplexAlgebraicNumber<C>>();
395        for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : SF.entrySet()) {
396            GenPolynomial<Complex<C>> sp = me.getKey();
397            List<Rectangle<C>> iv = cr.complexRoots(sp);
398            for (Rectangle<C> I : iv) {
399                Rectangle<C> Iv = I;
400                try {
401                    Iv = cr.complexRootRefinement(I, sp, eps);
402                } catch (InvalidBoundaryException e) {
403                    e.printStackTrace();
404                }
405                ComplexAlgebraicRing<C> car = new ComplexAlgebraicRing<C>(sp, Iv);
406                car.setEps(eps);
407                ComplexAlgebraicNumber<C> cn = car.getGenerator();
408                long mult = me.getValue();
409                for (int i = 0; i < mult; i++) {
410                    list.add(cn);
411                }
412            }
413        }
414        return list;
415    }
416
417
418    /**
419     * Complex algebraic numbers.
420     * @param f univariate (rational) polynomial.
421     * @return a list of different complex algebraic numbers.
422     */
423    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbers(
424                    GenPolynomial<C> f) {
425        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
426        return complexAlgebraicNumbersComplex(fc);
427    }
428
429
430    /**
431     * Complex algebraic numbers.
432     * @param f univariate (rational) polynomial.
433     * @param eps rational precision.
434     * @return a list of different complex algebraic numbers.
435     */
436    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbers(
437                    GenPolynomial<C> f, BigRational eps) {
438        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
439        return complexAlgebraicNumbersComplex(fc, eps);
440    }
441
442
443    /**
444     * Filter real roots from complex roots.
445     * @param f univariate polynomial.
446     * @param c list of complex algebraic numbers.
447     * @param r list of real algebraic numbers.
448     * @return c minus the real roots from r
449     */
450    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> filterOutRealRoots(
451                    GenPolynomial<C> f, List<ComplexAlgebraicNumber<C>> c, List<RealAlgebraicNumber<C>> r) {
452        if (c.isEmpty()) {
453            return c;
454        }
455        if (r.isEmpty()) {
456            return c;
457        }
458        List<ComplexAlgebraicNumber<C>> cmr = new ArrayList<ComplexAlgebraicNumber<C>>();
459        if (c.size() == r.size() /*&& r.size() == f.degree()*/ ) {
460            return cmr;
461        }
462        List<RealAlgebraicNumber<C>> rl = new LinkedList<RealAlgebraicNumber<C>>(r);
463        for (ComplexAlgebraicNumber<C> cn : c) {
464            RealAlgebraicNumber<C> rm = null; // ~boolean
465            for (RealAlgebraicNumber<C> rn : rl) {
466                if (isRealRoot(f, cn, rn)) {
467                    //System.out.println("filterOutRealRoots, cn = " + cn + ", rn = " + rn);
468                    rm = rn;
469                    break; // remove from r
470                }
471            }
472            if (rm == null) {
473                cmr.add(cn);
474            } else {
475                rl.remove(rm);
476            }
477        }
478        return cmr;
479    }
480
481
482    /**
483     * Filter real roots from complex roots.
484     * @param f univariate polynomial.
485     * @param c list of complex decimal numbers.
486     * @param r list of real decimal numbers.
487     * @param eps desired precision.
488     * @return c minus the real roots from r
489     */
490    public static <C extends GcdRingElem<C> & Rational> List<Complex<BigDecimal>> filterOutRealRoots(
491                    GenPolynomial<C> f, List<Complex<BigDecimal>> c, List<BigDecimal> r, BigRational eps) {
492        if (c.isEmpty()) {
493            return c;
494        }
495        if (r.isEmpty()) {
496            return c;
497        }
498        List<Complex<BigDecimal>> cmr = new ArrayList<Complex<BigDecimal>>();
499        if (c.size() == r.size() /*&& r.size() == f.degree()*/ ) {
500            return cmr;
501        }
502        List<BigDecimal> rl = new LinkedList<BigDecimal>(r);
503        for (Complex<BigDecimal> cn : c) {
504            BigDecimal rm = null; // ~boolean
505            for (BigDecimal rn : rl) {
506                if (isRealRoot(f, cn, rn, eps)) {
507                    //System.out.println("filterOutRealRoots, cn = " + cn + ", rn = " + rn);
508                    rm = rn;
509                    break; // remove from r
510                }
511            }
512            if (rm == null) {
513                cmr.add(cn);
514            } else {
515                rl.remove(rm);
516            }
517        }
518        return cmr;
519    }
520
521
522    /**
523     * Roots as real and complex algebraic numbers.
524     * @param f univariate polynomial.
525     * @return container of real and complex algebraic numbers.
526     */
527    public static <C extends GcdRingElem<C> & Rational> AlgebraicRoots<C> algebraicRoots(GenPolynomial<C> f) {
528        List<RealAlgebraicNumber<C>> rl = realAlgebraicNumbers(f);
529        List<ComplexAlgebraicNumber<C>> cl = complexAlgebraicNumbers(f);
530        GenPolynomial<Complex<C>> cf = null;
531        if (!cl.isEmpty()) {
532            cf = cl.get(0).ring.algebraic.modul;
533        }
534        cl = filterOutRealRoots(f, cl, rl);
535        AlgebraicRoots<C> ar = new AlgebraicRoots<C>(f, cf, rl, cl);
536        return ar;
537    }
538
539
540    /**
541     * Roots of unity of real and complex algebraic numbers.
542     * @param ar container of real and complex algebraic numbers.
543     * @return container of real and complex algebraic numbers which are roots
544     *         of unity.
545     */
546    public static <C extends GcdRingElem<C> & Rational> AlgebraicRoots<C> rootsOfUnity(AlgebraicRoots<C> ar) {
547        List<RealAlgebraicNumber<C>> rl = new ArrayList<RealAlgebraicNumber<C>>();
548        for (RealAlgebraicNumber<C> r : ar.real) {
549            if (r.isRootOfUnity()) {
550                rl.add(r);
551            }
552        }
553        List<ComplexAlgebraicNumber<C>> cl = new ArrayList<ComplexAlgebraicNumber<C>>();
554        for (ComplexAlgebraicNumber<C> c : ar.complex) {
555            if (c.isRootOfUnity()) {
556                cl.add(c);
557            }
558        }
559        AlgebraicRoots<C> ur = new AlgebraicRoots<C>(ar.p, ar.cp, rl, cl);
560        return ur;
561    }
562
563
564    /**
565     * Root refinement of real and complex algebraic numbers.
566     * @param a container of real and complex algebraic numbers.
567     * @param eps desired precision for root intervals and rectangles.
568     */
569    public static <C extends GcdRingElem<C> & Rational> void rootRefine(AlgebraicRoots<C> a,
570                    BigRational eps) {
571        for (RealAlgebraicNumber<C> r : a.real) {
572            r.ring.refineRoot(eps);
573        }
574        for (ComplexAlgebraicNumber<C> c : a.complex) {
575            c.ring.refineRoot(eps);
576        }
577        return; // a or void?
578    }
579
580
581    /**
582     * Roots as real and complex decimal numbers.
583     * @param f univariate polynomial.
584     * @param eps desired precision.
585     * @return container of real and complex decimal numbers.
586     */
587    public static <C extends GcdRingElem<C> & Rational> DecimalRoots<C> decimalRoots(GenPolynomial<C> f,
588                    BigRational eps) {
589        RealRootsAbstract<C> rengine = new RealRootsSturm<C>();
590        List<BigDecimal> rl = rengine.approximateRoots(f, eps);
591
592        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
593        ComplexRootsAbstract<C> cengine = new ComplexRootsSturm<C>(fc.ring.coFac);
594        List<Complex<BigDecimal>> cl = cengine.approximateRoots(fc, eps);
595
596        cl = filterOutRealRoots(f, cl, rl, eps);
597        DecimalRoots<C> ar = new DecimalRoots<C>(f, fc, rl, cl);
598        return ar;
599    }
600
601
602    /**
603     * Roots as real and complex decimal numbers.
604     * @param ar container for real and complex algebraic roots.
605     * @param eps desired precision.
606     * @return container of real and complex decimal numbers.
607     */
608    public static <C extends GcdRingElem<C> & Rational> DecimalRoots<C> decimalRoots(AlgebraicRoots<C> ar,
609                    BigRational eps) {
610        //no: rootRefine(ar, eps);
611        RealRootsAbstract<C> rengine = new RealRootsSturm<C>();
612        List<BigDecimal> rl = new ArrayList<BigDecimal>(ar.real.size());
613        for (RealAlgebraicNumber<C> r : ar.real) {
614            try {
615                BigDecimal d = rengine.approximateRoot(r.ring.root, ar.p, eps);
616                rl.add(d);
617            } catch (NoConvergenceException e) {
618                System.out.println("should not happen: " + e);
619            }
620        }
621        ComplexRootsAbstract<C> cengine = new ComplexRootsSturm<C>(ar.cp.ring.coFac);
622        List<Complex<BigDecimal>> cl = new ArrayList<Complex<BigDecimal>>(ar.complex.size());
623        for (ComplexAlgebraicNumber<C> c : ar.complex) {
624            try {
625                Complex<BigDecimal> d = cengine.approximateRoot(c.ring.root, ar.cp, eps);
626                cl.add(d);
627            } catch (NoConvergenceException e) {
628                System.out.println("should not happen: " + e);
629            }
630        }
631        DecimalRoots<C> dr = new DecimalRoots<C>(ar.p, ar.cp, rl, cl);
632        return dr;
633    }
634
635}