001/*
002 * $Id: PolyUtilApp.java 6007 2020-03-29 13:34:49Z kredel $
003 */
004
005package edu.jas.application;
006
007
008import java.util.ArrayList;
009import java.util.Arrays;
010import java.util.List;
011import java.util.Map;
012import java.util.SortedMap;
013import java.util.TreeMap;
014
015import org.apache.logging.log4j.Logger;
016import org.apache.logging.log4j.LogManager; 
017
018import edu.jas.arith.BigDecimal;
019import edu.jas.arith.BigRational;
020import edu.jas.arith.Product;
021import edu.jas.arith.ProductRing;
022import edu.jas.arith.Rational;
023import edu.jas.poly.AlgebraicNumber;
024import edu.jas.poly.AlgebraicNumberRing;
025import edu.jas.poly.Complex;
026import edu.jas.poly.ComplexRing;
027import edu.jas.poly.ExpVector;
028import edu.jas.poly.GenPolynomial;
029import edu.jas.poly.GenPolynomialRing;
030import edu.jas.poly.PolyUtil;
031import edu.jas.poly.PolynomialList;
032import edu.jas.poly.TermOrder;
033import edu.jas.root.ComplexRoots;
034import edu.jas.root.ComplexRootsAbstract;
035import edu.jas.root.ComplexRootsSturm;
036import edu.jas.root.Interval;
037import edu.jas.root.InvalidBoundaryException;
038import edu.jas.root.RealAlgebraicNumber;
039import edu.jas.root.RealAlgebraicRing;
040import edu.jas.root.RealRootTuple;
041import edu.jas.root.RealRootsAbstract;
042import edu.jas.root.RealRootsSturm;
043import edu.jas.root.Rectangle;
044import edu.jas.root.RootFactory;
045import edu.jas.structure.GcdRingElem;
046import edu.jas.structure.RingElem;
047import edu.jas.structure.RingFactory;
048import edu.jas.structure.UnaryFunctor;
049import edu.jas.util.ListUtil;
050
051
052/**
053 * Polynomial utilities for applications, for example conversion ExpVector to
054 * Product or zero dimensional ideal root computation.
055 * @param <C> coefficient type
056 * @author Heinz Kredel
057 */
058public class PolyUtilApp<C extends RingElem<C>> {
059
060
061    private static final Logger logger = LogManager.getLogger(PolyUtilApp.class);
062
063
064    private static final boolean debug = logger.isDebugEnabled();
065
066
067    /**
068     * Product representation.
069     * @param <C> coefficient type.
070     * @param pfac polynomial ring factory.
071     * @param L list of polynomials to be represented.
072     * @return Product represenation of L in the polynomial ring pfac.
073     */
074    public static <C extends GcdRingElem<C>> List<GenPolynomial<Product<Residue<C>>>> toProductRes(
075                    GenPolynomialRing<Product<Residue<C>>> pfac, List<GenPolynomial<GenPolynomial<C>>> L) {
076
077        List<GenPolynomial<Product<Residue<C>>>> list = new ArrayList<GenPolynomial<Product<Residue<C>>>>();
078        if (L == null || L.size() == 0) {
079            return list;
080        }
081        GenPolynomial<Product<Residue<C>>> b;
082        for (GenPolynomial<GenPolynomial<C>> a : L) {
083            b = toProductRes(pfac, a);
084            list.add(b);
085        }
086        return list;
087    }
088
089
090    /**
091     * Product representation.
092     * @param <C> coefficient type.
093     * @param pfac polynomial ring factory.
094     * @param A polynomial to be represented.
095     * @return Product represenation of A in the polynomial ring pfac.
096     */
097    public static <C extends GcdRingElem<C>> GenPolynomial<Product<Residue<C>>> toProductRes(
098                    GenPolynomialRing<Product<Residue<C>>> pfac, GenPolynomial<GenPolynomial<C>> A) {
099
100        GenPolynomial<Product<Residue<C>>> P = pfac.getZERO().copy();
101        if (A == null || A.isZERO()) {
102            return P;
103        }
104        RingFactory<Product<Residue<C>>> rpfac = pfac.coFac;
105        ProductRing<Residue<C>> fac = (ProductRing<Residue<C>>) rpfac;
106        Product<Residue<C>> p;
107        for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.getMap().entrySet()) {
108            ExpVector e = y.getKey();
109            GenPolynomial<C> a = y.getValue();
110            p = toProductRes(fac, a);
111            if (!p.isZERO()) {
112                P.doPutToMap(e, p);
113            }
114        }
115        return P;
116    }
117
118
119    /**
120     * Product representation.
121     * @param <C> coefficient type.
122     * @param pfac product ring factory.
123     * @param c coefficient to be represented.
124     * @return Product represenation of c in the ring pfac.
125     */
126    public static <C extends GcdRingElem<C>> Product<Residue<C>> toProductRes(ProductRing<Residue<C>> pfac,
127                    GenPolynomial<C> c) {
128
129        SortedMap<Integer, Residue<C>> elem = new TreeMap<Integer, Residue<C>>();
130        for (int i = 0; i < pfac.length(); i++) {
131            RingFactory<Residue<C>> rfac = pfac.getFactory(i);
132            ResidueRing<C> fac = (ResidueRing<C>) rfac;
133            Residue<C> u = new Residue<C>(fac, c);
134            //fac.fromInteger( c.getVal() );
135            if (!u.isZERO()) {
136                elem.put(i, u);
137            }
138        }
139        return new Product<Residue<C>>(pfac, elem);
140    }
141
142
143    /**
144     * Product residue representation.
145     * @param <C> coefficient type.
146     * @param CS list of ColoredSystems from comprehensive GB system.
147     * @return Product residue represenation of CS.
148     */
149    public static <C extends GcdRingElem<C>> List<GenPolynomial<Product<Residue<C>>>> toProductRes(
150                    List<ColoredSystem<C>> CS) {
151
152        List<GenPolynomial<Product<Residue<C>>>> list = new ArrayList<GenPolynomial<Product<Residue<C>>>>();
153        if (CS == null || CS.isEmpty()) {
154            return list;
155        }
156        GenPolynomialRing<GenPolynomial<C>> pr = null;
157        List<RingFactory<Residue<C>>> rrl = new ArrayList<RingFactory<Residue<C>>>(CS.size());
158        for (ColoredSystem<C> cs : CS) {
159            Ideal<C> id = cs.condition.zero;
160            ResidueRing<C> r = new ResidueRing<C>(id);
161            if (!rrl.contains(r)) {
162                rrl.add(r);
163            }
164            if (pr == null) {
165                if (cs.list.size() > 0) {
166                    pr = cs.list.get(0).green.ring;
167                }
168            }
169        }
170        if (pr == null) {
171            throw new IllegalArgumentException("no polynomial ring found");
172        }
173        ProductRing<Residue<C>> pfac;
174        pfac = new ProductRing<Residue<C>>(rrl);
175        //System.out.println("pfac = " + pfac);
176        GenPolynomialRing<Product<Residue<C>>> rf = new GenPolynomialRing<Product<Residue<C>>>(pfac, pr.nvar,
177                        pr.tord, pr.getVars());
178        GroebnerSystem<C> gs = new GroebnerSystem<C>(CS);
179        List<GenPolynomial<GenPolynomial<C>>> F = gs.getCGB();
180        list = PolyUtilApp.<C> toProductRes(rf, F);
181        return list;
182    }
183
184
185    /**
186     * Residue coefficient representation.
187     * @param pfac polynomial ring factory.
188     * @param L list of polynomials to be represented.
189     * @return Represenation of L in the polynomial ring pfac.
190     */
191    public static <C extends GcdRingElem<C>> List<GenPolynomial<Residue<C>>> toResidue(
192                    GenPolynomialRing<Residue<C>> pfac, List<GenPolynomial<GenPolynomial<C>>> L) {
193        List<GenPolynomial<Residue<C>>> list = new ArrayList<GenPolynomial<Residue<C>>>();
194        if (L == null || L.size() == 0) {
195            return list;
196        }
197        GenPolynomial<Residue<C>> b;
198        for (GenPolynomial<GenPolynomial<C>> a : L) {
199            b = toResidue(pfac, a);
200            if (!b.isZERO()) {
201                list.add(b);
202            }
203        }
204        return list;
205    }
206
207
208    /**
209     * Residue coefficient representation.
210     * @param pfac polynomial ring factory.
211     * @param A polynomial to be represented.
212     * @return Represenation of A in the polynomial ring pfac.
213     */
214    public static <C extends GcdRingElem<C>> GenPolynomial<Residue<C>> toResidue(
215                    GenPolynomialRing<Residue<C>> pfac, GenPolynomial<GenPolynomial<C>> A) {
216        GenPolynomial<Residue<C>> P = pfac.getZERO().copy();
217        if (A == null || A.isZERO()) {
218            return P;
219        }
220        RingFactory<Residue<C>> rpfac = pfac.coFac;
221        ResidueRing<C> fac = (ResidueRing<C>) rpfac;
222        Residue<C> p;
223        for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.getMap().entrySet()) {
224            ExpVector e = y.getKey();
225            GenPolynomial<C> a = y.getValue();
226            p = new Residue<C>(fac, a);
227            if (!p.isZERO()) {
228                P.doPutToMap(e, p);
229            }
230        }
231        return P;
232    }
233
234
235    /**
236     * Product slice.
237     * @param <C> coefficient type.
238     * @param L list of polynomials with product coefficients.
239     * @return Slices represenation of L.
240     */
241    public static <C extends GcdRingElem<C>> Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> productSlice(
242                    PolynomialList<Product<Residue<C>>> L) {
243
244        Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> map;
245        RingFactory<Product<Residue<C>>> fpr = L.ring.coFac;
246        ProductRing<Residue<C>> pr = (ProductRing<Residue<C>>) fpr;
247        int s = pr.length();
248        map = new TreeMap<Ideal<C>, PolynomialList<GenPolynomial<C>>>();
249        List<GenPolynomial<GenPolynomial<C>>> slist;
250
251        List<GenPolynomial<Product<Residue<C>>>> plist = L.list;
252        PolynomialList<GenPolynomial<C>> spl;
253
254        for (int i = 0; i < s; i++) {
255            RingFactory<Residue<C>> r = pr.getFactory(i);
256            ResidueRing<C> rr = (ResidueRing<C>) r;
257            Ideal<C> id = rr.ideal;
258            GenPolynomialRing<C> cof = rr.ring;
259            GenPolynomialRing<GenPolynomial<C>> pfc;
260            pfc = new GenPolynomialRing<GenPolynomial<C>>(cof, L.ring);
261            slist = fromProduct(pfc, plist, i);
262            spl = new PolynomialList<GenPolynomial<C>>(pfc, slist);
263            PolynomialList<GenPolynomial<C>> d = map.get(id);
264            if (d != null) {
265                throw new RuntimeException("ideal exists twice " + id);
266            }
267            map.put(id, spl);
268        }
269        return map;
270    }
271
272
273    /**
274     * Product slice at i.
275     * @param <C> coefficient type.
276     * @param L list of polynomials with product coeffients.
277     * @param i index of slice.
278     * @return Slice of of L at i.
279     */
280    public static <C extends GcdRingElem<C>> PolynomialList<GenPolynomial<C>> productSlice(
281                    PolynomialList<Product<Residue<C>>> L, int i) {
282
283        RingFactory<Product<Residue<C>>> fpr = L.ring.coFac;
284        ProductRing<Residue<C>> pr = (ProductRing<Residue<C>>) fpr;
285        List<GenPolynomial<GenPolynomial<C>>> slist;
286
287        List<GenPolynomial<Product<Residue<C>>>> plist = L.list;
288        PolynomialList<GenPolynomial<C>> spl;
289
290        RingFactory<Residue<C>> r = pr.getFactory(i);
291        ResidueRing<C> rr = (ResidueRing<C>) r;
292        GenPolynomialRing<C> cof = rr.ring;
293        GenPolynomialRing<GenPolynomial<C>> pfc;
294        pfc = new GenPolynomialRing<GenPolynomial<C>>(cof, L.ring);
295        slist = fromProduct(pfc, plist, i);
296        spl = new PolynomialList<GenPolynomial<C>>(pfc, slist);
297        return spl;
298    }
299
300
301    /**
302     * From product representation.
303     * @param <C> coefficient type.
304     * @param pfac polynomial ring factory.
305     * @param L list of polynomials to be converted from product representation.
306     * @param i index of product representation to be taken.
307     * @return Represenation of i-slice of L in the polynomial ring pfac.
308     */
309    public static <C extends GcdRingElem<C>> List<GenPolynomial<GenPolynomial<C>>> fromProduct(
310                    GenPolynomialRing<GenPolynomial<C>> pfac, List<GenPolynomial<Product<Residue<C>>>> L,
311                    int i) {
312
313        List<GenPolynomial<GenPolynomial<C>>> list = new ArrayList<GenPolynomial<GenPolynomial<C>>>();
314
315        if (L == null || L.size() == 0) {
316            return list;
317        }
318        GenPolynomial<GenPolynomial<C>> b;
319        for (GenPolynomial<Product<Residue<C>>> a : L) {
320            b = fromProduct(pfac, a, i);
321            if (!b.isZERO()) {
322                b = b.abs();
323                if (!list.contains(b)) {
324                    list.add(b);
325                }
326            }
327        }
328        return list;
329    }
330
331
332    /**
333     * From product representation.
334     * @param <C> coefficient type.
335     * @param pfac polynomial ring factory.
336     * @param P polynomial to be converted from product representation.
337     * @param i index of product representation to be taken.
338     * @return Represenation of i-slice of P in the polynomial ring pfac.
339     */
340    public static <C extends GcdRingElem<C>> GenPolynomial<GenPolynomial<C>> fromProduct(
341                    GenPolynomialRing<GenPolynomial<C>> pfac, GenPolynomial<Product<Residue<C>>> P, int i) {
342
343        GenPolynomial<GenPolynomial<C>> b = pfac.getZERO().copy();
344        if (P == null || P.isZERO()) {
345            return b;
346        }
347
348        for (Map.Entry<ExpVector, Product<Residue<C>>> y : P.getMap().entrySet()) {
349            ExpVector e = y.getKey();
350            Product<Residue<C>> a = y.getValue();
351            Residue<C> r = a.get(i);
352            if (r != null && !r.isZERO()) {
353                GenPolynomial<C> p = r.val;
354                if (!p.isZERO()) {
355                    b.doPutToMap(e, p);
356                }
357            }
358        }
359        return b;
360    }
361
362
363    /**
364     * Product slice to String.
365     * @param <C> coefficient type.
366     * @param L list of polynomials with to be represented.
367     * @return Product represenation of L in the polynomial ring pfac.
368     */
369    public static <C extends GcdRingElem<C>> String productSliceToString(
370                    Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> L) {
371        //Set<GenPolynomial<GenPolynomial<C>>> sl = new TreeSet<GenPolynomial<GenPolynomial<C>>>();
372        PolynomialList<GenPolynomial<C>> pl = null;
373        StringBuffer sb = new StringBuffer(); //"\nproductSlice ----------------- begin");
374        for (Map.Entry<Ideal<C>, PolynomialList<GenPolynomial<C>>> en : L.entrySet()) {
375            sb.append("\n\ncondition == 0:\n");
376            sb.append(en.getKey().list.toScript());
377            pl = en.getValue(); //L.get(id);
378            //sl.addAll(pl.list);
379            sb.append("\ncorresponding ideal:\n");
380            sb.append(pl.toScript());
381        }
382        //List<GenPolynomial<GenPolynomial<C>>> sll 
383        //   = new ArrayList<GenPolynomial<GenPolynomial<C>>>( sl );
384        //pl = new PolynomialList<GenPolynomial<C>>(pl.ring,sll);
385        // sb.append("\nunion = " + pl.toString());
386        //sb.append("\nproductSlice ------------------------- end\n");
387        return sb.toString();
388    }
389
390
391    /**
392     * Product slice to String.
393     * @param <C> coefficient type.
394     * @param L list of polynomials with product coefficients.
395     * @return string represenation of slices of L.
396     */
397    public static <C extends GcdRingElem<C>> String productToString(PolynomialList<Product<Residue<C>>> L) {
398        Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> M;
399        M = productSlice(L);
400        String s = productSliceToString(M);
401        return s;
402    }
403
404
405    /**
406     * Construct superset of complex roots for zero dimensional ideal(G).
407     * @param I zero dimensional ideal.
408     * @param eps desired precision.
409     * @return list of coordinates of complex roots for ideal(G)
410     */
411    public static <D extends GcdRingElem<D> & Rational> List<List<Complex<BigDecimal>>> complexRootTuples(
412                    Ideal<D> I, BigRational eps) {
413        List<GenPolynomial<D>> univs = I.constructUnivariate();
414        if (logger.isInfoEnabled()) {
415            logger.info("univs = " + univs);
416        }
417        return complexRoots(I, univs, eps);
418    }
419
420
421    /**
422     * Construct superset of complex roots for zero dimensional ideal(G).
423     * @param I zero dimensional ideal.
424     * @param univs list of univariate polynomials.
425     * @param eps desired precision.
426     * @return list of coordinates of complex roots for ideal(G)
427     */
428    public static <D extends GcdRingElem<D> & Rational> List<List<Complex<BigDecimal>>> complexRoots(
429                    Ideal<D> I, List<GenPolynomial<D>> univs, BigRational eps) {
430        List<List<Complex<BigDecimal>>> croots = new ArrayList<List<Complex<BigDecimal>>>();
431        RingFactory<D> cf = I.list.ring.coFac;
432        ComplexRing<D> cr = new ComplexRing<D>(cf);
433        ComplexRootsAbstract<D> cra = new ComplexRootsSturm<D>(cr);
434        List<GenPolynomial<Complex<D>>> cunivs = new ArrayList<GenPolynomial<Complex<D>>>();
435        for (GenPolynomial<D> p : univs) {
436            GenPolynomialRing<Complex<D>> pfac = new GenPolynomialRing<Complex<D>>(cr, p.ring);
437            //System.out.println("pfac = " + pfac.toScript());
438            GenPolynomial<Complex<D>> cp = PolyUtil.<D> toComplex(pfac, p);
439            cunivs.add(cp);
440            //System.out.println("cp = " + cp);
441        }
442        for (int i = 0; i < I.list.ring.nvar; i++) {
443            List<Complex<BigDecimal>> cri = cra.approximateRoots(cunivs.get(i), eps);
444            //System.out.println("cri = " + cri);
445            croots.add(cri);
446        }
447        croots = ListUtil.<Complex<BigDecimal>> tupleFromList(croots);
448        return croots;
449    }
450
451
452    /**
453     * Construct superset of complex roots for zero dimensional ideal(G).
454     * @param Il list of zero dimensional ideals with univariate polynomials.
455     * @param eps desired precision.
456     * @return list of coordinates of complex roots for ideal(cap_i(G_i))
457     */
458    public static <D extends GcdRingElem<D> & Rational> List<List<Complex<BigDecimal>>> complexRootTuples(
459                    List<IdealWithUniv<D>> Il, BigRational eps) {
460        List<List<Complex<BigDecimal>>> croots = new ArrayList<List<Complex<BigDecimal>>>();
461        for (IdealWithUniv<D> I : Il) {
462            List<List<Complex<BigDecimal>>> cr = complexRoots(I.ideal, I.upolys, eps);
463            croots.addAll(cr);
464        }
465        return croots;
466    }
467
468
469    /**
470     * Construct superset of complex roots for zero dimensional ideal(G).
471     * @param Il list of zero dimensional ideals with univariate polynomials.
472     * @param eps desired precision.
473     * @return list of ideals with coordinates of complex roots for
474     *         ideal(cap_i(G_i))
475     */
476    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexRoots<D>> complexRoots(
477                    List<IdealWithUniv<D>> Il, BigRational eps) {
478        List<IdealWithComplexRoots<D>> Ic = new ArrayList<IdealWithComplexRoots<D>>(Il.size());
479        for (IdealWithUniv<D> I : Il) {
480            List<List<Complex<BigDecimal>>> cr = complexRoots(I.ideal, I.upolys, eps);
481            IdealWithComplexRoots<D> ic = new IdealWithComplexRoots<D>(I, cr);
482            Ic.add(ic);
483        }
484        return Ic;
485    }
486
487
488    /**
489     * Construct superset of complex roots for zero dimensional ideal(G).
490     * @param G list of polynomials of a of zero dimensional ideal.
491     * @param eps desired precision.
492     * @return list of ideals with coordinates of complex roots for ideal(G)
493     */
494    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexRoots<D>> complexRoots(
495                    Ideal<D> G, BigRational eps) {
496        List<IdealWithUniv<D>> Il = G.zeroDimDecomposition();
497        return complexRoots(Il, eps);
498    }
499
500
501    /**
502     * Construct superset of real roots for zero dimensional ideal(G).
503     * @param I zero dimensional ideal.
504     * @param eps desired precision.
505     * @return list of coordinates of real roots for ideal(G)
506     */
507    public static <D extends GcdRingElem<D> & Rational> List<List<BigDecimal>> realRootTuples(Ideal<D> I,
508                    BigRational eps) {
509        List<GenPolynomial<D>> univs = I.constructUnivariate();
510        if (logger.isInfoEnabled()) {
511            logger.info("univs = " + univs);
512        }
513        return realRoots(I, univs, eps);
514    }
515
516
517    /**
518     * Construct superset of real roots for zero dimensional ideal(G).
519     * @param I zero dimensional ideal.
520     * @param univs list of univariate polynomials.
521     * @param eps desired precision.
522     * @return list of coordinates of real roots for ideal(G)
523     */
524    public static <D extends GcdRingElem<D> & Rational> List<List<BigDecimal>> realRoots(Ideal<D> I,
525                    List<GenPolynomial<D>> univs, BigRational eps) {
526        List<List<BigDecimal>> roots = new ArrayList<List<BigDecimal>>();
527        //RingFactory<D> cf = (RingFactory<D>) I.list.ring.coFac;
528        RealRootsAbstract<D> rra = new RealRootsSturm<D>();
529        for (int i = 0; i < I.list.ring.nvar; i++) {
530            List<BigDecimal> rri = rra.approximateRoots(univs.get(i), eps);
531            //System.out.println("rri = " + rri);
532            roots.add(rri);
533        }
534        //System.out.println("roots-1 = " + roots);
535        roots = ListUtil.<BigDecimal> tupleFromList(roots);
536        //System.out.println("roots-2 = " + roots);
537        return roots;
538    }
539
540
541    /**
542     * Construct superset of real roots for zero dimensional ideal(G).
543     * @param Il list of zero dimensional ideals with univariate polynomials.
544     * @param eps desired precision.
545     * @return list of coordinates of real roots for ideal(cap_i(G_i))
546     */
547    public static <D extends GcdRingElem<D> & Rational> List<List<BigDecimal>> realRootTuples(
548                    List<IdealWithUniv<D>> Il, BigRational eps) {
549        List<List<BigDecimal>> rroots = new ArrayList<List<BigDecimal>>();
550        for (IdealWithUniv<D> I : Il) {
551            List<List<BigDecimal>> rr = realRoots(I.ideal, I.upolys, eps);
552            rroots.addAll(rr);
553        }
554        return rroots;
555    }
556
557
558    /**
559     * Construct superset of real roots for zero dimensional ideal(G).
560     * @param Il list of zero dimensional ideals with univariate polynomials.
561     * @param eps desired precision.
562     * @return list of ideals with coordinates of real roots for
563     *         ideal(cap_i(G_i))
564     */
565    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealRoots<D>> realRoots(
566                    List<IdealWithUniv<D>> Il, BigRational eps) {
567        List<IdealWithRealRoots<D>> Ir = new ArrayList<IdealWithRealRoots<D>>(Il.size());
568        for (IdealWithUniv<D> I : Il) {
569            List<List<BigDecimal>> rr = realRoots(I.ideal, I.upolys, eps);
570            IdealWithRealRoots<D> ir = new IdealWithRealRoots<D>(I, rr);
571            Ir.add(ir);
572        }
573        return Ir;
574    }
575
576
577    /**
578     * Construct superset of real roots for zero dimensional ideal(G).
579     * @param G list of polynomials of a of zero dimensional ideal.
580     * @param eps desired precision.
581     * @return list of ideals with coordinates of real roots for ideal(G)
582     */
583    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealRoots<D>> realRoots(Ideal<D> G,
584                    BigRational eps) {
585        List<IdealWithUniv<D>> Il = G.zeroDimDecomposition();
586        return realRoots(Il, eps);
587    }
588
589
590    /**
591     * Test for real roots of zero dimensional ideal(L).
592     * @param L list of polynomials.
593     * @param roots list of real roots for ideal(G).
594     * @param eps desired precision.
595     * @return true if root is a list of coordinates of real roots for ideal(L)
596     */
597    public static boolean isRealRoots(List<GenPolynomial<BigDecimal>> L, List<List<BigDecimal>> roots,
598                    BigDecimal eps) {
599        if (L == null || L.size() == 0) {
600            return true;
601        }
602        // polynomials with decimal coefficients
603        BigDecimal dc = BigDecimal.ONE;
604        //GenPolynomialRing<BigDecimal> dfac = L.get(0).ring;
605        //System.out.println("dfac = " + dfac);
606        for (GenPolynomial<BigDecimal> dp : L) {
607            //System.out.println("dp = " + dp);
608            for (List<BigDecimal> r : roots) {
609                //System.out.println("r = " + r);
610                BigDecimal ev = PolyUtil.<BigDecimal> evaluateAll(dc, dp, r);
611                if (ev.abs().compareTo(eps) > 0) {
612                    System.out.println("ev = " + ev);
613                    return false;
614                }
615            }
616        }
617        return true;
618    }
619
620
621    /**
622     * Test for complex roots of zero dimensional ideal(L).
623     * @param L list of polynomials.
624     * @param roots list of real roots for ideal(G).
625     * @param eps desired precision.
626     * @return true if root is a list of coordinates of complex roots for
627     *         ideal(L)
628     */
629    public static boolean isComplexRoots(List<GenPolynomial<Complex<BigDecimal>>> L,
630                    List<List<Complex<BigDecimal>>> roots, BigDecimal eps) {
631        if (L == null || L.size() == 0) {
632            return true;
633        }
634        // polynomials with decimal coefficients
635        BigDecimal dc = BigDecimal.ONE;
636        ComplexRing<BigDecimal> dcc = new ComplexRing<BigDecimal>(dc);
637        //GenPolynomialRing<Complex<BigDecimal>> dfac = L.get(0).ring;
638        //System.out.println("dfac = " + dfac);
639        for (GenPolynomial<Complex<BigDecimal>> dp : L) {
640            //System.out.println("dp = " + dp);
641            for (List<Complex<BigDecimal>> r : roots) {
642                //System.out.println("r = " + r);
643                Complex<BigDecimal> ev = PolyUtil.<Complex<BigDecimal>> evaluateAll(dcc, dp, r);
644                if (ev.norm().getRe().compareTo(eps) > 0) {
645                    System.out.println("ev = " + ev);
646                    return false;
647                }
648            }
649        }
650        return true;
651    }
652
653
654    /**
655     * Construct real roots for zero dimensional ideal(G).
656     * @param I zero dimensional ideal with univariate irreducible polynomials
657     *            and bi-variate polynomials.
658     * @return real algebraic roots for ideal(G)
659     */
660    public static <D extends GcdRingElem<D> & Rational> IdealWithRealAlgebraicRoots<D> realAlgebraicRoots(
661                    IdealWithUniv<D> I) {
662        List<List<RealAlgebraicNumber<D>>> ran = new ArrayList<List<RealAlgebraicNumber<D>>>();
663        if (I == null) {
664            throw new IllegalArgumentException("null ideal not permitted");
665        }
666        if (I.ideal == null || I.upolys == null) {
667            throw new IllegalArgumentException("null ideal components not permitted " + I);
668        }
669        if (I.ideal.isZERO() || I.upolys.size() == 0) {
670            return new IdealWithRealAlgebraicRoots<D>(I, ran);
671        }
672        GenPolynomialRing<D> fac = I.ideal.list.ring;
673        // case i == 0:
674        GenPolynomial<D> p0 = I.upolys.get(0);
675        GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
676        if (p0p == null) {
677            throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of  " + I.ideal);
678        }
679        //System.out.println("p0  = " + p0);
680        if (logger.isInfoEnabled()) {
681            logger.info("p0p = " + p0p);
682        }
683        int[] dep0 = p0p.degreeVector().dependencyOnVariables();
684        //System.out.println("dep0 = " + Arrays.toString(dep0));
685        if (dep0.length != 1) {
686            throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
687        }
688        List<RealAlgebraicNumber<D>> rra = RootFactory.<D> realAlgebraicNumbersIrred(p0);
689        if (logger.isInfoEnabled()) {
690            List<Interval<D>> il = new ArrayList<Interval<D>>();
691            for (RealAlgebraicNumber<D> rr : rra) {
692                il.add(rr.ring.getRoot());
693            }
694            logger.info("roots(p0) = " + il);
695        }
696        for (RealAlgebraicNumber<D> rr : rra) {
697            List<RealAlgebraicNumber<D>> rl = new ArrayList<RealAlgebraicNumber<D>>();
698            rl.add(rr);
699            ran.add(rl);
700        }
701        // case i > 0:
702        for (int i = 1; i < I.upolys.size(); i++) {
703            List<List<RealAlgebraicNumber<D>>> rn = new ArrayList<List<RealAlgebraicNumber<D>>>();
704            GenPolynomial<D> pi = I.upolys.get(i);
705            GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
706            if (pip == null) {
707                throw new RuntimeException(
708                                "no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
709            }
710            //System.out.println("i   = " + i);
711            //System.out.println("pi  = " + pi);
712            if (logger.isInfoEnabled()) {
713                logger.info("pi  = " + pi);
714                logger.info("pip = " + pip);
715            }
716            int[] depi = pip.degreeVector().dependencyOnVariables();
717            //System.out.println("depi = " + Arrays.toString(depi));
718            if (depi.length < 1 || depi.length > 2) {
719                throw new RuntimeException("wrong number of variables " + Arrays.toString(depi));
720            }
721            rra = RootFactory.<D> realAlgebraicNumbersIrred(pi);
722            if (logger.isInfoEnabled()) {
723                List<Interval<D>> il = new ArrayList<Interval<D>>();
724                for (RealAlgebraicNumber<D> rr : rra) {
725                    il.add(rr.ring.getRoot());
726                }
727                logger.info("roots(pi) = " + il);
728            }
729            if (depi.length == 1) {
730                // all combinations are roots of the ideal I
731                for (RealAlgebraicNumber<D> rr : rra) {
732                    //System.out.println("rr.ring = " + rr.ring);
733                    for (List<RealAlgebraicNumber<D>> rx : ran) {
734                        //System.out.println("rx = " + rx);
735                        List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
736                        ry.addAll(rx);
737                        ry.add(rr);
738                        rn.add(ry);
739                    }
740                }
741            } else { // depi.length == 2
742                // select roots of the ideal I
743                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
744                //System.out.println("pip2 = " + pip2.ring);
745                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
746                TermOrder to = new TermOrder(TermOrder.INVLEX);
747                GenPolynomialRing<GenPolynomial<D>> rfac = new GenPolynomialRing<GenPolynomial<D>>(ufac, 1,
748                                                                                                   to); // new vars
749                GenPolynomial<GenPolynomial<D>> pip2r = PolyUtil.<D> recursive(rfac, pip2);
750                int ix = fac.nvar - 1 - depi[depi.length - 1];
751                //System.out.println("ix = " + ix);
752                for (RealAlgebraicNumber<D> rr : rra) {
753                    //System.out.println("rr.ring = " + rr.ring);
754                    Interval<D> rroot = rr.ring.getRoot();
755                    GenPolynomial<D> pip2el = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.left);
756                    GenPolynomial<D> pip2er = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.right);
757                    GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
758                    GenPolynomial<D> pip2elc = convert(upfac, pip2el);
759                    GenPolynomial<D> pip2erc = convert(upfac, pip2er);
760                    //System.out.println("pip2elc = " + pip2elc);
761                    //System.out.println("pip2erc = " + pip2erc);
762                    for (List<RealAlgebraicNumber<D>> rx : ran) {
763                        //System.out.println("rx = " + rx);
764                        RealAlgebraicRing<D> rar = rx.get(ix).ring;
765                        //System.out.println("rar = " + rar.toScript());
766                        RealAlgebraicNumber<D> rel = new RealAlgebraicNumber<D>(rar, pip2elc);
767                        RealAlgebraicNumber<D> rer = new RealAlgebraicNumber<D>(rar, pip2erc);
768                        int sl = rel.signum();
769                        int sr = rer.signum();
770                        //System.out.println("sl = " + sl + ", sr = " + sr + ", sl*sr = " + (sl*sr));
771                        if (sl * sr <= 0) {
772                            //System.out.println("sl * sr <= 0: rar = " + rar.toScript());
773                            List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
774                            ry.addAll(rx);
775                            ry.add(rr);
776                            rn.add(ry);
777                        }
778                    }
779                }
780            }
781            ran = rn;
782        }
783        if (logger.isInfoEnabled()) {
784            for (List<RealAlgebraicNumber<D>> rz : ran) {
785                List<Interval<D>> il = new ArrayList<Interval<D>>();
786                for (RealAlgebraicNumber<D> rr : rz) {
787                    il.add(rr.ring.getRoot());
788                }
789                logger.info("root-tuple = " + il);
790            }
791        }
792        IdealWithRealAlgebraicRoots<D> Ir = new IdealWithRealAlgebraicRoots<D>(I, ran);
793        return Ir;
794    }
795
796
797    /**
798     * Construct real roots for zero dimensional ideal(G).
799     * @param I list of zero dimensional ideal with univariate irreducible
800     *            polynomials and bi-variate polynomials.
801     * @return list of real algebraic roots for all ideal(I_i)
802     */
803    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots(
804                    List<IdealWithUniv<D>> I) {
805        List<IdealWithRealAlgebraicRoots<D>> lir = new ArrayList<IdealWithRealAlgebraicRoots<D>>(I.size());
806        for (IdealWithUniv<D> iu : I) {
807            IdealWithRealAlgebraicRoots<D> iur = PolyUtilApp.<D> realAlgebraicRoots(iu);
808            //System.out.println("iur = " + iur);
809            lir.add(iur);
810        }
811        return lir;
812    }
813
814
815    /**
816     * Construct complex roots for zero dimensional ideal(G).
817     * @param I zero dimensional ideal with univariate irreducible polynomials
818     *            and bi-variate polynomials.
819     * @return complex algebraic roots for ideal(G) <b>Note:</b> implementation
820     *         contains errors, do not use.
821     */
822    public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRootsWrong( // Wrong
823                    IdealWithUniv<D> I) {
824        List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can;
825        can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
826        if (I == null) {
827            throw new IllegalArgumentException("null ideal not permitted");
828        }
829        if (I.ideal == null || I.upolys == null) {
830            throw new IllegalArgumentException("null ideal components not permitted " + I);
831        }
832        if (I.ideal.isZERO() || I.upolys.size() == 0) {
833            return new IdealWithComplexAlgebraicRoots<D>(I, can);
834        }
835        GenPolynomialRing<D> fac = I.ideal.list.ring;
836        if (fac.nvar == 0) {
837            return new IdealWithComplexAlgebraicRoots<D>(I, can);
838        }
839        if (fac.nvar != I.upolys.size()) {
840            throw new IllegalArgumentException("ideal not zero dimnsional: " + I);
841        }
842        // case i == 0:
843        GenPolynomial<D> p0 = I.upolys.get(0);
844        GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
845        if (p0p == null) {
846            throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of  " + I.ideal);
847        }
848        if (logger.isInfoEnabled()) {
849            logger.info("p0  = " + p0);
850            logger.info("p0p = " + p0p);
851        }
852        int[] dep0 = p0p.degreeVector().dependencyOnVariables();
853        //System.out.println("dep0 = " + Arrays.toString(dep0));
854        if (dep0.length != 1) {
855            throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
856        }
857        RingFactory<D> cfac = p0.ring.coFac;
858        ComplexRing<D> ccfac = new ComplexRing<D>(cfac);
859        GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring);
860        GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0);
861        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra;
862        cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(p0c);
863        logger.info("#roots(p0c) = " + cra.size());
864        if (debug) {
865            boolean t = edu.jas.application.RootFactoryApp.<D> isRoot(p0c, cra);
866            if (!t) {
867                throw new RuntimeException("no roots of " + p0c);
868            }
869        }
870        for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
871            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl;
872            cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
873            cl.add(cr);
874            can.add(cl);
875        }
876        if (fac.nvar == 1) {
877            return new IdealWithComplexAlgebraicRoots<D>(I, can);
878        }
879        // case i > 0:
880        for (int i = 1; i < I.upolys.size(); i++) {
881            List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn;
882            cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
883            GenPolynomial<D> pi = I.upolys.get(i);
884            GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
885            if (pip == null) {
886                throw new RuntimeException(
887                                "no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
888            }
889            if (logger.isInfoEnabled()) {
890                logger.info("pi(" + i + ") = " + pi);
891                logger.info("pip  = " + pip);
892            }
893            facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring);
894            GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi);
895            int[] depi = pip.degreeVector().dependencyOnVariables();
896            //System.out.println("depi = " + Arrays.toString(depi));
897            if (depi.length < 1 || depi.length > 2) {
898                throw new RuntimeException(
899                                "wrong number of variables " + Arrays.toString(depi) + " for " + pip);
900            }
901            cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(pic);
902            logger.info("#roots(pic) = " + cra.size());
903            if (debug) {
904                boolean t = edu.jas.application.RootFactoryApp.<D> isRoot(pic, cra);
905                if (!t) {
906                    throw new RuntimeException("no roots of " + pic);
907                }
908            }
909            if (depi.length == 1) {
910                // all combinations are roots of the ideal I
911                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
912                    //System.out.println("cr.ring = " + cr.ring);
913                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
914                        //System.out.println("cx = " + cx);
915                        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
916                        cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
917                        cy.addAll(cx);
918                        cy.add(cr);
919                        cn.add(cy);
920                    }
921                }
922            } else { // depi.length == 2
923                // select roots of the ideal I
924                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
925                GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1);
926                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
927                GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac);
928                GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring);
929                GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2);
930                //System.out.println("pip2c = " + pip2c);
931                GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac;
932                rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac);
933                GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac,
934                                pip2c);
935                //System.out.println("pip2cr = " + pip2cr);
936
937                int ix = fac.nvar - 1 - depi[depi.length - 1];
938                //System.out.println("ix = " + ix);
939                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
940                    System.out.println("cr = " + toString(cr)); // <----------------------------------
941                    edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring;
942                    RealRootTuple<D> rroot = cring.getRoot();
943                    List<RealAlgebraicNumber<D>> rlist = rroot.tuple;
944                    Interval<D> vr = rlist.get(0).ring.getRoot();
945                    Interval<D> vi = rlist.get(1).ring.getRoot();
946                    logger.info("vr = " + vr + ", vi = " + vi);
947                    if (vr.length().isZERO()) {
948                        D e = vr.left.factory().parse("1/2");
949                        D m = vr.left; //middle();
950                        vr = new Interval<D>(m.subtract(e), m.sum(e));
951                        logger.info("|vr| == 0: " + vr);
952                    }
953                    if (vi.length().isZERO()) {
954                        D e = vi.left.factory().parse("1/2");
955                        D m = vi.left; //middle();
956                        vi = new Interval<D>(m.subtract(e), m.sum(e));
957                        logger.info("|vi| == 0: " + vi);
958                    }
959                    Complex<D> sw = new Complex<D>(ccfac, vr.left, vi.left);
960                    Complex<D> ne = new Complex<D>(ccfac, vr.right, vi.right);
961                    logger.info("sw   = " + toString1(sw) + ", ne   = " + toString1(ne));
962                    GenPolynomial<Complex<D>> pip2cesw, pip2cene;
963                    pip2cesw = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, sw);
964                    pip2cene = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, ne);
965                    GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
966                    GenPolynomialRing<Complex<D>> upcfac = new GenPolynomialRing<Complex<D>>(ccfac, upfac);
967                    //System.out.println("upfac = " + upfac);
968                    //System.out.println("upcfac = " + upcfac);
969                    GenPolynomial<Complex<D>> pip2eswc = convertComplexComplex(upcfac, pip2cesw);
970                    GenPolynomial<Complex<D>> pip2enec = convertComplexComplex(upcfac, pip2cene);
971                    //System.out.println("pip2eswc = " + pip2eswc);
972                    //System.out.println("pip2enec = " + pip2enec);
973                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
974                        //System.out.println("cxi = " + toString(cx.get(ix)));
975                        Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix);
976                        ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring;
977                        edu.jas.application.RealAlgebraicRing<D> rar = (edu.jas.application.RealAlgebraicRing<D>) car.ring;
978                        //System.out.println("car = " + car);
979                        //System.out.println("rar = " + rar);
980                        TermOrder to = new TermOrder(TermOrder.INVLEX);
981                        String vvr = rar.algebraic.ring.getVars()[0];
982                        String vvi = rar.algebraic.ring.getVars()[1];
983                        String[] vars = new String[] { vvr, vvi };
984                        GenPolynomialRing<Complex<D>> tfac = new GenPolynomialRing<Complex<D>>(ccfac, to,
985                                        vars);
986                        GenPolynomial<Complex<D>> t = tfac.univariate(1, 1L)
987                                        .sum(tfac.univariate(0, 1L).multiply(ccfac.getIMAG()));
988                        //System.out.println("t  = " + t); // t = x + i y
989                        GenPolynomialRing<D> rtfac = new GenPolynomialRing<D>(cfac, tfac);
990                        GenPolynomial<Complex<D>> su;
991                        GenPolynomial<D> re, im;
992                        su = PolyUtil.<Complex<D>> substituteUnivariate(pip2eswc, t);
993                        //su = su.monic(); not here
994                        re = PolyUtil.<D> realPartFromComplex(rtfac, su);
995                        im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
996                        //System.out.println("re = " + re);
997                        //System.out.println("im = " + im);
998                        edu.jas.application.RealAlgebraicNumber<D> resw, imsw, rene, imne;
999                        resw = new edu.jas.application.RealAlgebraicNumber<D>(rar, re);
1000                        //System.out.println("resw = " + resw);
1001                        int sswr = resw.signum();
1002                        imsw = new edu.jas.application.RealAlgebraicNumber<D>(rar, im);
1003                        //System.out.println("imsw = " + imsw);
1004                        int sswi = imsw.signum();
1005                        su = PolyUtil.<Complex<D>> substituteUnivariate(pip2enec, t);
1006                        //su = su.monic(); not here
1007                        re = PolyUtil.<D> realPartFromComplex(rtfac, su);
1008                        im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
1009                        //System.out.println("re = " + re);
1010                        //System.out.println("im = " + im);
1011                        rene = new edu.jas.application.RealAlgebraicNumber<D>(rar, re);
1012                        //System.out.println("rene = " + rene);
1013                        int sner = rene.signum();
1014                        imne = new edu.jas.application.RealAlgebraicNumber<D>(rar, im);
1015                        //System.out.println("imne = " + imne);
1016                        int snei = imne.signum();
1017                        //System.out.println("sswr = " + sswr + ", sswi = " + sswi);
1018                        //System.out.println("sner = " + sner + ", snei = " + snei);
1019                        if ((sswr * sner <= 0 && sswi * snei <= 0)) { // wrong !
1020                            logger.info("   hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1021                            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1022                            cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1023                            cy.addAll(cx);
1024                            cy.add(cr);
1025                            cn.add(cy);
1026                        } else {
1027                            logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1028                        }
1029                    }
1030                }
1031            }
1032            can = cn;
1033        }
1034        IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can);
1035        return Ic;
1036    }
1037
1038
1039    /**
1040     * Construct complex roots for zero dimensional ideal(G).
1041     * @param I zero dimensional ideal with univariate irreducible polynomials
1042     *            and bi-variate polynomials.
1043     * @return complex algebraic roots for ideal(G) <b>Note:</b> not jet
1044     *         completed oin all cases.
1045     */
1046    public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRoots(
1047                    IdealWithUniv<D> I) {
1048        List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can;
1049        can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
1050        if (I == null) {
1051            throw new IllegalArgumentException("null ideal not permitted");
1052        }
1053        if (I.ideal == null || I.upolys == null) {
1054            throw new IllegalArgumentException("null ideal components not permitted " + I);
1055        }
1056        if (I.ideal.isZERO() || I.upolys.size() == 0) {
1057            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1058        }
1059        GenPolynomialRing<D> fac = I.ideal.list.ring;
1060        if (fac.nvar == 0) {
1061            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1062        }
1063        if (fac.nvar != I.upolys.size()) {
1064            throw new IllegalArgumentException("ideal not zero dimnsional: " + I);
1065        }
1066        // case i == 0:
1067        GenPolynomial<D> p0 = I.upolys.get(0);
1068        GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
1069        if (p0p == null) {
1070            throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of  " + I.ideal);
1071        }
1072        if (logger.isInfoEnabled()) {
1073            logger.info("p0  = " + p0);
1074            logger.info("p0p = " + p0p);
1075        }
1076        int[] dep0 = p0p.degreeVector().dependencyOnVariables();
1077        //System.out.println("dep0 = " + Arrays.toString(dep0));
1078        if (dep0.length != 1) {
1079            throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
1080        }
1081        RingFactory<D> cfac = p0.ring.coFac;
1082        ComplexRing<D> ccfac = new ComplexRing<D>(cfac);
1083        GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring);
1084        GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0);
1085        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra;
1086        cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(p0c);
1087        logger.info("#roots(p0c) = " + cra.size());
1088        for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1089            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl;
1090            cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1091            cl.add(cr);
1092            can.add(cl);
1093        }
1094        if (fac.nvar == 1) {
1095            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1096        }
1097        // case i > 0:
1098        for (int i = 1; i < I.upolys.size(); i++) {
1099            List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn;
1100            cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
1101            GenPolynomial<D> pi = I.upolys.get(i);
1102            GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
1103            if (pip == null) {
1104                throw new RuntimeException(
1105                                "no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
1106            }
1107            if (logger.isInfoEnabled()) {
1108                logger.info("pi(" + i + ") = " + pi);
1109                logger.info("pip  = " + pip);
1110            }
1111            facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring);
1112            GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi);
1113            int[] depi = pip.degreeVector().dependencyOnVariables();
1114            //System.out.println("depi = " + Arrays.toString(depi));
1115            if (depi.length < 1 || depi.length > 2) {
1116                throw new RuntimeException(
1117                                "wrong number of variables " + Arrays.toString(depi) + " for " + pip);
1118            }
1119            cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(pic);
1120            logger.info("#roots(pic) = " + cra.size());
1121            if (depi.length == 1) {
1122                // all combinations are roots of the ideal I
1123                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1124                    //System.out.println("cr.ring = " + cr.ring);
1125                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
1126                        //System.out.println("cx = " + cx);
1127                        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1128                        cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1129                        cy.addAll(cx);
1130                        cy.add(cr);
1131                        cn.add(cy);
1132                    }
1133                }
1134            } else { // depi.length == 2
1135                // select roots of the ideal I
1136                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
1137                pip2 = PolyUtil.<D> removeUnusedLowerVariables(pip2);
1138                pip2 = PolyUtil.<D> removeUnusedMiddleVariables(pip2);
1139                GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1);
1140                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
1141                GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac);
1142                GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring);
1143                GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2);
1144                GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac;
1145                rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac);
1146                GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac,
1147                                pip2c);
1148                //System.out.println("pip2cr = " + pip2cr);
1149                int ix = fac.nvar - 1 - depi[depi.length - 1];
1150                //System.out.println("ix = " + ix);
1151                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1152                    //System.out.println("cr = " + toString(cr)); 
1153                    edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring;
1154                    RealRootTuple<D> rroot = cring.getRoot();
1155                    List<RealAlgebraicNumber<D>> rlist = rroot.tuple;
1156                    //System.out.println("rlist = " + rlist);
1157                    Interval<D> vr = rlist.get(0).ring.getRoot();
1158                    Interval<D> vi = rlist.get(1).ring.getRoot();
1159                    //logger.info("vr = " + vr + ", vi = " + vi);
1160                    edu.jas.application.RealAlgebraicNumber<D> vrl, vil, vrr, vir;
1161                    vrl = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.left);
1162                    vil = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.left);
1163                    vrr = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.right);
1164                    vir = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.right);
1165                    ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> crr;
1166                    crr = new ComplexRing<edu.jas.application.RealAlgebraicNumber<D>>(cring);
1167                    Complex<edu.jas.application.RealAlgebraicNumber<D>> csw, cne;
1168                    csw = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrl, vil);
1169                    cne = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrr, vir);
1170                    //logger.info("csw  = " + toString(csw)   + ", cne  = " + toString(cne));
1171                    Rectangle<edu.jas.application.RealAlgebraicNumber<D>> rec;
1172                    rec = new Rectangle<edu.jas.application.RealAlgebraicNumber<D>>(csw, cne);
1173                    //System.out.println("rec = " + rec);
1174                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
1175                        Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix);
1176                        //System.out.println("cax = " + toString(cax));
1177                        ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring;
1178                        //System.out.println("car = " + car);
1179                        GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcrfac;
1180                        pcrfac = new GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(
1181                                        car, rcfac);
1182                        GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcr;
1183                        pcr = evaluateToComplexRealCoefficients(pcrfac, pip2cr, cax);
1184                        //System.out.println("pcr = " + pcr);
1185                        ComplexRoots<edu.jas.application.RealAlgebraicNumber<D>> rengine;
1186                        rengine = new ComplexRootsSturm<edu.jas.application.RealAlgebraicNumber<D>>(car);
1187                        long nr = 0;
1188                        try {
1189                            nr = rengine.complexRootCount(rec, pcr);
1190                            //logger.info("rootCount = " + nr);
1191                        } catch (InvalidBoundaryException e) {
1192                            e.printStackTrace();
1193                        }
1194                        if (nr == 1) { // one root
1195                            logger.info("   hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1196                            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1197                            cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1198                            cy.addAll(cx);
1199                            cy.add(cr);
1200                            cn.add(cy);
1201                        } else if (nr > 1) {
1202                            logger.error("to many roots, cxi = " + toString(cx.get(ix)) + ", cr = "
1203                                            + toString(cr));
1204                        } else { // no root
1205                            logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1206                        }
1207                    }
1208                }
1209            }
1210            can = cn;
1211        }
1212        IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can);
1213        return Ic;
1214    }
1215
1216
1217    /**
1218     * String representation of a deximal approximation of a complex number.
1219     * @param c compelx number.
1220     * @return String representation of c
1221     */
1222    public static <D extends GcdRingElem<D> & Rational> String toString(
1223                    Complex<edu.jas.application.RealAlgebraicNumber<D>> c) {
1224        edu.jas.application.RealAlgebraicNumber<D> re = c.getRe();
1225        edu.jas.application.RealAlgebraicNumber<D> im = c.getIm();
1226        String s = re.decimalMagnitude().toString();
1227        if (!im.isZERO()) {
1228            s = s + "i" + im.decimalMagnitude();
1229        }
1230        return s;
1231    }
1232
1233
1234    /**
1235     * String representation of a deximal approximation of a complex number.
1236     * @param c compelx number.
1237     * @return String representation of c
1238     */
1239    public static <D extends GcdRingElem<D> & Rational> String toString1(Complex<D> c) {
1240        D re = c.getRe();
1241        D im = c.getIm();
1242        String s = new BigDecimal(re.getRational()).toString();
1243        if (!im.isZERO()) {
1244            s = s + "i" + new BigDecimal(im.getRational());
1245        }
1246        return s;
1247    }
1248
1249
1250    /**
1251     * Construct complex roots for zero dimensional ideal(G).
1252     * @param I list of zero dimensional ideal with univariate irreducible
1253     *            polynomials and bi-variate polynomials.
1254     * @return list of complex algebraic roots for ideal(G)
1255     */
1256    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots(
1257                    List<IdealWithUniv<D>> I) {
1258        List<IdealWithComplexAlgebraicRoots<D>> lic = new ArrayList<IdealWithComplexAlgebraicRoots<D>>();
1259        for (IdealWithUniv<D> iu : I) {
1260            IdealWithComplexAlgebraicRoots<D> iuc = PolyUtilApp.<D> complexAlgebraicRoots(iu);
1261            //System.out.println("iuc = " + iuc);
1262            lic.add(iuc);
1263        }
1264        return lic;
1265    }
1266
1267
1268    /**
1269     * Construct exact set of complex roots for zero dimensional ideal(G).
1270     * @param I zero dimensional ideal.
1271     * @return list of coordinates of complex roots for ideal(G)
1272     */
1273    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots(
1274                    Ideal<D> I) {
1275        List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition();
1276        //System.out.println("Ir = " + Ir);
1277        List<IdealWithComplexAlgebraicRoots<D>> roots = PolyUtilApp.<D> complexAlgebraicRoots(Ir);
1278        return roots;
1279    }
1280
1281
1282    /*
1283     * Convert to a polynomial in given ring.
1284     * @param fac result polynomial ring.
1285     * @param p polynomial.
1286     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1287     *         in fac then the results are unpredictable.
1288     */
1289    static <C extends RingElem<C>> GenPolynomial<C> convert(GenPolynomialRing<C> fac, GenPolynomial<C> p) {
1290        if (fac.equals(p.factory())) {
1291            return p;
1292        }
1293        GenPolynomial<C> q = fac.parse(p.toString());
1294        if (!q.toString().equals(p.toString())) {
1295            throw new RuntimeException("convert(" + p + ") = " + q);
1296        }
1297        return q;
1298    }
1299
1300
1301    /*
1302     * Convert to a polynomial in given ring.
1303     * @param fac result polynomial ring.
1304     * @param p polynomial.
1305     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1306     *         in fac then the results are unpredictable.
1307     */
1308    static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplex(GenPolynomialRing<Complex<C>> fac,
1309                    GenPolynomial<C> p) {
1310        GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1311        if (!q.toString().equals(p.toString())) {
1312            throw new RuntimeException("convert(" + p + ") = " + q);
1313        }
1314        return q;
1315    }
1316
1317
1318    /*
1319     * Convert to a polynomial in given ring.
1320     * @param fac result polynomial ring.
1321     * @param p complex polynomial.
1322     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1323     *         in fac then the results are unpredictable.
1324     */
1325    static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplexComplex(
1326                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<Complex<C>> p) {
1327        if (fac.equals(p.factory())) {
1328            return p;
1329        }
1330        GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1331        if (!q.toString().equals(p.toString())) {
1332            throw new RuntimeException("convert(" + p + ") = " + q);
1333        }
1334        return q;
1335    }
1336
1337
1338    /**
1339     * Construct exact set of real roots for zero dimensional ideal(G).
1340     * @param I zero dimensional ideal.
1341     * @return list of coordinates of real roots for ideal(G)
1342     */
1343    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots(
1344                    Ideal<D> I) {
1345        List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition();
1346        //System.out.println("Ir = " + Ir);
1347        List<IdealWithRealAlgebraicRoots<D>> roots = PolyUtilApp.<D> realAlgebraicRoots(Ir);
1348        return roots;
1349    }
1350
1351
1352    /**
1353     * Construct primitive element for double field extension.
1354     * @param a algebraic number ring with squarefree monic minimal polynomial
1355     * @param b algebraic number ring with squarefree monic minimal polynomial
1356     * @return primitive element container with algebraic number ring c, with
1357     *         Q(c) = Q(a,b)
1358     */
1359    public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(AlgebraicNumberRing<C> a,
1360                    AlgebraicNumberRing<C> b) {
1361        GenPolynomial<C> ap = a.modul;
1362        GenPolynomial<C> bp = b.modul;
1363
1364        // setup bivariate polynomial ring
1365        String[] cv = new String[2];
1366        cv[0] = ap.ring.getVars()[0];
1367        cv[1] = bp.ring.getVars()[0];
1368        TermOrder to = new TermOrder(TermOrder.INVLEX);
1369        GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1370        GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1371        GenPolynomial<C> bs = bp.extendUnivariate(cfac, 1);
1372        List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1373        L.add(as);
1374        L.add(bs);
1375        List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1376
1377        Ideal<C> id = new Ideal<C>(cfac, L);
1378        //System.out.println("id = " + id);
1379        IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1380        //System.out.println("iu = " + iu);
1381
1382        // extract result polynomials
1383        List<GenPolynomial<C>> Np = iu.ideal.getList();
1384        //System.out.println("Np = " + Np);
1385        as = PolyUtil.<C> selectWithVariable(Np, 1);
1386        bs = PolyUtil.<C> selectWithVariable(Np, 0);
1387        GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1388        //System.out.println("as = " + as);
1389        //System.out.println("bs = " + bs);
1390        //System.out.println("cs = " + cs);
1391        String[] ev = new String[] { cs.ring.getVars()[0] };
1392        GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1393        //System.out.println("efac = " + efac);
1394        cs = cs.contractCoeff(efac);
1395        //System.out.println("cs = " + cs);
1396        as = as.reductum().contractCoeff(efac);
1397        as = as.negate();
1398        //System.out.println("as = " + as);
1399        bs = bs.reductum().contractCoeff(efac);
1400        bs = bs.negate();
1401        //System.out.println("bs = " + bs);
1402        AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1403        AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1404        AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1405        PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb, a, b);
1406        if (logger.isInfoEnabled()) {
1407            logger.info("primitive element = " + c);
1408        }
1409        return pe;
1410    }
1411
1412
1413    /**
1414     * Convert to primitive element ring.
1415     * @param cfac primitive element ring.
1416     * @param A algebraic number representing the generating element of a in the
1417     *            new ring.
1418     * @param a algebraic number to convert.
1419     * @return a converted to the primitive element ring
1420     */
1421    public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1422                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> a) {
1423        GenPolynomialRing<C> aufac = a.ring.ring;
1424        GenPolynomialRing<AlgebraicNumber<C>> ar = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, aufac);
1425        GenPolynomial<AlgebraicNumber<C>> aps = PolyUtil.<C> convertToAlgebraicCoefficients(ar, a.val);
1426        AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, A);
1427        return ac;
1428    }
1429
1430
1431    /**
1432     * Convert coefficients to primitive element ring.
1433     * @param cfac primitive element ring.
1434     * @param A algebraic number representing the generating element of a in the
1435     *            new ring.
1436     * @param a polynomial with coefficients algebraic number to convert.
1437     * @return a with coefficients converted to the primitive element ring
1438     */
1439    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1440                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, GenPolynomial<AlgebraicNumber<C>> a) {
1441        GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1442        return PolyUtil.<AlgebraicNumber<C>, AlgebraicNumber<C>> map(cr, a, new CoeffConvertAlg<C>(cfac, A));
1443    }
1444
1445
1446    /**
1447     * Convert to primitive element ring.
1448     * @param cfac primitive element ring.
1449     * @param A algebraic number representing the generating element of a in the
1450     *            new ring.
1451     * @param a recursive algebraic number to convert.
1452     * @return a converted to the primitive element ring
1453     */
1454    public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1455                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1456                    AlgebraicNumber<AlgebraicNumber<C>> a) {
1457        GenPolynomial<AlgebraicNumber<C>> aps = PolyUtilApp.<C> convertToPrimitiveElem(cfac, A, a.val);
1458        AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, B);
1459        return ac;
1460    }
1461
1462
1463    /**
1464     * Construct primitive element for double field extension.
1465     * @param b algebraic number ring with squarefree monic minimal polynomial
1466     *            over Q(a)
1467     * @return primitive element container with algebraic number ring c, with
1468     *         Q(c) = Q(a)(b)
1469     */
1470    public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(
1471                    AlgebraicNumberRing<AlgebraicNumber<C>> b) {
1472        GenPolynomial<AlgebraicNumber<C>> bp = b.modul;
1473        AlgebraicNumberRing<C> a = (AlgebraicNumberRing<C>) b.ring.coFac;
1474        GenPolynomial<C> ap = a.modul;
1475
1476        // setup bivariate polynomial ring
1477        String[] cv = new String[2];
1478        cv[0] = ap.ring.getVars()[0];
1479        cv[1] = bp.ring.getVars()[0];
1480        TermOrder to = new TermOrder(TermOrder.INVLEX);
1481        GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1482        GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(a.ring, 1,
1483                        bp.ring.getVars());
1484        GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1485        GenPolynomial<GenPolynomial<C>> bss = PolyUtil.<C> fromAlgebraicCoefficients(rfac, bp);
1486        GenPolynomial<C> bs = PolyUtil.<C> distribute(cfac, bss);
1487        List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1488        L.add(as);
1489        L.add(bs);
1490        List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1491
1492        Ideal<C> id = new Ideal<C>(cfac, L);
1493        //System.out.println("id = " + id);
1494        IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1495        //System.out.println("iu = " + iu);
1496
1497        // extract result polynomials
1498        List<GenPolynomial<C>> Np = iu.ideal.getList();
1499        as = PolyUtil.<C> selectWithVariable(Np, 1);
1500        bs = PolyUtil.<C> selectWithVariable(Np, 0);
1501        GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1502        //System.out.println("as = " + as);
1503        //System.out.println("bs = " + bs);
1504        //System.out.println("cs = " + cs);
1505        String[] ev = new String[] { cs.ring.getVars()[0] };
1506        GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1507        // System.out.println("efac = " + efac);
1508        cs = cs.contractCoeff(efac);
1509        // System.out.println("cs = " + cs);
1510        as = as.reductum().contractCoeff(efac);
1511        as = as.negate();
1512        // System.out.println("as = " + as);
1513        bs = bs.reductum().contractCoeff(efac);
1514        bs = bs.negate();
1515        //System.out.println("bs = " + bs);
1516        AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1517        AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1518        AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1519        PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb); // missing ,a,b);
1520        if (logger.isInfoEnabled()) {
1521            logger.info("primitive element = " + pe);
1522        }
1523        return pe;
1524    }
1525
1526
1527    /**
1528     * Convert to primitive element ring.
1529     * @param cfac primitive element ring.
1530     * @param A algebraic number representing the generating element of a in the
1531     *            new ring.
1532     * @param a polynomial with recursive algebraic number coefficients to
1533     *            convert.
1534     * @return a converted to the primitive element ring
1535     */
1536    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1537                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1538                    GenPolynomial<AlgebraicNumber<AlgebraicNumber<C>>> a) {
1539        GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1540        return PolyUtil.<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> map(cr, a,
1541                        new CoeffRecConvertAlg<C>(cfac, A, B));
1542    }
1543
1544
1545    /**
1546     * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1547     * RealAlgebraicNumber<C> coefficients from package
1548     * <code>edu.jas.root.</code>
1549     * @param afac result polynomial factory.
1550     * @param A polynomial with RealAlgebraicNumber&lt;C&gt; coefficients to be
1551     *            converted.
1552     * @return polynomial with RealAlgebraicNumber&lt;C&gt; coefficients.
1553     */
1554    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> realAlgFromRealCoefficients(
1555                    GenPolynomialRing<edu.jas.root.RealAlgebraicNumber<C>> afac,
1556                    GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> A) {
1557        edu.jas.root.RealAlgebraicRing<C> cfac = (edu.jas.root.RealAlgebraicRing<C>) afac.coFac;
1558        return PolyUtil.<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> map(
1559                        afac, A, new ReAlgFromRealCoeff<C>(cfac));
1560    }
1561
1562
1563    /**
1564     * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1565     * RealAlgebraicNumber<C> coefficients from package
1566     * 
1567     * <pre>
1568     * edu.jas.application
1569     * </pre>
1570     * 
1571     * .
1572     * @param rfac result polynomial factory.
1573     * @param A polynomial with RealAlgebraicNumber&lt;C&gt; coefficients to be
1574     *            converted.
1575     * @return polynomial with RealAlgebraicNumber&lt;C&gt; coefficients.
1576     */
1577    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> realFromRealAlgCoefficients(
1578                    GenPolynomialRing<edu.jas.application.RealAlgebraicNumber<C>> rfac,
1579                    GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> A) {
1580        edu.jas.application.RealAlgebraicRing<C> cfac = (edu.jas.application.RealAlgebraicRing<C>) rfac.coFac;
1581        return PolyUtil.<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> map(
1582                        rfac, A, new RealFromReAlgCoeff<C>(cfac));
1583    }
1584
1585
1586    /**
1587     * Convert to Complex&lt;RealAlgebraicNumber&gt; coefficients. Represent as
1588     * polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients, C is
1589     * e.g. BigRational.
1590     * @param pfac result polynomial factory.
1591     * @param A polynomial with Complex coefficients to be converted.
1592     * @return polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients.
1593     */
1594    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> convertToComplexRealCoefficients(
1595                    GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac,
1596                    GenPolynomial<Complex<C>> A) {
1597        ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> afac;
1598        afac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) pfac.coFac;
1599        return PolyUtil.<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(pfac, A,
1600                        new CoeffToComplexReal<C>(afac));
1601    }
1602
1603
1604    /**
1605     * Evaluate to Complex&lt;RealAlgebraicNumber&gt; coefficients. Represent as
1606     * polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients, C is
1607     * e.g. BigRational.
1608     * @param pfac result polynomial factory.
1609     * @param A = A(x,Y) a recursive polynomial with
1610     *            GenPolynomial&lt;Complex&gt; coefficients to be converted.
1611     * @param r Complex&lt;RealAlgebraicNumber&gt; to be evaluated at.
1612     * @return A(r,Y), a polynomial with Complex&lt;RealAlgebraicNumber&gt;
1613     *         coefficients.
1614     */
1615    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateToComplexRealCoefficients(
1616                    GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac,
1617                    GenPolynomial<GenPolynomial<Complex<C>>> A,
1618                    Complex<edu.jas.application.RealAlgebraicNumber<C>> r) {
1619        return PolyUtil.<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(
1620                        pfac, A, new EvaluateToComplexReal<C>(pfac, r));
1621    }
1622
1623
1624}
1625
1626
1627/**
1628 * Coefficient to convert algebriac functor.
1629 */
1630class CoeffConvertAlg<C extends GcdRingElem<C>>
1631                implements UnaryFunctor<AlgebraicNumber<C>, AlgebraicNumber<C>> {
1632
1633
1634    final protected AlgebraicNumberRing<C> afac;
1635
1636
1637    final protected AlgebraicNumber<C> A;
1638
1639
1640    public CoeffConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a) {
1641        if (fac == null || a == null) {
1642            throw new IllegalArgumentException("fac and a must not be null");
1643        }
1644        afac = fac;
1645        A = a;
1646    }
1647
1648
1649    public AlgebraicNumber<C> eval(AlgebraicNumber<C> c) {
1650        if (c == null) {
1651            return afac.getZERO();
1652        }
1653        return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, c);
1654    }
1655}
1656
1657
1658/**
1659 * Coefficient recursive to convert algebriac functor.
1660 */
1661class CoeffRecConvertAlg<C extends GcdRingElem<C>>
1662                implements UnaryFunctor<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> {
1663
1664
1665    final protected AlgebraicNumberRing<C> afac;
1666
1667
1668    final protected AlgebraicNumber<C> A;
1669
1670
1671    final protected AlgebraicNumber<C> B;
1672
1673
1674    public CoeffRecConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a, AlgebraicNumber<C> b) {
1675        if (fac == null || a == null || b == null) {
1676            throw new IllegalArgumentException("fac, a and b must not be null");
1677        }
1678        afac = fac;
1679        A = a;
1680        B = b;
1681    }
1682
1683
1684    public AlgebraicNumber<C> eval(AlgebraicNumber<AlgebraicNumber<C>> c) {
1685        if (c == null) {
1686            return afac.getZERO();
1687        }
1688        return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, B, c);
1689    }
1690}
1691
1692
1693/**
1694 * Coefficient to real algebriac from real algebraic functor.
1695 */
1696class ReAlgFromRealCoeff<C extends GcdRingElem<C> & Rational> implements
1697                UnaryFunctor<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> {
1698
1699
1700    final protected edu.jas.root.RealAlgebraicRing<C> afac;
1701
1702
1703    public ReAlgFromRealCoeff(edu.jas.root.RealAlgebraicRing<C> fac) {
1704        if (fac == null) {
1705            throw new IllegalArgumentException("fac must not be null");
1706        }
1707        afac = fac;
1708    }
1709
1710
1711    @SuppressWarnings("unchecked")
1712    public edu.jas.root.RealAlgebraicNumber<C> eval(edu.jas.application.RealAlgebraicNumber<C> c) {
1713        if (c == null) {
1714            return afac.getZERO();
1715        }
1716        return (edu.jas.root.RealAlgebraicNumber<C>) (Object) c.number; // force ignore recursion
1717    }
1718}
1719
1720
1721/**
1722 * Coefficient to real algebriac from algebraic functor.
1723 */
1724class RealFromReAlgCoeff<C extends GcdRingElem<C> & Rational> implements
1725                UnaryFunctor<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> {
1726
1727
1728    final protected edu.jas.application.RealAlgebraicRing<C> rfac;
1729
1730
1731    public RealFromReAlgCoeff(edu.jas.application.RealAlgebraicRing<C> fac) {
1732        if (fac == null) {
1733            throw new IllegalArgumentException("fac must not be null");
1734        }
1735        rfac = fac;
1736    }
1737
1738
1739    @SuppressWarnings("unchecked")
1740    public edu.jas.application.RealAlgebraicNumber<C> eval(edu.jas.root.RealAlgebraicNumber<C> c) {
1741        if (c == null) {
1742            return rfac.getZERO();
1743        }
1744        edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>> rrc = (edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>>) (Object) c; // force resurrect recursion
1745        return new edu.jas.application.RealAlgebraicNumber<C>(rfac, rrc);
1746    }
1747}
1748
1749
1750/**
1751 * Coefficient to complex real algebriac functor.
1752 */
1753class CoeffToComplexReal<C extends GcdRingElem<C> & Rational>
1754                implements UnaryFunctor<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> {
1755
1756
1757    final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac;
1758
1759
1760    final edu.jas.application.RealAlgebraicRing<C> afac;
1761
1762
1763    final GenPolynomialRing<C> pfac;
1764
1765
1766    public CoeffToComplexReal(ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> fac) {
1767        if (fac == null) {
1768            throw new IllegalArgumentException("fac must not be null");
1769        }
1770        cfac = fac;
1771        afac = (edu.jas.application.RealAlgebraicRing<C>) cfac.ring;
1772        pfac = afac.univs.ideal.getRing();
1773    }
1774
1775
1776    public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(Complex<C> c) {
1777        if (c == null) {
1778            return cfac.getZERO();
1779        }
1780        GenPolynomial<C> pr, pi;
1781        pr = new GenPolynomial<C>(pfac, c.getRe());
1782        pi = new GenPolynomial<C>(pfac, c.getIm());
1783        //System.out.println("pr = " + pr);
1784        //System.out.println("pi = " + pi);
1785        edu.jas.application.RealAlgebraicNumber<C> re, im;
1786        re = new edu.jas.application.RealAlgebraicNumber<C>(afac, pr);
1787        im = new edu.jas.application.RealAlgebraicNumber<C>(afac, pi);
1788        //System.out.println("re = " + re);
1789        //System.out.println("im = " + im);
1790        return new Complex<edu.jas.application.RealAlgebraicNumber<C>>(cfac, re, im);
1791    }
1792}
1793
1794
1795/**
1796 * Polynomial coefficient to complex real algebriac evaluation functor.
1797 */
1798class EvaluateToComplexReal<C extends GcdRingElem<C> & Rational> implements
1799                UnaryFunctor<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> {
1800
1801
1802    final protected GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac;
1803
1804
1805    final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac;
1806
1807
1808    final protected Complex<edu.jas.application.RealAlgebraicNumber<C>> root;
1809
1810
1811    public EvaluateToComplexReal(GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> fac,
1812                    Complex<edu.jas.application.RealAlgebraicNumber<C>> r) {
1813        if (fac == null) {
1814            throw new IllegalArgumentException("fac must not be null");
1815        }
1816        if (r == null) {
1817            throw new IllegalArgumentException("r must not be null");
1818        }
1819        pfac = fac;
1820        cfac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) fac.coFac;
1821        root = r;
1822        //System.out.println("cfac  = " + cfac);
1823        //System.out.println("root  = " + root);
1824    }
1825
1826
1827    public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(GenPolynomial<Complex<C>> c) {
1828        if (c == null) {
1829            return cfac.getZERO();
1830        }
1831        //System.out.println("c  = " + c);
1832        GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> cp;
1833        cp = PolyUtilApp.<C> convertToComplexRealCoefficients(pfac, c);
1834        Complex<edu.jas.application.RealAlgebraicNumber<C>> cr;
1835        cr = PolyUtil.<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateMain(cfac, cp, root);
1836        return cr;
1837    }
1838}