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