001/*
002 * $Id: PolyUtilApp.java 5313 2015-10-26 22:46:38Z 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.Set;
013import java.util.SortedMap;
014import java.util.TreeMap;
015import java.util.TreeSet;
016
017import org.apache.log4j.Logger;
018
019import edu.jas.arith.BigDecimal;
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 = Logger.getLogger(PolyUtilApp.class);
062
063
064    private static 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, D 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, D 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, D 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, D 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, D 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                    D 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, D 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, D 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, D 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                    D 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("no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
708            }
709            //System.out.println("i   = " + i);
710            //System.out.println("pi  = " + pi);
711            if (logger.isInfoEnabled()) {
712                logger.info("pi  = " + pi);
713                logger.info("pip = " + pip);
714            }
715            int[] depi = pip.degreeVector().dependencyOnVariables();
716            //System.out.println("depi = " + Arrays.toString(depi));
717            if (depi.length < 1 || depi.length > 2) {
718                throw new RuntimeException("wrong number of variables " + Arrays.toString(depi));
719            }
720            rra = RootFactory.<D> realAlgebraicNumbersIrred(pi);
721            if (logger.isInfoEnabled()) {
722                List<Interval<D>> il = new ArrayList<Interval<D>>();
723                for (RealAlgebraicNumber<D> rr : rra) {
724                    il.add(rr.ring.getRoot());
725                }
726                logger.info("roots(pi) = " + il);
727            }
728            if (depi.length == 1) {
729                // all combinations are roots of the ideal I
730                for (RealAlgebraicNumber<D> rr : rra) {
731                    //System.out.println("rr.ring = " + rr.ring);
732                    for (List<RealAlgebraicNumber<D>> rx : ran) {
733                        //System.out.println("rx = " + rx);
734                        List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
735                        ry.addAll(rx);
736                        ry.add(rr);
737                        rn.add(ry);
738                    }
739                }
740            } else { // depi.length == 2
741                // select roots of the ideal I
742                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
743                //System.out.println("pip2 = " + pip2.ring);
744                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
745                TermOrder to = new TermOrder(TermOrder.INVLEX);
746                GenPolynomialRing<GenPolynomial<D>> rfac = new GenPolynomialRing<GenPolynomial<D>>(ufac, 1,
747                                to);
748                GenPolynomial<GenPolynomial<D>> pip2r = PolyUtil.<D> recursive(rfac, pip2);
749                int ix = fac.nvar - 1 - depi[depi.length - 1];
750                //System.out.println("ix = " + ix);
751                for (RealAlgebraicNumber<D> rr : rra) {
752                    //System.out.println("rr.ring = " + rr.ring);
753                    Interval<D> rroot = rr.ring.getRoot();
754                    GenPolynomial<D> pip2el = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.left);
755                    GenPolynomial<D> pip2er = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.right);
756                    GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
757                    GenPolynomial<D> pip2elc = convert(upfac, pip2el);
758                    GenPolynomial<D> pip2erc = convert(upfac, pip2er);
759                    //System.out.println("pip2elc = " + pip2elc);
760                    //System.out.println("pip2erc = " + pip2erc);
761                    for (List<RealAlgebraicNumber<D>> rx : ran) {
762                        //System.out.println("rx = " + rx);
763                        RealAlgebraicRing<D> rar = rx.get(ix).ring;
764                        //System.out.println("rar = " + rar.toScript());
765                        RealAlgebraicNumber<D> rel = new RealAlgebraicNumber<D>(rar, pip2elc);
766                        RealAlgebraicNumber<D> rer = new RealAlgebraicNumber<D>(rar, pip2erc);
767                        int sl = rel.signum();
768                        int sr = rer.signum();
769                        //System.out.println("sl = " + sl + ", sr = " + sr + ", sl*sr = " + (sl*sr));
770                        if (sl * sr <= 0) {
771                            //System.out.println("sl * sr <= 0: rar = " + rar.toScript());
772                            List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
773                            ry.addAll(rx);
774                            ry.add(rr);
775                            rn.add(ry);
776                        }
777                    }
778                }
779            }
780            ran = rn;
781        }
782        if (logger.isInfoEnabled()) {
783            for (List<RealAlgebraicNumber<D>> rz : ran) {
784                List<Interval<D>> il = new ArrayList<Interval<D>>();
785                for (RealAlgebraicNumber<D> rr : rz) {
786                    il.add(rr.ring.getRoot());
787                }
788                logger.info("root-tuple = " + il);
789            }
790        }
791        IdealWithRealAlgebraicRoots<D> Ir = new IdealWithRealAlgebraicRoots<D>(I, ran);
792        return Ir;
793    }
794
795
796    /**
797     * Construct real roots for zero dimensional ideal(G).
798     * @param I list of zero dimensional ideal with univariate irreducible
799     *            polynomials and bi-variate polynomials.
800     * @return list of real algebraic roots for all ideal(I_i)
801     */
802    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots(
803                    List<IdealWithUniv<D>> I) {
804        List<IdealWithRealAlgebraicRoots<D>> lir = new ArrayList<IdealWithRealAlgebraicRoots<D>>(I.size());
805        for (IdealWithUniv<D> iu : I) {
806            IdealWithRealAlgebraicRoots<D> iur = PolyUtilApp.<D> realAlgebraicRoots(iu);
807            //System.out.println("iur = " + iur);
808            lir.add(iur);
809        }
810        return lir;
811    }
812
813
814    /**
815     * Construct complex roots for zero dimensional ideal(G).
816     * @param I zero dimensional ideal with univariate irreducible polynomials
817     *            and bi-variate polynomials.
818     * @return complex algebraic roots for ideal(G) <b>Note:</b> implementation
819     *         contains errors, do not use.
820     */
821    public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRootsWrong( // Wrong
822                    IdealWithUniv<D> I) {
823        List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can;
824        can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
825        if (I == null) {
826            throw new IllegalArgumentException("null ideal not permitted");
827        }
828        if (I.ideal == null || I.upolys == null) {
829            throw new IllegalArgumentException("null ideal components not permitted " + I);
830        }
831        if (I.ideal.isZERO() || I.upolys.size() == 0) {
832            return new IdealWithComplexAlgebraicRoots<D>(I, can);
833        }
834        GenPolynomialRing<D> fac = I.ideal.list.ring;
835        if (fac.nvar == 0) {
836            return new IdealWithComplexAlgebraicRoots<D>(I, can);
837        }
838        if (fac.nvar != I.upolys.size()) {
839            throw new IllegalArgumentException("ideal not zero dimnsional: " + I);
840        }
841        // case i == 0:
842        GenPolynomial<D> p0 = I.upolys.get(0);
843        GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
844        if (p0p == null) {
845            throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of  " + I.ideal);
846        }
847        if (logger.isInfoEnabled()) {
848            logger.info("p0  = " + p0);
849            logger.info("p0p = " + p0p);
850        }
851        int[] dep0 = p0p.degreeVector().dependencyOnVariables();
852        //System.out.println("dep0 = " + Arrays.toString(dep0));
853        if (dep0.length != 1) {
854            throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
855        }
856        RingFactory<D> cfac = p0.ring.coFac;
857        ComplexRing<D> ccfac = new ComplexRing<D>(cfac);
858        GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring);
859        GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0);
860        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra;
861        cra = edu.jas.application.RootFactory.<D> complexAlgebraicNumbersSquarefree(p0c);
862        logger.info("#roots(p0c) = " + cra.size());
863        if (debug) {
864            boolean t = edu.jas.application.RootFactory.<D> isRoot(p0c, cra);
865            if (!t) {
866                throw new RuntimeException("no roots of " + p0c);
867            }
868        }
869        for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
870            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl;
871            cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
872            cl.add(cr);
873            can.add(cl);
874        }
875        if (fac.nvar == 1) {
876            return new IdealWithComplexAlgebraicRoots<D>(I, can);
877        }
878        // case i > 0:
879        for (int i = 1; i < I.upolys.size(); i++) {
880            List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn;
881            cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
882            GenPolynomial<D> pi = I.upolys.get(i);
883            GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
884            if (pip == null) {
885                throw new RuntimeException("no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
886            }
887            if (logger.isInfoEnabled()) {
888                logger.info("pi(" + i + ") = " + pi);
889                logger.info("pip  = " + pip);
890            }
891            facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring);
892            GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi);
893            int[] depi = pip.degreeVector().dependencyOnVariables();
894            //System.out.println("depi = " + Arrays.toString(depi));
895            if (depi.length < 1 || depi.length > 2) {
896                throw new RuntimeException("wrong number of variables " + Arrays.toString(depi) + " for "
897                                + pip);
898            }
899            cra = edu.jas.application.RootFactory.<D> complexAlgebraicNumbersSquarefree(pic);
900            logger.info("#roots(pic) = " + cra.size());
901            if (debug) {
902                boolean t = edu.jas.application.RootFactory.<D> isRoot(pic, cra);
903                if (!t) {
904                    throw new RuntimeException("no roots of " + pic);
905                }
906            }
907            if (depi.length == 1) {
908                // all combinations are roots of the ideal I
909                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
910                    //System.out.println("cr.ring = " + cr.ring);
911                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
912                        //System.out.println("cx = " + cx);
913                        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
914                        cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
915                        cy.addAll(cx);
916                        cy.add(cr);
917                        cn.add(cy);
918                    }
919                }
920            } else { // depi.length == 2
921                // select roots of the ideal I
922                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
923                GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1);
924                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
925                GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac);
926                GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring);
927                GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2);
928                //System.out.println("pip2c = " + pip2c);
929                GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac;
930                rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac);
931                GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac,
932                                pip2c);
933                //System.out.println("pip2cr = " + pip2cr);
934
935                int ix = fac.nvar - 1 - depi[depi.length - 1];
936                //System.out.println("ix = " + ix);
937                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
938                    System.out.println("cr = " + toString(cr)); // <----------------------------------
939                    edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring;
940                    RealRootTuple<D> rroot = cring.getRoot();
941                    List<RealAlgebraicNumber<D>> rlist = rroot.tuple;
942                    Interval<D> vr = rlist.get(0).ring.getRoot();
943                    Interval<D> vi = rlist.get(1).ring.getRoot();
944                    logger.info("vr = " + vr + ", vi = " + vi);
945                    if (vr.length().isZERO()) {
946                        D e = vr.left.factory().parse("1/2");
947                        D m = vr.left; //middle();
948                        vr = new Interval<D>(m.subtract(e), m.sum(e));
949                        logger.info("|vr| == 0: " + vr);
950                    }
951                    if (vi.length().isZERO()) {
952                        D e = vi.left.factory().parse("1/2");
953                        D m = vi.left; //middle();
954                        vi = new Interval<D>(m.subtract(e), m.sum(e));
955                        logger.info("|vi| == 0: " + vi);
956                    }
957                    Complex<D> sw = new Complex<D>(ccfac, vr.left, vi.left);
958                    Complex<D> ne = new Complex<D>(ccfac, vr.right, vi.right);
959                    logger.info("sw   = " + toString1(sw) + ", ne   = " + toString1(ne));
960                    GenPolynomial<Complex<D>> pip2cesw, pip2cene;
961                    pip2cesw = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, sw);
962                    pip2cene = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, ne);
963                    GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
964                    GenPolynomialRing<Complex<D>> upcfac = new GenPolynomialRing<Complex<D>>(ccfac, upfac);
965                    //System.out.println("upfac = " + upfac);
966                    //System.out.println("upcfac = " + upcfac);
967                    GenPolynomial<Complex<D>> pip2eswc = convertComplexComplex(upcfac, pip2cesw);
968                    GenPolynomial<Complex<D>> pip2enec = convertComplexComplex(upcfac, pip2cene);
969                    //System.out.println("pip2eswc = " + pip2eswc);
970                    //System.out.println("pip2enec = " + pip2enec);
971                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
972                        //System.out.println("cxi = " + toString(cx.get(ix)));
973                        Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix);
974                        ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring;
975                        edu.jas.application.RealAlgebraicRing<D> rar = (edu.jas.application.RealAlgebraicRing<D>) car.ring;
976                        //System.out.println("car = " + car);
977                        //System.out.println("rar = " + rar);
978                        TermOrder to = new TermOrder(TermOrder.INVLEX);
979                        String vvr = rar.algebraic.ring.getVars()[0];
980                        String vvi = rar.algebraic.ring.getVars()[1];
981                        String[] vars = new String[] { vvr, vvi };
982                        GenPolynomialRing<Complex<D>> tfac = new GenPolynomialRing<Complex<D>>(ccfac, to,
983                                        vars);
984                        GenPolynomial<Complex<D>> t = tfac.univariate(1, 1L).sum(
985                                        tfac.univariate(0, 1L).multiply(ccfac.getIMAG()));
986                        //System.out.println("t  = " + t); // t = x + i y
987                        GenPolynomialRing<D> rtfac = new GenPolynomialRing<D>(cfac, tfac);
988                        GenPolynomial<Complex<D>> su;
989                        GenPolynomial<D> re, im;
990                        su = PolyUtil.<Complex<D>> substituteUnivariate(pip2eswc, t);
991                        //su = su.monic(); not here
992                        re = PolyUtil.<D> realPartFromComplex(rtfac, su);
993                        im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
994                        //System.out.println("re = " + re);
995                        //System.out.println("im = " + im);
996                        edu.jas.application.RealAlgebraicNumber<D> resw, imsw, rene, imne;
997                        resw = new edu.jas.application.RealAlgebraicNumber<D>(rar, re);
998                        //System.out.println("resw = " + resw);
999                        int sswr = resw.signum();
1000                        imsw = new edu.jas.application.RealAlgebraicNumber<D>(rar, im);
1001                        //System.out.println("imsw = " + imsw);
1002                        int sswi = imsw.signum();
1003                        su = PolyUtil.<Complex<D>> substituteUnivariate(pip2enec, t);
1004                        //su = su.monic(); not here
1005                        re = PolyUtil.<D> realPartFromComplex(rtfac, su);
1006                        im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
1007                        //System.out.println("re = " + re);
1008                        //System.out.println("im = " + im);
1009                        rene = new edu.jas.application.RealAlgebraicNumber<D>(rar, re);
1010                        //System.out.println("rene = " + rene);
1011                        int sner = rene.signum();
1012                        imne = new edu.jas.application.RealAlgebraicNumber<D>(rar, im);
1013                        //System.out.println("imne = " + imne);
1014                        int snei = imne.signum();
1015                        //System.out.println("sswr = " + sswr + ", sswi = " + sswi);
1016                        //System.out.println("sner = " + sner + ", snei = " + snei);
1017                        if ((sswr * sner <= 0 && sswi * snei <= 0)) { // wrong !
1018                            logger.info("   hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1019                            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1020                            cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1021                            cy.addAll(cx);
1022                            cy.add(cr);
1023                            cn.add(cy);
1024                        } else {
1025                            logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1026                        }
1027                    }
1028                }
1029            }
1030            can = cn;
1031        }
1032        IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can);
1033        return Ic;
1034    }
1035
1036
1037    /**
1038     * Construct complex roots for zero dimensional ideal(G).
1039     * @param I zero dimensional ideal with univariate irreducible polynomials
1040     *            and bi-variate polynomials.
1041     * @return complex algebraic roots for ideal(G) <b>Note:</b> not jet
1042     *         completed oin all cases.
1043     */
1044    public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRoots(
1045                    IdealWithUniv<D> I) {
1046        List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can;
1047        can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
1048        if (I == null) {
1049            throw new IllegalArgumentException("null ideal not permitted");
1050        }
1051        if (I.ideal == null || I.upolys == null) {
1052            throw new IllegalArgumentException("null ideal components not permitted " + I);
1053        }
1054        if (I.ideal.isZERO() || I.upolys.size() == 0) {
1055            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1056        }
1057        GenPolynomialRing<D> fac = I.ideal.list.ring;
1058        if (fac.nvar == 0) {
1059            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1060        }
1061        if (fac.nvar != I.upolys.size()) {
1062            throw new IllegalArgumentException("ideal not zero dimnsional: " + I);
1063        }
1064        // case i == 0:
1065        GenPolynomial<D> p0 = I.upolys.get(0);
1066        GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
1067        if (p0p == null) {
1068            throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of  " + I.ideal);
1069        }
1070        if (logger.isInfoEnabled()) {
1071            logger.info("p0  = " + p0);
1072            logger.info("p0p = " + p0p);
1073        }
1074        int[] dep0 = p0p.degreeVector().dependencyOnVariables();
1075        //System.out.println("dep0 = " + Arrays.toString(dep0));
1076        if (dep0.length != 1) {
1077            throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
1078        }
1079        RingFactory<D> cfac = p0.ring.coFac;
1080        ComplexRing<D> ccfac = new ComplexRing<D>(cfac);
1081        GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring);
1082        GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0);
1083        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra;
1084        cra = edu.jas.application.RootFactory.<D> complexAlgebraicNumbersSquarefree(p0c);
1085        logger.info("#roots(p0c) = " + cra.size());
1086        for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1087            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl;
1088            cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1089            cl.add(cr);
1090            can.add(cl);
1091        }
1092        if (fac.nvar == 1) {
1093            return new IdealWithComplexAlgebraicRoots<D>(I, can);
1094        }
1095        // case i > 0:
1096        for (int i = 1; i < I.upolys.size(); i++) {
1097            List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn;
1098            cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
1099            GenPolynomial<D> pi = I.upolys.get(i);
1100            GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
1101            if (pip == null) {
1102                throw new RuntimeException("no polynomial found in " + (fac.nvar - 1 - i) + " of  " + I.ideal);
1103            }
1104            if (logger.isInfoEnabled()) {
1105                logger.info("pi(" + i + ") = " + pi);
1106                logger.info("pip  = " + pip);
1107            }
1108            facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring);
1109            GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi);
1110            int[] depi = pip.degreeVector().dependencyOnVariables();
1111            //System.out.println("depi = " + Arrays.toString(depi));
1112            if (depi.length < 1 || depi.length > 2) {
1113                throw new RuntimeException("wrong number of variables " + Arrays.toString(depi) + " for "
1114                                + pip);
1115            }
1116            cra = edu.jas.application.RootFactory.<D> complexAlgebraicNumbersSquarefree(pic);
1117            logger.info("#roots(pic) = " + cra.size());
1118            if (depi.length == 1) {
1119                // all combinations are roots of the ideal I
1120                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1121                    //System.out.println("cr.ring = " + cr.ring);
1122                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
1123                        //System.out.println("cx = " + cx);
1124                        List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1125                        cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1126                        cy.addAll(cx);
1127                        cy.add(cr);
1128                        cn.add(cy);
1129                    }
1130                }
1131            } else { // depi.length == 2
1132                // select roots of the ideal I
1133                GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip);
1134                pip2 = PolyUtil.<D> removeUnusedLowerVariables(pip2);
1135                pip2 = PolyUtil.<D> removeUnusedMiddleVariables(pip2);
1136                GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1);
1137                GenPolynomialRing<D> ufac = pip2.ring.contract(1);
1138                GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac);
1139                GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring);
1140                GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2);
1141                GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac;
1142                rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac);
1143                GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac,
1144                                pip2c);
1145                //System.out.println("pip2cr = " + pip2cr);
1146                int ix = fac.nvar - 1 - depi[depi.length - 1];
1147                //System.out.println("ix = " + ix);
1148                for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
1149                    //System.out.println("cr = " + toString(cr)); 
1150                    edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring;
1151                    RealRootTuple<D> rroot = cring.getRoot();
1152                    List<RealAlgebraicNumber<D>> rlist = rroot.tuple;
1153                    //System.out.println("rlist = " + rlist);
1154                    Interval<D> vr = rlist.get(0).ring.getRoot();
1155                    Interval<D> vi = rlist.get(1).ring.getRoot();
1156                    //logger.info("vr = " + vr + ", vi = " + vi);
1157                    edu.jas.application.RealAlgebraicNumber<D> vrl, vil, vrr, vir;
1158                    vrl = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.left);
1159                    vil = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.left);
1160                    vrr = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.right);
1161                    vir = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.right);
1162                    ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> crr;
1163                    crr = new ComplexRing<edu.jas.application.RealAlgebraicNumber<D>>(cring);
1164                    Complex<edu.jas.application.RealAlgebraicNumber<D>> csw, cne;
1165                    csw = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrl, vil);
1166                    cne = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrr, vir);
1167                    //logger.info("csw  = " + toString(csw)   + ", cne  = " + toString(cne));
1168                    Rectangle<edu.jas.application.RealAlgebraicNumber<D>> rec;
1169                    rec = new Rectangle<edu.jas.application.RealAlgebraicNumber<D>>(csw, cne);
1170                    //System.out.println("rec = " + rec);
1171                    for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
1172                        Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix);
1173                        //System.out.println("cax = " + toString(cax));
1174                        ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring;
1175                        //System.out.println("car = " + car);
1176                        GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcrfac;
1177                        pcrfac = new GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(
1178                                        car, rcfac);
1179                        GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcr;
1180                        pcr = evaluateToComplexRealCoefficients(pcrfac, pip2cr, cax);
1181                        //System.out.println("pcr = " + pcr);
1182                        ComplexRoots<edu.jas.application.RealAlgebraicNumber<D>> rengine;
1183                        rengine = new ComplexRootsSturm<edu.jas.application.RealAlgebraicNumber<D>>(car);
1184                        long nr = 0;
1185                        try {
1186                            nr = rengine.complexRootCount(rec, pcr);
1187                            //logger.info("rootCount = " + nr);
1188                        } catch (InvalidBoundaryException e) {
1189                            e.printStackTrace();
1190                        }
1191                        if (nr == 1) { // one root
1192                            logger.info("   hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1193                            List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy;
1194                            cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
1195                            cy.addAll(cx);
1196                            cy.add(cr);
1197                            cn.add(cy);
1198                        } else if (nr > 1) {
1199                            logger.error("to many roots, cxi = " + toString(cx.get(ix)) + ", cr = "
1200                                            + toString(cr));
1201                        } else { // no root
1202                            logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr));
1203                        }
1204                    }
1205                }
1206            }
1207            can = cn;
1208        }
1209        IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can);
1210        return Ic;
1211    }
1212
1213
1214    /**
1215     * String representation of a deximal approximation of a complex number.
1216     * @param c compelx number.
1217     * @return String representation of c
1218     */
1219    public static <D extends GcdRingElem<D> & Rational> String toString(
1220                    Complex<edu.jas.application.RealAlgebraicNumber<D>> c) {
1221        edu.jas.application.RealAlgebraicNumber<D> re = c.getRe();
1222        edu.jas.application.RealAlgebraicNumber<D> im = c.getIm();
1223        String s = re.decimalMagnitude().toString();
1224        if (!im.isZERO()) {
1225            s = s + "i" + im.decimalMagnitude();
1226        }
1227        return s;
1228    }
1229
1230
1231    /**
1232     * String representation of a deximal approximation of a complex number.
1233     * @param c compelx number.
1234     * @return String representation of c
1235     */
1236    public static <D extends GcdRingElem<D> & Rational> String toString1(Complex<D> c) {
1237        D re = c.getRe();
1238        D im = c.getIm();
1239        String s = new BigDecimal(re.getRational()).toString();
1240        if (!im.isZERO()) {
1241            s = s + "i" + new BigDecimal(im.getRational());
1242        }
1243        return s;
1244    }
1245
1246
1247    /**
1248     * Construct complex roots for zero dimensional ideal(G).
1249     * @param I list of zero dimensional ideal with univariate irreducible
1250     *            polynomials and bi-variate polynomials.
1251     * @return list of complex algebraic roots for ideal(G)
1252     */
1253    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots(
1254                    List<IdealWithUniv<D>> I) {
1255        List<IdealWithComplexAlgebraicRoots<D>> lic = new ArrayList<IdealWithComplexAlgebraicRoots<D>>();
1256        for (IdealWithUniv<D> iu : I) {
1257            IdealWithComplexAlgebraicRoots<D> iuc = PolyUtilApp.<D> complexAlgebraicRoots(iu);
1258            //System.out.println("iuc = " + iuc);
1259            lic.add(iuc);
1260        }
1261        return lic;
1262    }
1263
1264
1265    /**
1266     * Construct exact set of complex roots for zero dimensional ideal(G).
1267     * @param I zero dimensional ideal.
1268     * @return list of coordinates of complex roots for ideal(G)
1269     */
1270    public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots(
1271                    Ideal<D> I) {
1272        List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition();
1273        //System.out.println("Ir = " + Ir);
1274        List<IdealWithComplexAlgebraicRoots<D>> roots = PolyUtilApp.<D> complexAlgebraicRoots(Ir);
1275        return roots;
1276    }
1277
1278
1279    /*
1280     * Convert to a polynomial in given ring.
1281     * @param fac result polynomial ring.
1282     * @param p polynomial.
1283     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1284     *         in fac then the results are unpredictable.
1285     */
1286    static <C extends RingElem<C>> GenPolynomial<C> convert(GenPolynomialRing<C> fac, GenPolynomial<C> p) {
1287        if (fac.equals(p.factory())) {
1288            return p;
1289        }
1290        GenPolynomial<C> q = fac.parse(p.toString());
1291        if (!q.toString().equals(p.toString())) {
1292            throw new RuntimeException("convert(" + p + ") = " + q);
1293        }
1294        return q;
1295    }
1296
1297
1298    /*
1299     * Convert to a polynomial in given ring.
1300     * @param fac result polynomial ring.
1301     * @param p polynomial.
1302     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1303     *         in fac then the results are unpredictable.
1304     */
1305    static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplex(
1306                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<C> p) {
1307        GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1308        if (!q.toString().equals(p.toString())) {
1309            throw new RuntimeException("convert(" + p + ") = " + q);
1310        }
1311        return q;
1312    }
1313
1314
1315    /*
1316     * Convert to a polynomial in given ring.
1317     * @param fac result polynomial ring.
1318     * @param p complex polynomial.
1319     * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1320     *         in fac then the results are unpredictable.
1321     */
1322    static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplexComplex(
1323                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<Complex<C>> p) {
1324        if (fac.equals(p.factory())) {
1325            return p;
1326        }
1327        GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1328        if (!q.toString().equals(p.toString())) {
1329            throw new RuntimeException("convert(" + p + ") = " + q);
1330        }
1331        return q;
1332    }
1333
1334
1335    /**
1336     * Construct exact set of real roots for zero dimensional ideal(G).
1337     * @param I zero dimensional ideal.
1338     * @return list of coordinates of real roots for ideal(G)
1339     */
1340    public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots(
1341                    Ideal<D> I) {
1342        List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition();
1343        //System.out.println("Ir = " + Ir);
1344        List<IdealWithRealAlgebraicRoots<D>> roots = PolyUtilApp.<D> realAlgebraicRoots(Ir);
1345        return roots;
1346    }
1347
1348
1349    /**
1350     * Construct primitive element for double field extension.
1351     * @param a algebraic number ring with squarefree monic minimal polynomial
1352     * @param b algebraic number ring with squarefree monic minimal polynomial
1353     * @return primitive element container with algebraic number ring c, with
1354     *         Q(c) = Q(a,b)
1355     */
1356    public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(AlgebraicNumberRing<C> a,
1357                    AlgebraicNumberRing<C> b) {
1358        GenPolynomial<C> ap = a.modul;
1359        GenPolynomial<C> bp = b.modul;
1360
1361        // setup bivariate polynomial ring
1362        String[] cv = new String[2];
1363        cv[0] = ap.ring.getVars()[0];
1364        cv[1] = bp.ring.getVars()[0];
1365        TermOrder to = new TermOrder(TermOrder.INVLEX);
1366        GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1367        GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1368        GenPolynomial<C> bs = bp.extendUnivariate(cfac, 1);
1369        List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1370        L.add(as);
1371        L.add(bs);
1372        List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1373
1374        Ideal<C> id = new Ideal<C>(cfac, L);
1375        //System.out.println("id = " + id);
1376        IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1377        //System.out.println("iu = " + iu);
1378
1379        // extract result polynomials
1380        List<GenPolynomial<C>> Np = iu.ideal.getList();
1381        //System.out.println("Np = " + Np);
1382        as = PolyUtil.<C> selectWithVariable(Np, 1);
1383        bs = PolyUtil.<C> selectWithVariable(Np, 0);
1384        GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1385        //System.out.println("as = " + as);
1386        //System.out.println("bs = " + bs);
1387        //System.out.println("cs = " + cs);
1388        String[] ev = new String[] { cs.ring.getVars()[0] };
1389        GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1390        //System.out.println("efac = " + efac);
1391        cs = cs.contractCoeff(efac);
1392        //System.out.println("cs = " + cs);
1393        as = as.reductum().contractCoeff(efac);
1394        as = as.negate();
1395        //System.out.println("as = " + as);
1396        bs = bs.reductum().contractCoeff(efac);
1397        bs = bs.negate();
1398        //System.out.println("bs = " + bs);
1399        AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1400        AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1401        AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1402        PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb, a, b);
1403        if (logger.isInfoEnabled()) {
1404            logger.info("primitive element = " + c);
1405        }
1406        return pe;
1407    }
1408
1409
1410    /**
1411     * Convert to primitive element ring.
1412     * @param cfac primitive element ring.
1413     * @param A algebraic number representing the generating element of a in the
1414     *            new ring.
1415     * @param a algebraic number to convert.
1416     * @return a converted to the primitive element ring
1417     */
1418    public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1419                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> a) {
1420        GenPolynomialRing<C> aufac = a.ring.ring;
1421        GenPolynomialRing<AlgebraicNumber<C>> ar = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, aufac);
1422        GenPolynomial<AlgebraicNumber<C>> aps = PolyUtil.<C> convertToAlgebraicCoefficients(ar, a.val);
1423        AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, A);
1424        return ac;
1425    }
1426
1427
1428    /**
1429     * Convert coefficients to primitive element ring.
1430     * @param cfac primitive element ring.
1431     * @param A algebraic number representing the generating element of a in the
1432     *            new ring.
1433     * @param a polynomial with coefficients algebraic number to convert.
1434     * @return a with coefficients converted to the primitive element ring
1435     */
1436    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1437                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, GenPolynomial<AlgebraicNumber<C>> a) {
1438        GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1439        return PolyUtil.<AlgebraicNumber<C>, AlgebraicNumber<C>> map(cr, a, new CoeffConvertAlg<C>(cfac, A));
1440    }
1441
1442
1443    /**
1444     * Convert to primitive element ring.
1445     * @param cfac primitive element ring.
1446     * @param A algebraic number representing the generating element of a in the
1447     *            new ring.
1448     * @param a recursive algebraic number to convert.
1449     * @return a converted to the primitive element ring
1450     */
1451    public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1452                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1453                    AlgebraicNumber<AlgebraicNumber<C>> a) {
1454        GenPolynomial<AlgebraicNumber<C>> aps = PolyUtilApp.<C> convertToPrimitiveElem(cfac, A, a.val);
1455        AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, B);
1456        return ac;
1457    }
1458
1459
1460    /**
1461     * Construct primitive element for double field extension.
1462     * @param b algebraic number ring with squarefree monic minimal polynomial
1463     *            over Q(a)
1464     * @return primitive element container with algebraic number ring c, with
1465     *         Q(c) = Q(a)(b)
1466     */
1467    public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(
1468                    AlgebraicNumberRing<AlgebraicNumber<C>> b) {
1469        GenPolynomial<AlgebraicNumber<C>> bp = b.modul;
1470        AlgebraicNumberRing<C> a = (AlgebraicNumberRing<C>) b.ring.coFac;
1471        GenPolynomial<C> ap = a.modul;
1472
1473        // setup bivariate polynomial ring
1474        String[] cv = new String[2];
1475        cv[0] = ap.ring.getVars()[0];
1476        cv[1] = bp.ring.getVars()[0];
1477        TermOrder to = new TermOrder(TermOrder.INVLEX);
1478        GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1479        GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(a.ring, 1,
1480                        bp.ring.getVars());
1481        GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1482        GenPolynomial<GenPolynomial<C>> bss = PolyUtil.<C> fromAlgebraicCoefficients(rfac, bp);
1483        GenPolynomial<C> bs = PolyUtil.<C> distribute(cfac, bss);
1484        List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1485        L.add(as);
1486        L.add(bs);
1487        List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1488
1489        Ideal<C> id = new Ideal<C>(cfac, L);
1490        //System.out.println("id = " + id);
1491        IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1492        //System.out.println("iu = " + iu);
1493
1494        // extract result polynomials
1495        List<GenPolynomial<C>> Np = iu.ideal.getList();
1496        as = PolyUtil.<C> selectWithVariable(Np, 1);
1497        bs = PolyUtil.<C> selectWithVariable(Np, 0);
1498        GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1499        //System.out.println("as = " + as);
1500        //System.out.println("bs = " + bs);
1501        //System.out.println("cs = " + cs);
1502        String[] ev = new String[] { cs.ring.getVars()[0] };
1503        GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1504        // System.out.println("efac = " + efac);
1505        cs = cs.contractCoeff(efac);
1506        // System.out.println("cs = " + cs);
1507        as = as.reductum().contractCoeff(efac);
1508        as = as.negate();
1509        // System.out.println("as = " + as);
1510        bs = bs.reductum().contractCoeff(efac);
1511        bs = bs.negate();
1512        //System.out.println("bs = " + bs);
1513        AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1514        AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1515        AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1516        PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb); // missing ,a,b);
1517        if (logger.isInfoEnabled()) {
1518            logger.info("primitive element = " + pe);
1519        }
1520        return pe;
1521    }
1522
1523
1524    /**
1525     * Convert to primitive element ring.
1526     * @param cfac primitive element ring.
1527     * @param A algebraic number representing the generating element of a in the
1528     *            new ring.
1529     * @param a polynomial with recursive algebraic number coefficients to
1530     *            convert.
1531     * @return a converted to the primitive element ring
1532     */
1533    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1534                    AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1535                    GenPolynomial<AlgebraicNumber<AlgebraicNumber<C>>> a) {
1536        GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1537        return PolyUtil.<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> map(cr, a,
1538                        new CoeffRecConvertAlg<C>(cfac, A, B));
1539    }
1540
1541
1542    /**
1543     * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1544     * RealAlgebraicNumber<C> coefficients from package
1545     * 
1546     * <pre>
1547     * edu.jas.root
1548     * </pre>
1549     * 
1550     * .
1551     * @param afac result polynomial factory.
1552     * @param A polynomial with RealAlgebraicNumber&lt;C&gt; coefficients to be
1553     *            converted.
1554     * @return polynomial with RealAlgebraicNumber&lt;C&gt; coefficients.
1555     */
1556    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> realAlgFromRealCoefficients(
1557                    GenPolynomialRing<edu.jas.root.RealAlgebraicNumber<C>> afac,
1558                    GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> A) {
1559        edu.jas.root.RealAlgebraicRing<C> cfac = (edu.jas.root.RealAlgebraicRing<C>) afac.coFac;
1560        return PolyUtil.<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> map(
1561                        afac, A, new ReAlgFromRealCoeff<C>(cfac));
1562    }
1563
1564
1565    /**
1566     * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1567     * RealAlgebraicNumber<C> coefficients from package
1568     * 
1569     * <pre>
1570     * edu.jas.application
1571     * </pre>
1572     * 
1573     * .
1574     * @param rfac result polynomial factory.
1575     * @param A polynomial with RealAlgebraicNumber&lt;C&gt; coefficients to be
1576     *            converted.
1577     * @return polynomial with RealAlgebraicNumber&lt;C&gt; coefficients.
1578     */
1579    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> realFromRealAlgCoefficients(
1580                    GenPolynomialRing<edu.jas.application.RealAlgebraicNumber<C>> rfac,
1581                    GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> A) {
1582        edu.jas.application.RealAlgebraicRing<C> cfac = (edu.jas.application.RealAlgebraicRing<C>) rfac.coFac;
1583        return PolyUtil.<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> map(
1584                        rfac, A, new RealFromReAlgCoeff<C>(cfac));
1585    }
1586
1587
1588    /**
1589     * Convert to Complex&lt;RealAlgebraicNumber&gt; coefficients. Represent as
1590     * polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients, C is
1591     * e.g. BigRational.
1592     * @param pfac result polynomial factory.
1593     * @param A polynomial with Complex coefficients to be converted.
1594     * @return polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients.
1595     */
1596    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> convertToComplexRealCoefficients(
1597                    GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac,
1598                    GenPolynomial<Complex<C>> A) {
1599        ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> afac;
1600        afac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) pfac.coFac;
1601        return PolyUtil.<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(pfac, A,
1602                        new CoeffToComplexReal<C>(afac));
1603    }
1604
1605
1606    /**
1607     * Evaluate to Complex&lt;RealAlgebraicNumber&gt; coefficients. Represent as
1608     * polynomial with Complex&lt;RealAlgebraicNumber&gt; coefficients, C is
1609     * e.g. BigRational.
1610     * @param pfac result polynomial factory.
1611     * @param A = A(x,Y) a recursive polynomial with
1612     *            GenPolynomial&lt;Complex&gt; coefficients to be converted.
1613     * @param r Complex&lt;RealAlgebraicNumber&gt; to be evaluated at.
1614     * @return A(r,Y), a polynomial with Complex&lt;RealAlgebraicNumber&gt;
1615     *         coefficients.
1616     */
1617    public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateToComplexRealCoefficients(
1618                    GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac,
1619                    GenPolynomial<GenPolynomial<Complex<C>>> A,
1620                    Complex<edu.jas.application.RealAlgebraicNumber<C>> r) {
1621        return PolyUtil.<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(
1622                        pfac, A, new EvaluateToComplexReal<C>(pfac, r));
1623    }
1624
1625
1626}
1627
1628
1629/**
1630 * Coefficient to convert algebriac functor.
1631 */
1632class CoeffConvertAlg<C extends GcdRingElem<C>> implements
1633                UnaryFunctor<AlgebraicNumber<C>, AlgebraicNumber<C>> {
1634
1635
1636    final protected AlgebraicNumberRing<C> afac;
1637
1638
1639    final protected AlgebraicNumber<C> A;
1640
1641
1642    public CoeffConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a) {
1643        if (fac == null || a == null) {
1644            throw new IllegalArgumentException("fac and a must not be null");
1645        }
1646        afac = fac;
1647        A = a;
1648    }
1649
1650
1651    public AlgebraicNumber<C> eval(AlgebraicNumber<C> c) {
1652        if (c == null) {
1653            return afac.getZERO();
1654        }
1655        return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, c);
1656    }
1657}
1658
1659
1660/**
1661 * Coefficient recursive to convert algebriac functor.
1662 */
1663class CoeffRecConvertAlg<C extends GcdRingElem<C>> implements
1664                UnaryFunctor<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> {
1665
1666
1667    final protected AlgebraicNumberRing<C> afac;
1668
1669
1670    final protected AlgebraicNumber<C> A;
1671
1672
1673    final protected AlgebraicNumber<C> B;
1674
1675
1676    public CoeffRecConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a, AlgebraicNumber<C> b) {
1677        if (fac == null || a == null || b == null) {
1678            throw new IllegalArgumentException("fac, a and b must not be null");
1679        }
1680        afac = fac;
1681        A = a;
1682        B = b;
1683    }
1684
1685
1686    public AlgebraicNumber<C> eval(AlgebraicNumber<AlgebraicNumber<C>> c) {
1687        if (c == null) {
1688            return afac.getZERO();
1689        }
1690        return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, B, c);
1691    }
1692}
1693
1694
1695/**
1696 * Coefficient to real algebriac from real algebraic functor.
1697 */
1698class ReAlgFromRealCoeff<C extends GcdRingElem<C> & Rational> implements
1699                UnaryFunctor<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> {
1700
1701
1702    final protected edu.jas.root.RealAlgebraicRing<C> afac;
1703
1704
1705    public ReAlgFromRealCoeff(edu.jas.root.RealAlgebraicRing<C> fac) {
1706        if (fac == null) {
1707            throw new IllegalArgumentException("fac must not be null");
1708        }
1709        afac = fac;
1710    }
1711
1712
1713    @SuppressWarnings("cast")
1714    public edu.jas.root.RealAlgebraicNumber<C> eval(edu.jas.application.RealAlgebraicNumber<C> c) {
1715        if (c == null) {
1716            return afac.getZERO();
1717        }
1718        return (edu.jas.root.RealAlgebraicNumber<C>) (Object) c.number; // force ignore recursion
1719    }
1720}
1721
1722
1723/**
1724 * Coefficient to real algebriac from algebraic functor.
1725 */
1726class RealFromReAlgCoeff<C extends GcdRingElem<C> & Rational> implements
1727                UnaryFunctor<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> {
1728
1729
1730    final protected edu.jas.application.RealAlgebraicRing<C> rfac;
1731
1732
1733    public RealFromReAlgCoeff(edu.jas.application.RealAlgebraicRing<C> fac) {
1734        if (fac == null) {
1735            throw new IllegalArgumentException("fac must not be null");
1736        }
1737        rfac = fac;
1738    }
1739
1740
1741    @SuppressWarnings("cast")
1742    public edu.jas.application.RealAlgebraicNumber<C> eval(edu.jas.root.RealAlgebraicNumber<C> c) {
1743        if (c == null) {
1744            return rfac.getZERO();
1745        }
1746        edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>> rrc = (edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>>) (Object) c; // force resurrect recursion
1747        return new edu.jas.application.RealAlgebraicNumber<C>(rfac, rrc);
1748    }
1749}
1750
1751
1752/**
1753 * Coefficient to complex real algebriac functor.
1754 */
1755class CoeffToComplexReal<C extends GcdRingElem<C> & Rational> implements
1756                UnaryFunctor<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> {
1757
1758
1759    final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac;
1760
1761
1762    final edu.jas.application.RealAlgebraicRing<C> afac;
1763
1764
1765    final GenPolynomialRing<C> pfac;
1766
1767
1768    public CoeffToComplexReal(ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> fac) {
1769        if (fac == null) {
1770            throw new IllegalArgumentException("fac must not be null");
1771        }
1772        cfac = fac;
1773        afac = (edu.jas.application.RealAlgebraicRing<C>) cfac.ring;
1774        pfac = afac.univs.ideal.getRing();
1775    }
1776
1777
1778    public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(Complex<C> c) {
1779        if (c == null) {
1780            return cfac.getZERO();
1781        }
1782        GenPolynomial<C> pr, pi;
1783        pr = new GenPolynomial<C>(pfac, c.getRe());
1784        pi = new GenPolynomial<C>(pfac, c.getIm());
1785        //System.out.println("pr = " + pr);
1786        //System.out.println("pi = " + pi);
1787        edu.jas.application.RealAlgebraicNumber<C> re, im;
1788        re = new edu.jas.application.RealAlgebraicNumber<C>(afac, pr);
1789        im = new edu.jas.application.RealAlgebraicNumber<C>(afac, pi);
1790        //System.out.println("re = " + re);
1791        //System.out.println("im = " + im);
1792        return new Complex<edu.jas.application.RealAlgebraicNumber<C>>(cfac, re, im);
1793    }
1794}
1795
1796
1797/**
1798 * Polynomial coefficient to complex real algebriac evaluation functor.
1799 */
1800class EvaluateToComplexReal<C extends GcdRingElem<C> & Rational> implements
1801                UnaryFunctor<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> {
1802
1803
1804    final protected GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac;
1805
1806
1807    final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac;
1808
1809
1810    final protected Complex<edu.jas.application.RealAlgebraicNumber<C>> root;
1811
1812
1813    public EvaluateToComplexReal(GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> fac,
1814                    Complex<edu.jas.application.RealAlgebraicNumber<C>> r) {
1815        if (fac == null) {
1816            throw new IllegalArgumentException("fac must not be null");
1817        }
1818        if (r == null) {
1819            throw new IllegalArgumentException("r must not be null");
1820        }
1821        pfac = fac;
1822        cfac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) fac.coFac;
1823        root = r;
1824        //System.out.println("cfac  = " + cfac);
1825        //System.out.println("root  = " + root);
1826    }
1827
1828
1829    public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(GenPolynomial<Complex<C>> c) {
1830        if (c == null) {
1831            return cfac.getZERO();
1832        }
1833        //System.out.println("c  = " + c);
1834        GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> cp;
1835        cp = PolyUtilApp.<C> convertToComplexRealCoefficients(pfac, c);
1836        Complex<edu.jas.application.RealAlgebraicNumber<C>> cr;
1837        cr = PolyUtil.<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateMain(cfac, cp, root);
1838        return cr;
1839    }
1840}