001/*
002 * $Id$
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        if (z.getIm().abs().compareTo(e) <= 0 && z.getRe().abs().compareTo(e) <= 0) {
348            return true;
349        }
350        System.out.println("z != 0: " + z);
351        return false;
352    }
353
354
355    /**
356     * Complex algebraic numbers.
357     * @param f univariate polynomial.
358     * @return a list of different complex algebraic numbers.
359     */
360    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(
361                    GenPolynomial<Complex<C>> f) {
362        ComplexRoots<C> cr = new ComplexRootsSturm<C>(f.ring.coFac);
363        SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory
364                        .<Complex<C>> getImplementation(f.ring.coFac);
365        Map<GenPolynomial<Complex<C>>, Long> SF = engine.squarefreeFactors(f);
366        //Set<GenPolynomial<Complex<C>>> S = SF.keySet();
367        List<ComplexAlgebraicNumber<C>> list = new ArrayList<ComplexAlgebraicNumber<C>>();
368        for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : SF.entrySet()) {
369            GenPolynomial<Complex<C>> sp = me.getKey();
370            List<Rectangle<C>> iv = cr.complexRoots(sp);
371            for (Rectangle<C> I : iv) {
372                ComplexAlgebraicRing<C> car = new ComplexAlgebraicRing<C>(sp, I);
373                ComplexAlgebraicNumber<C> cn = car.getGenerator();
374                long mult = me.getValue();
375                for (int i = 0; i < mult; i++) {
376                    list.add(cn);
377                }
378            }
379        }
380        return list;
381    }
382
383
384    /**
385     * Complex algebraic numbers.
386     * @param f univariate polynomial.
387     * @param eps rational precision.
388     * @return a list of different complex algebraic numbers.
389     */
390    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(
391                    GenPolynomial<Complex<C>> f, BigRational eps) {
392        ComplexRoots<C> cr = new ComplexRootsSturm<C>(f.ring.coFac);
393        SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory
394                        .<Complex<C>> getImplementation(f.ring.coFac);
395        Map<GenPolynomial<Complex<C>>, Long> SF = engine.squarefreeFactors(f);
396        //Set<GenPolynomial<Complex<C>>> S = SF.keySet();
397        List<ComplexAlgebraicNumber<C>> list = new ArrayList<ComplexAlgebraicNumber<C>>();
398        for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : SF.entrySet()) {
399            GenPolynomial<Complex<C>> sp = me.getKey();
400            List<Rectangle<C>> iv = cr.complexRoots(sp);
401            for (Rectangle<C> I : iv) {
402                Rectangle<C> Iv = I;
403                try {
404                    Iv = cr.complexRootRefinement(I, sp, eps);
405                } catch (InvalidBoundaryException e) {
406                    e.printStackTrace();
407                }
408                ComplexAlgebraicRing<C> car = new ComplexAlgebraicRing<C>(sp, Iv);
409                car.setEps(eps);
410                ComplexAlgebraicNumber<C> cn = car.getGenerator();
411                long mult = me.getValue();
412                for (int i = 0; i < mult; i++) {
413                    list.add(cn);
414                }
415            }
416        }
417        return list;
418    }
419
420
421    /**
422     * Complex algebraic numbers.
423     * @param f univariate (rational) polynomial.
424     * @return a list of different complex algebraic numbers.
425     */
426    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbers(
427                    GenPolynomial<C> f) {
428        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
429        return complexAlgebraicNumbersComplex(fc);
430    }
431
432
433    /**
434     * Complex algebraic numbers.
435     * @param f univariate (rational) polynomial.
436     * @param eps rational precision.
437     * @return a list of different complex algebraic numbers.
438     */
439    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbers(
440                    GenPolynomial<C> f, BigRational eps) {
441        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
442        return complexAlgebraicNumbersComplex(fc, eps);
443    }
444
445
446    /**
447     * Filter real roots from complex roots.
448     * @param f univariate polynomial.
449     * @param c list of complex algebraic numbers.
450     * @param r list of real algebraic numbers.
451     * @return c minus the real roots from r
452     */
453    public static <C extends GcdRingElem<C> & Rational> List<ComplexAlgebraicNumber<C>> filterOutRealRoots(
454                    GenPolynomial<C> f, List<ComplexAlgebraicNumber<C>> c, List<RealAlgebraicNumber<C>> r) {
455        if (c.isEmpty()) {
456            return c;
457        }
458        if (r.isEmpty()) {
459            return c;
460        }
461        List<ComplexAlgebraicNumber<C>> cmr = new ArrayList<ComplexAlgebraicNumber<C>>();
462        if (c.size() == r.size() /*&& r.size() == f.degree()*/ ) {
463            return cmr;
464        }
465        List<RealAlgebraicNumber<C>> rl = new LinkedList<RealAlgebraicNumber<C>>(r);
466        for (ComplexAlgebraicNumber<C> cn : c) {
467            RealAlgebraicNumber<C> rm = null; // ~boolean
468            for (RealAlgebraicNumber<C> rn : rl) {
469                if (isRealRoot(f, cn, rn)) {
470                    //System.out.println("filterOutRealRoots, cn = " + cn + ", rn = " + rn);
471                    rm = rn;
472                    break; // remove from r
473                }
474            }
475            if (rm == null) {
476                cmr.add(cn);
477            } else {
478                rl.remove(rm);
479            }
480        }
481        return cmr;
482    }
483
484
485    /**
486     * Filter real roots from complex roots.
487     * @param f univariate polynomial.
488     * @param c list of complex decimal numbers.
489     * @param r list of real decimal numbers.
490     * @param eps desired precision.
491     * @return c minus the real roots from r
492     */
493    public static <C extends GcdRingElem<C> & Rational> List<Complex<BigDecimal>> filterOutRealRoots(
494                    GenPolynomial<C> f, List<Complex<BigDecimal>> c, List<BigDecimal> r, BigRational eps) {
495        if (c.isEmpty()) {
496            return c;
497        }
498        if (r.isEmpty()) {
499            return c;
500        }
501        List<Complex<BigDecimal>> cmr = new ArrayList<Complex<BigDecimal>>();
502        if (c.size() == r.size() /*&& r.size() == f.degree()*/ ) {
503            return cmr;
504        }
505        List<BigDecimal> rl = new LinkedList<BigDecimal>(r);
506        for (Complex<BigDecimal> cn : c) {
507            BigDecimal rm = null; // ~boolean
508            for (BigDecimal rn : rl) {
509                if (isRealRoot(f, cn, rn, eps)) {
510                    //System.out.println("filterOutRealRoots, cn = " + cn + ", rn = " + rn);
511                    rm = rn;
512                    break; // remove from r
513                }
514            }
515            if (rm == null) {
516                cmr.add(cn);
517            } else {
518                rl.remove(rm);
519            }
520        }
521        return cmr;
522    }
523
524
525    /**
526     * Roots as real and complex algebraic numbers.
527     * @param f univariate polynomial.
528     * @return container of real and complex algebraic numbers.
529     */
530    public static <C extends GcdRingElem<C> & Rational> AlgebraicRoots<C> algebraicRoots(GenPolynomial<C> f) {
531        List<RealAlgebraicNumber<C>> rl = realAlgebraicNumbers(f);
532        List<ComplexAlgebraicNumber<C>> cl = complexAlgebraicNumbers(f);
533        GenPolynomial<Complex<C>> cf = null;
534        if (!cl.isEmpty()) {
535            cf = cl.get(0).ring.algebraic.modul;
536        }
537        cl = filterOutRealRoots(f, cl, rl);
538        AlgebraicRoots<C> ar = new AlgebraicRoots<C>(f, cf, rl, cl);
539        return ar;
540    }
541
542
543    /**
544     * Roots of unity of real and complex algebraic numbers.
545     * @param ar container of real and complex algebraic numbers.
546     * @return container of real and complex algebraic numbers which are roots
547     *         of unity.
548     */
549    public static <C extends GcdRingElem<C> & Rational> AlgebraicRoots<C> rootsOfUnity(AlgebraicRoots<C> ar) {
550        List<RealAlgebraicNumber<C>> rl = new ArrayList<RealAlgebraicNumber<C>>();
551        for (RealAlgebraicNumber<C> r : ar.real) {
552            if (r.isRootOfUnity()) {
553                rl.add(r);
554            }
555        }
556        List<ComplexAlgebraicNumber<C>> cl = new ArrayList<ComplexAlgebraicNumber<C>>();
557        for (ComplexAlgebraicNumber<C> c : ar.complex) {
558            if (c.isRootOfUnity()) {
559                cl.add(c);
560            }
561        }
562        AlgebraicRoots<C> ur = new AlgebraicRoots<C>(ar.p, ar.cp, rl, cl);
563        return ur;
564    }
565
566
567    /**
568     * Root refinement of real and complex algebraic numbers.
569     * @param a container of real and complex algebraic numbers.
570     * @param eps desired precision for root intervals and rectangles.
571     */
572    public static <C extends GcdRingElem<C> & Rational> void rootRefine(AlgebraicRoots<C> a,
573                    BigRational eps) {
574        for (RealAlgebraicNumber<C> r : a.real) {
575            r.ring.refineRoot(eps);
576        }
577        for (ComplexAlgebraicNumber<C> c : a.complex) {
578            c.ring.refineRoot(eps);
579        }
580        return; // a or void?
581    }
582
583
584    /**
585     * Roots as real and complex decimal numbers.
586     * @param f univariate polynomial.
587     * @param eps desired precision.
588     * @return container of real and complex decimal numbers.
589     */
590    public static <C extends GcdRingElem<C> & Rational> DecimalRoots<C> decimalRoots(GenPolynomial<C> f,
591                    BigRational eps) {
592        RealRootsAbstract<C> rengine = new RealRootsSturm<C>();
593        List<BigDecimal> rl = rengine.approximateRoots(f, eps);
594
595        GenPolynomial<Complex<C>> fc = PolyUtilRoot.<C> complexFromAny(f);
596        ComplexRootsAbstract<C> cengine = new ComplexRootsSturm<C>(fc.ring.coFac);
597        List<Complex<BigDecimal>> cl = cengine.approximateRoots(fc, eps);
598
599        cl = filterOutRealRoots(f, cl, rl, eps);
600        DecimalRoots<C> ar = new DecimalRoots<C>(f, fc, rl, cl);
601        return ar;
602    }
603
604
605    /**
606     * Roots as real and complex decimal numbers.
607     * @param ar container for real and complex algebraic roots.
608     * @param eps desired precision.
609     * @return container of real and complex decimal numbers.
610     */
611    public static <C extends GcdRingElem<C> & Rational> DecimalRoots<C> decimalRoots(AlgebraicRoots<C> ar,
612                    BigRational eps) {
613        //no: rootRefine(ar, eps);
614        RealRootsAbstract<C> rengine = new RealRootsSturm<C>();
615        List<BigDecimal> rl = new ArrayList<BigDecimal>(ar.real.size());
616        for (RealAlgebraicNumber<C> r : ar.real) {
617            try {
618                BigDecimal d = rengine.approximateRoot(r.ring.root, ar.p, eps);
619                rl.add(d);
620            } catch (NoConvergenceException e) {
621                System.out.println("should not happen: " + e);
622            }
623        }
624        ComplexRootsAbstract<C> cengine = new ComplexRootsSturm<C>(ar.cp.ring.coFac);
625        List<Complex<BigDecimal>> cl = new ArrayList<Complex<BigDecimal>>(ar.complex.size());
626        for (ComplexAlgebraicNumber<C> c : ar.complex) {
627            try {
628                Complex<BigDecimal> d = cengine.approximateRoot(c.ring.root, ar.cp, eps);
629                cl.add(d);
630            } catch (NoConvergenceException e) {
631                System.out.println("should not happen: " + e);
632            }
633        }
634        DecimalRoots<C> dr = new DecimalRoots<C>(ar.p, ar.cp, rl, cl);
635        return dr;
636    }
637
638}