001/*
002 * $Id: PolyUtil.java 6010 2020-04-01 10:39:15Z kredel $
003 */
004
005package edu.jas.poly;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.Map;
011import java.util.SortedMap;
012import java.util.TreeMap;
013
014import org.apache.logging.log4j.LogManager;
015import org.apache.logging.log4j.Logger;
016
017import edu.jas.arith.BigComplex;
018import edu.jas.arith.BigDecimal;
019import edu.jas.arith.BigDecimalComplex;
020import edu.jas.arith.BigInteger;
021import edu.jas.arith.BigRational;
022import edu.jas.arith.ModInteger;
023import edu.jas.arith.ModIntegerRing;
024import edu.jas.arith.Modular;
025import edu.jas.arith.ModularRingFactory;
026import edu.jas.arith.Product;
027import edu.jas.arith.ProductRing;
028import edu.jas.arith.Rational;
029import edu.jas.arith.Roots;
030import edu.jas.structure.Element;
031import edu.jas.structure.GcdRingElem;
032import edu.jas.structure.RingElem;
033import edu.jas.structure.RingFactory;
034import edu.jas.structure.StarRingElem;
035import edu.jas.structure.UnaryFunctor;
036import edu.jas.util.ListUtil;
037
038
039/**
040 * Polynomial utilities, for example conversion between different
041 * representations, evaluation and interpolation.
042 * @author Heinz Kredel
043 */
044
045public class PolyUtil {
046
047
048    private static final Logger logger = LogManager.getLogger(PolyUtil.class);
049
050
051    private static final boolean debug = logger.isDebugEnabled();
052
053
054    /**
055     * Recursive representation. Represent as polynomial in i variables with
056     * coefficients in n-i variables. Works for arbitrary term orders.
057     * @param <C> coefficient type.
058     * @param rfac recursive polynomial ring factory.
059     * @param A polynomial to be converted.
060     * @return Recursive represenations of this in the ring rfac.
061     */
062    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursive(
063                    GenPolynomialRing<GenPolynomial<C>> rfac, GenPolynomial<C> A) {
064
065        GenPolynomial<GenPolynomial<C>> B = rfac.getZERO().copy();
066        if (A.isZERO()) {
067            return B;
068        }
069        int i = rfac.nvar;
070        GenPolynomial<C> zero = rfac.getZEROCoefficient();
071        Map<ExpVector, GenPolynomial<C>> Bv = B.val; //getMap();
072        for (Map.Entry<ExpVector, C> y : A.getMap().entrySet()) {
073            ExpVector e = y.getKey();
074            C a = y.getValue();
075            ExpVector f = e.contract(0, i);
076            ExpVector g = e.contract(i, e.length() - i);
077            GenPolynomial<C> p = Bv.get(f);
078            if (p == null) {
079                p = zero;
080            }
081            p = p.sum(a, g);
082            Bv.put(f, p);
083        }
084        return B;
085    }
086
087
088    /**
089     * Distribute a recursive polynomial to a generic polynomial. Works for
090     * arbitrary term orders.
091     * @param <C> coefficient type.
092     * @param dfac combined polynomial ring factory of coefficients and this.
093     * @param B polynomial to be converted.
094     * @return distributed polynomial.
095     */
096    public static <C extends RingElem<C>> GenPolynomial<C> distribute(GenPolynomialRing<C> dfac,
097                    GenPolynomial<GenPolynomial<C>> B) {
098        GenPolynomial<C> C = dfac.getZERO().copy();
099        if (B.isZERO()) {
100            return C;
101        }
102        Map<ExpVector, C> Cm = C.val; //getMap();
103        for (Map.Entry<ExpVector, GenPolynomial<C>> y : B.getMap().entrySet()) {
104            ExpVector e = y.getKey();
105            GenPolynomial<C> A = y.getValue();
106            for (Map.Entry<ExpVector, C> x : A.val.entrySet()) {
107                ExpVector f = x.getKey();
108                C b = x.getValue();
109                ExpVector g = e.combine(f);
110                assert (Cm.get(g) == null);
111                Cm.put(g, b);
112            }
113        }
114        return C;
115    }
116
117
118    /**
119     * Recursive representation. Represent as polynomials in i variables with
120     * coefficients in n-i variables. Works for arbitrary term orders.
121     * @param <C> coefficient type.
122     * @param rfac recursive polynomial ring factory.
123     * @param L list of polynomials to be converted.
124     * @return Recursive represenations of the list in the ring rfac.
125     */
126    public static <C extends RingElem<C>> List<GenPolynomial<GenPolynomial<C>>> recursive(
127                    GenPolynomialRing<GenPolynomial<C>> rfac, List<GenPolynomial<C>> L) {
128        return ListUtil.<GenPolynomial<C>, GenPolynomial<GenPolynomial<C>>> map(L, new DistToRec<C>(rfac));
129    }
130
131
132    /**
133     * Distribute a recursive polynomial list to a generic polynomial list.
134     * Works for arbitrary term orders.
135     * @param <C> coefficient type.
136     * @param dfac combined polynomial ring factory of coefficients and this.
137     * @param L list of polynomials to be converted.
138     * @return distributed polynomial list.
139     */
140    public static <C extends RingElem<C>> List<GenPolynomial<C>> distribute(GenPolynomialRing<C> dfac,
141                    List<GenPolynomial<GenPolynomial<C>>> L) {
142        return ListUtil.<GenPolynomial<GenPolynomial<C>>, GenPolynomial<C>> map(L, new RecToDist<C>(dfac));
143    }
144
145
146    /**
147     * BigInteger from ModInteger coefficients, symmetric. Represent as
148     * polynomial with BigInteger coefficients by removing the modules and
149     * making coefficients symmetric to 0.
150     * @param fac result polynomial factory.
151     * @param A polynomial with ModInteger coefficients to be converted.
152     * @return polynomial with BigInteger coefficients.
153     */
154    public static <C extends RingElem<C> & Modular> GenPolynomial<BigInteger> integerFromModularCoefficients(
155                    GenPolynomialRing<BigInteger> fac, GenPolynomial<C> A) {
156        return PolyUtil.<C, BigInteger> map(fac, A, new ModSymToInt<C>());
157    }
158
159
160    /**
161     * BigInteger from ModInteger coefficients, symmetric. Represent as
162     * polynomial with BigInteger coefficients by removing the modules and
163     * making coefficients symmetric to 0.
164     * @param fac result polynomial factory.
165     * @param L list of polynomials with ModInteger coefficients to be
166     *            converted.
167     * @return list of polynomials with BigInteger coefficients.
168     */
169    public static <C extends RingElem<C> & Modular> List<GenPolynomial<BigInteger>> integerFromModularCoefficients(
170                    final GenPolynomialRing<BigInteger> fac, List<GenPolynomial<C>> L) {
171        return ListUtil.<GenPolynomial<C>, GenPolynomial<BigInteger>> map(L,
172                        new UnaryFunctor<GenPolynomial<C>, GenPolynomial<BigInteger>>() {
173
174
175                            public GenPolynomial<BigInteger> eval(GenPolynomial<C> c) {
176                                return PolyUtil.<C> integerFromModularCoefficients(fac, c);
177                            }
178                        });
179    }
180
181
182    /**
183     * BigInteger from ModInteger coefficients, positive. Represent as
184     * polynomial with BigInteger coefficients by removing the modules.
185     * @param fac result polynomial factory.
186     * @param A polynomial with ModInteger coefficients to be converted.
187     * @return polynomial with BigInteger coefficients.
188     */
189    public static <C extends RingElem<C> & Modular> GenPolynomial<BigInteger> integerFromModularCoefficientsPositive(
190                    GenPolynomialRing<BigInteger> fac, GenPolynomial<C> A) {
191        return PolyUtil.<C, BigInteger> map(fac, A, new ModToInt<C>());
192    }
193
194
195    /**
196     * BigInteger from BigRational coefficients. Represent as polynomial with
197     * BigInteger coefficients by multiplication with the lcm of the numerators
198     * of the BigRational coefficients.
199     * @param fac result polynomial factory.
200     * @param A polynomial with BigRational coefficients to be converted.
201     * @return polynomial with BigInteger coefficients.
202     */
203    public static GenPolynomial<BigInteger> integerFromRationalCoefficients(GenPolynomialRing<BigInteger> fac,
204                    GenPolynomial<BigRational> A) {
205        if (A == null || A.isZERO()) {
206            return fac.getZERO();
207        }
208        java.math.BigInteger c = null;
209        int s = 0;
210        // lcm of denominators
211        for (BigRational y : A.val.values()) {
212            java.math.BigInteger x = y.denominator();
213            // c = lcm(c,x)
214            if (c == null) {
215                c = x;
216                s = x.signum();
217            } else {
218                java.math.BigInteger d = c.gcd(x);
219                c = c.multiply(x.divide(d));
220            }
221        }
222        if (s < 0) {
223            c = c.negate();
224        }
225        return PolyUtil.<BigRational, BigInteger> map(fac, A, new RatToInt(c));
226    }
227
228
229    /**
230     * BigInteger from BigRational coefficients. Represent as polynomial with
231     * BigInteger coefficients by multiplication with the gcd of the numerators
232     * and the lcm of the denominators of the BigRational coefficients. <br />
233     * <b>Author:</b> Axel Kramer
234     * @param fac result polynomial factory.
235     * @param A polynomial with BigRational coefficients to be converted.
236     * @return Object[] with 3 entries: [0]->gcd [1]->lcm and [2]->polynomial
237     *         with BigInteger coefficients.
238     */
239    public static Object[] integerFromRationalCoefficientsFactor(GenPolynomialRing<BigInteger> fac,
240                    GenPolynomial<BigRational> A) {
241        Object[] result = new Object[3];
242        if (A == null || A.isZERO()) {
243            result[0] = java.math.BigInteger.ONE;
244            result[1] = java.math.BigInteger.ZERO;
245            result[2] = fac.getZERO();
246            return result;
247        }
248        java.math.BigInteger gcd = null;
249        java.math.BigInteger lcm = null;
250        int sLCM = 0;
251        int sGCD = 0;
252        // lcm of denominators
253        for (BigRational y : A.val.values()) {
254            java.math.BigInteger numerator = y.numerator();
255            java.math.BigInteger denominator = y.denominator();
256            // lcm = lcm(lcm,x)
257            if (lcm == null) {
258                lcm = denominator;
259                sLCM = denominator.signum();
260            } else {
261                java.math.BigInteger d = lcm.gcd(denominator);
262                lcm = lcm.multiply(denominator.divide(d));
263            }
264            // gcd = gcd(gcd,x)
265            if (gcd == null) {
266                gcd = numerator;
267                sGCD = numerator.signum();
268            } else {
269                gcd = gcd.gcd(numerator);
270            }
271        }
272        if (sLCM < 0) {
273            lcm = lcm.negate();
274        }
275        if (sGCD < 0) {
276            gcd = gcd.negate();
277        }
278        //System.out.println("gcd* = " + gcd + ", lcm = " + lcm);
279        result[0] = gcd;
280        result[1] = lcm;
281        result[2] = PolyUtil.<BigRational, BigInteger> map(fac, A, new RatToIntFactor(gcd, lcm));
282        return result;
283    }
284
285
286    /**
287     * BigInteger from BigRational coefficients. Represent as polynomial with
288     * BigInteger coefficients by multiplication with the gcd of the numerators
289     * and the lcm of the denominators of the BigRational coefficients. <br />
290     * @param fac result polynomial factory.
291     * @param gcd of rational coefficient numerators.
292     * @param lcm of rational coefficient denominators.
293     * @param A polynomial with BigRational coefficients to be converted.
294     * @return polynomial with BigInteger coefficients.
295     */
296    public static GenPolynomial<BigInteger> integerFromRationalCoefficients(GenPolynomialRing<BigInteger> fac,
297                    java.math.BigInteger gcd, java.math.BigInteger lcm, GenPolynomial<BigRational> A) {
298        //System.out.println("gcd = " + gcd + ", lcm = " + lcm);
299        GenPolynomial<BigInteger> Ai = PolyUtil.<BigRational, BigInteger> map(fac, A,
300                        new RatToIntFactor(gcd, lcm));
301        return Ai;
302    }
303
304
305    /**
306     * BigInteger from BigRational coefficients. Represent as list of
307     * polynomials with BigInteger coefficients by multiplication with the lcm
308     * of the numerators of the BigRational coefficients of each polynomial.
309     * @param fac result polynomial factory.
310     * @param L list of polynomials with BigRational coefficients to be
311     *            converted.
312     * @return polynomial list with BigInteger coefficients.
313     */
314    public static List<GenPolynomial<BigInteger>> integerFromRationalCoefficients(
315                    GenPolynomialRing<BigInteger> fac, List<GenPolynomial<BigRational>> L) {
316        return ListUtil.<GenPolynomial<BigRational>, GenPolynomial<BigInteger>> map(L, new RatToIntPoly(fac));
317    }
318
319
320    /**
321     * From BigInteger coefficients. Represent as polynomial with type C
322     * coefficients, e.g. ModInteger or BigRational.
323     * @param <C> coefficient type.
324     * @param fac result polynomial factory.
325     * @param A polynomial with BigInteger coefficients to be converted.
326     * @return polynomial with type C coefficients.
327     */
328    public static <C extends RingElem<C>> GenPolynomial<C> fromIntegerCoefficients(GenPolynomialRing<C> fac,
329                    GenPolynomial<BigInteger> A) {
330        return PolyUtil.<BigInteger, C> map(fac, A, new FromInteger<C>(fac.coFac));
331    }
332
333
334    /**
335     * From BigInteger coefficients. Represent as list of polynomials with type
336     * C coefficients, e.g. ModInteger or BigRational.
337     * @param <C> coefficient type.
338     * @param fac result polynomial factory.
339     * @param L list of polynomials with BigInteger coefficients to be
340     *            converted.
341     * @return list of polynomials with type C coefficients.
342     */
343    public static <C extends RingElem<C>> List<GenPolynomial<C>> fromIntegerCoefficients(
344                    GenPolynomialRing<C> fac, List<GenPolynomial<BigInteger>> L) {
345        return ListUtil.<GenPolynomial<BigInteger>, GenPolynomial<C>> map(L, new FromIntegerPoly<C>(fac));
346    }
347
348
349    /**
350     * Convert to decimal coefficients.
351     * @param fac result polynomial factory.
352     * @param A polynomial with Rational coefficients to be converted.
353     * @return polynomial with BigDecimal coefficients.
354     */
355    public static <C extends RingElem<C> & Rational> GenPolynomial<BigDecimal> decimalFromRational(
356                    GenPolynomialRing<BigDecimal> fac, GenPolynomial<C> A) {
357        return PolyUtil.<C, BigDecimal> map(fac, A, new RatToDec<C>());
358    }
359
360
361    /**
362     * Convert to complex decimal coefficients.
363     * @param fac result polynomial factory.
364     * @param A polynomial with complex Rational coefficients to be converted.
365     * @return polynomial with Complex BigDecimal coefficients.
366     */
367    public static <C extends RingElem<C> & Rational> GenPolynomial<Complex<BigDecimal>> complexDecimalFromRational(
368                    GenPolynomialRing<Complex<BigDecimal>> fac, GenPolynomial<Complex<C>> A) {
369        return PolyUtil.<Complex<C>, Complex<BigDecimal>> map(fac, A, new CompRatToDec<C>(fac.coFac));
370    }
371
372
373    /**
374     * Real part.
375     * @param fac result polynomial factory.
376     * @param A polynomial with BigComplex coefficients to be converted.
377     * @return polynomial with real part of the coefficients.
378     */
379    public static GenPolynomial<BigRational> realPart(GenPolynomialRing<BigRational> fac,
380                    GenPolynomial<BigComplex> A) {
381        return PolyUtil.<BigComplex, BigRational> map(fac, A, new RealPart());
382    }
383
384
385    /**
386     * Imaginary part.
387     * @param fac result polynomial factory.
388     * @param A polynomial with BigComplex coefficients to be converted.
389     * @return polynomial with imaginary part of coefficients.
390     */
391    public static GenPolynomial<BigRational> imaginaryPart(GenPolynomialRing<BigRational> fac,
392                    GenPolynomial<BigComplex> A) {
393        return PolyUtil.<BigComplex, BigRational> map(fac, A, new ImagPart());
394    }
395
396
397    /**
398     * Real part.
399     * @param fac result polynomial factory.
400     * @param A polynomial with BigComplex coefficients to be converted.
401     * @return polynomial with real part of the coefficients.
402     */
403    public static <C extends RingElem<C>> GenPolynomial<C> realPartFromComplex(GenPolynomialRing<C> fac,
404                    GenPolynomial<Complex<C>> A) {
405        return PolyUtil.<Complex<C>, C> map(fac, A, new RealPartComplex<C>());
406    }
407
408
409    /**
410     * Imaginary part.
411     * @param fac result polynomial factory.
412     * @param A polynomial with BigComplex coefficients to be converted.
413     * @return polynomial with imaginary part of coefficients.
414     */
415    public static <C extends RingElem<C>> GenPolynomial<C> imaginaryPartFromComplex(GenPolynomialRing<C> fac,
416                    GenPolynomial<Complex<C>> A) {
417        return PolyUtil.<Complex<C>, C> map(fac, A, new ImagPartComplex<C>());
418    }
419
420
421    /**
422     * Complex from real polynomial.
423     * @param fac result polynomial factory.
424     * @param A polynomial with C coefficients to be converted.
425     * @return polynomial with Complex<C> coefficients.
426     */
427    public static <C extends RingElem<C>> GenPolynomial<Complex<C>> toComplex(
428                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<C> A) {
429        return PolyUtil.<C, Complex<C>> map(fac, A, new ToComplex<C>(fac.coFac));
430    }
431
432
433    /**
434     * Complex from rational coefficients.
435     * @param fac result polynomial factory.
436     * @param A polynomial with BigRational coefficients to be converted.
437     * @return polynomial with BigComplex coefficients.
438     */
439    public static GenPolynomial<BigComplex> complexFromRational(GenPolynomialRing<BigComplex> fac,
440                    GenPolynomial<BigRational> A) {
441        return PolyUtil.<BigRational, BigComplex> map(fac, A, new RatToCompl());
442    }
443
444
445    /**
446     * Complex from ring element coefficients.
447     * @param fac result polynomial factory.
448     * @param A polynomial with RingElem coefficients to be converted.
449     * @return polynomial with Complex coefficients.
450     */
451    public static <C extends GcdRingElem<C>> GenPolynomial<Complex<C>> complexFromAny(
452                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<C> A) {
453        ComplexRing<C> cr = (ComplexRing<C>) fac.coFac;
454        return PolyUtil.<C, Complex<C>> map(fac, A, new AnyToComplex<C>(cr));
455    }
456
457
458    /**
459     * From AlgebraicNumber coefficients. Represent as polynomial with type
460     * GenPolynomial&lt;C&gt; coefficients, e.g. ModInteger or BigRational.
461     * @param rfac result polynomial factory.
462     * @param A polynomial with AlgebraicNumber coefficients to be converted.
463     * @return polynomial with type GenPolynomial&lt;C&gt; coefficients.
464     */
465    public static <C extends GcdRingElem<C>> GenPolynomial<GenPolynomial<C>> fromAlgebraicCoefficients(
466                    GenPolynomialRing<GenPolynomial<C>> rfac, GenPolynomial<AlgebraicNumber<C>> A) {
467        return PolyUtil.<AlgebraicNumber<C>, GenPolynomial<C>> map(rfac, A, new AlgToPoly<C>());
468    }
469
470
471    /**
472     * Convert to AlgebraicNumber coefficients. Represent as polynomial with
473     * AlgebraicNumber<C> coefficients, C is e.g. ModInteger or BigRational.
474     * @param pfac result polynomial factory.
475     * @param A polynomial with C coefficients to be converted.
476     * @return polynomial with AlgebraicNumber&lt;C&gt; coefficients.
477     */
478    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToAlgebraicCoefficients(
479                    GenPolynomialRing<AlgebraicNumber<C>> pfac, GenPolynomial<C> A) {
480        AlgebraicNumberRing<C> afac = (AlgebraicNumberRing<C>) pfac.coFac;
481        return PolyUtil.<C, AlgebraicNumber<C>> map(pfac, A, new CoeffToAlg<C>(afac));
482    }
483
484
485    /**
486     * Convert to recursive AlgebraicNumber coefficients. Represent as
487     * polynomial with recursive AlgebraicNumber<C> coefficients, C is e.g.
488     * ModInteger or BigRational.
489     * @param depth recursion depth of AlgebraicNumber coefficients.
490     * @param pfac result polynomial factory.
491     * @param A polynomial with C coefficients to be converted.
492     * @return polynomial with AlgebraicNumber&lt;C&gt; coefficients.
493     */
494    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToRecAlgebraicCoefficients(
495                    int depth, GenPolynomialRing<AlgebraicNumber<C>> pfac, GenPolynomial<C> A) {
496        AlgebraicNumberRing<C> afac = (AlgebraicNumberRing<C>) pfac.coFac;
497        return PolyUtil.<C, AlgebraicNumber<C>> map(pfac, A, new CoeffToRecAlg<C>(depth, afac));
498    }
499
500
501    /**
502     * Convert to AlgebraicNumber coefficients. Represent as polynomial with
503     * AlgebraicNumber<C> coefficients, C is e.g. ModInteger or BigRational.
504     * @param pfac result polynomial factory.
505     * @param A recursive polynomial with GenPolynomial&lt;BigInteger&gt;
506     *            coefficients to be converted.
507     * @return polynomial with AlgebraicNumber&lt;C&gt; coefficients.
508     */
509    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertRecursiveToAlgebraicCoefficients(
510                    GenPolynomialRing<AlgebraicNumber<C>> pfac, GenPolynomial<GenPolynomial<C>> A) {
511        AlgebraicNumberRing<C> afac = (AlgebraicNumberRing<C>) pfac.coFac;
512        return PolyUtil.<GenPolynomial<C>, AlgebraicNumber<C>> map(pfac, A, new PolyToAlg<C>(afac));
513    }
514
515
516    /**
517     * Complex from algebraic coefficients.
518     * @param fac result polynomial factory.
519     * @param A polynomial with AlgebraicNumber coefficients Q(i) to be
520     *            converted.
521     * @return polynomial with Complex coefficients.
522     */
523    public static <C extends GcdRingElem<C>> GenPolynomial<Complex<C>> complexFromAlgebraic(
524                    GenPolynomialRing<Complex<C>> fac, GenPolynomial<AlgebraicNumber<C>> A) {
525        ComplexRing<C> cfac = (ComplexRing<C>) fac.coFac;
526        return PolyUtil.<AlgebraicNumber<C>, Complex<C>> map(fac, A, new AlgebToCompl<C>(cfac));
527    }
528
529
530    /**
531     * AlgebraicNumber from complex coefficients.
532     * @param fac result polynomial factory over Q(i).
533     * @param A polynomial with Complex coefficients to be converted.
534     * @return polynomial with AlgebraicNumber coefficients.
535     */
536    public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> algebraicFromComplex(
537                    GenPolynomialRing<AlgebraicNumber<C>> fac, GenPolynomial<Complex<C>> A) {
538        AlgebraicNumberRing<C> afac = (AlgebraicNumberRing<C>) fac.coFac;
539        return PolyUtil.<Complex<C>, AlgebraicNumber<C>> map(fac, A, new ComplToAlgeb<C>(afac));
540    }
541
542
543    /**
544     * ModInteger chinese remainder algorithm on coefficients.
545     * @param fac GenPolynomial&lt;ModInteger&gt; result factory with
546     *            A.coFac.modul*B.coFac.modul = C.coFac.modul.
547     * @param A GenPolynomial&lt;ModInteger&gt;.
548     * @param B other GenPolynomial&lt;ModInteger&gt;.
549     * @param mi inverse of A.coFac.modul in ring B.coFac.
550     * @return S = cra(A,B), with S mod A.coFac.modul == A and S mod
551     *         B.coFac.modul == B.
552     */
553    public static <C extends RingElem<C> & Modular> GenPolynomial<C> chineseRemainder(
554                    GenPolynomialRing<C> fac, GenPolynomial<C> A, C mi, GenPolynomial<C> B) {
555        ModularRingFactory<C> cfac = (ModularRingFactory<C>) fac.coFac; // get RingFactory
556        GenPolynomial<C> S = fac.getZERO().copy();
557        GenPolynomial<C> Ap = A.copy();
558        SortedMap<ExpVector, C> av = Ap.val; //getMap();
559        SortedMap<ExpVector, C> bv = B.getMap();
560        SortedMap<ExpVector, C> sv = S.val; //getMap();
561        C c = null;
562        for (Map.Entry<ExpVector, C> me : bv.entrySet()) {
563            ExpVector e = me.getKey();
564            C y = me.getValue(); //bv.get(e); // assert y != null
565            C x = av.get(e);
566            if (x != null) {
567                av.remove(e);
568                c = cfac.chineseRemainder(x, mi, y);
569                if (!c.isZERO()) { // 0 cannot happen
570                    sv.put(e, c);
571                }
572            } else {
573                //c = cfac.fromInteger( y.getVal() );
574                c = cfac.chineseRemainder(A.ring.coFac.getZERO(), mi, y);
575                if (!c.isZERO()) { // 0 cannot happen
576                    sv.put(e, c); // c != null
577                }
578            }
579        }
580        // assert bv is empty = done
581        for (Map.Entry<ExpVector, C> me : av.entrySet()) { // rest of av
582            ExpVector e = me.getKey();
583            C x = me.getValue(); // av.get(e); // assert x != null
584            //c = cfac.fromInteger( x.getVal() );
585            c = cfac.chineseRemainder(x, mi, B.ring.coFac.getZERO());
586            if (!c.isZERO()) { // 0 cannot happen
587                sv.put(e, c); // c != null
588            }
589        }
590        return S;
591    }
592
593
594    /**
595     * GenPolynomial monic, i.e. leadingBaseCoefficient == 1. If
596     * leadingBaseCoefficient is not invertible returns this unmodified.
597     * @param <C> coefficient type.
598     * @param p recursive GenPolynomial<GenPolynomial<C>>.
599     * @return monic(p).
600     */
601    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> monic(
602                    GenPolynomial<GenPolynomial<C>> p) {
603        if (p == null || p.isZERO()) {
604            return p;
605        }
606        C lc = p.leadingBaseCoefficient().leadingBaseCoefficient();
607        if (!lc.isUnit()) {
608            return p;
609        }
610        C lm = lc.inverse();
611        GenPolynomial<C> L = p.ring.coFac.getONE();
612        L = L.multiply(lm);
613        return p.multiplyLeft(L);
614    }
615
616
617    /**
618     * GenSolvablePolynomial monic, i.e. leadingBaseCoefficient == 1. If
619     * leadingBaseCoefficient is not invertible returns this unmodified.
620     * @param <C> coefficient type.
621     * @param p recursive GenSolvablePolynomial<GenPolynomial<C>>.
622     * @return monic(p).
623     */
624    public static <C extends RingElem<C>> GenSolvablePolynomial<GenPolynomial<C>> monic(
625                    GenSolvablePolynomial<GenPolynomial<C>> p) {
626        if (p == null || p.isZERO()) {
627            return p;
628        }
629        C lc = p.leadingBaseCoefficient().leadingBaseCoefficient();
630        if (!lc.isUnit()) {
631            return p;
632        }
633        C lm = lc.inverse();
634        GenPolynomial<C> L = p.ring.coFac.getONE();
635        L = L.multiply(lm);
636        return p.multiplyLeft(L);
637    }
638
639
640    /**
641     * Polynomial list monic.
642     * @param <C> coefficient type.
643     * @param L list of polynomials with field coefficients.
644     * @return list of polynomials with leading coefficient 1.
645     */
646    public static <C extends RingElem<C>> List<GenPolynomial<C>> monic(List<GenPolynomial<C>> L) {
647        return ListUtil.<GenPolynomial<C>, GenPolynomial<C>> map(L,
648                        new UnaryFunctor<GenPolynomial<C>, GenPolynomial<C>>() {
649
650
651                            public GenPolynomial<C> eval(GenPolynomial<C> c) {
652                                if (c == null) {
653                                    return null;
654                                }
655                                return c.monic();
656                            }
657                        });
658    }
659
660
661    /**
662     * Word polynomial list monic.
663     * @param <C> coefficient type.
664     * @param L list of word polynomials with field coefficients.
665     * @return list of word polynomials with leading coefficient 1.
666     */
667    public static <C extends RingElem<C>> List<GenWordPolynomial<C>> wordMonic(List<GenWordPolynomial<C>> L) {
668        return ListUtil.<GenWordPolynomial<C>, GenWordPolynomial<C>> map(L,
669                        new UnaryFunctor<GenWordPolynomial<C>, GenWordPolynomial<C>>() {
670
671
672                            public GenWordPolynomial<C> eval(GenWordPolynomial<C> c) {
673                                if (c == null) {
674                                    return null;
675                                }
676                                return c.monic();
677                            }
678                        });
679    }
680
681
682    /**
683     * Recursive polynomial list monic.
684     * @param <C> coefficient type.
685     * @param L list of recursive polynomials with field coefficients.
686     * @return list of polynomials with leading base coefficient 1.
687     */
688    public static <C extends RingElem<C>> List<GenPolynomial<GenPolynomial<C>>> monicRec(
689                    List<GenPolynomial<GenPolynomial<C>>> L) {
690        return ListUtil.<GenPolynomial<GenPolynomial<C>>, GenPolynomial<GenPolynomial<C>>> map(L,
691                        new UnaryFunctor<GenPolynomial<GenPolynomial<C>>, GenPolynomial<GenPolynomial<C>>>() {
692
693
694                            public GenPolynomial<GenPolynomial<C>> eval(GenPolynomial<GenPolynomial<C>> c) {
695                                if (c == null) {
696                                    return null;
697                                }
698                                return PolyUtil.<C> monic(c);
699                            }
700                        });
701    }
702
703
704    /**
705     * Polynomial list leading exponent vectors.
706     * @param <C> coefficient type.
707     * @param L list of polynomials.
708     * @return list of leading exponent vectors.
709     */
710    public static <C extends RingElem<C>> List<ExpVector> leadingExpVector(List<GenPolynomial<C>> L) {
711        return ListUtil.<GenPolynomial<C>, ExpVector> map(L, new UnaryFunctor<GenPolynomial<C>, ExpVector>() {
712
713
714            public ExpVector eval(GenPolynomial<C> c) {
715                if (c == null) {
716                    return null;
717                }
718                return c.leadingExpVector();
719            }
720        });
721    }
722
723
724    /**
725     * Extend coefficient variables. Extend all coefficient ExpVectors by i
726     * elements and multiply by x_j^k.
727     * @param pfac extended polynomial ring factory (by i variables in the
728     *            coefficients).
729     * @param j index of variable to be used for multiplication.
730     * @param k exponent for x_j.
731     * @return extended polynomial.
732     */
733    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> extendCoefficients(
734                    GenPolynomialRing<GenPolynomial<C>> pfac, GenPolynomial<GenPolynomial<C>> A, int j,
735                    long k) {
736        GenPolynomial<GenPolynomial<C>> Cp = pfac.getZERO().copy();
737        if (A.isZERO()) {
738            return Cp;
739        }
740        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) pfac.coFac;
741        //GenPolynomialRing<C> acfac = (GenPolynomialRing<C>) A.ring.coFac;
742        //int i = cfac.nvar - acfac.nvar;
743        Map<ExpVector, GenPolynomial<C>> CC = Cp.val; //getMap();
744        for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.val.entrySet()) {
745            ExpVector e = y.getKey();
746            GenPolynomial<C> a = y.getValue();
747            GenPolynomial<C> f = a.extend(cfac, j, k);
748            CC.put(e, f);
749        }
750        return Cp;
751    }
752
753
754    /**
755     * Extend coefficient variables. Extend all coefficient ExpVectors by i
756     * elements and multiply by x_j^k.
757     * @param pfac extended polynomial ring factory (by i variables in the
758     *            coefficients).
759     * @param j index of variable to be used for multiplication.
760     * @param k exponent for x_j.
761     * @return extended polynomial.
762     */
763    public static <C extends RingElem<C>> GenSolvablePolynomial<GenPolynomial<C>> extendCoefficients(
764                    GenSolvablePolynomialRing<GenPolynomial<C>> pfac,
765                    GenSolvablePolynomial<GenPolynomial<C>> A, int j, long k) {
766        GenSolvablePolynomial<GenPolynomial<C>> Cp = pfac.getZERO().copy();
767        if (A.isZERO()) {
768            return Cp;
769        }
770        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) pfac.coFac;
771        //GenPolynomialRing<C> acfac = (GenPolynomialRing<C>) A.ring.coFac;
772        //int i = cfac.nvar - acfac.nvar;
773        Map<ExpVector, GenPolynomial<C>> CC = Cp.val; //getMap();
774        for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.val.entrySet()) {
775            ExpVector e = y.getKey();
776            GenPolynomial<C> a = y.getValue();
777            GenPolynomial<C> f = a.extend(cfac, j, k);
778            CC.put(e, f);
779        }
780        return Cp;
781    }
782
783
784    /**
785     * To recursive representation. Represent as polynomial in i+r variables
786     * with coefficients in i variables. Works for arbitrary term orders.
787     * @param <C> coefficient type.
788     * @param rfac recursive polynomial ring factory.
789     * @param A polynomial to be converted.
790     * @return Recursive represenations of A in the ring rfac.
791     */
792    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> toRecursive(
793                    GenPolynomialRing<GenPolynomial<C>> rfac, GenPolynomial<C> A) {
794
795        GenPolynomial<GenPolynomial<C>> B = rfac.getZERO().copy();
796        if (A.isZERO()) {
797            return B;
798        }
799        //int i = rfac.nvar;
800        //GenPolynomial<C> zero = rfac.getZEROCoefficient();
801        GenPolynomial<C> one = rfac.getONECoefficient();
802        Map<ExpVector, GenPolynomial<C>> Bv = B.val; //getMap();
803        for (Monomial<C> m : A) {
804            ExpVector e = m.e;
805            C a = m.c;
806            GenPolynomial<C> p = one.multiply(a);
807            Bv.put(e, p);
808        }
809        return B;
810    }
811
812
813    /**
814     * To recursive representation. Represent as solvable polynomial in i+r
815     * variables with coefficients in i variables. Works for arbitrary term
816     * orders.
817     * @param <C> coefficient type.
818     * @param rfac recursive solvable polynomial ring factory.
819     * @param A solvable polynomial to be converted.
820     * @return Recursive represenations of A in the ring rfac.
821     */
822    public static <C extends RingElem<C>> GenSolvablePolynomial<GenPolynomial<C>> toRecursive(
823                    GenSolvablePolynomialRing<GenPolynomial<C>> rfac, GenSolvablePolynomial<C> A) {
824
825        GenSolvablePolynomial<GenPolynomial<C>> B = rfac.getZERO().copy();
826        if (A.isZERO()) {
827            return B;
828        }
829        //int i = rfac.nvar;
830        //GenPolynomial<C> zero = rfac.getZEROCoefficient();
831        GenPolynomial<C> one = rfac.getONECoefficient();
832        Map<ExpVector, GenPolynomial<C>> Bv = B.val; //getMap();
833        for (Monomial<C> m : A) {
834            ExpVector e = m.e;
835            C a = m.c;
836            GenPolynomial<C> p = one.multiply(a);
837            Bv.put(e, p);
838        }
839        return B;
840    }
841
842
843    /**
844     * GenPolynomial coefficient wise remainder.
845     * @param <C> coefficient type.
846     * @param P GenPolynomial.
847     * @param s nonzero coefficient.
848     * @return coefficient wise remainder.
849     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
850     */
851    public static <C extends RingElem<C>> GenPolynomial<C> baseRemainderPoly(GenPolynomial<C> P, C s) {
852        if (s == null || s.isZERO()) {
853            throw new ArithmeticException(P + " division by zero " + s);
854        }
855        GenPolynomial<C> h = P.ring.getZERO().copy();
856        Map<ExpVector, C> hm = h.val; //getMap();
857        for (Map.Entry<ExpVector, C> m : P.getMap().entrySet()) {
858            ExpVector f = m.getKey();
859            C a = m.getValue();
860            C x = a.remainder(s);
861            hm.put(f, x);
862        }
863        return h;
864    }
865
866
867    /**
868     * GenPolynomial sparse pseudo remainder. For univariate polynomials.
869     * @param <C> coefficient type.
870     * @param P GenPolynomial.
871     * @param S nonzero GenPolynomial.
872     * @return remainder with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
873     *         m' &le; deg(P)-deg(S)
874     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
875     * @deprecated Use
876     *             {@link #baseSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)}
877     *             instead
878     */
879    @Deprecated
880    public static <C extends RingElem<C>> GenPolynomial<C> basePseudoRemainder(GenPolynomial<C> P,
881                    GenPolynomial<C> S) {
882        return baseSparsePseudoRemainder(P, S);
883    }
884
885
886    /**
887     * GenPolynomial sparse pseudo remainder. For univariate polynomials.
888     * @param <C> coefficient type.
889     * @param P GenPolynomial.
890     * @param S nonzero GenPolynomial.
891     * @return remainder with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
892     *         m' &le; deg(P)-deg(S)
893     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
894     */
895    public static <C extends RingElem<C>> GenPolynomial<C> baseSparsePseudoRemainder(GenPolynomial<C> P,
896                    GenPolynomial<C> S) {
897        if (S == null || S.isZERO()) {
898            throw new ArithmeticException(P.toString() + " division by zero " + S);
899        }
900        if (P.isZERO()) {
901            return P;
902        }
903        if (S.isConstant()) {
904            return P.ring.getZERO();
905        }
906        C c = S.leadingBaseCoefficient();
907        ExpVector e = S.leadingExpVector();
908        GenPolynomial<C> h;
909        GenPolynomial<C> r = P;
910        while (!r.isZERO()) {
911            ExpVector f = r.leadingExpVector();
912            if (f.multipleOf(e)) {
913                C a = r.leadingBaseCoefficient();
914                ExpVector fe = f.subtract(e);
915                C x = a.remainder(c);
916                if (x.isZERO()) {
917                    C y = a.divide(c);
918                    h = S.multiply(y, fe); // coeff a
919                } else {
920                    r = r.multiply(c); // coeff ac
921                    h = S.multiply(a, fe); // coeff ac
922                }
923                r = r.subtract(h);
924            } else {
925                break;
926            }
927        }
928        return r;
929    }
930
931
932    /**
933     * GenPolynomial dense pseudo remainder. For univariate polynomials.
934     * @param P GenPolynomial.
935     * @param S nonzero GenPolynomial.
936     * @return remainder with ldcf(S)<sup>m</sup> P = quotient * S + remainder.
937     *         m == deg(P)-deg(S)
938     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
939     */
940    public static <C extends RingElem<C>> GenPolynomial<C> baseDensePseudoRemainder(GenPolynomial<C> P,
941                    GenPolynomial<C> S) {
942        if (S == null || S.isZERO()) {
943            throw new ArithmeticException(P + " division by zero " + S);
944        }
945        if (P.isZERO()) {
946            return P;
947        }
948        if (S.isConstant()) {
949            return P.ring.getZERO();
950        }
951        long m = P.degree(0);
952        long n = S.degree(0);
953        C c = S.leadingBaseCoefficient();
954        ExpVector e = S.leadingExpVector();
955        GenPolynomial<C> h;
956        GenPolynomial<C> r = P;
957        for (long i = m; i >= n; i--) {
958            if (r.isZERO()) {
959                return r;
960            }
961            long k = r.degree(0);
962            if (i == k) {
963                ExpVector f = r.leadingExpVector();
964                C a = r.leadingBaseCoefficient();
965                f = f.subtract(e); // EVDIF( f, e );
966                //System.out.println("red div = " + f);
967                r = r.multiply(c); // coeff ac
968                h = S.multiply(a, f); // coeff ac
969                r = r.subtract(h);
970            } else {
971                r = r.multiply(c);
972            }
973        }
974        return r;
975    }
976
977
978    /**
979     * GenPolynomial dense pseudo quotient. For univariate polynomials.
980     * @param P GenPolynomial.
981     * @param S nonzero GenPolynomial.
982     * @return quotient with ldcf(S)<sup>m</sup> P = quotient * S + remainder. m
983     *         == deg(P)-deg(S)
984     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
985     */
986    public static <C extends RingElem<C>> GenPolynomial<C> baseDensePseudoQuotient(GenPolynomial<C> P,
987                    GenPolynomial<C> S) {
988        if (S == null || S.isZERO()) {
989            throw new ArithmeticException(P + " division by zero " + S);
990        }
991        if (P.isZERO()) {
992            return P;
993        }
994        //if (S.degree() <= 0) {
995        //    return l^n P; //P.ring.getZERO();
996        //}
997        long m = P.degree(0);
998        long n = S.degree(0);
999        C c = S.leadingBaseCoefficient();
1000        ExpVector e = S.leadingExpVector();
1001        GenPolynomial<C> q = P.ring.getZERO();
1002        GenPolynomial<C> h;
1003        GenPolynomial<C> r = P;
1004        for (long i = m; i >= n; i--) {
1005            if (r.isZERO()) {
1006                return q;
1007            }
1008            long k = r.degree(0);
1009            if (i == k) {
1010                ExpVector f = r.leadingExpVector();
1011                C a = r.leadingBaseCoefficient();
1012                f = f.subtract(e); // EVDIF( f, e );
1013                //System.out.println("red div = " + f);
1014                r = r.multiply(c); // coeff ac
1015                h = S.multiply(a, f); // coeff ac
1016                r = r.subtract(h);
1017                q = q.multiply(c);
1018                q = q.sum(a, f);
1019            } else {
1020                q = q.multiply(c);
1021                r = r.multiply(c);
1022            }
1023        }
1024        return q;
1025    }
1026
1027
1028    /**
1029     * GenPolynomial sparse pseudo divide. For univariate polynomials or exact
1030     * division.
1031     * @param <C> coefficient type.
1032     * @param P GenPolynomial.
1033     * @param S nonzero GenPolynomial.
1034     * @return quotient with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
1035     *         m' &le; deg(P)-deg(S)
1036     * @see edu.jas.poly.GenPolynomial#divide(edu.jas.poly.GenPolynomial).
1037     */
1038    public static <C extends RingElem<C>> GenPolynomial<C> basePseudoDivide(GenPolynomial<C> P,
1039                    GenPolynomial<C> S) {
1040        if (S == null || S.isZERO()) {
1041            throw new ArithmeticException(P.toString() + " division by zero " + S);
1042        }
1043        //if (S.ring.nvar != 1) {
1044        // ok if exact division
1045        // throw new RuntimeException(this.getClass().getName()
1046        //                            + " univariate polynomials only");
1047        //}
1048        if (P.isZERO() || S.isONE()) {
1049            return P;
1050        }
1051        C c = S.leadingBaseCoefficient();
1052        ExpVector e = S.leadingExpVector();
1053        GenPolynomial<C> h;
1054        GenPolynomial<C> r = P;
1055        GenPolynomial<C> q = S.ring.getZERO().copy();
1056
1057        while (!r.isZERO()) {
1058            ExpVector f = r.leadingExpVector();
1059            if (f.multipleOf(e)) {
1060                C a = r.leadingBaseCoefficient();
1061                f = f.subtract(e);
1062                C x = a.remainder(c);
1063                if (x.isZERO()) {
1064                    C y = a.divide(c);
1065                    q = q.sum(y, f);
1066                    h = S.multiply(y, f); // coeff a
1067                } else {
1068                    q = q.multiply(c);
1069                    q = q.sum(a, f);
1070                    r = r.multiply(c); // coeff ac
1071                    h = S.multiply(a, f); // coeff ac
1072                }
1073                r = r.subtract(h);
1074            } else {
1075                break;
1076            }
1077        }
1078        return q;
1079    }
1080
1081
1082    /**
1083     * GenPolynomial sparse pseudo quotient and remainder. For univariate
1084     * polynomials or exact division.
1085     * @param <C> coefficient type.
1086     * @param P GenPolynomial.
1087     * @param S nonzero GenPolynomial.
1088     * @return [ quotient, remainder ] with ldcf(S)<sup>m'</sup> P = quotient *
1089     *         S + remainder. m' &le; deg(P)-deg(S)
1090     * @see edu.jas.poly.GenPolynomial#divide(edu.jas.poly.GenPolynomial).
1091     */
1092    @SuppressWarnings("unchecked")
1093    public static <C extends RingElem<C>> GenPolynomial<C>[] basePseudoQuotientRemainder(GenPolynomial<C> P,
1094                    GenPolynomial<C> S) {
1095        if (S == null || S.isZERO()) {
1096            throw new ArithmeticException(P.toString() + " division by zero " + S);
1097        }
1098        //if (S.ring.nvar != 1) {
1099        // ok if exact division
1100        // throw new RuntimeException(this.getClass().getName()
1101        //                            + " univariate polynomials only");
1102        //}
1103        GenPolynomial<C>[] ret = new GenPolynomial[2];
1104        ret[0] = null;
1105        ret[1] = null;
1106        if (P.isZERO() || S.isONE()) {
1107            ret[0] = P;
1108            ret[1] = S.ring.getZERO();
1109            return ret;
1110        }
1111        C c = S.leadingBaseCoefficient();
1112        ExpVector e = S.leadingExpVector();
1113        GenPolynomial<C> h;
1114        GenPolynomial<C> r = P;
1115        GenPolynomial<C> q = S.ring.getZERO().copy();
1116
1117        while (!r.isZERO()) {
1118            ExpVector f = r.leadingExpVector();
1119            if (f.multipleOf(e)) {
1120                C a = r.leadingBaseCoefficient();
1121                f = f.subtract(e);
1122                C x = a.remainder(c);
1123                if (x.isZERO()) {
1124                    C y = a.divide(c);
1125                    q = q.sum(y, f);
1126                    h = S.multiply(y, f); // coeff a
1127                } else {
1128                    q = q.multiply(c);
1129                    q = q.sum(a, f);
1130                    r = r.multiply(c); // coeff a c
1131                    h = S.multiply(a, f); // coeff c a
1132                }
1133                r = r.subtract(h);
1134            } else {
1135                break;
1136            }
1137        }
1138        //GenPolynomial<C> rhs = q.multiply(S).sum(r);
1139        //GenPolynomial<C> lhs = P;
1140        ret[0] = q;
1141        ret[1] = r;
1142        return ret;
1143    }
1144
1145
1146    /**
1147     * Is GenPolynomial pseudo quotient and remainder. For univariate
1148     * polynomials.
1149     * @param <C> coefficient type.
1150     * @param P base GenPolynomial.
1151     * @param S nonzero base GenPolynomial.
1152     * @return true, if P = q * S + r, else false.
1153     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1154     *      <b>Note:</b> not always meaningful and working
1155     */
1156    public static <C extends RingElem<C>> boolean isBasePseudoQuotientRemainder(GenPolynomial<C> P,
1157                    GenPolynomial<C> S, GenPolynomial<C> q, GenPolynomial<C> r) {
1158        GenPolynomial<C> rhs = q.multiply(S).sum(r);
1159        //System.out.println("rhs,1 = " + rhs);
1160        GenPolynomial<C> lhs = P;
1161        C ldcf = S.leadingBaseCoefficient();
1162        long d = P.degree(0) - S.degree(0) + 1;
1163        d = (d > 0 ? d : -d + 2);
1164        for (long i = 0; i <= d; i++) {
1165            //System.out.println("lhs-rhs = " + lhs.subtract(rhs));
1166            if (lhs.equals(rhs) || lhs.negate().equals(rhs)) {
1167                //System.out.println("lhs,1 = " + lhs);
1168                return true;
1169            }
1170            lhs = lhs.multiply(ldcf);
1171        }
1172        GenPolynomial<C> Pp = P;
1173        rhs = q.multiply(S);
1174        //System.out.println("rhs,2 = " + rhs);
1175        for (long i = 0; i <= d; i++) {
1176            lhs = Pp.subtract(r);
1177            //System.out.println("lhs-rhs = " + lhs.subtract(rhs));
1178            if (lhs.equals(rhs) || lhs.negate().equals(rhs)) {
1179                //System.out.println("lhs,2 = " + lhs);
1180                return true;
1181            }
1182            Pp = Pp.multiply(ldcf);
1183        }
1184        C a = P.leadingBaseCoefficient();
1185        rhs = q.multiply(S).sum(r);
1186        C b = rhs.leadingBaseCoefficient();
1187        C gcd = a.gcd(b);
1188        C p = a.multiply(b);
1189        C lcm = p.divide(gcd);
1190        C ap = lcm.divide(a);
1191        C bp = lcm.divide(b);
1192        if (P.multiply(ap).equals(rhs.multiply(bp))) {
1193            return true;
1194        }
1195        return false;
1196    }
1197
1198
1199    /**
1200     * GenPolynomial divide. For recursive polynomials. Division by coefficient
1201     * ring element.
1202     * @param <C> coefficient type.
1203     * @param P recursive GenPolynomial.
1204     * @param s GenPolynomial.
1205     * @return this/s.
1206     */
1207    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursiveDivide(
1208                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<C> s) {
1209        if (s == null || s.isZERO()) {
1210            throw new ArithmeticException("division by zero " + P + ", " + s);
1211        }
1212        if (P.isZERO()) {
1213            return P;
1214        }
1215        if (s.isONE()) {
1216            return P;
1217        }
1218        GenPolynomial<GenPolynomial<C>> p = P.ring.getZERO().copy();
1219        SortedMap<ExpVector, GenPolynomial<C>> pv = p.val; //getMap();
1220        for (Map.Entry<ExpVector, GenPolynomial<C>> m1 : P.getMap().entrySet()) {
1221            GenPolynomial<C> c1 = m1.getValue();
1222            ExpVector e1 = m1.getKey();
1223            GenPolynomial<C> c = PolyUtil.<C> basePseudoDivide(c1, s);
1224            if (!c.isZERO()) {
1225                pv.put(e1, c); // or m1.setValue( c )
1226            } else {
1227                logger.warn("rDiv, P  = " + P);
1228                logger.warn("rDiv, c1 = " + c1);
1229                logger.warn("rDiv, s  = " + s);
1230                logger.warn("rDiv, c  = " + c);
1231                throw new RuntimeException("something is wrong");
1232            }
1233        }
1234        return p;
1235    }
1236
1237
1238    /**
1239     * GenPolynomial divide. For recursive polynomials. Division by coefficient
1240     * ring element.
1241     * @param <C> coefficient type.
1242     * @param P recursive GenPolynomial.
1243     * @param s GenPolynomial.
1244     * @return this/s.
1245     */
1246    public static <C extends RingElem<C>> GenWordPolynomial<GenPolynomial<C>> recursiveDivide(
1247                    GenWordPolynomial<GenPolynomial<C>> P, GenPolynomial<C> s) {
1248        if (s == null || s.isZERO()) {
1249            throw new ArithmeticException("division by zero " + P + ", " + s);
1250        }
1251        if (P.isZERO()) {
1252            return P;
1253        }
1254        if (s.isONE()) {
1255            return P;
1256        }
1257        GenWordPolynomial<GenPolynomial<C>> p = P.ring.getZERO().copy();
1258        SortedMap<Word, GenPolynomial<C>> pv = p.val; //getMap();
1259        for (Map.Entry<Word, GenPolynomial<C>> m1 : P.getMap().entrySet()) {
1260            GenPolynomial<C> c1 = m1.getValue();
1261            Word e1 = m1.getKey();
1262            GenPolynomial<C> c = PolyUtil.<C> basePseudoDivide(c1, s);
1263            if (!c.isZERO()) {
1264                pv.put(e1, c); // or m1.setValue( c )
1265            } else {
1266                logger.warn("rDiv, P  = " + P);
1267                logger.warn("rDiv, c1 = " + c1);
1268                logger.warn("rDiv, s  = " + s);
1269                logger.warn("rDiv, c  = " + c);
1270                throw new RuntimeException("something is wrong");
1271            }
1272        }
1273        return p;
1274    }
1275
1276
1277    /**
1278     * GenPolynomial base divide. For recursive polynomials. Division by
1279     * coefficient ring element.
1280     * @param <C> coefficient type.
1281     * @param P recursive GenPolynomial.
1282     * @param s coefficient.
1283     * @return this/s.
1284     */
1285    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> baseRecursiveDivide(
1286                    GenPolynomial<GenPolynomial<C>> P, C s) {
1287        if (s == null || s.isZERO()) {
1288            throw new ArithmeticException("division by zero " + P + ", " + s);
1289        }
1290        if (P.isZERO()) {
1291            return P;
1292        }
1293        if (s.isONE()) {
1294            return P;
1295        }
1296        GenPolynomial<GenPolynomial<C>> p = P.ring.getZERO().copy();
1297        SortedMap<ExpVector, GenPolynomial<C>> pv = p.val; //getMap();
1298        for (Map.Entry<ExpVector, GenPolynomial<C>> m1 : P.getMap().entrySet()) {
1299            GenPolynomial<C> c1 = m1.getValue();
1300            ExpVector e1 = m1.getKey();
1301            GenPolynomial<C> c = PolyUtil.<C> coefficientBasePseudoDivide(c1, s);
1302            if (!c.isZERO()) {
1303                pv.put(e1, c); // or m1.setValue( c )
1304            } else {
1305                logger.warn("pu, c1 = " + c1);
1306                logger.warn("pu, s  = " + s);
1307                logger.warn("pu, c  = " + c);
1308                throw new RuntimeException("something is wrong");
1309            }
1310        }
1311        return p;
1312    }
1313
1314
1315    /**
1316     * GenPolynomial sparse pseudo remainder. For recursive polynomials.
1317     * @param <C> coefficient type.
1318     * @param P recursive GenPolynomial.
1319     * @param S nonzero recursive GenPolynomial.
1320     * @return remainder with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
1321     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1322     * @deprecated Use
1323     *             {@link #recursiveSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)}
1324     *             instead
1325     */
1326    @Deprecated
1327    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursivePseudoRemainder(
1328                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
1329        return recursiveSparsePseudoRemainder(P, S);
1330    }
1331
1332
1333    /**
1334     * GenPolynomial sparse pseudo remainder. For recursive polynomials.
1335     * @param <C> coefficient type.
1336     * @param P recursive GenPolynomial.
1337     * @param S nonzero recursive GenPolynomial.
1338     * @return remainder with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
1339     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1340     */
1341    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursiveSparsePseudoRemainder(
1342                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
1343        if (S == null || S.isZERO()) {
1344            throw new ArithmeticException(P + " division by zero " + S);
1345        }
1346        if (P == null || P.isZERO()) {
1347            return P;
1348        }
1349        if (S.isConstant()) {
1350            return P.ring.getZERO();
1351        }
1352        GenPolynomial<C> c = S.leadingBaseCoefficient();
1353        ExpVector e = S.leadingExpVector();
1354        GenPolynomial<GenPolynomial<C>> h;
1355        GenPolynomial<GenPolynomial<C>> r = P;
1356        while (!r.isZERO()) {
1357            ExpVector f = r.leadingExpVector();
1358            if (f.multipleOf(e)) {
1359                GenPolynomial<C> a = r.leadingBaseCoefficient();
1360                f = f.subtract(e);
1361                GenPolynomial<C> x = c; //test basePseudoRemainder(a,c);
1362                if (x.isZERO()) {
1363                    GenPolynomial<C> y = PolyUtil.<C> basePseudoDivide(a, c);
1364                    h = S.multiply(y, f); // coeff a
1365                } else {
1366                    r = r.multiply(c); // coeff a c
1367                    h = S.multiply(a, f); // coeff c a
1368                }
1369                r = r.subtract(h);
1370            } else {
1371                break;
1372            }
1373        }
1374        return r;
1375    }
1376
1377
1378    /**
1379     * GenPolynomial dense pseudo remainder. For recursive polynomials.
1380     * @param P recursive GenPolynomial.
1381     * @param S nonzero recursive GenPolynomial.
1382     * @return remainder with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
1383     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1384     */
1385    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursiveDensePseudoRemainder(
1386                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
1387        if (S == null || S.isZERO()) {
1388            throw new ArithmeticException(P + " division by zero " + S);
1389        }
1390        if (P == null || P.isZERO()) {
1391            return P;
1392        }
1393        if (S.isConstant()) {
1394            return P.ring.getZERO();
1395        }
1396        long m = P.degree(0);
1397        long n = S.degree(0);
1398        GenPolynomial<C> c = S.leadingBaseCoefficient();
1399        ExpVector e = S.leadingExpVector();
1400        GenPolynomial<GenPolynomial<C>> h;
1401        GenPolynomial<GenPolynomial<C>> r = P;
1402        for (long i = m; i >= n; i--) {
1403            if (r.isZERO()) {
1404                return r;
1405            }
1406            long k = r.degree(0);
1407            if (i == k) {
1408                ExpVector f = r.leadingExpVector();
1409                GenPolynomial<C> a = r.leadingBaseCoefficient();
1410                f = f.subtract(e); //EVDIF( f, e );
1411                //System.out.println("red div = " + f);
1412                r = r.multiply(c); // coeff ac
1413                h = S.multiply(a, f); // coeff ac
1414                r = r.subtract(h);
1415            } else {
1416                r = r.multiply(c);
1417            }
1418        }
1419        return r;
1420    }
1421
1422
1423    /**
1424     * GenPolynomial recursive pseudo divide. For recursive polynomials.
1425     * @param <C> coefficient type.
1426     * @param P recursive GenPolynomial.
1427     * @param S nonzero recursive GenPolynomial.
1428     * @return quotient with ldcf(S)<sup>m'</sup> P = quotient * S + remainder.
1429     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1430     */
1431    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursivePseudoDivide(
1432                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
1433        if (S == null || S.isZERO()) {
1434            throw new ArithmeticException(P + " division by zero " + S);
1435        }
1436        //if (S.ring.nvar != 1) {
1437        // ok if exact division
1438        // throw new RuntimeException(this.getClass().getName()
1439        //                            + " univariate polynomials only");
1440        //}
1441        if (P == null || P.isZERO()) {
1442            return P;
1443        }
1444        if (S.isONE()) {
1445            return P;
1446        }
1447        GenPolynomial<C> c = S.leadingBaseCoefficient();
1448        ExpVector e = S.leadingExpVector();
1449        GenPolynomial<GenPolynomial<C>> h;
1450        GenPolynomial<GenPolynomial<C>> r = P;
1451        GenPolynomial<GenPolynomial<C>> q = S.ring.getZERO().copy();
1452        while (!r.isZERO()) {
1453            ExpVector f = r.leadingExpVector();
1454            if (f.multipleOf(e)) {
1455                GenPolynomial<C> a = r.leadingBaseCoefficient();
1456                f = f.subtract(e);
1457                GenPolynomial<C> x = PolyUtil.<C> baseSparsePseudoRemainder(a, c);
1458                if (x.isZERO() && !c.isConstant()) {
1459                    GenPolynomial<C> y = PolyUtil.<C> basePseudoDivide(a, c);
1460                    q = q.sum(y, f);
1461                    h = S.multiply(y, f); // coeff a
1462                } else {
1463                    q = q.multiply(c);
1464                    q = q.sum(a, f);
1465                    r = r.multiply(c); // coeff ac
1466                    h = S.multiply(a, f); // coeff ac
1467                }
1468                r = r.subtract(h);
1469            } else {
1470                break;
1471            }
1472        }
1473        return q;
1474    }
1475
1476
1477    /**
1478     * Is recursive GenPolynomial pseudo quotient and remainder. For recursive
1479     * polynomials.
1480     * @param <C> coefficient type.
1481     * @param P recursive GenPolynomial.
1482     * @param S nonzero recursive GenPolynomial.
1483     * @return true, if P ~= q * S + r, else false.
1484     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1485     *      <b>Note:</b> not always meaningful and working
1486     */
1487    public static <C extends RingElem<C>> boolean isRecursivePseudoQuotientRemainder(
1488                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S,
1489                    GenPolynomial<GenPolynomial<C>> q, GenPolynomial<GenPolynomial<C>> r) {
1490        GenPolynomial<GenPolynomial<C>> rhs = q.multiply(S).sum(r);
1491        GenPolynomial<GenPolynomial<C>> lhs = P;
1492        GenPolynomial<C> ldcf = S.leadingBaseCoefficient();
1493        long d = P.degree(0) - S.degree(0) + 1;
1494        d = (d > 0 ? d : -d + 2);
1495        for (long i = 0; i <= d; i++) {
1496            //System.out.println("lhs = " + lhs);
1497            //System.out.println("rhs = " + rhs);
1498            //System.out.println("lhs-rhs = " + lhs.subtract(rhs));
1499            if (lhs.equals(rhs)) {
1500                return true;
1501            }
1502            lhs = lhs.multiply(ldcf);
1503        }
1504        GenPolynomial<GenPolynomial<C>> Pp = P;
1505        rhs = q.multiply(S);
1506        //System.out.println("rhs,2 = " + rhs);
1507        for (long i = 0; i <= d; i++) {
1508            lhs = Pp.subtract(r);
1509            //System.out.println("lhs-rhs = " + lhs.subtract(rhs));
1510            if (lhs.equals(rhs)) {
1511                //System.out.println("lhs,2 = " + lhs);
1512                return true;
1513            }
1514            Pp = Pp.multiply(ldcf);
1515        }
1516        GenPolynomial<C> a = P.leadingBaseCoefficient();
1517        rhs = q.multiply(S).sum(r);
1518        GenPolynomial<C> b = rhs.leadingBaseCoefficient();
1519        if (P.multiply(b).equals(rhs.multiply(a))) {
1520            return true;
1521        }
1522        return false;
1523    }
1524
1525
1526    /**
1527     * GenPolynomial pseudo divide. For recursive polynomials.
1528     * @param <C> coefficient type.
1529     * @param P recursive GenPolynomial.
1530     * @param s nonzero GenPolynomial.
1531     * @return quotient with ldcf(s)<sup>m</sup> P = quotient * s + remainder.
1532     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1533     */
1534    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> coefficientPseudoDivide(
1535                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<C> s) {
1536        if (s == null || s.isZERO()) {
1537            throw new ArithmeticException(P + " division by zero " + s);
1538        }
1539        if (P.isZERO()) {
1540            return P;
1541        }
1542        GenPolynomial<GenPolynomial<C>> p = P.ring.getZERO().copy();
1543        SortedMap<ExpVector, GenPolynomial<C>> pv = p.val;
1544        for (Map.Entry<ExpVector, GenPolynomial<C>> m : P.getMap().entrySet()) {
1545            ExpVector e = m.getKey();
1546            GenPolynomial<C> c1 = m.getValue();
1547            GenPolynomial<C> c = basePseudoDivide(c1, s);
1548            if (debug) {
1549                GenPolynomial<C> x = c1.remainder(s);
1550                if (!x.isZERO()) {
1551                    logger.info("divide x = " + x);
1552                    throw new ArithmeticException("no exact division: " + c1 + "/" + s);
1553                }
1554            }
1555            if (c.isZERO()) {
1556                //logger.warn("no exact division: " + c1 + "/" + s);
1557                throw new ArithmeticException("no exact division: " + c1 + "/" + s);
1558            }
1559            pv.put(e, c); // or m1.setValue( c )
1560        }
1561        return p;
1562    }
1563
1564
1565    /**
1566     * GenPolynomial pseudo divide. For polynomials.
1567     * @param <C> coefficient type.
1568     * @param P GenPolynomial.
1569     * @param s nonzero coefficient.
1570     * @return quotient with ldcf(s)<sup>m</sup> P = quotient * s + remainder.
1571     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
1572     */
1573    public static <C extends RingElem<C>> GenPolynomial<C> coefficientBasePseudoDivide(GenPolynomial<C> P,
1574                    C s) {
1575        if (s == null || s.isZERO()) {
1576            throw new ArithmeticException(P + " division by zero " + s);
1577        }
1578        if (P.isZERO()) {
1579            return P;
1580        }
1581        GenPolynomial<C> p = P.ring.getZERO().copy();
1582        SortedMap<ExpVector, C> pv = p.val;
1583        for (Map.Entry<ExpVector, C> m : P.getMap().entrySet()) {
1584            ExpVector e = m.getKey();
1585            C c1 = m.getValue();
1586            C c = c1.divide(s);
1587            if (debug) {
1588                C x = c1.remainder(s);
1589                if (!x.isZERO()) {
1590                    logger.info("divide x = " + x);
1591                    throw new ArithmeticException("no exact division: " + c1 + "/" + s);
1592                }
1593            }
1594            if (c.isZERO()) {
1595                //logger.warn("no exact division: " + c1 + "/" + s);
1596                throw new ArithmeticException("no exact division: " + c1 + "/" + s);
1597            }
1598            pv.put(e, c); // or m1.setValue( c )
1599        }
1600        return p;
1601    }
1602
1603
1604    /**
1605     * GenPolynomial polynomial derivative main variable.
1606     * @param <C> coefficient type.
1607     * @param P GenPolynomial.
1608     * @return deriviative(P).
1609     */
1610    public static <C extends RingElem<C>> GenPolynomial<C> baseDeriviative(GenPolynomial<C> P) {
1611        if (P == null || P.isZERO()) {
1612            return P;
1613        }
1614        GenPolynomialRing<C> pfac = P.ring;
1615        if (pfac.nvar == 0) {
1616            return pfac.getZERO();
1617        }
1618        if (pfac.nvar > 1) {
1619            // baseContent not possible by return type
1620            throw new IllegalArgumentException(P.getClass().getName() + " only for univariate polynomials");
1621        }
1622        RingFactory<C> rf = pfac.coFac;
1623        GenPolynomial<C> d = pfac.getZERO().copy();
1624        Map<ExpVector, C> dm = d.val; //getMap();
1625        for (Map.Entry<ExpVector, C> m : P.getMap().entrySet()) {
1626            ExpVector f = m.getKey();
1627            long fl = f.getVal(0);
1628            if (fl > 0) {
1629                C cf = rf.fromInteger(fl);
1630                C a = m.getValue();
1631                C x = a.multiply(cf);
1632                if (x != null && !x.isZERO()) {
1633                    ExpVector e = ExpVector.create(1, 0, fl - 1L);
1634                    dm.put(e, x);
1635                }
1636            }
1637        }
1638        return d;
1639    }
1640
1641
1642    /**
1643     * GenPolynomial polynomial partial derivative variable r.
1644     * @param <C> coefficient type.
1645     * @param P GenPolynomial.
1646     * @param r variable for partial deriviate.
1647     * @return deriviative(P,r).
1648     */
1649    public static <C extends RingElem<C>> GenPolynomial<C> baseDeriviative(GenPolynomial<C> P, int r) {
1650        if (P == null || P.isZERO()) {
1651            return P;
1652        }
1653        GenPolynomialRing<C> pfac = P.ring;
1654        if (r < 0 || pfac.nvar <= r) {
1655            throw new IllegalArgumentException(
1656                            P.getClass().getName() + " deriviative variable out of bound " + r);
1657        }
1658        int rp = pfac.nvar - 1 - r;
1659        RingFactory<C> rf = pfac.coFac;
1660        GenPolynomial<C> d = pfac.getZERO().copy();
1661        Map<ExpVector, C> dm = d.val; //getMap();
1662        for (Map.Entry<ExpVector, C> m : P.getMap().entrySet()) {
1663            ExpVector f = m.getKey();
1664            long fl = f.getVal(rp);
1665            if (fl > 0) {
1666                C cf = rf.fromInteger(fl);
1667                C a = m.getValue();
1668                C x = a.multiply(cf);
1669                if (x != null && !x.isZERO()) {
1670                    ExpVector e = f.subst(rp, fl - 1L);
1671                    dm.put(e, x);
1672                }
1673            }
1674        }
1675        return d;
1676    }
1677
1678
1679    /**
1680     * GenPolynomial polynomial integral main variable.
1681     * @param <C> coefficient type.
1682     * @param P GenPolynomial.
1683     * @return integral(P).
1684     */
1685    public static <C extends RingElem<C>> GenPolynomial<C> baseIntegral(GenPolynomial<C> P) {
1686        if (P == null || P.isZERO()) {
1687            return P;
1688        }
1689        GenPolynomialRing<C> pfac = P.ring;
1690        if (pfac.nvar == 0) {
1691            return pfac.getONE();
1692        }
1693        if (pfac.nvar > 1) {
1694            // baseContent not possible by return type
1695            throw new IllegalArgumentException(P.getClass().getName() + " only for univariate polynomials");
1696        }
1697        RingFactory<C> rf = pfac.coFac;
1698        GenPolynomial<C> d = pfac.getZERO().copy();
1699        Map<ExpVector, C> dm = d.val; //getMap();
1700        for (Map.Entry<ExpVector, C> m : P.getMap().entrySet()) {
1701            ExpVector f = m.getKey();
1702            long fl = f.getVal(0);
1703            fl++;
1704            C cf = rf.fromInteger(fl);
1705            C a = m.getValue();
1706            C x = a.divide(cf);
1707            if (x != null && !x.isZERO()) {
1708                ExpVector e = ExpVector.create(1, 0, fl);
1709                dm.put(e, x);
1710            }
1711        }
1712        return d;
1713    }
1714
1715
1716    /**
1717     * GenPolynomial recursive polynomial derivative main variable.
1718     * @param <C> coefficient type.
1719     * @param P recursive GenPolynomial.
1720     * @return deriviative(P).
1721     */
1722    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> recursiveDeriviative(
1723                    GenPolynomial<GenPolynomial<C>> P) {
1724        if (P == null || P.isZERO()) {
1725            return P;
1726        }
1727        GenPolynomialRing<GenPolynomial<C>> pfac = P.ring;
1728        if (pfac.nvar == 0) {
1729            return pfac.getZERO();
1730        }
1731        if (pfac.nvar > 1) {
1732            // baseContent not possible by return type
1733            throw new IllegalArgumentException(P.getClass().getName() + " only for univariate polynomials");
1734        }
1735        GenPolynomialRing<C> pr = (GenPolynomialRing<C>) pfac.coFac;
1736        RingFactory<C> rf = pr.coFac;
1737        GenPolynomial<GenPolynomial<C>> d = pfac.getZERO().copy();
1738        Map<ExpVector, GenPolynomial<C>> dm = d.val; //getMap();
1739        for (Map.Entry<ExpVector, GenPolynomial<C>> m : P.getMap().entrySet()) {
1740            ExpVector f = m.getKey();
1741            long fl = f.getVal(0);
1742            if (fl > 0) {
1743                C cf = rf.fromInteger(fl);
1744                GenPolynomial<C> a = m.getValue();
1745                GenPolynomial<C> x = a.multiply(cf);
1746                if (x != null && !x.isZERO()) {
1747                    ExpVector e = ExpVector.create(1, 0, fl - 1L);
1748                    dm.put(e, x);
1749                }
1750            }
1751        }
1752        return d;
1753    }
1754
1755
1756    /**
1757     * Factor coefficient bound. See SACIPOL.IPFCB: the product of all maxNorms
1758     * of potential factors is less than or equal to 2**b times the maxNorm of
1759     * A.
1760     * @param e degree vector of a GenPolynomial A.
1761     * @return 2**b.
1762     */
1763    public static BigInteger factorBound(ExpVector e) {
1764        int n = 0;
1765        java.math.BigInteger p = java.math.BigInteger.ONE;
1766        java.math.BigInteger v;
1767        if (e == null || e.isZERO()) {
1768            return BigInteger.ONE;
1769        }
1770        for (int i = 0; i < e.length(); i++) {
1771            if (e.getVal(i) > 0) {
1772                n += (2 * e.getVal(i) - 1);
1773                v = new java.math.BigInteger("" + (e.getVal(i) - 1));
1774                p = p.multiply(v);
1775            }
1776        }
1777        n += (p.bitCount() + 1); // log2(p)
1778        n /= 2;
1779        v = new java.math.BigInteger("" + 2);
1780        v = v.shiftLeft(n);
1781        BigInteger N = new BigInteger(v);
1782        return N;
1783    }
1784
1785
1786    /**
1787     * Absoulte norm. Square root of the sum of the squared coefficients.
1788     * @param p GenPolynomial
1789     * @return sqrt( sum<sub>i</sub> |c<sub>i</sub>|<sup>2</sup> ).
1790     */
1791    @SuppressWarnings("unchecked")
1792    public static <C extends RingElem<C>> C absNorm(GenPolynomial<C> p) {
1793        if (p == null) {
1794            return null;
1795        }
1796        C a = p.ring.getZEROCoefficient();
1797        if (a instanceof StarRingElem) {
1798            //System.out.println("StarRingElem case");
1799            for (C c : p.val.values()) {
1800                @SuppressWarnings("unchecked")
1801                C n = (C) ((StarRingElem) c).norm();
1802                a = a.sum(n);
1803            }
1804        } else {
1805            for (C c : p.val.values()) {
1806                C n = c.multiply(c);
1807                a = a.sum(n);
1808            }
1809        }
1810        // compute square root if possible
1811        // todo refactor for sqrt in RingElem (?)
1812        if (a instanceof BigRational) {
1813            BigRational b = (BigRational) a;
1814            a = (C) Roots.sqrt(b);
1815        } else if (a instanceof BigComplex) {
1816            BigComplex b = (BigComplex) a;
1817            a = (C) Roots.sqrt(b);
1818        } else if (a instanceof BigInteger) {
1819            BigInteger b = (BigInteger) a;
1820            a = (C) Roots.sqrt(b);
1821        } else if (a instanceof BigDecimal) {
1822            BigDecimal b = (BigDecimal) a;
1823            a = (C) Roots.sqrt(b);
1824        } else if (a instanceof BigDecimalComplex) {
1825            BigDecimalComplex b = (BigDecimalComplex) a;
1826            a = (C) Roots.sqrt(b);
1827        } else {
1828            logger.error("no square root implemented for " + a.toScriptFactory());
1829        }
1830        return a;
1831    }
1832
1833
1834    /**
1835     * Evaluate at main variable.
1836     * @param <C> coefficient type.
1837     * @param cfac coefficent polynomial ring factory.
1838     * @param A recursive polynomial to be evaluated.
1839     * @param a value to evaluate at.
1840     * @return A( x_1, ..., x_{n-1}, a ).
1841     */
1842    public static <C extends RingElem<C>> GenPolynomial<C> evaluateMainRecursive(GenPolynomialRing<C> cfac,
1843                    GenPolynomial<GenPolynomial<C>> A, C a) {
1844        if (A == null || A.isZERO()) {
1845            return cfac.getZERO();
1846        }
1847        if (A.ring.nvar != 1) {
1848            throw new IllegalArgumentException("evaluateMain no univariate polynomial");
1849        }
1850        if (a == null || a.isZERO()) {
1851            return A.trailingBaseCoefficient();
1852        }
1853        // assert descending exponents, i.e. compatible term order
1854        Map<ExpVector, GenPolynomial<C>> val = A.getMap();
1855        GenPolynomial<C> B = null;
1856        long el1 = -1; // undefined
1857        long el2 = -1;
1858        for (Map.Entry<ExpVector, GenPolynomial<C>> me : val.entrySet()) {
1859            ExpVector e = me.getKey();
1860            el2 = e.getVal(0);
1861            if (B == null /*el1 < 0*/) { // first turn
1862                B = me.getValue();
1863            } else {
1864                for (long i = el2; i < el1; i++) {
1865                    B = B.multiply(a);
1866                }
1867                B = B.sum(me.getValue());
1868            }
1869            el1 = el2;
1870        }
1871        for (long i = 0; i < el2; i++) {
1872            B = B.multiply(a);
1873        }
1874        return B;
1875    }
1876
1877
1878    /**
1879     * Evaluate at main variable.
1880     * @param <C> coefficient type.
1881     * @param cfac coefficent polynomial ring factory.
1882     * @param A distributed polynomial to be evaluated.
1883     * @param a value to evaluate at.
1884     * @return A( x_1, ..., x_{n-1}, a ).
1885     */
1886    public static <C extends RingElem<C>> GenPolynomial<C> evaluateMain(GenPolynomialRing<C> cfac,
1887                    GenPolynomial<C> A, C a) {
1888        if (A == null || A.isZERO()) {
1889            return cfac.getZERO();
1890        }
1891        GenPolynomialRing<GenPolynomial<C>> rfac = A.ring.recursive(1);
1892        if (rfac.nvar + cfac.nvar != A.ring.nvar) {
1893            throw new IllegalArgumentException("evaluateMain number of variabes mismatch");
1894        }
1895        GenPolynomial<GenPolynomial<C>> Ap = recursive(rfac, A);
1896        return PolyUtil.<C> evaluateMainRecursive(cfac, Ap, a);
1897    }
1898
1899
1900    /**
1901     * Evaluate at main variable.
1902     * @param <C> coefficient type.
1903     * @param cfac coefficent ring factory.
1904     * @param L list of univariate polynomials to be evaluated.
1905     * @param a value to evaluate at.
1906     * @return list( A( x_1, ..., x_{n-1}, a ) ) for A in L.
1907     */
1908    public static <C extends RingElem<C>> List<GenPolynomial<C>> evaluateMain(GenPolynomialRing<C> cfac,
1909                    List<GenPolynomial<C>> L, C a) {
1910        return ListUtil.<GenPolynomial<C>, GenPolynomial<C>> map(L, new EvalMainPol<C>(cfac, a));
1911    }
1912
1913
1914    /**
1915     * Evaluate at main variable.
1916     * @param <C> coefficient type.
1917     * @param cfac coefficent ring factory.
1918     * @param A univariate polynomial to be evaluated.
1919     * @param a value to evaluate at.
1920     * @return A( a ).
1921     */
1922    public static <C extends RingElem<C>> C evaluateMain(RingFactory<C> cfac, GenPolynomial<C> A, C a) {
1923        if (A == null || A.isZERO()) {
1924            return cfac.getZERO();
1925        }
1926        if (A.ring.nvar != 1) {
1927            throw new IllegalArgumentException("evaluateMain no univariate polynomial");
1928        }
1929        if (a == null || a.isZERO()) {
1930            return A.trailingBaseCoefficient();
1931        }
1932        // assert decreasing exponents, i.e. compatible term order
1933        Map<ExpVector, C> val = A.getMap();
1934        C B = null;
1935        long el1 = -1; // undefined
1936        long el2 = -1;
1937        for (Map.Entry<ExpVector, C> me : val.entrySet()) {
1938            ExpVector e = me.getKey();
1939            el2 = e.getVal(0);
1940            if (B == null /*el1 < 0*/) { // first turn
1941                B = me.getValue();
1942            } else {
1943                for (long i = el2; i < el1; i++) {
1944                    B = B.multiply(a);
1945                }
1946                B = B.sum(me.getValue());
1947            }
1948            el1 = el2;
1949        }
1950        for (long i = 0; i < el2; i++) {
1951            B = B.multiply(a);
1952        }
1953        return B;
1954    }
1955
1956
1957    /**
1958     * Evaluate at main variable.
1959     * @param <C> coefficient type.
1960     * @param cfac coefficent ring factory.
1961     * @param L list of univariate polynomial to be evaluated.
1962     * @param a value to evaluate at.
1963     * @return list( A( a ) ) for A in L.
1964     */
1965    public static <C extends RingElem<C>> List<C> evaluateMain(RingFactory<C> cfac, List<GenPolynomial<C>> L,
1966                    C a) {
1967        return ListUtil.<GenPolynomial<C>, C> map(L, new EvalMain<C>(cfac, a));
1968    }
1969
1970
1971    /**
1972     * Evaluate at k-th variable.
1973     * @param <C> coefficient type.
1974     * @param cfac coefficient polynomial ring in k variables C[x_1, ..., x_k]
1975     *            factory.
1976     * @param rfac coefficient polynomial ring C[x_1, ..., x_{k-1}] [x_k]
1977     *            factory, a recursive polynomial ring in 1 variable with
1978     *            coefficients in k-1 variables.
1979     * @param nfac polynomial ring in n-1 varaibles C[x_1, ..., x_{k-1}]
1980     *            [x_{k+1}, ..., x_n] factory, a recursive polynomial ring in
1981     *            n-k+1 variables with coefficients in k-1 variables.
1982     * @param dfac polynomial ring in n-1 variables. C[x_1, ..., x_{k-1},
1983     *            x_{k+1}, ..., x_n] factory.
1984     * @param A polynomial to be evaluated.
1985     * @param a value to evaluate at.
1986     * @return A( x_1, ..., x_{k-1}, a, x_{k+1}, ..., x_n).
1987     */
1988    public static <C extends RingElem<C>> GenPolynomial<C> evaluate(GenPolynomialRing<C> cfac,
1989                    GenPolynomialRing<GenPolynomial<C>> rfac, GenPolynomialRing<GenPolynomial<C>> nfac,
1990                    GenPolynomialRing<C> dfac, GenPolynomial<C> A, C a) {
1991        if (rfac.nvar != 1) {
1992            throw new IllegalArgumentException("evaluate coefficient ring not univariate");
1993        }
1994        if (A == null || A.isZERO()) {
1995            return cfac.getZERO();
1996        }
1997        Map<ExpVector, GenPolynomial<C>> Ap = A.contract(cfac);
1998        GenPolynomialRing<C> rcf = (GenPolynomialRing<C>) rfac.coFac;
1999        GenPolynomial<GenPolynomial<C>> Ev = nfac.getZERO().copy();
2000        Map<ExpVector, GenPolynomial<C>> Evm = Ev.val; //getMap();
2001        for (Map.Entry<ExpVector, GenPolynomial<C>> m : Ap.entrySet()) {
2002            ExpVector e = m.getKey();
2003            GenPolynomial<C> b = m.getValue();
2004            GenPolynomial<GenPolynomial<C>> c = recursive(rfac, b);
2005            GenPolynomial<C> d = evaluateMainRecursive(rcf, c, a);
2006            if (d != null && !d.isZERO()) {
2007                Evm.put(e, d);
2008            }
2009        }
2010        GenPolynomial<C> B = distribute(dfac, Ev);
2011        return B;
2012    }
2013
2014
2015    /**
2016     * Evaluate at first (lowest) variable.
2017     * @param <C> coefficient type.
2018     * @param cfac coefficient polynomial ring in first variable C[x_1] factory.
2019     * @param dfac polynomial ring in n-1 variables. C[x_2, ..., x_n] factory.
2020     * @param A polynomial to be evaluated.
2021     * @param a value to evaluate at.
2022     * @return A( a, x_2, ..., x_n).
2023     */
2024    public static <C extends RingElem<C>> GenPolynomial<C> evaluateFirst(GenPolynomialRing<C> cfac,
2025                    GenPolynomialRing<C> dfac, GenPolynomial<C> A, C a) {
2026        if (A == null || A.isZERO()) {
2027            return dfac.getZERO();
2028        }
2029        Map<ExpVector, GenPolynomial<C>> Ap = A.contract(cfac);
2030        //RingFactory<C> rcf = cfac.coFac; // == dfac.coFac
2031
2032        GenPolynomial<C> B = dfac.getZERO().copy();
2033        Map<ExpVector, C> Bm = B.val; //getMap();
2034
2035        for (Map.Entry<ExpVector, GenPolynomial<C>> m : Ap.entrySet()) {
2036            ExpVector e = m.getKey();
2037            GenPolynomial<C> b = m.getValue();
2038            C d = evaluateMain(cfac.coFac, b, a);
2039            if (d != null && !d.isZERO()) {
2040                Bm.put(e, d);
2041            }
2042        }
2043        return B;
2044    }
2045
2046
2047    /**
2048     * Evaluate at first (lowest) variable. Could also be called
2049     * <code>evaluateFirst()</code>, but type erasure of parameter
2050     * <code>A</code> does not allow the same name.
2051     * @param <C> coefficient type.
2052     * @param cfac coefficient polynomial ring in first variable C[x_1] factory.
2053     * @param dfac polynomial ring in n-1 variables. C[x_2, ..., x_n] factory.
2054     * @param A recursive polynomial to be evaluated.
2055     * @param a value to evaluate at.
2056     * @return A( a, x_2, ..., x_n).
2057     */
2058    public static <C extends RingElem<C>> GenPolynomial<C> evaluateFirstRec(GenPolynomialRing<C> cfac,
2059                    GenPolynomialRing<C> dfac, GenPolynomial<GenPolynomial<C>> A, C a) {
2060        if (A == null || A.isZERO()) {
2061            return dfac.getZERO();
2062        }
2063        Map<ExpVector, GenPolynomial<C>> Ap = A.getMap();
2064        GenPolynomial<C> B = dfac.getZERO().copy();
2065        Map<ExpVector, C> Bm = B.val; //getMap();
2066        for (Map.Entry<ExpVector, GenPolynomial<C>> m : Ap.entrySet()) {
2067            ExpVector e = m.getKey();
2068            GenPolynomial<C> b = m.getValue();
2069            C d = evaluateMain(cfac.coFac, b, a);
2070            if (d != null && !d.isZERO()) {
2071                Bm.put(e, d);
2072            }
2073        }
2074        return B;
2075    }
2076
2077
2078    /**
2079     * Evaluate all variables.
2080     * @param <C> coefficient type.
2081     * @param cfac coefficient ring factory.
2082     * @param A polynomial to be evaluated.
2083     * @param a = (a_1, a_2, ..., a_n) a tuple of values to evaluate at.
2084     * @return A(a_1, a_2, ..., a_n).
2085     */
2086    public static <C extends RingElem<C>> C evaluateAll(RingFactory<C> cfac, GenPolynomial<C> A, List<C> a) {
2087        if (A == null || A.isZERO()) {
2088            return cfac.getZERO();
2089        }
2090        GenPolynomialRing<C> dfac = A.ring;
2091        if (a == null || a.size() != dfac.nvar) {
2092            throw new IllegalArgumentException("evaluate tuple size not equal to number of variables");
2093        }
2094        if (dfac.nvar == 0) {
2095            return A.trailingBaseCoefficient();
2096        }
2097        if (dfac.nvar == 1) {
2098            return evaluateMain(cfac, A, a.get(0));
2099        }
2100        C b = cfac.getZERO();
2101        GenPolynomial<C> Ap = A;
2102        for (int k = 0; k < dfac.nvar - 1; k++) {
2103            C ap = a.get(k);
2104            GenPolynomialRing<C> c1fac = new GenPolynomialRing<C>(cfac, 1); // no vars
2105            GenPolynomialRing<C> cnfac = new GenPolynomialRing<C>(cfac, dfac.nvar - 1 - k); // no vars
2106            GenPolynomial<C> Bp = evaluateFirst(c1fac, cnfac, Ap, ap);
2107            if (Bp.isZERO()) {
2108                return b;
2109            }
2110            Ap = Bp;
2111            //System.out.println("Ap = " + Ap);
2112        }
2113        C ap = a.get(dfac.nvar - 1);
2114        b = evaluateMain(cfac, Ap, ap);
2115        return b;
2116    }
2117
2118
2119    /**
2120     * Substitute main variable.
2121     * @param A univariate polynomial.
2122     * @param s polynomial for substitution.
2123     * @return polynomial A(x <- s).
2124     */
2125    public static <C extends RingElem<C>> GenPolynomial<C> substituteMain(GenPolynomial<C> A,
2126                    GenPolynomial<C> s) {
2127        return substituteUnivariate(A, s);
2128    }
2129
2130
2131    /**
2132     * Substitute univariate polynomial.
2133     * @param f univariate polynomial.
2134     * @param t polynomial for substitution.
2135     * @return polynomial f(x <- t).
2136     */
2137    public static <C extends RingElem<C>> GenPolynomial<C> substituteUnivariate(GenPolynomial<C> f,
2138                    GenPolynomial<C> t) {
2139        if (f == null || t == null) {
2140            return null;
2141        }
2142        GenPolynomialRing<C> fac = f.ring;
2143        if (fac.nvar > 1) {
2144            throw new IllegalArgumentException("only for univariate polynomial f");
2145        }
2146        if (f.isZERO() || f.isConstant()) {
2147            return f;
2148        }
2149        if (t.ring.nvar > 1) {
2150            fac = t.ring;
2151        }
2152        // assert decending exponents, i.e. compatible term order
2153        Map<ExpVector, C> val = f.getMap();
2154        GenPolynomial<C> s = null;
2155        long el1 = -1; // undefined
2156        long el2 = -1;
2157        for (Map.Entry<ExpVector, C> me : val.entrySet()) {
2158            ExpVector e = me.getKey();
2159            el2 = e.getVal(0);
2160            if (s == null /*el1 < 0*/) { // first turn
2161                s = fac.getZERO().sum(me.getValue());
2162            } else {
2163                for (long i = el2; i < el1; i++) {
2164                    s = s.multiply(t);
2165                }
2166                s = s.sum(me.getValue());
2167            }
2168            el1 = el2;
2169        }
2170        for (long i = 0; i < el2; i++) {
2171            s = s.multiply(t);
2172        }
2173        //System.out.println("s = " + s);
2174        return s;
2175    }
2176
2177
2178    /**
2179     * Substitute univariate polynomial with multivariate coefficients.
2180     * @param f univariate polynomial with multivariate coefficients.
2181     * @param t polynomial for substitution.
2182     * @return polynomial f(x <- t).
2183     */
2184    public static <C extends RingElem<C>> GenPolynomial<C> substituteUnivariateMult(GenPolynomial<C> f,
2185                    GenPolynomial<C> t) {
2186        if (f == null || t == null) {
2187            return null;
2188        }
2189        GenPolynomialRing<C> fac = f.ring;
2190        if (fac.nvar == 1) {
2191            return substituteUnivariate(f, t);
2192        }
2193        GenPolynomialRing<GenPolynomial<C>> rfac = fac.recursive(1);
2194        GenPolynomial<GenPolynomial<C>> fr = PolyUtil.<C> recursive(rfac, f);
2195        GenPolynomial<GenPolynomial<C>> tr = PolyUtil.<C> recursive(rfac, t);
2196        //System.out.println("fr = " + fr);
2197        //System.out.println("tr = " + tr);
2198        GenPolynomial<GenPolynomial<C>> sr = PolyUtil.<GenPolynomial<C>> substituteUnivariate(fr, tr);
2199        //System.out.println("sr = " + sr);
2200        GenPolynomial<C> s = PolyUtil.<C> distribute(fac, sr);
2201        return s;
2202    }
2203
2204
2205    /**
2206     * Taylor series for polynomial.
2207     * @param f univariate polynomial.
2208     * @param a expansion point.
2209     * @return Taylor series (a polynomial) of f at a.
2210     */
2211    public static <C extends RingElem<C>> GenPolynomial<C> seriesOfTaylor(GenPolynomial<C> f, C a) {
2212        if (f == null) {
2213            return null;
2214        }
2215        GenPolynomialRing<C> fac = f.ring;
2216        if (fac.nvar > 1) {
2217            throw new IllegalArgumentException("only for univariate polynomials");
2218        }
2219        if (f.isZERO() || f.isConstant()) {
2220            return f;
2221        }
2222        GenPolynomial<C> s = fac.getZERO();
2223        C fa = PolyUtil.<C> evaluateMain(fac.coFac, f, a);
2224        s = s.sum(fa);
2225        long n = 1;
2226        long i = 0;
2227        GenPolynomial<C> g = PolyUtil.<C> baseDeriviative(f);
2228        //GenPolynomial<C> p = fac.getONE();
2229        while (!g.isZERO()) {
2230            i++;
2231            n *= i;
2232            fa = PolyUtil.<C> evaluateMain(fac.coFac, g, a);
2233            GenPolynomial<C> q = fac.univariate(0, i); //p;
2234            q = q.multiply(fa);
2235            q = q.divide(fac.fromInteger(n));
2236            s = s.sum(q);
2237            g = PolyUtil.<C> baseDeriviative(g);
2238        }
2239        //System.out.println("s = " + s);
2240        return s;
2241    }
2242
2243
2244    /**
2245     * ModInteger interpolate on first variable.
2246     * @param <C> coefficient type.
2247     * @param fac GenPolynomial<C> result factory.
2248     * @param A GenPolynomial.
2249     * @param M GenPolynomial interpolation modul of A.
2250     * @param mi inverse of M(am) in ring fac.coFac.
2251     * @param B evaluation of other GenPolynomial.
2252     * @param am evaluation point (interpolation modul) of B, i.e. P(am) = B.
2253     * @return S, with S mod M == A and S(am) == B.
2254     */
2255    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> interpolate(
2256                    GenPolynomialRing<GenPolynomial<C>> fac, GenPolynomial<GenPolynomial<C>> A,
2257                    GenPolynomial<C> M, C mi, GenPolynomial<C> B, C am) {
2258        GenPolynomial<GenPolynomial<C>> S = fac.getZERO().copy();
2259        GenPolynomial<GenPolynomial<C>> Ap = A.copy();
2260        SortedMap<ExpVector, GenPolynomial<C>> av = Ap.val; //getMap();
2261        SortedMap<ExpVector, C> bv = B.getMap();
2262        SortedMap<ExpVector, GenPolynomial<C>> sv = S.val; //getMap();
2263        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) fac.coFac;
2264        RingFactory<C> bfac = cfac.coFac;
2265        GenPolynomial<C> c = null;
2266        for (Map.Entry<ExpVector, C> me : bv.entrySet()) {
2267            ExpVector e = me.getKey();
2268            C y = me.getValue(); //bv.get(e); // assert y != null
2269            GenPolynomial<C> x = av.get(e);
2270            if (x != null) {
2271                av.remove(e);
2272                c = PolyUtil.<C> interpolate(cfac, x, M, mi, y, am);
2273                if (!c.isZERO()) { // 0 cannot happen
2274                    sv.put(e, c);
2275                }
2276            } else {
2277                c = PolyUtil.<C> interpolate(cfac, cfac.getZERO(), M, mi, y, am);
2278                if (!c.isZERO()) { // 0 cannot happen
2279                    sv.put(e, c); // c != null
2280                }
2281            }
2282        }
2283        // assert bv is empty = done
2284        for (Map.Entry<ExpVector, GenPolynomial<C>> me : av.entrySet()) { // rest of av
2285            ExpVector e = me.getKey();
2286            GenPolynomial<C> x = me.getValue(); //av.get(e); // assert x != null
2287            c = PolyUtil.<C> interpolate(cfac, x, M, mi, bfac.getZERO(), am);
2288            if (!c.isZERO()) { // 0 cannot happen
2289                sv.put(e, c); // c != null
2290            }
2291        }
2292        return S;
2293    }
2294
2295
2296    /**
2297     * Univariate polynomial interpolation.
2298     * @param <C> coefficient type.
2299     * @param fac GenPolynomial<C> result factory.
2300     * @param A GenPolynomial.
2301     * @param M GenPolynomial interpolation modul of A.
2302     * @param mi inverse of M(am) in ring fac.coFac.
2303     * @param a evaluation of other GenPolynomial.
2304     * @param am evaluation point (interpolation modul) of a, i.e. P(am) = a.
2305     * @return S, with S mod M == A and S(am) == a.
2306     */
2307    public static <C extends RingElem<C>> GenPolynomial<C> interpolate(GenPolynomialRing<C> fac,
2308                    GenPolynomial<C> A, GenPolynomial<C> M, C mi, C a, C am) {
2309        GenPolynomial<C> s;
2310        C b = PolyUtil.<C> evaluateMain(fac.coFac, A, am);
2311        // A mod a.modul
2312        C d = a.subtract(b); // a-A mod a.modul
2313        if (d.isZERO()) {
2314            return A;
2315        }
2316        b = d.multiply(mi); // b = (a-A)*mi mod a.modul
2317        // (M*b)+A mod M = A mod M = 
2318        // (M*mi*(a-A)+A) mod a.modul = a mod a.modul
2319        s = M.multiply(b);
2320        s = s.sum(A);
2321        return s;
2322    }
2323
2324
2325    /**
2326     * Recursive GenPolynomial switch varaible blocks.
2327     * @param <C> coefficient type.
2328     * @param P recursive GenPolynomial in R[X,Y].
2329     * @return this in R[Y,X].
2330     */
2331    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> switchVariables(
2332                    GenPolynomial<GenPolynomial<C>> P) {
2333        if (P == null) {
2334            throw new IllegalArgumentException("P == null");
2335        }
2336        GenPolynomialRing<GenPolynomial<C>> rfac1 = P.ring;
2337        GenPolynomialRing<C> cfac1 = (GenPolynomialRing<C>) rfac1.coFac;
2338        GenPolynomialRing<C> cfac2 = new GenPolynomialRing<C>(cfac1.coFac, rfac1);
2339        GenPolynomial<C> zero = cfac2.getZERO();
2340        GenPolynomialRing<GenPolynomial<C>> rfac2 = new GenPolynomialRing<GenPolynomial<C>>(cfac2, cfac1);
2341        GenPolynomial<GenPolynomial<C>> B = rfac2.getZERO().copy();
2342        if (P.isZERO()) {
2343            return B;
2344        }
2345        for (Monomial<GenPolynomial<C>> mr : P) {
2346            GenPolynomial<C> cr = mr.c;
2347            for (Monomial<C> mc : cr) {
2348                GenPolynomial<C> c = zero.sum(mc.c, mr.e);
2349                B = B.sum(c, mc.e);
2350            }
2351        }
2352        return B;
2353    }
2354
2355
2356    /**
2357     * Maximal degree of leading terms of a polynomial list.
2358     * @return maximum degree of the leading terms of a polynomial list.
2359     */
2360    public static <C extends RingElem<C>> long totalDegreeLeadingTerm(List<GenPolynomial<C>> P) {
2361        long degree = 0;
2362        for (GenPolynomial<C> g : P) {
2363            long total = g.leadingExpVector().totalDeg();
2364            if (degree < total) {
2365                degree = total;
2366            }
2367        }
2368        return degree;
2369    }
2370
2371
2372    /**
2373     * Total degree of polynomial list.
2374     * @return total degree of the polynomial list.
2375     */
2376    public static <C extends RingElem<C>> long totalDegree(List<GenPolynomial<C>> P) {
2377        long degree = 0;
2378        for (GenPolynomial<C> g : P) {
2379            long total = g.totalDegree();
2380            if (degree < total) {
2381                degree = total;
2382            }
2383        }
2384        return degree;
2385    }
2386
2387
2388    /**
2389     * Maximal degree of polynomial list.
2390     * @return maximal degree of the polynomial list.
2391     */
2392    public static <C extends RingElem<C>> long maxDegree(List<GenPolynomial<C>> P) {
2393        long degree = 0;
2394        for (GenPolynomial<C> g : P) {
2395            long total = g.degree();
2396            if (degree < total) {
2397                degree = total;
2398            }
2399        }
2400        return degree;
2401    }
2402
2403
2404    /**
2405     * Maximal degree in the coefficient polynomials.
2406     * @param <C> coefficient type.
2407     * @return maximal degree in the coefficients.
2408     */
2409    public static <C extends RingElem<C>> long coeffMaxDegree(GenPolynomial<GenPolynomial<C>> A) {
2410        if (A.isZERO()) {
2411            return 0; // 0 or -1 ?;
2412        }
2413        long deg = 0;
2414        for (GenPolynomial<C> a : A.getMap().values()) {
2415            long d = a.degree();
2416            if (d > deg) {
2417                deg = d;
2418            }
2419        }
2420        return deg;
2421    }
2422
2423
2424    /**
2425     * Map a unary function to the coefficients.
2426     * @param ring result polynomial ring factory.
2427     * @param p polynomial.
2428     * @param f evaluation functor.
2429     * @return new polynomial with coefficients f(p(e)).
2430     */
2431    public static <C extends RingElem<C>, D extends RingElem<D>> GenPolynomial<D> map(
2432                    GenPolynomialRing<D> ring, GenPolynomial<C> p, UnaryFunctor<C, D> f) {
2433        GenPolynomial<D> n = ring.getZERO().copy();
2434        SortedMap<ExpVector, D> nv = n.val;
2435        for (Monomial<C> m : p) {
2436            D c = f.eval(m.c);
2437            if (c != null && !c.isZERO()) {
2438                nv.put(m.e, c);
2439            }
2440        }
2441        return n;
2442    }
2443
2444
2445    /**
2446     * Product representation.
2447     * @param <C> coefficient type.
2448     * @param pfac polynomial ring factory.
2449     * @param L list of polynomials to be represented.
2450     * @return Product represenation of L in the polynomial ring pfac.
2451     */
2452    public static <C extends GcdRingElem<C>> List<GenPolynomial<Product<C>>> toProductGen(
2453                    GenPolynomialRing<Product<C>> pfac, List<GenPolynomial<C>> L) {
2454
2455        List<GenPolynomial<Product<C>>> list = new ArrayList<GenPolynomial<Product<C>>>();
2456        if (L == null || L.size() == 0) {
2457            return list;
2458        }
2459        for (GenPolynomial<C> a : L) {
2460            GenPolynomial<Product<C>> b = toProductGen(pfac, a);
2461            list.add(b);
2462        }
2463        return list;
2464    }
2465
2466
2467    /**
2468     * Product representation.
2469     * @param <C> coefficient type.
2470     * @param pfac polynomial ring factory.
2471     * @param A polynomial to be represented.
2472     * @return Product represenation of A in the polynomial ring pfac.
2473     */
2474    public static <C extends GcdRingElem<C>> GenPolynomial<Product<C>> toProductGen(
2475                    GenPolynomialRing<Product<C>> pfac, GenPolynomial<C> A) {
2476
2477        GenPolynomial<Product<C>> P = pfac.getZERO().copy();
2478        if (A == null || A.isZERO()) {
2479            return P;
2480        }
2481        RingFactory<Product<C>> rpfac = pfac.coFac;
2482        ProductRing<C> rfac = (ProductRing<C>) rpfac;
2483        for (Map.Entry<ExpVector, C> y : A.getMap().entrySet()) {
2484            ExpVector e = y.getKey();
2485            C a = y.getValue();
2486            Product<C> p = toProductGen(rfac, a);
2487            if (!p.isZERO()) {
2488                P.doPutToMap(e, p);
2489            }
2490        }
2491        return P;
2492    }
2493
2494
2495    /**
2496     * Product representation.
2497     * @param <C> coefficient type.
2498     * @param pfac product ring factory.
2499     * @param c coefficient to be represented.
2500     * @return Product represenation of c in the ring pfac.
2501     */
2502    public static <C extends GcdRingElem<C>> Product<C> toProductGen(ProductRing<C> pfac, C c) {
2503
2504        SortedMap<Integer, C> elem = new TreeMap<Integer, C>();
2505        for (int i = 0; i < pfac.length(); i++) {
2506            RingFactory<C> rfac = pfac.getFactory(i);
2507            C u = rfac.copy(c);
2508            if (u != null && !u.isZERO()) {
2509                elem.put(i, u);
2510            }
2511        }
2512        return new Product<C>(pfac, elem);
2513    }
2514
2515
2516    /**
2517     * Product representation.
2518     * @param <C> coefficient type.
2519     * @param pfac product polynomial ring factory.
2520     * @param c coefficient to be used.
2521     * @param e exponent vector.
2522     * @return Product represenation of c X^e in the ring pfac.
2523     */
2524    public static <C extends RingElem<C>> Product<GenPolynomial<C>> toProduct(
2525                    ProductRing<GenPolynomial<C>> pfac, C c, ExpVector e) {
2526        SortedMap<Integer, GenPolynomial<C>> elem = new TreeMap<Integer, GenPolynomial<C>>();
2527        for (int i = 0; i < e.length(); i++) {
2528            RingFactory<GenPolynomial<C>> rfac = pfac.getFactory(i);
2529            GenPolynomialRing<C> fac = (GenPolynomialRing<C>) rfac;
2530            //GenPolynomialRing<C> cfac = fac.ring;
2531            long a = e.getVal(i);
2532            GenPolynomial<C> u;
2533            if (a == 0) {
2534                u = fac.getONE();
2535            } else {
2536                u = fac.univariate(0, a);
2537            }
2538            u = u.multiply(c);
2539            elem.put(i, u);
2540        }
2541        return new Product<GenPolynomial<C>>(pfac, elem);
2542    }
2543
2544
2545    /**
2546     * Product representation.
2547     * @param <C> coefficient type.
2548     * @param pfac product polynomial ring factory.
2549     * @param A polynomial.
2550     * @return Product represenation of the terms of A in the ring pfac.
2551     */
2552    public static <C extends RingElem<C>> Product<GenPolynomial<C>> toProduct(
2553                    ProductRing<GenPolynomial<C>> pfac, GenPolynomial<C> A) {
2554        Product<GenPolynomial<C>> P = pfac.getZERO();
2555        if (A == null || A.isZERO()) {
2556            return P;
2557        }
2558        for (Map.Entry<ExpVector, C> y : A.getMap().entrySet()) {
2559            ExpVector e = y.getKey();
2560            C a = y.getValue();
2561            Product<GenPolynomial<C>> p = toProduct(pfac, a, e);
2562            P = P.sum(p);
2563        }
2564        return P;
2565    }
2566
2567
2568    /**
2569     * Product representation.
2570     * @param pfac product ring factory.
2571     * @param c coefficient to be represented.
2572     * @return Product represenation of c in the ring pfac.
2573     */
2574    public static Product<ModInteger> toProduct(ProductRing<ModInteger> pfac, BigInteger c) {
2575
2576        SortedMap<Integer, ModInteger> elem = new TreeMap<Integer, ModInteger>();
2577        for (int i = 0; i < pfac.length(); i++) {
2578            RingFactory<ModInteger> rfac = pfac.getFactory(i);
2579            ModIntegerRing fac = (ModIntegerRing) rfac;
2580            ModInteger u = fac.fromInteger(c.getVal());
2581            if (!u.isZERO()) {
2582                elem.put(i, u);
2583            }
2584        }
2585        return new Product<ModInteger>(pfac, elem);
2586    }
2587
2588
2589    /**
2590     * Product representation.
2591     * @param pfac polynomial ring factory.
2592     * @param A polynomial to be represented.
2593     * @return Product represenation of A in the polynomial ring pfac.
2594     */
2595    public static GenPolynomial<Product<ModInteger>> toProduct(GenPolynomialRing<Product<ModInteger>> pfac,
2596                    GenPolynomial<BigInteger> A) {
2597
2598        GenPolynomial<Product<ModInteger>> P = pfac.getZERO().copy();
2599        if (A == null || A.isZERO()) {
2600            return P;
2601        }
2602        RingFactory<Product<ModInteger>> rpfac = pfac.coFac;
2603        ProductRing<ModInteger> fac = (ProductRing<ModInteger>) rpfac;
2604        for (Map.Entry<ExpVector, BigInteger> y : A.getMap().entrySet()) {
2605            ExpVector e = y.getKey();
2606            BigInteger a = y.getValue();
2607            Product<ModInteger> p = toProduct(fac, a);
2608            if (!p.isZERO()) {
2609                P.doPutToMap(e, p);
2610            }
2611        }
2612        return P;
2613    }
2614
2615
2616    /**
2617     * Product representation.
2618     * @param pfac polynomial ring factory.
2619     * @param L list of polynomials to be represented.
2620     * @return Product represenation of L in the polynomial ring pfac.
2621     */
2622    public static List<GenPolynomial<Product<ModInteger>>> toProduct(
2623                    GenPolynomialRing<Product<ModInteger>> pfac, List<GenPolynomial<BigInteger>> L) {
2624
2625        List<GenPolynomial<Product<ModInteger>>> list = new ArrayList<GenPolynomial<Product<ModInteger>>>();
2626        if (L == null || L.size() == 0) {
2627            return list;
2628        }
2629        for (GenPolynomial<BigInteger> a : L) {
2630            GenPolynomial<Product<ModInteger>> b = toProduct(pfac, a);
2631            list.add(b);
2632        }
2633        return list;
2634    }
2635
2636
2637    /**
2638     * Intersection. Intersection of a list of polynomials with a polynomial
2639     * ring. The polynomial ring must be a contraction of the polynomial ring of
2640     * the list of polynomials and the TermOrder must be an elimination order.
2641     * @param R polynomial ring
2642     * @param F list of polynomials
2643     * @return R \cap F
2644     */
2645    public static <C extends RingElem<C>> List<GenPolynomial<C>> intersect(GenPolynomialRing<C> R,
2646                    List<GenPolynomial<C>> F) {
2647        if (F == null || F.isEmpty()) {
2648            return F;
2649        }
2650        GenPolynomialRing<C> pfac = F.get(0).ring;
2651        int d = pfac.nvar - R.nvar;
2652        if (d <= 0) {
2653            return F;
2654        }
2655        List<GenPolynomial<C>> H = new ArrayList<GenPolynomial<C>>(F.size());
2656        for (GenPolynomial<C> p : F) {
2657            Map<ExpVector, GenPolynomial<C>> m = null;
2658            m = p.contract(R);
2659            if (logger.isDebugEnabled()) {
2660                logger.debug("intersect contract m = " + m);
2661            }
2662            if (m.size() == 1) { // contains one power of variables
2663                for (Map.Entry<ExpVector, GenPolynomial<C>> me : m.entrySet()) {
2664                    ExpVector e = me.getKey();
2665                    if (e.isZERO()) {
2666                        H.add(me.getValue());
2667                    }
2668                }
2669            }
2670        }
2671        GenPolynomialRing<C> tfac = pfac.contract(d);
2672        if (tfac.equals(R)) { // check 
2673            return H;
2674        }
2675        logger.warn("tfac != R: tfac = " + tfac.toScript() + ", R = " + R.toScript() + ", pfac = "
2676                        + pfac.toScript());
2677        // throw new RuntimeException("contract(pfac) != R");
2678        return H;
2679    }
2680
2681
2682    /**
2683     * Intersection. Intersection of a list of solvable polynomials with a
2684     * solvable polynomial ring. The solvable polynomial ring must be a
2685     * contraction of the solvable polynomial ring of the list of polynomials
2686     * and the TermOrder must be an elimination order.
2687     * @param R solvable polynomial ring
2688     * @param F list of solvable polynomials
2689     * @return R \cap F
2690     */
2691    @SuppressWarnings("cast")
2692    public static <C extends RingElem<C>> List<GenSolvablePolynomial<C>> intersect(
2693                    GenSolvablePolynomialRing<C> R, List<GenSolvablePolynomial<C>> F) {
2694        List<GenPolynomial<C>> Fp = PolynomialList.<C> castToList(F);
2695        GenPolynomialRing<C> Rp = (GenPolynomialRing<C>) R;
2696        List<GenPolynomial<C>> H = intersect(Rp, Fp);
2697        return PolynomialList.<C> castToSolvableList(H);
2698    }
2699
2700
2701    /**
2702     * Intersection. Intersection of a list of word polynomials with a word
2703     * polynomial ring. The polynomial ring must be a contraction of the
2704     * polynomial ring of the list of polynomials,
2705     * @param R word polynomial ring
2706     * @param F list of word polynomials
2707     * @return R \cap F
2708     */
2709    public static <C extends RingElem<C>> List<GenWordPolynomial<C>> intersect(GenWordPolynomialRing<C> R,
2710                    List<GenWordPolynomial<C>> F) {
2711        if (F == null || F.isEmpty()) {
2712            return F;
2713        }
2714        GenWordPolynomialRing<C> pfac = F.get(0).ring;
2715        assert pfac.alphabet.isSubFactory(R.alphabet) : "pfac=" + pfac.alphabet + ", R=" + R.alphabet;
2716        List<GenWordPolynomial<C>> H = new ArrayList<GenWordPolynomial<C>>(F.size());
2717        for (GenWordPolynomial<C> p : F) {
2718            if (p == null || p.isZERO()) {
2719                continue;
2720            }
2721            GenWordPolynomial<C> m = p.contract(R);
2722            if (logger.isDebugEnabled()) {
2723                logger.debug("intersect contract m = " + m);
2724            }
2725            if (!m.isZERO()) {
2726                H.add(m);
2727            }
2728        }
2729        // throw new RuntimeException("contract(pfac) != R");
2730        return H;
2731    }
2732
2733
2734    /**
2735     * Remove all upper variables which do not occur in polynomial.
2736     * @param p polynomial.
2737     * @return polynomial with removed variables
2738     */
2739    public static <C extends RingElem<C>> GenPolynomial<C> removeUnusedUpperVariables(GenPolynomial<C> p) {
2740        GenPolynomialRing<C> fac = p.ring;
2741        if (fac.nvar <= 1) { // univariate
2742            return p;
2743        }
2744        int[] dep = p.degreeVector().dependencyOnVariables();
2745        if (fac.nvar == dep.length) { // all variables appear
2746            return p;
2747        }
2748        if (dep.length == 0) { // no variables
2749            GenPolynomialRing<C> fac0 = new GenPolynomialRing<C>(fac.coFac, 0);
2750            GenPolynomial<C> p0 = new GenPolynomial<C>(fac0, p.leadingBaseCoefficient());
2751            return p0;
2752        }
2753        int l = dep[0]; // higher variable
2754        int r = dep[dep.length - 1]; // lower variable
2755        if (l == 0 /*|| l == fac.nvar-1*/) { // upper variable appears
2756            return p;
2757        }
2758        int n = l;
2759        GenPolynomialRing<C> facr = fac.contract(n);
2760        Map<ExpVector, GenPolynomial<C>> mpr = p.contract(facr);
2761        if (mpr.size() != 1) {
2762            logger.warn("upper ex, l = " + l + ", r = " + r + ", p = " + p + ", fac = " + fac.toScript());
2763            throw new RuntimeException("this should not happen " + mpr);
2764        }
2765        GenPolynomial<C> pr = mpr.values().iterator().next();
2766        n = fac.nvar - 1 - r;
2767        if (n == 0) {
2768            return pr;
2769        } // else case not implemented
2770        return pr;
2771    }
2772
2773
2774    /**
2775     * Remove all lower variables which do not occur in polynomial.
2776     * @param p polynomial.
2777     * @return polynomial with removed variables
2778     */
2779    public static <C extends RingElem<C>> GenPolynomial<C> removeUnusedLowerVariables(GenPolynomial<C> p) {
2780        GenPolynomialRing<C> fac = p.ring;
2781        if (fac.nvar <= 1) { // univariate
2782            return p;
2783        }
2784        int[] dep = p.degreeVector().dependencyOnVariables();
2785        if (fac.nvar == dep.length) { // all variables appear
2786            return p;
2787        }
2788        if (dep.length == 0) { // no variables
2789            GenPolynomialRing<C> fac0 = new GenPolynomialRing<C>(fac.coFac, 0);
2790            GenPolynomial<C> p0 = new GenPolynomial<C>(fac0, p.leadingBaseCoefficient());
2791            return p0;
2792        }
2793        int l = dep[0]; // higher variable
2794        int r = dep[dep.length - 1]; // lower variable
2795        if (r == fac.nvar - 1) { // lower variable appears
2796            return p;
2797        }
2798        int n = r + 1;
2799        GenPolynomialRing<GenPolynomial<C>> rfac = fac.recursive(n);
2800        GenPolynomial<GenPolynomial<C>> mpr = recursive(rfac, p);
2801        if (mpr.length() != p.length()) {
2802            logger.warn("lower ex, l = " + l + ", r = " + r + ", p = " + p + ", fac = " + fac.toScript());
2803            throw new RuntimeException("this should not happen " + mpr);
2804        }
2805        RingFactory<C> cf = fac.coFac;
2806        GenPolynomialRing<C> facl = new GenPolynomialRing<C>(cf, rfac);
2807        GenPolynomial<C> pr = facl.getZERO().copy();
2808        for (Monomial<GenPolynomial<C>> m : mpr) {
2809            ExpVector e = m.e;
2810            GenPolynomial<C> a = m.c;
2811            if (!a.isConstant()) {
2812                throw new RuntimeException("this can not happen " + a);
2813            }
2814            C c = a.leadingBaseCoefficient();
2815            pr.doPutToMap(e, c);
2816        }
2817        return pr;
2818    }
2819
2820
2821    /**
2822     * Remove upper block of middle variables which do not occur in polynomial.
2823     * @param p polynomial.
2824     * @return polynomial with removed variables
2825     */
2826    public static <C extends RingElem<C>> GenPolynomial<C> removeUnusedMiddleVariables(GenPolynomial<C> p) {
2827        GenPolynomialRing<C> fac = p.ring;
2828        if (fac.nvar <= 2) { // univariate or bi-variate
2829            return p;
2830        }
2831        int[] dep = p.degreeVector().dependencyOnVariables();
2832        if (fac.nvar == dep.length) { // all variables appear
2833            return p;
2834        }
2835        if (dep.length == 0) { // no variables
2836            GenPolynomialRing<C> fac0 = new GenPolynomialRing<C>(fac.coFac, 0);
2837            GenPolynomial<C> p0 = new GenPolynomial<C>(fac0, p.leadingBaseCoefficient());
2838            return p0;
2839        }
2840        ExpVector e1 = p.leadingExpVector();
2841        if (dep.length == 1) { // one variable
2842            TermOrder to = new TermOrder(fac.tord.getEvord());
2843            int i = dep[0];
2844            String v1 = e1.indexVarName(i, fac.getVars());
2845            String[] vars = new String[] { v1 };
2846            GenPolynomialRing<C> fac1 = new GenPolynomialRing<C>(fac.coFac, to, vars);
2847            GenPolynomial<C> p1 = fac1.getZERO().copy();
2848            for (Monomial<C> m : p) {
2849                ExpVector e = m.e;
2850                ExpVector f = ExpVector.create(1, 0, e.getVal(i));
2851                p1.doPutToMap(f, m.c);
2852            }
2853            return p1;
2854        }
2855        GenPolynomialRing<GenPolynomial<C>> rfac = fac.recursive(1);
2856        GenPolynomial<GenPolynomial<C>> mpr = recursive(rfac, p);
2857
2858        int l = dep[0]; // higher variable
2859        int r = fac.nvar - dep[1]; // next variable
2860        //System.out.println("l  = " + l);
2861        //System.out.println("r  = " + r);
2862
2863        TermOrder to = new TermOrder(fac.tord.getEvord());
2864        String[] vs = fac.getVars();
2865        String[] vars = new String[r + 1];
2866        for (int i = 0; i < r; i++) {
2867            vars[i] = vs[i];
2868        }
2869        vars[r] = e1.indexVarName(l, vs);
2870        //System.out.println("fac  = " + fac);
2871        GenPolynomialRing<C> dfac = new GenPolynomialRing<C>(fac.coFac, to, vars);
2872        //System.out.println("dfac = " + dfac);
2873        GenPolynomialRing<GenPolynomial<C>> fac2 = dfac.recursive(1);
2874        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) fac2.coFac;
2875        GenPolynomial<GenPolynomial<C>> p2r = fac2.getZERO().copy();
2876        for (Monomial<GenPolynomial<C>> m : mpr) {
2877            ExpVector e = m.e;
2878            GenPolynomial<C> a = m.c;
2879            Map<ExpVector, GenPolynomial<C>> cc = a.contract(cfac);
2880            for (Map.Entry<ExpVector, GenPolynomial<C>> me : cc.entrySet()) {
2881                ExpVector f = me.getKey();
2882                if (f.isZERO()) {
2883                    GenPolynomial<C> c = me.getValue(); //cc.get(f);
2884                    p2r.doPutToMap(e, c);
2885                } else {
2886                    throw new RuntimeException("this should not happen " + cc);
2887                }
2888            }
2889        }
2890        GenPolynomial<C> p2 = distribute(dfac, p2r);
2891        return p2;
2892    }
2893
2894
2895    /**
2896     * Select polynomial with univariate leading term in variable i.
2897     * @param i variable index.
2898     * @return polynomial with head term in variable i
2899     */
2900    public static <C extends RingElem<C>> GenPolynomial<C> selectWithVariable(List<GenPolynomial<C>> P,
2901                    int i) {
2902        for (GenPolynomial<C> p : P) {
2903            int[] dep = p.leadingExpVector().dependencyOnVariables();
2904            if (dep.length == 1 && dep[0] == i) {
2905                return p;
2906            }
2907        }
2908        return null; // not found       
2909    }
2910
2911}
2912
2913
2914/**
2915 * Conversion of distributive to recursive representation.
2916 */
2917class DistToRec<C extends RingElem<C>>
2918                implements UnaryFunctor<GenPolynomial<C>, GenPolynomial<GenPolynomial<C>>> {
2919
2920
2921    GenPolynomialRing<GenPolynomial<C>> fac;
2922
2923
2924    public DistToRec(GenPolynomialRing<GenPolynomial<C>> fac) {
2925        this.fac = fac;
2926    }
2927
2928
2929    public GenPolynomial<GenPolynomial<C>> eval(GenPolynomial<C> c) {
2930        if (c == null) {
2931            return fac.getZERO();
2932        }
2933        return PolyUtil.<C> recursive(fac, c);
2934    }
2935}
2936
2937
2938/**
2939 * Conversion of recursive to distributive representation.
2940 */
2941class RecToDist<C extends RingElem<C>>
2942                implements UnaryFunctor<GenPolynomial<GenPolynomial<C>>, GenPolynomial<C>> {
2943
2944
2945    GenPolynomialRing<C> fac;
2946
2947
2948    public RecToDist(GenPolynomialRing<C> fac) {
2949        this.fac = fac;
2950    }
2951
2952
2953    public GenPolynomial<C> eval(GenPolynomial<GenPolynomial<C>> c) {
2954        if (c == null) {
2955            return fac.getZERO();
2956        }
2957        return PolyUtil.<C> distribute(fac, c);
2958    }
2959}
2960
2961
2962/**
2963 * BigRational numerator functor.
2964 */
2965class RatNumer implements UnaryFunctor<BigRational, BigInteger> {
2966
2967
2968    public BigInteger eval(BigRational c) {
2969        if (c == null) {
2970            return new BigInteger();
2971        }
2972        return new BigInteger(c.numerator());
2973    }
2974}
2975
2976
2977/**
2978 * Conversion of symmetric ModInteger to BigInteger functor.
2979 */
2980class ModSymToInt<C extends RingElem<C> & Modular> implements UnaryFunctor<C, BigInteger> {
2981
2982
2983    public BigInteger eval(C c) {
2984        if (c == null) {
2985            return new BigInteger();
2986        }
2987        return c.getSymmetricInteger();
2988    }
2989}
2990
2991
2992/**
2993 * Conversion of ModInteger to BigInteger functor.
2994 */
2995class ModToInt<C extends RingElem<C> & Modular> implements UnaryFunctor<C, BigInteger> {
2996
2997
2998    public BigInteger eval(C c) {
2999        if (c == null) {
3000            return new BigInteger();
3001        }
3002        return c.getInteger();
3003    }
3004}
3005
3006
3007/**
3008 * Conversion of BigRational to BigInteger with division by lcm functor. result
3009 * = num*(lcm/denom).
3010 */
3011class RatToInt implements UnaryFunctor<BigRational, BigInteger> {
3012
3013
3014    java.math.BigInteger lcm;
3015
3016
3017    public RatToInt(java.math.BigInteger lcm) {
3018        this.lcm = lcm; //.getVal();
3019    }
3020
3021
3022    public BigInteger eval(BigRational c) {
3023        if (c == null) {
3024            return new BigInteger();
3025        }
3026        // p = num*(lcm/denom)
3027        java.math.BigInteger b = lcm.divide(c.denominator());
3028        return new BigInteger(c.numerator().multiply(b));
3029    }
3030}
3031
3032
3033/**
3034 * Conversion of BigRational to BigInteger. result = (num/gcd)*(lcm/denom).
3035 */
3036class RatToIntFactor implements UnaryFunctor<BigRational, BigInteger> {
3037
3038
3039    final java.math.BigInteger lcm;
3040
3041
3042    final java.math.BigInteger gcd;
3043
3044
3045    public RatToIntFactor(java.math.BigInteger gcd, java.math.BigInteger lcm) {
3046        this.gcd = gcd;
3047        this.lcm = lcm; // .getVal();
3048    }
3049
3050
3051    public BigInteger eval(BigRational c) {
3052        if (c == null) {
3053            return new BigInteger();
3054        }
3055        if (gcd.equals(java.math.BigInteger.ONE)) {
3056            // p = num*(lcm/denom)
3057            java.math.BigInteger b = lcm.divide(c.denominator());
3058            return new BigInteger(c.numerator().multiply(b));
3059        }
3060        // p = (num/gcd)*(lcm/denom)
3061        java.math.BigInteger a = c.numerator().divide(gcd);
3062        java.math.BigInteger b = lcm.divide(c.denominator());
3063        return new BigInteger(a.multiply(b));
3064    }
3065}
3066
3067
3068/**
3069 * Conversion of Rational to BigDecimal. result = decimal(r).
3070 */
3071class RatToDec<C extends Element<C> & Rational> implements UnaryFunctor<C, BigDecimal> {
3072
3073
3074    public BigDecimal eval(C c) {
3075        if (c == null) {
3076            return new BigDecimal();
3077        }
3078        return new BigDecimal(c.getRational());
3079    }
3080}
3081
3082
3083/**
3084 * Conversion of Complex Rational to Complex BigDecimal. result = decimal(r).
3085 */
3086class CompRatToDec<C extends RingElem<C> & Rational>
3087                implements UnaryFunctor<Complex<C>, Complex<BigDecimal>> {
3088
3089
3090    ComplexRing<BigDecimal> ring;
3091
3092
3093    public CompRatToDec(RingFactory<Complex<BigDecimal>> ring) {
3094        this.ring = (ComplexRing<BigDecimal>) ring;
3095    }
3096
3097
3098    public Complex<BigDecimal> eval(Complex<C> c) {
3099        if (c == null) {
3100            return ring.getZERO();
3101        }
3102        BigDecimal r = new BigDecimal(c.getRe().getRational());
3103        BigDecimal i = new BigDecimal(c.getIm().getRational());
3104        return new Complex<BigDecimal>(ring, r, i);
3105    }
3106}
3107
3108
3109/**
3110 * Conversion from BigInteger functor.
3111 */
3112class FromInteger<D extends RingElem<D>> implements UnaryFunctor<BigInteger, D> {
3113
3114
3115    RingFactory<D> ring;
3116
3117
3118    public FromInteger(RingFactory<D> ring) {
3119        this.ring = ring;
3120    }
3121
3122
3123    public D eval(BigInteger c) {
3124        if (c == null) {
3125            return ring.getZERO();
3126        }
3127        return ring.fromInteger(c.getVal());
3128    }
3129}
3130
3131
3132/**
3133 * Conversion from GenPolynomial<BigInteger> functor.
3134 */
3135class FromIntegerPoly<D extends RingElem<D>>
3136                implements UnaryFunctor<GenPolynomial<BigInteger>, GenPolynomial<D>> {
3137
3138
3139    GenPolynomialRing<D> ring;
3140
3141
3142    FromInteger<D> fi;
3143
3144
3145    public FromIntegerPoly(GenPolynomialRing<D> ring) {
3146        if (ring == null) {
3147            throw new IllegalArgumentException("ring must not be null");
3148        }
3149        this.ring = ring;
3150        fi = new FromInteger<D>(ring.coFac);
3151    }
3152
3153
3154    public GenPolynomial<D> eval(GenPolynomial<BigInteger> c) {
3155        if (c == null) {
3156            return ring.getZERO();
3157        }
3158        return PolyUtil.<BigInteger, D> map(ring, c, fi);
3159    }
3160}
3161
3162
3163/**
3164 * Conversion from GenPolynomial<BigRational> to GenPolynomial <BigInteger>
3165 * functor.
3166 */
3167class RatToIntPoly implements UnaryFunctor<GenPolynomial<BigRational>, GenPolynomial<BigInteger>> {
3168
3169
3170    GenPolynomialRing<BigInteger> ring;
3171
3172
3173    public RatToIntPoly(GenPolynomialRing<BigInteger> ring) {
3174        if (ring == null) {
3175            throw new IllegalArgumentException("ring must not be null");
3176        }
3177        this.ring = ring;
3178    }
3179
3180
3181    public GenPolynomial<BigInteger> eval(GenPolynomial<BigRational> c) {
3182        if (c == null) {
3183            return ring.getZERO();
3184        }
3185        return PolyUtil.integerFromRationalCoefficients(ring, c);
3186    }
3187}
3188
3189
3190/**
3191 * Real part functor.
3192 */
3193class RealPart implements UnaryFunctor<BigComplex, BigRational> {
3194
3195
3196    public BigRational eval(BigComplex c) {
3197        if (c == null) {
3198            return new BigRational();
3199        }
3200        return c.getRe();
3201    }
3202}
3203
3204
3205/**
3206 * Imaginary part functor.
3207 */
3208class ImagPart implements UnaryFunctor<BigComplex, BigRational> {
3209
3210
3211    public BigRational eval(BigComplex c) {
3212        if (c == null) {
3213            return new BigRational();
3214        }
3215        return c.getIm();
3216    }
3217}
3218
3219
3220/**
3221 * Real part functor.
3222 */
3223class RealPartComplex<C extends RingElem<C>> implements UnaryFunctor<Complex<C>, C> {
3224
3225
3226    public C eval(Complex<C> c) {
3227        if (c == null) {
3228            return null;
3229        }
3230        return c.getRe();
3231    }
3232}
3233
3234
3235/**
3236 * Imaginary part functor.
3237 */
3238class ImagPartComplex<C extends RingElem<C>> implements UnaryFunctor<Complex<C>, C> {
3239
3240
3241    public C eval(Complex<C> c) {
3242        if (c == null) {
3243            return null;
3244        }
3245        return c.getIm();
3246    }
3247}
3248
3249
3250/**
3251 * Rational to complex functor.
3252 */
3253class ToComplex<C extends RingElem<C>> implements UnaryFunctor<C, Complex<C>> {
3254
3255
3256    final protected ComplexRing<C> cfac;
3257
3258
3259    @SuppressWarnings("unchecked")
3260    public ToComplex(RingFactory<Complex<C>> fac) {
3261        if (fac == null) {
3262            throw new IllegalArgumentException("fac must not be null");
3263        }
3264        cfac = (ComplexRing<C>) fac;
3265    }
3266
3267
3268    public Complex<C> eval(C c) {
3269        if (c == null) {
3270            return cfac.getZERO();
3271        }
3272        return new Complex<C>(cfac, c);
3273    }
3274}
3275
3276
3277/**
3278 * Rational to complex functor.
3279 */
3280class RatToCompl implements UnaryFunctor<BigRational, BigComplex> {
3281
3282
3283    public BigComplex eval(BigRational c) {
3284        if (c == null) {
3285            return new BigComplex();
3286        }
3287        return new BigComplex(c);
3288    }
3289}
3290
3291
3292/**
3293 * Any ring element to generic complex functor.
3294 */
3295class AnyToComplex<C extends GcdRingElem<C>> implements UnaryFunctor<C, Complex<C>> {
3296
3297
3298    final protected ComplexRing<C> cfac;
3299
3300
3301    public AnyToComplex(ComplexRing<C> fac) {
3302        if (fac == null) {
3303            throw new IllegalArgumentException("fac must not be null");
3304        }
3305        cfac = fac;
3306    }
3307
3308
3309    public AnyToComplex(RingFactory<C> fac) {
3310        this(new ComplexRing<C>(fac));
3311    }
3312
3313
3314    public Complex<C> eval(C a) {
3315        if (a == null || a.isZERO()) { // should not happen
3316            return cfac.getZERO();
3317        } else if (a.isONE()) {
3318            return cfac.getONE();
3319        } else {
3320            return new Complex<C>(cfac, a);
3321        }
3322    }
3323}
3324
3325
3326/**
3327 * Algebraic to generic complex functor.
3328 */
3329class AlgebToCompl<C extends GcdRingElem<C>> implements UnaryFunctor<AlgebraicNumber<C>, Complex<C>> {
3330
3331
3332    final protected ComplexRing<C> cfac;
3333
3334
3335    public AlgebToCompl(ComplexRing<C> fac) {
3336        if (fac == null) {
3337            throw new IllegalArgumentException("fac must not be null");
3338        }
3339        cfac = fac;
3340    }
3341
3342
3343    public Complex<C> eval(AlgebraicNumber<C> a) {
3344        if (a == null || a.isZERO()) { // should not happen
3345            return cfac.getZERO();
3346        } else if (a.isONE()) {
3347            return cfac.getONE();
3348        } else {
3349            GenPolynomial<C> p = a.getVal();
3350            C real = cfac.ring.getZERO();
3351            C imag = cfac.ring.getZERO();
3352            for (Monomial<C> m : p) {
3353                if (m.exponent().getVal(0) == 1L) {
3354                    imag = m.coefficient();
3355                } else if (m.exponent().getVal(0) == 0L) {
3356                    real = m.coefficient();
3357                } else {
3358                    throw new IllegalArgumentException("unexpected monomial " + m);
3359                }
3360            }
3361            //Complex<C> c = new Complex<C>(cfac,real,imag);
3362            return new Complex<C>(cfac, real, imag);
3363        }
3364    }
3365}
3366
3367
3368/**
3369 * Ceneric complex to algebraic number functor.
3370 */
3371class ComplToAlgeb<C extends GcdRingElem<C>> implements UnaryFunctor<Complex<C>, AlgebraicNumber<C>> {
3372
3373
3374    final protected AlgebraicNumberRing<C> afac;
3375
3376
3377    final protected AlgebraicNumber<C> I;
3378
3379
3380    public ComplToAlgeb(AlgebraicNumberRing<C> fac) {
3381        if (fac == null) {
3382            throw new IllegalArgumentException("fac must not be null");
3383        }
3384        afac = fac;
3385        I = afac.getGenerator();
3386    }
3387
3388
3389    public AlgebraicNumber<C> eval(Complex<C> c) {
3390        if (c == null || c.isZERO()) { // should not happen
3391            return afac.getZERO();
3392        } else if (c.isONE()) {
3393            return afac.getONE();
3394        } else if (c.isIMAG()) {
3395            return I;
3396        } else {
3397            return I.multiply(c.getIm()).sum(c.getRe());
3398        }
3399    }
3400}
3401
3402
3403/**
3404 * Algebraic to polynomial functor.
3405 */
3406class AlgToPoly<C extends GcdRingElem<C>> implements UnaryFunctor<AlgebraicNumber<C>, GenPolynomial<C>> {
3407
3408
3409    public GenPolynomial<C> eval(AlgebraicNumber<C> c) {
3410        if (c == null) {
3411            return null;
3412        }
3413        return c.val;
3414    }
3415}
3416
3417
3418/**
3419 * Polynomial to algebriac functor.
3420 */
3421class PolyToAlg<C extends GcdRingElem<C>> implements UnaryFunctor<GenPolynomial<C>, AlgebraicNumber<C>> {
3422
3423
3424    final protected AlgebraicNumberRing<C> afac;
3425
3426
3427    public PolyToAlg(AlgebraicNumberRing<C> fac) {
3428        if (fac == null) {
3429            throw new IllegalArgumentException("fac must not be null");
3430        }
3431        afac = fac;
3432    }
3433
3434
3435    public AlgebraicNumber<C> eval(GenPolynomial<C> c) {
3436        if (c == null) {
3437            return afac.getZERO();
3438        }
3439        return new AlgebraicNumber<C>(afac, c);
3440    }
3441}
3442
3443
3444/**
3445 * Coefficient to algebriac functor.
3446 */
3447class CoeffToAlg<C extends GcdRingElem<C>> implements UnaryFunctor<C, AlgebraicNumber<C>> {
3448
3449
3450    final protected AlgebraicNumberRing<C> afac;
3451
3452
3453    final protected GenPolynomial<C> zero;
3454
3455
3456    public CoeffToAlg(AlgebraicNumberRing<C> fac) {
3457        if (fac == null) {
3458            throw new IllegalArgumentException("fac must not be null");
3459        }
3460        afac = fac;
3461        GenPolynomialRing<C> pfac = afac.ring;
3462        zero = pfac.getZERO();
3463    }
3464
3465
3466    public AlgebraicNumber<C> eval(C c) {
3467        if (c == null) {
3468            return afac.getZERO();
3469        }
3470        return new AlgebraicNumber<C>(afac, zero.sum(c));
3471    }
3472}
3473
3474
3475/**
3476 * Coefficient to recursive algebriac functor.
3477 */
3478class CoeffToRecAlg<C extends GcdRingElem<C>> implements UnaryFunctor<C, AlgebraicNumber<C>> {
3479
3480
3481    final protected List<AlgebraicNumberRing<C>> lfac;
3482
3483
3484    final int depth;
3485
3486
3487    @SuppressWarnings({ "unchecked", "cast" })
3488    public CoeffToRecAlg(int depth, AlgebraicNumberRing<C> fac) {
3489        if (fac == null) {
3490            throw new IllegalArgumentException("fac must not be null");
3491        }
3492        AlgebraicNumberRing<C> afac = fac;
3493        this.depth = depth;
3494        lfac = new ArrayList<AlgebraicNumberRing<C>>(this.depth);
3495        lfac.add(fac);
3496        for (int i = 1; i < this.depth; i++) {
3497            RingFactory<C> rf = afac.ring.coFac;
3498            if (!(rf instanceof AlgebraicNumberRing)) {
3499                throw new IllegalArgumentException("fac depth to low");
3500            }
3501            afac = (AlgebraicNumberRing<C>) (Object) rf;
3502            lfac.add(afac);
3503        }
3504    }
3505
3506
3507    @SuppressWarnings({ "unchecked" })
3508    public AlgebraicNumber<C> eval(C c) {
3509        if (c == null) {
3510            return lfac.get(0).getZERO();
3511        }
3512        C ac = c;
3513        AlgebraicNumberRing<C> af = lfac.get(lfac.size() - 1);
3514        GenPolynomial<C> zero = af.ring.getZERO();
3515        AlgebraicNumber<C> an = new AlgebraicNumber<C>(af, zero.sum(ac));
3516        for (int i = lfac.size() - 2; i >= 0; i--) {
3517            af = lfac.get(i);
3518            zero = af.ring.getZERO();
3519            ac = (C) (Object) an;
3520            an = new AlgebraicNumber<C>(af, zero.sum(ac));
3521        }
3522        return an;
3523    }
3524}
3525
3526
3527/**
3528 * Evaluate main variable functor.
3529 */
3530class EvalMain<C extends RingElem<C>> implements UnaryFunctor<GenPolynomial<C>, C> {
3531
3532
3533    final RingFactory<C> cfac;
3534
3535
3536    final C a;
3537
3538
3539    public EvalMain(RingFactory<C> cfac, C a) {
3540        this.cfac = cfac;
3541        this.a = a;
3542    }
3543
3544
3545    public C eval(GenPolynomial<C> c) {
3546        if (c == null) {
3547            return cfac.getZERO();
3548        }
3549        return PolyUtil.<C> evaluateMain(cfac, c, a);
3550    }
3551}
3552
3553
3554/**
3555 * Evaluate main variable functor.
3556 */
3557class EvalMainPol<C extends RingElem<C>> implements UnaryFunctor<GenPolynomial<C>, GenPolynomial<C>> {
3558
3559
3560    final GenPolynomialRing<C> cfac;
3561
3562
3563    final C a;
3564
3565
3566    public EvalMainPol(GenPolynomialRing<C> cfac, C a) {
3567        this.cfac = cfac;
3568        this.a = a;
3569    }
3570
3571
3572    public GenPolynomial<C> eval(GenPolynomial<C> c) {
3573        if (c == null) {
3574            return cfac.getZERO();
3575        }
3576        return PolyUtil.<C> evaluateMain(cfac, c, a);
3577    }
3578}