001/*
002 * $Id: GenPolynomial.java 4291 2012-11-04 18:16:27Z kredel $
003 */
004
005package edu.jas.poly;
006
007
008import java.util.Collections;
009import java.util.Iterator;
010import java.util.Map;
011import java.util.SortedMap;
012import java.util.TreeMap;
013
014import org.apache.log4j.Logger;
015
016import edu.jas.kern.PreemptingException;
017import edu.jas.kern.PrettyPrint;
018import edu.jas.structure.NotInvertibleException;
019import edu.jas.structure.RingElem;
020import edu.jas.structure.UnaryFunctor;
021
022
023/**
024 * GenPolynomial generic polynomials implementing RingElem. n-variate ordered
025 * polynomials over C. Objects of this class are intended to be immutable. The
026 * implementation is based on TreeMap respectively SortedMap from exponents to
027 * coefficients. Only the coefficients are modeled with generic types, the
028 * exponents are fixed to ExpVector with long entries (this will eventually be
029 * changed in the future). C can also be a non integral domain, e.g. a
030 * ModInteger, i.e. it may contain zero divisors, since multiply() does now
031 * check for zeros. <b>Note:</b> multiply() now checks for wrong method dispatch
032 * for GenSolvablePolynomial.
033 * @param <C> coefficient type
034 * @author Heinz Kredel
035 */
036public class GenPolynomial<C extends RingElem<C>> implements RingElem<GenPolynomial<C>>, /* not yet Polynomial<C> */
037Iterable<Monomial<C>> {
038
039
040    /**
041     * The factory for the polynomial ring.
042     */
043    public final GenPolynomialRing<C> ring;
044
045
046    /**
047     * The data structure for polynomials.
048     */
049    protected final SortedMap<ExpVector, C> val; // do not change to TreeMap
050
051
052    private static final Logger logger = Logger.getLogger(GenPolynomial.class);
053
054
055    private final boolean debug = logger.isDebugEnabled();
056
057
058    // protected GenPolynomial() { ring = null; val = null; } // don't use
059
060
061    /**
062     * Private constructor for GenPolynomial.
063     * @param r polynomial ring factory.
064     * @param t TreeMap with correct ordering.
065     */
066    private GenPolynomial(GenPolynomialRing<C> r, TreeMap<ExpVector, C> t) {
067        ring = r;
068        val = t;
069        if (ring.checkPreempt) {
070            if (Thread.currentThread().isInterrupted()) {
071                logger.debug("throw PreemptingException");
072                throw new PreemptingException();
073            }
074        }
075    }
076
077
078    /**
079     * Constructor for zero GenPolynomial.
080     * @param r polynomial ring factory.
081     */
082    public GenPolynomial(GenPolynomialRing<C> r) {
083        this(r, new TreeMap<ExpVector, C>(r.tord.getDescendComparator()));
084    }
085
086
087    /**
088     * Constructor for GenPolynomial c * x<sup>e</sup>.
089     * @param r polynomial ring factory.
090     * @param c coefficient.
091     * @param e exponent.
092     */
093    public GenPolynomial(GenPolynomialRing<C> r, C c, ExpVector e) {
094        this(r);
095        if (!c.isZERO()) {
096            val.put(e, c);
097        }
098    }
099
100
101    /**
102     * Constructor for GenPolynomial c * x<sup>0</sup>.
103     * @param r polynomial ring factory.
104     * @param c coefficient.
105     */
106    public GenPolynomial(GenPolynomialRing<C> r, C c) {
107        this(r, c, r.evzero);
108    }
109
110
111    /**
112     * Constructor for GenPolynomial x<sup>e</sup>.
113     * @param r polynomial ring factory.
114     * @param e exponent.
115     */
116    public GenPolynomial(GenPolynomialRing<C> r, ExpVector e) {
117        this(r, r.coFac.getONE(), e);
118    }
119
120
121    /**
122     * Constructor for GenPolynomial.
123     * @param r polynomial ring factory.
124     * @param v the SortedMap of some other polynomial.
125     */
126    protected GenPolynomial(GenPolynomialRing<C> r, SortedMap<ExpVector, C> v) {
127        this(r);
128        val.putAll(v); // assume no zero coefficients
129    }
130
131
132    /**
133     * Get the corresponding element factory.
134     * @return factory for this Element.
135     * @see edu.jas.structure.Element#factory()
136     */
137    public GenPolynomialRing<C> factory() {
138        return ring;
139    }
140
141
142    /**
143     * Copy this GenPolynomial.
144     * @return copy of this.
145     */
146    public GenPolynomial<C> copy() {
147        return new GenPolynomial<C>(ring, this.val);
148    }
149
150
151    /**
152     * Length of GenPolynomial.
153     * @return number of coefficients of this GenPolynomial.
154     */
155    public int length() {
156        return val.size();
157    }
158
159
160    /**
161     * ExpVector to coefficient map of GenPolynomial.
162     * @return val as unmodifiable SortedMap.
163     */
164    public SortedMap<ExpVector, C> getMap() {
165        // return val;
166        return Collections.<ExpVector, C> unmodifiableSortedMap(val);
167    }
168
169
170    /**
171     * Put an ExpVector to coefficient entry into the internal map of this
172     * GenPolynomial. <b>Note:</b> Do not use this method unless you are
173     * constructing a new polynomial. this is modified and breaks the
174     * immutability promise of this class.
175     * @param c coefficient.
176     * @param e exponent.
177     */
178    public void doPutToMap(ExpVector e, C c) {
179        if (debug) {
180            C a = val.get(e);
181            if (a != null) {
182                logger.error("map entry exists " + e + " to " + a + " new " + c);
183            }
184        }
185        if (!c.isZERO()) {
186            val.put(e, c);
187        }
188    }
189
190
191    /**
192     * Remove an ExpVector to coefficient entry from the internal map of this
193     * GenPolynomial. <b>Note:</b> Do not use this method unless you are
194     * constructing a new polynomial. this is modified and breaks the
195     * immutability promise of this class.
196     * @param e exponent.
197     * @param c expected coefficient, null for ignore.
198     */
199    public void doRemoveFromMap(ExpVector e, C c) {
200        C b = val.remove(e);
201        if (debug) {
202            if ( c == null ) { // ignore b
203                return;
204            }
205            if ( ! c.equals(b) ) {
206                logger.error("map entry wrong " + e + " to " + c + " old " + b);
207            }
208        }
209    }
210
211
212    /**
213     * Put an a sorted map of exponents to coefficients into the internal map of
214     * this GenPolynomial. <b>Note:</b> Do not use this method unless you are
215     * constructing a new polynomial. this is modified and breaks the
216     * immutability promise of this class.
217     * @param vals sorted map of exponents and coefficients.
218     */
219    public void doPutToMap(SortedMap<ExpVector, C> vals) {
220        for (Map.Entry<ExpVector, C> me : vals.entrySet()) {
221            ExpVector e = me.getKey();
222            if (debug) {
223                C a = val.get(e);
224                if (a != null) {
225                    logger.error("map entry exists " + e + " to " + a + " new " + me.getValue());
226                }
227            }
228            C c = me.getValue();
229            if (!c.isZERO()) {
230                val.put(e, c);
231            }
232        }
233    }
234
235
236    /**
237     * String representation of GenPolynomial.
238     * @see java.lang.Object#toString()
239     */
240    @Override
241    public String toString() {
242        if (ring.vars != null) {
243            return toString(ring.vars);
244        }
245        StringBuffer s = new StringBuffer();
246        s.append(this.getClass().getSimpleName() + ":");
247        s.append(ring.coFac.getClass().getSimpleName());
248        if (ring.coFac.characteristic().signum() != 0) {
249            s.append("(" + ring.coFac.characteristic() + ")");
250        }
251        s.append("[ ");
252        boolean first = true;
253        for (Map.Entry<ExpVector, C> m : val.entrySet()) {
254            if (first) {
255                first = false;
256            } else {
257                s.append(", ");
258            }
259            s.append(m.getValue().toString());
260            s.append(" ");
261            s.append(m.getKey().toString());
262        }
263        s.append(" ] "); // no not use: ring.toString() );
264        return s.toString();
265    }
266
267
268    /**
269     * String representation of GenPolynomial.
270     * @param v names for variables.
271     * @see java.lang.Object#toString()
272     */
273    public String toString(String[] v) {
274        StringBuffer s = new StringBuffer();
275        if (PrettyPrint.isTrue()) {
276            if (val.size() == 0) {
277                s.append("0");
278            } else {
279                // s.append( "( " );
280                boolean first = true;
281                for (Map.Entry<ExpVector, C> m : val.entrySet()) {
282                    C c = m.getValue();
283                    if (first) {
284                        first = false;
285                    } else {
286                        if (c.signum() < 0) {
287                            s.append(" - ");
288                            c = c.negate();
289                        } else {
290                            s.append(" + ");
291                        }
292                    }
293                    ExpVector e = m.getKey();
294                    if (!c.isONE() || e.isZERO()) {
295                        String cs = c.toString();
296                        if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
297                            s.append("( ");
298                            s.append(cs);
299                            s.append(" )");
300                        } else {
301                            s.append(cs);
302                        }
303                        s.append(" ");
304                    }
305                    if (e != null && v != null) {
306                        s.append(e.toString(v));
307                    } else {
308                        s.append(e);
309                    }
310                }
311                //s.append(" )");
312            }
313        } else {
314            s.append(this.getClass().getSimpleName() + "[ ");
315            if (val.size() == 0) {
316                s.append("0");
317            } else {
318                boolean first = true;
319                for (Map.Entry<ExpVector, C> m : val.entrySet()) {
320                    C c = m.getValue();
321                    if (first) {
322                        first = false;
323                    } else {
324                        if (c.signum() < 0) {
325                            s.append(" - ");
326                            c = c.negate();
327                        } else {
328                            s.append(" + ");
329                        }
330                    }
331                    ExpVector e = m.getKey();
332                    if (!c.isONE() || e.isZERO()) {
333                        s.append(c.toString());
334                        s.append(" ");
335                    }
336                    s.append(e.toString(v));
337                }
338            }
339            s.append(" ] "); // no not use: ring.toString() );
340        }
341        return s.toString();
342    }
343
344
345    /**
346     * Get a scripting compatible string representation.
347     * @return script compatible representation for this Element.
348     * @see edu.jas.structure.Element#toScript()
349     */
350    //JAVA6only: @Override
351    public String toScript() {
352        // Python case
353        if (isZERO()) {
354            return "0";
355        }
356        StringBuffer s = new StringBuffer();
357        if (val.size() > 1) {
358            s.append("( ");
359        }
360        String[] v = ring.vars;
361        if (v == null) {
362            v = GenPolynomialRing.newVars("x", ring.nvar);
363        }
364        boolean parenthesis = false;
365        //if (ring.coFac instanceof GenPolynomialRing || ring.coFac instanceof AlgebraicNumberRing
366        //    || ring.coFac instanceof RealAlgebraicRing
367        //) {
368            // inactive: parenthesis = true;
369        //}
370        boolean first = true;
371        for (Map.Entry<ExpVector, C> m : val.entrySet()) {
372            C c = m.getValue();
373            if (first) {
374                first = false;
375            } else {
376                if (c.signum() < 0) {
377                    s.append(" - ");
378                    c = c.negate();
379                } else {
380                    s.append(" + ");
381                }
382            }
383            ExpVector e = m.getKey();
384            if (!c.isONE() || e.isZERO()) {
385                if (parenthesis) {
386                    s.append("( ");
387                }
388                s.append(c.toScript());
389                if (parenthesis) {
390                    s.append(" )");
391                }
392                if (!e.isZERO()) {
393                    s.append(" * ");
394                }
395            }
396            s.append(e.toScript(v));
397        }
398        if (val.size() > 1) {
399            s.append(" )");
400        }
401        return s.toString();
402    }
403
404
405    /**
406     * Get a scripting compatible string representation of the factory.
407     * @return script compatible representation for this ElemFactory.
408     * @see edu.jas.structure.Element#toScriptFactory()
409     */
410    //JAVA6only: @Override
411    public String toScriptFactory() {
412        // Python case
413        return factory().toScript();
414    }
415
416
417    /**
418     * Is GenPolynomial&lt;C&gt; zero.
419     * @return If this is 0 then true is returned, else false.
420     * @see edu.jas.structure.RingElem#isZERO()
421     */
422    public boolean isZERO() {
423        return (val.size() == 0);
424    }
425
426
427    /**
428     * Is GenPolynomial&lt;C&gt; one.
429     * @return If this is 1 then true is returned, else false.
430     * @see edu.jas.structure.RingElem#isONE()
431     */
432    public boolean isONE() {
433        if (val.size() != 1) {
434            return false;
435        }
436        C c = val.get(ring.evzero);
437        if (c == null) {
438            return false;
439        }
440        return c.isONE();
441    }
442
443
444    /**
445     * Is GenPolynomial&lt;C&gt; a unit.
446     * @return If this is a unit then true is returned, else false.
447     * @see edu.jas.structure.RingElem#isUnit()
448     */
449    public boolean isUnit() {
450        if (val.size() != 1) {
451            return false;
452        }
453        C c = val.get(ring.evzero);
454        if (c == null) {
455            return false;
456        }
457        return c.isUnit();
458    }
459
460
461    /**
462     * Is GenPolynomial&lt;C&gt; a constant.
463     * @return If this is a constant polynomial then true is returned, else
464     *         false.
465     */
466    public boolean isConstant() {
467        if (val.size() != 1) {
468            return false;
469        }
470        C c = val.get(ring.evzero);
471        if (c == null) {
472            return false;
473        }
474        return true;
475    }
476
477
478    /**
479     * Is GenPolynomial&lt;C&gt; homogeneous.
480     * @return true, if this is homogeneous, else false.
481     */
482    public boolean isHomogeneous() {
483        if (val.size() <= 1) {
484            return true;
485        }
486        long deg = -1;
487        for (ExpVector e : val.keySet()) {
488            if (deg < 0) {
489                deg = e.totalDeg();
490            } else if (deg != e.totalDeg()) {
491                return false;
492            }
493        }
494        return true;
495    }
496
497
498    /**
499     * Comparison with any other object.
500     * @see java.lang.Object#equals(java.lang.Object)
501     */
502    @Override
503    @SuppressWarnings("unchecked")
504    public boolean equals(Object B) {
505        if (!(B instanceof GenPolynomial)) {
506            return false;
507        }
508        GenPolynomial<C> a = null;
509        try {
510            a = (GenPolynomial<C>) B;
511        } catch (ClassCastException ignored) {
512        }
513        if (a == null) {
514            return false;
515        }
516        return this.compareTo(a) == 0;
517    }
518
519
520    /**
521     * Hash code for this polynomial.
522     * @see java.lang.Object#hashCode()
523     */
524    @Override
525    public int hashCode() {
526        int h;
527        h = (ring.hashCode() << 27);
528        h += val.hashCode();
529        return h;
530    }
531
532
533    /**
534     * GenPolynomial comparison.
535     * @param b GenPolynomial.
536     * @return sign(this-b).
537     */
538    public int compareTo(GenPolynomial<C> b) {
539        if (b == null) {
540            return 1;
541        }
542        SortedMap<ExpVector, C> av = this.val;
543        SortedMap<ExpVector, C> bv = b.val;
544        Iterator<Map.Entry<ExpVector, C>> ai = av.entrySet().iterator();
545        Iterator<Map.Entry<ExpVector, C>> bi = bv.entrySet().iterator();
546        int s = 0;
547        int c = 0;
548        while (ai.hasNext() && bi.hasNext()) {
549            Map.Entry<ExpVector, C> aie = ai.next();
550            Map.Entry<ExpVector, C> bie = bi.next();
551            ExpVector ae = aie.getKey();
552            ExpVector be = bie.getKey();
553            s = ae.compareTo(be);
554            //System.out.println("s = " + s + ", " + ae + ", " +be);
555            if (s != 0) {
556                return s;
557            }
558            if (c == 0) {
559                C ac = aie.getValue(); //av.get(ae);
560                C bc = bie.getValue(); //bv.get(be);
561                c = ac.compareTo(bc);
562            }
563        }
564        if (ai.hasNext()) {
565            return 1;
566        }
567        if (bi.hasNext()) {
568            return -1;
569        }
570        // now all keys are equal
571        return c;
572    }
573
574
575    /**
576     * GenPolynomial signum.
577     * @return sign(ldcf(this)).
578     */
579    public int signum() {
580        if (this.isZERO()) {
581            return 0;
582        }
583        ExpVector t = val.firstKey();
584        C c = val.get(t);
585        return c.signum();
586    }
587
588
589    /**
590     * Number of variables.
591     * @return ring.nvar.
592     */
593    public int numberOfVariables() {
594        return ring.nvar;
595    }
596
597
598    /**
599     * Leading monomial.
600     * @return first map entry.
601     */
602    public Map.Entry<ExpVector, C> leadingMonomial() {
603        if (val.size() == 0)
604            return null;
605        Iterator<Map.Entry<ExpVector, C>> ai = val.entrySet().iterator();
606        return ai.next();
607    }
608
609
610    /**
611     * Leading exponent vector.
612     * @return first exponent.
613     */
614    public ExpVector leadingExpVector() {
615        if (val.size() == 0) {
616            return null; // ring.evzero? needs many changes 
617        }
618        return val.firstKey();
619    }
620
621
622    /**
623     * Trailing exponent vector.
624     * @return last exponent.
625     */
626    public ExpVector trailingExpVector() {
627        if (val.size() == 0) {
628            return ring.evzero; // or null ?;
629        }
630        return val.lastKey();
631    }
632
633
634    /**
635     * Leading base coefficient.
636     * @return first coefficient.
637     */
638    public C leadingBaseCoefficient() {
639        if (val.size() == 0) {
640            return ring.coFac.getZERO();
641        }
642        return val.get(val.firstKey());
643    }
644
645
646    /**
647     * Trailing base coefficient.
648     * @return coefficient of constant term.
649     */
650    public C trailingBaseCoefficient() {
651        C c = val.get(ring.evzero);
652        if (c == null) {
653            return ring.coFac.getZERO();
654        }
655        return c;
656    }
657
658
659    /**
660     * Coefficient.
661     * @param e exponent.
662     * @return coefficient for given exponent.
663     */
664    public C coefficient(ExpVector e) {
665        C c = val.get(e);
666        if (c == null) {
667            c = ring.coFac.getZERO();
668        }
669        return c;
670    }
671
672
673    /**
674     * Reductum.
675     * @return this - leading monomial.
676     */
677    public GenPolynomial<C> reductum() {
678        if (val.size() <= 1) {
679            return ring.getZERO();
680        }
681        Iterator<ExpVector> ai = val.keySet().iterator();
682        ExpVector lt = ai.next();
683        lt = ai.next(); // size > 1
684        SortedMap<ExpVector, C> red = val.tailMap(lt);
685        return new GenPolynomial<C>(ring, red);
686    }
687
688
689    /**
690     * Degree in variable i.
691     * @return maximal degree in the variable i.
692     */
693    public long degree(int i) {
694        if (val.size() == 0) {
695            return 0; // 0 or -1 ?;
696        }
697        int j;
698        if (i >= 0) {
699            j = ring.nvar - 1 - i;
700        } else { // python like -1 means main variable
701            j = ring.nvar + i;
702        }
703        long deg = 0;
704        for (ExpVector e : val.keySet()) {
705            long d = e.getVal(j);
706            if (d > deg) {
707                deg = d;
708            }
709        }
710        return deg;
711    }
712
713
714    /**
715     * Maximal degree.
716     * @return maximal degree in any variables.
717     */
718    public long degree() {
719        if (val.size() == 0) {
720            return 0; // 0 or -1 ?;
721        }
722        long deg = 0;
723        for (ExpVector e : val.keySet()) {
724            long d = e.maxDeg();
725            if (d > deg) {
726                deg = d;
727            }
728        }
729        return deg;
730    }
731
732
733    /**
734     * Total degree.
735     * @return total degree in any variables.
736     */
737    public long totalDegree() {
738        if (val.size() == 0) {
739            return 0; // 0 or -1 ?;
740        }
741        long deg = 0;
742        for (ExpVector e : val.keySet()) {
743            long d = e.totalDeg();
744            if (d > deg) {
745                deg = d;
746            }
747        }
748        return deg;
749    }
750
751
752    /**
753     * Maximal degree vector.
754     * @return maximal degree vector of all variables.
755     */
756    public ExpVector degreeVector() {
757        ExpVector deg = ring.evzero;
758        if (val.size() == 0) {
759            return deg;
760        }
761        for (ExpVector e : val.keySet()) {
762            deg = deg.lcm(e);
763        }
764        return deg;
765    }
766
767
768    /**
769     * GenPolynomial maximum norm.
770     * @return ||this||.
771     */
772    public C maxNorm() {
773        C n = ring.getZEROCoefficient();
774        for (C c : val.values()) {
775            C x = c.abs();
776            if (n.compareTo(x) < 0) {
777                n = x;
778            }
779        }
780        return n;
781    }
782
783
784    /**
785     * GenPolynomial sum norm.
786     * @return sum of all absolute values of coefficients.
787     */
788    public C sumNorm() {
789        C n = ring.getZEROCoefficient();
790        for (C c : val.values()) {
791            C x = c.abs();
792            n = n.sum(x);
793        }
794        return n;
795    }
796
797
798    /**
799     * GenPolynomial summation.
800     * @param S GenPolynomial.
801     * @return this+S.
802     */
803    //public <T extends GenPolynomial<C>> T sum(T /*GenPolynomial<C>*/ S) {
804    public GenPolynomial<C> sum(GenPolynomial<C> S) {
805        if (S == null) {
806            return this;
807        }
808        if (S.isZERO()) {
809            return this;
810        }
811        if (this.isZERO()) {
812            return S;
813        }
814        assert (ring.nvar == S.ring.nvar);
815        GenPolynomial<C> n = this.copy(); //new GenPolynomial<C>(ring, val); 
816        SortedMap<ExpVector, C> nv = n.val;
817        SortedMap<ExpVector, C> sv = S.val;
818        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
819            ExpVector e = me.getKey();
820            C y = me.getValue(); //sv.get(e); // assert y != null
821            C x = nv.get(e);
822            if (x != null) {
823                x = x.sum(y);
824                if (!x.isZERO()) {
825                    nv.put(e, x);
826                } else {
827                    nv.remove(e);
828                }
829            } else {
830                nv.put(e, y);
831            }
832        }
833        return n;
834    }
835
836
837    /**
838     * GenPolynomial addition. This method is not very efficient, since this is
839     * copied.
840     * @param a coefficient.
841     * @param e exponent.
842     * @return this + a x<sup>e</sup>.
843     */
844    public GenPolynomial<C> sum(C a, ExpVector e) {
845        if (a == null) {
846            return this;
847        }
848        if (a.isZERO()) {
849            return this;
850        }
851        GenPolynomial<C> n = this.copy(); //new GenPolynomial<C>(ring, val); 
852        SortedMap<ExpVector, C> nv = n.val;
853        //if ( nv.size() == 0 ) { nv.put(e,a); return n; }
854        C x = nv.get(e);
855        if (x != null) {
856            x = x.sum(a);
857            if (!x.isZERO()) {
858                nv.put(e, x);
859            } else {
860                nv.remove(e);
861            }
862        } else {
863            nv.put(e, a);
864        }
865        return n;
866    }
867
868
869    /**
870     * GenPolynomial addition. This method is not very efficient, since this is
871     * copied.
872     * @param a coefficient.
873     * @return this + a x<sup>0</sup>.
874     */
875    public GenPolynomial<C> sum(C a) {
876        return sum(a, ring.evzero);
877    }
878
879
880    /**
881     * GenPolynomial subtraction.
882     * @param S GenPolynomial.
883     * @return this-S.
884     */
885    public GenPolynomial<C> subtract(GenPolynomial<C> S) {
886        if (S == null) {
887            return this;
888        }
889        if (S.isZERO()) {
890            return this;
891        }
892        if (this.isZERO()) {
893            return S.negate();
894        }
895        assert (ring.nvar == S.ring.nvar);
896        GenPolynomial<C> n = this.copy(); //new GenPolynomial<C>(ring, val); 
897        SortedMap<ExpVector, C> nv = n.val;
898        SortedMap<ExpVector, C> sv = S.val;
899        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
900            ExpVector e = me.getKey();
901            C y = me.getValue(); //sv.get(e); // assert y != null
902            C x = nv.get(e);
903            if (x != null) {
904                x = x.subtract(y);
905                if (!x.isZERO()) {
906                    nv.put(e, x);
907                } else {
908                    nv.remove(e);
909                }
910            } else {
911                nv.put(e, y.negate());
912            }
913        }
914        return n;
915    }
916
917
918    /**
919     * GenPolynomial subtraction. This method is not very efficient, since this
920     * is copied.
921     * @param a coefficient.
922     * @param e exponent.
923     * @return this - a x<sup>e</sup>.
924     */
925    public GenPolynomial<C> subtract(C a, ExpVector e) {
926        if (a == null) {
927            return this;
928        }
929        if (a.isZERO()) {
930            return this;
931        }
932        GenPolynomial<C> n = this.copy(); 
933        SortedMap<ExpVector, C> nv = n.val;
934        C x = nv.get(e);
935        if (x != null) {
936            x = x.subtract(a);
937            if (!x.isZERO()) {
938                nv.put(e, x);
939            } else {
940                nv.remove(e);
941            }
942        } else {
943            nv.put(e, a.negate());
944        }
945        return n;
946    }
947
948
949    /**
950     * GenPolynomial subtract. This method is not very efficient, since this is
951     * copied.
952     * @param a coefficient.
953     * @return this + a x<sup>0</sup>.
954     */
955    public GenPolynomial<C> subtract(C a) {
956        return subtract(a, ring.evzero);
957    }
958
959
960    /**
961     * GenPolynomial subtract a multiple.
962     * @param a coefficient.
963     * @param S GenPolynomial.
964     * @return this - a x<sup>e</sup> S.
965     */
966    public GenPolynomial<C> subtractMultiple(C a, GenPolynomial<C> S) {
967        if (a == null) {
968            return this;
969        }
970        if (a.isZERO()) {
971            return this;
972        }
973        if (S == null) {
974            return this;
975        }
976        if (S.isZERO()) {
977            return this;
978        }
979        if (this.isZERO()) {
980            return S.multiply(a.negate());
981        }
982        assert (ring.nvar == S.ring.nvar);
983        GenPolynomial<C> n = this.copy(); 
984        SortedMap<ExpVector, C> nv = n.val;
985        SortedMap<ExpVector, C> sv = S.val;
986        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
987            ExpVector f = me.getKey();
988            C y = me.getValue(); // assert y != null
989            y = a.multiply(y);
990            C x = nv.get(f);
991            if (x != null) {
992                x = x.subtract(y);
993                if (!x.isZERO()) {
994                    nv.put(f, x);
995                } else {
996                    nv.remove(f);
997                }
998            } else if ( !y.isZERO() ) {
999                nv.put(f, y.negate());
1000            }
1001        }
1002        return n;
1003    }
1004
1005
1006    /**
1007     * GenPolynomial subtract a multiple.
1008     * @param a coefficient.
1009     * @param e exponent.
1010     * @param S GenPolynomial.
1011     * @return this - a x<sup>e</sup> S.
1012     */
1013    public GenPolynomial<C> subtractMultiple(C a, ExpVector e, GenPolynomial<C> S) {
1014        if (a == null) {
1015            return this;
1016        }
1017        if (a.isZERO()) {
1018            return this;
1019        }
1020        if (S == null) {
1021            return this;
1022        }
1023        if (S.isZERO()) {
1024            return this;
1025        }
1026        if (this.isZERO()) {
1027            return S.multiply(a.negate(),e);
1028        }
1029        assert (ring.nvar == S.ring.nvar);
1030        GenPolynomial<C> n = this.copy(); 
1031        SortedMap<ExpVector, C> nv = n.val;
1032        SortedMap<ExpVector, C> sv = S.val;
1033        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
1034            ExpVector f = me.getKey();
1035            f = e.sum(f);
1036            C y = me.getValue(); // assert y != null
1037            y = a.multiply(y);
1038            C x = nv.get(f);
1039            if (x != null) {
1040                x = x.subtract(y);
1041                if (!x.isZERO()) {
1042                    nv.put(f, x);
1043                } else {
1044                    nv.remove(f);
1045                }
1046            } else if ( !y.isZERO() ) {
1047                nv.put(f, y.negate());
1048            }
1049        }
1050        return n;
1051    }
1052
1053
1054    /**
1055     * GenPolynomial scale and subtract a multiple.
1056     * @param b scale factor.
1057     * @param a coefficient.
1058     * @param e exponent.
1059     * @param S GenPolynomial.
1060     * @return this * b - a x<sup>e</sup> S.
1061     */
1062    public GenPolynomial<C> scaleSubtractMultiple(C b, C a, ExpVector e, GenPolynomial<C> S) {
1063        if (a == null || S == null) {
1064            return this.multiply(b);
1065        }
1066        if (a.isZERO() || S.isZERO()) {
1067            return this.multiply(b);
1068        }
1069        if (this.isZERO() || b == null || b.isZERO()) {
1070            return S.multiply(a.negate(),e);
1071        }
1072        if (b.isONE()) {
1073            return subtractMultiple(a,e,S);
1074        }
1075        assert (ring.nvar == S.ring.nvar);
1076        GenPolynomial<C> n = this.multiply(b);
1077        SortedMap<ExpVector, C> nv = n.val;
1078        SortedMap<ExpVector, C> sv = S.val;
1079        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
1080            ExpVector f = me.getKey();
1081            f = e.sum(f);
1082            C y = me.getValue(); // assert y != null
1083            y = a.multiply(y); // now y can be zero
1084            C x = nv.get(f);
1085            if (x != null) {
1086                x = x.subtract(y);
1087                if (!x.isZERO()) {
1088                    nv.put(f, x);
1089                } else {
1090                    nv.remove(f);
1091                }
1092            } else if ( !y.isZERO() ) {
1093                nv.put(f, y.negate());
1094            }
1095        }
1096        return n;
1097    }
1098
1099
1100    /**
1101     * GenPolynomial scale and subtract a multiple.
1102     * @param b scale factor.
1103     * @param g scale exponent.
1104     * @param a coefficient.
1105     * @param e exponent.
1106     * @param S GenPolynomial.
1107     * @return this * a x<sup>g</sup> - a x<sup>e</sup> S.
1108     */
1109    public GenPolynomial<C> scaleSubtractMultiple(C b, ExpVector g, C a, ExpVector e, GenPolynomial<C> S) {
1110        if (a == null || S == null) {
1111            return this.multiply(b,g);
1112        }
1113        if (a.isZERO() || S.isZERO()) {
1114            return this.multiply(b,g);
1115        }
1116        if (this.isZERO() || b == null || b.isZERO()) {
1117            return S.multiply(a.negate(),e);
1118        }
1119        if (b.isONE() && g.isZERO()) {
1120            return subtractMultiple(a,e,S);
1121        }
1122        assert (ring.nvar == S.ring.nvar);
1123        GenPolynomial<C> n = this.multiply(b,g);
1124        SortedMap<ExpVector, C> nv = n.val;
1125        SortedMap<ExpVector, C> sv = S.val;
1126        for (Map.Entry<ExpVector, C> me : sv.entrySet()) {
1127            ExpVector f = me.getKey();
1128            f = e.sum(f);
1129            C y = me.getValue(); // assert y != null
1130            y = a.multiply(y); // y can be zero now
1131            C x = nv.get(f);
1132            if (x != null) {
1133                x = x.subtract(y);
1134                if (!x.isZERO()) {
1135                    nv.put(f, x);
1136                } else {
1137                    nv.remove(f);
1138                }
1139            } else if ( !y.isZERO() ) {
1140                nv.put(f, y.negate());
1141            }
1142        }
1143        return n;
1144    }
1145
1146
1147    /**
1148     * GenPolynomial negation.
1149     * @return -this.
1150     */
1151    public GenPolynomial<C> negate() {
1152        GenPolynomial<C> n = ring.getZERO().copy();
1153        //new GenPolynomial<C>(ring, ring.getZERO().val);
1154        SortedMap<ExpVector, C> v = n.val;
1155        for (Map.Entry<ExpVector, C> m : val.entrySet()) {
1156            C x = m.getValue(); // != null, 0
1157            v.put(m.getKey(), x.negate());
1158            // or m.setValue( x.negate() ) if this cloned 
1159        }
1160        return n;
1161    }
1162
1163
1164    /**
1165     * GenPolynomial absolute value, i.e. leadingCoefficient &gt; 0.
1166     * @return abs(this).
1167     */
1168    public GenPolynomial<C> abs() {
1169        if (leadingBaseCoefficient().signum() < 0) {
1170            return this.negate();
1171        }
1172        return this;
1173    }
1174
1175
1176    /**
1177     * GenPolynomial multiplication.
1178     * @param S GenPolynomial.
1179     * @return this*S.
1180     */
1181    public GenPolynomial<C> multiply(GenPolynomial<C> S) {
1182        if (S == null) {
1183            return ring.getZERO();
1184        }
1185        if (S.isZERO()) {
1186            return ring.getZERO();
1187        }
1188        if (this.isZERO()) {
1189            return this;
1190        }
1191        assert (ring.nvar == S.ring.nvar);
1192        if (this instanceof GenSolvablePolynomial || S instanceof GenSolvablePolynomial) {
1193            //throw new RuntimeException("wrong method dispatch in JRE ");
1194            logger.warn("wrong method dispatch in JRE multiply(S) - trying to fix");
1195            GenSolvablePolynomial<C> T = (GenSolvablePolynomial<C>) this;
1196            GenSolvablePolynomial<C> Sp = (GenSolvablePolynomial<C>) S;
1197            return T.multiply(Sp);
1198        }
1199        GenPolynomial<C> p = ring.getZERO().copy();
1200        SortedMap<ExpVector, C> pv = p.val;
1201        for (Map.Entry<ExpVector, C> m1 : val.entrySet()) {
1202            C c1 = m1.getValue();
1203            ExpVector e1 = m1.getKey();
1204            for (Map.Entry<ExpVector, C> m2 : S.val.entrySet()) {
1205                C c2 = m2.getValue();
1206                ExpVector e2 = m2.getKey();
1207                C c = c1.multiply(c2); // check non zero if not domain
1208                if (!c.isZERO()) {
1209                    ExpVector e = e1.sum(e2);
1210                    C c0 = pv.get(e);
1211                    if (c0 == null) {
1212                        pv.put(e, c);
1213                    } else {
1214                        c0 = c0.sum(c);
1215                        if (!c0.isZERO()) {
1216                            pv.put(e, c0);
1217                        } else {
1218                            pv.remove(e);
1219                        }
1220                    }
1221                }
1222            }
1223        }
1224        return p;
1225    }
1226
1227
1228    /**
1229     * GenPolynomial multiplication. Product with coefficient ring element.
1230     * @param s coefficient.
1231     * @return this*s.
1232     */
1233    public GenPolynomial<C> multiply(C s) {
1234        if (s == null) {
1235            return ring.getZERO();
1236        }
1237        if (s.isZERO()) {
1238            return ring.getZERO();
1239        }
1240        if (this.isZERO()) {
1241            return this;
1242        }
1243        GenPolynomial<C> p = ring.getZERO().copy();
1244        SortedMap<ExpVector, C> pv = p.val;
1245        for (Map.Entry<ExpVector, C> m1 : val.entrySet()) {
1246            C c1 = m1.getValue();
1247            ExpVector e1 = m1.getKey();
1248            C c = c1.multiply(s); // check non zero if not domain
1249            if (!c.isZERO()) {
1250                pv.put(e1, c); // or m1.setValue( c )
1251            }
1252        }
1253        return p;
1254    }
1255
1256
1257    /**
1258     * GenPolynomial monic, i.e. leadingCoefficient == 1. If leadingCoefficient
1259     * is not invertible returns this unmodified.
1260     * @return monic(this).
1261     */
1262    public GenPolynomial<C> monic() {
1263        if (this.isZERO()) {
1264            return this;
1265        }
1266        C lc = leadingBaseCoefficient();
1267        if (!lc.isUnit()) {
1268            //System.out.println("lc = "+lc);
1269            return this;
1270        }
1271        C lm = lc.inverse();
1272        return multiply(lm);
1273    }
1274
1275
1276    /**
1277     * GenPolynomial multiplication. Product with ring element and exponent
1278     * vector.
1279     * @param s coefficient.
1280     * @param e exponent.
1281     * @return this * s x<sup>e</sup>.
1282     */
1283    public GenPolynomial<C> multiply(C s, ExpVector e) {
1284        if (s == null) {
1285            return ring.getZERO();
1286        }
1287        if (s.isZERO()) {
1288            return ring.getZERO();
1289        }
1290        if (this.isZERO()) {
1291            return this;
1292        }
1293        if (this instanceof GenSolvablePolynomial) {
1294            //throw new RuntimeException("wrong method dispatch in JRE ");
1295            logger.warn("wrong method dispatch in JRE multiply(s,e) - trying to fix");
1296            GenSolvablePolynomial<C> T = (GenSolvablePolynomial<C>) this;
1297            return T.multiply(s, e);
1298        }
1299        GenPolynomial<C> p = ring.getZERO().copy();
1300        SortedMap<ExpVector, C> pv = p.val;
1301        for (Map.Entry<ExpVector, C> m1 : val.entrySet()) {
1302            C c1 = m1.getValue();
1303            ExpVector e1 = m1.getKey();
1304            C c = c1.multiply(s); // check non zero if not domain
1305            if (!c.isZERO()) {
1306                ExpVector e2 = e1.sum(e);
1307                pv.put(e2, c);
1308            }
1309        }
1310        return p;
1311    }
1312
1313
1314    /**
1315     * GenPolynomial multiplication. Product with exponent vector.
1316     * @param e exponent (!= null).
1317     * @return this * x<sup>e</sup>.
1318     */
1319    public GenPolynomial<C> multiply(ExpVector e) {
1320        // assert e != null. This is never allowed.
1321        if (this.isZERO()) {
1322            return this;
1323        }
1324        if (this instanceof GenSolvablePolynomial) {
1325            //throw new RuntimeException("wrong method dispatch in JRE ");
1326            logger.warn("wrong method dispatch in JRE multiply(e) - trying to fix");
1327            GenSolvablePolynomial<C> T = (GenSolvablePolynomial<C>) this;
1328            return T.multiply(e);
1329        }
1330        GenPolynomial<C> p = ring.getZERO().copy();
1331        SortedMap<ExpVector, C> pv = p.val;
1332        for (Map.Entry<ExpVector, C> m1 : val.entrySet()) {
1333            C c1 = m1.getValue();
1334            ExpVector e1 = m1.getKey();
1335            ExpVector e2 = e1.sum(e);
1336            pv.put(e2, c1);
1337        }
1338        return p;
1339    }
1340
1341
1342    /**
1343     * GenPolynomial multiplication. Product with 'monomial'.
1344     * @param m 'monomial'.
1345     * @return this * m.
1346     */
1347    public GenPolynomial<C> multiply(Map.Entry<ExpVector, C> m) {
1348        if (m == null) {
1349            return ring.getZERO();
1350        }
1351        return multiply(m.getValue(), m.getKey());
1352    }
1353
1354
1355    /**
1356     * GenPolynomial division. Division by coefficient ring element. Fails, if
1357     * exact division is not possible.
1358     * @param s coefficient.
1359     * @return this/s.
1360     */
1361    public GenPolynomial<C> divide(C s) {
1362        if (s == null || s.isZERO()) {
1363            throw new ArithmeticException("division by zero");
1364        }
1365        if (this.isZERO()) {
1366            return this;
1367        }
1368        //C t = s.inverse();
1369        //return multiply(t);
1370        GenPolynomial<C> p = ring.getZERO().copy();
1371        SortedMap<ExpVector, C> pv = p.val;
1372        for (Map.Entry<ExpVector, C> m : val.entrySet()) {
1373            ExpVector e = m.getKey();
1374            C c1 = m.getValue();
1375            C c = c1.divide(s);
1376            if (debug) {
1377                C x = c1.remainder(s);
1378                if (!x.isZERO()) {
1379                    logger.info("divide x = " + x);
1380                    throw new ArithmeticException("no exact division: " + c1 + "/" + s);
1381                }
1382            }
1383            if (c.isZERO()) {
1384                throw new ArithmeticException("no exact division: " + c1 + "/" + s + ", in " + this);
1385            }
1386            pv.put(e, c); // or m1.setValue( c )
1387        }
1388        return p;
1389    }
1390
1391
1392    /**
1393     * GenPolynomial division with remainder. Fails, if exact division by
1394     * leading base coefficient is not possible. Meaningful only for univariate
1395     * polynomials over fields, but works in any case.
1396     * @param S nonzero GenPolynomial with invertible leading coefficient.
1397     * @return [ quotient , remainder ] with this = quotient * S + remainder and
1398     *         deg(remainder) &lt; deg(S) or remiander = 0.
1399     * @see edu.jas.poly.PolyUtil#baseSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)
1400     *      .
1401     */
1402    @SuppressWarnings("unchecked")
1403    public GenPolynomial<C>[] quotientRemainder(GenPolynomial<C> S) {
1404        if (S == null || S.isZERO()) {
1405            throw new ArithmeticException("division by zero");
1406        }
1407        C c = S.leadingBaseCoefficient();
1408        if (!c.isUnit()) {
1409            throw new ArithmeticException("lbcf not invertible " + c);
1410        }
1411        C ci = c.inverse();
1412        assert (ring.nvar == S.ring.nvar);
1413        ExpVector e = S.leadingExpVector();
1414        GenPolynomial<C> h;
1415        GenPolynomial<C> q = ring.getZERO().copy();
1416        GenPolynomial<C> r = this.copy();
1417        while (!r.isZERO()) {
1418            ExpVector f = r.leadingExpVector();
1419            if (f.multipleOf(e)) {
1420                C a = r.leadingBaseCoefficient();
1421                f = f.subtract(e);
1422                a = a.multiply(ci);
1423                q = q.sum(a, f);
1424                h = S.multiply(a, f);
1425                r = r.subtract(h);
1426            } else {
1427                break;
1428            }
1429        }
1430        GenPolynomial<C>[] ret = new GenPolynomial[2];
1431        ret[0] = q;
1432        ret[1] = r;
1433        return ret;
1434    }
1435
1436
1437    /**
1438     * GenPolynomial division with remainder. Fails, if exact division by
1439     * leading base coefficient is not possible. Meaningful only for univariate
1440     * polynomials over fields, but works in any case.
1441     * @param S nonzero GenPolynomial with invertible leading coefficient.
1442     * @return [ quotient , remainder ] with this = quotient * S + remainder and
1443     *         deg(remainder) &lt; deg(S) or remiander = 0.
1444     * @see edu.jas.poly.PolyUtil#baseSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)
1445     *      .
1446     * @deprecated use quotientRemainder()
1447     */
1448    @Deprecated
1449    public GenPolynomial<C>[] divideAndRemainder(GenPolynomial<C> S) {
1450        return quotientRemainder(S);
1451    }
1452
1453
1454    /**
1455     * GenPolynomial division. Fails, if exact division by leading base
1456     * coefficient is not possible. Meaningful only for univariate polynomials
1457     * over fields, but works in any case.
1458     * @param S nonzero GenPolynomial with invertible leading coefficient.
1459     * @return quotient with this = quotient * S + remainder.
1460     * @see edu.jas.poly.PolyUtil#baseSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)
1461     *      .
1462     */
1463    public GenPolynomial<C> divide(GenPolynomial<C> S) {
1464        return quotientRemainder(S)[0];
1465    }
1466
1467
1468    /**
1469     * GenPolynomial remainder. Fails, if exact division by leading base
1470     * coefficient is not possible. Meaningful only for univariate polynomials
1471     * over fields, but works in any case.
1472     * @param S nonzero GenPolynomial with invertible leading coefficient.
1473     * @return remainder with this = quotient * S + remainder.
1474     * @see edu.jas.poly.PolyUtil#baseSparsePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)
1475     *      .
1476     */
1477    public GenPolynomial<C> remainder(GenPolynomial<C> S) {
1478        if (S == null || S.isZERO()) {
1479            throw new ArithmeticException("division by zero");
1480        }
1481        C c = S.leadingBaseCoefficient();
1482        if (!c.isUnit()) {
1483            throw new ArithmeticException("lbc not invertible " + c);
1484        }
1485        C ci = c.inverse();
1486        assert (ring.nvar == S.ring.nvar);
1487        ExpVector e = S.leadingExpVector();
1488        GenPolynomial<C> h;
1489        GenPolynomial<C> r = this.copy();
1490        while (!r.isZERO()) {
1491            ExpVector f = r.leadingExpVector();
1492            if (f.multipleOf(e)) {
1493                C a = r.leadingBaseCoefficient();
1494                f = f.subtract(e);
1495                //logger.info("red div = " + e);
1496                a = a.multiply(ci);
1497                h = S.multiply(a, f);
1498                r = r.subtract(h);
1499            } else {
1500                break;
1501            }
1502        }
1503        return r;
1504    }
1505
1506
1507    /**
1508     * GenPolynomial greatest common divisor. Only for univariate polynomials
1509     * over fields.
1510     * @param S GenPolynomial.
1511     * @return gcd(this,S).
1512     */
1513    public GenPolynomial<C> gcd(GenPolynomial<C> S) {
1514        if (S == null || S.isZERO()) {
1515            return this;
1516        }
1517        if (this.isZERO()) {
1518            return S;
1519        }
1520        if (ring.nvar != 1) {
1521            throw new IllegalArgumentException("not univariate polynomials" + ring);
1522        }
1523        GenPolynomial<C> x;
1524        GenPolynomial<C> q = this;
1525        GenPolynomial<C> r = S;
1526        while (!r.isZERO()) {
1527            x = q.remainder(r);
1528            q = r;
1529            r = x;
1530        }
1531        return q.monic(); // normalize
1532    }
1533
1534
1535    /**
1536     * GenPolynomial extended greatest comon divisor. Only for univariate
1537     * polynomials over fields.
1538     * @param S GenPolynomial.
1539     * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
1540     */
1541    @SuppressWarnings("unchecked")
1542    public GenPolynomial<C>[] egcd(GenPolynomial<C> S) {
1543        GenPolynomial<C>[] ret = new GenPolynomial[3];
1544        ret[0] = null;
1545        ret[1] = null;
1546        ret[2] = null;
1547        if (S == null || S.isZERO()) {
1548            ret[0] = this;
1549            ret[1] = this.ring.getONE();
1550            ret[2] = this.ring.getZERO();
1551            return ret;
1552        }
1553        if (this.isZERO()) {
1554            ret[0] = S;
1555            ret[1] = this.ring.getZERO();
1556            ret[2] = this.ring.getONE();
1557            return ret;
1558        }
1559        if (ring.nvar != 1) {
1560            throw new IllegalArgumentException(this.getClass().getName() + " not univariate polynomials"
1561                            + ring);
1562        }
1563        if (this.isConstant() && S.isConstant()) {
1564            C t = this.leadingBaseCoefficient();
1565            C s = S.leadingBaseCoefficient();
1566            C[] gg = t.egcd(s);
1567            //System.out.println("coeff gcd = " + Arrays.toString(gg));
1568            GenPolynomial<C> z = this.ring.getZERO();
1569            ret[0] = z.sum(gg[0]);
1570            ret[1] = z.sum(gg[1]);
1571            ret[2] = z.sum(gg[2]);
1572            return ret;
1573        }
1574        GenPolynomial<C>[] qr;
1575        GenPolynomial<C> q = this;
1576        GenPolynomial<C> r = S;
1577        GenPolynomial<C> c1 = ring.getONE().copy();
1578        GenPolynomial<C> d1 = ring.getZERO().copy();
1579        GenPolynomial<C> c2 = ring.getZERO().copy();
1580        GenPolynomial<C> d2 = ring.getONE().copy();
1581        GenPolynomial<C> x1;
1582        GenPolynomial<C> x2;
1583        while (!r.isZERO()) {
1584            qr = q.quotientRemainder(r);
1585            q = qr[0];
1586            x1 = c1.subtract(q.multiply(d1));
1587            x2 = c2.subtract(q.multiply(d2));
1588            c1 = d1;
1589            c2 = d2;
1590            d1 = x1;
1591            d2 = x2;
1592            q = r;
1593            r = qr[1];
1594        }
1595        // normalize ldcf(q) to 1, i.e. make monic
1596        C g = q.leadingBaseCoefficient();
1597        if (g.isUnit()) {
1598            C h = g.inverse();
1599            q = q.multiply(h);
1600            c1 = c1.multiply(h);
1601            c2 = c2.multiply(h);
1602        }
1603        //assert ( ((c1.multiply(this)).sum( c2.multiply(S)).equals(q) )); 
1604        ret[0] = q;
1605        ret[1] = c1;
1606        ret[2] = c2;
1607        return ret;
1608    }
1609
1610
1611    /**
1612     * GenPolynomial half extended greatest comon divisor. Only for univariate
1613     * polynomials over fields.
1614     * @param S GenPolynomial.
1615     * @return [ gcd(this,S), a ] with a*this + b*S = gcd(this,S).
1616     */
1617    @SuppressWarnings("unchecked")
1618    public GenPolynomial<C>[] hegcd(GenPolynomial<C> S) {
1619        GenPolynomial<C>[] ret = new GenPolynomial[2];
1620        ret[0] = null;
1621        ret[1] = null;
1622        if (S == null || S.isZERO()) {
1623            ret[0] = this;
1624            ret[1] = this.ring.getONE();
1625            return ret;
1626        }
1627        if (this.isZERO()) {
1628            ret[0] = S;
1629            return ret;
1630        }
1631        if (ring.nvar != 1) {
1632            throw new IllegalArgumentException(this.getClass().getName() + " not univariate polynomials"
1633                            + ring);
1634        }
1635        GenPolynomial<C>[] qr;
1636        GenPolynomial<C> q = this;
1637        GenPolynomial<C> r = S;
1638        GenPolynomial<C> c1 = ring.getONE().copy();
1639        GenPolynomial<C> d1 = ring.getZERO().copy();
1640        GenPolynomial<C> x1;
1641        while (!r.isZERO()) {
1642            qr = q.quotientRemainder(r);
1643            q = qr[0];
1644            x1 = c1.subtract(q.multiply(d1));
1645            c1 = d1;
1646            d1 = x1;
1647            q = r;
1648            r = qr[1];
1649        }
1650        // normalize ldcf(q) to 1, i.e. make monic
1651        C g = q.leadingBaseCoefficient();
1652        if (g.isUnit()) {
1653            C h = g.inverse();
1654            q = q.multiply(h);
1655            c1 = c1.multiply(h);
1656        }
1657        //assert ( ((c1.multiply(this)).remainder(S).equals(q) )); 
1658        ret[0] = q;
1659        ret[1] = c1;
1660        return ret;
1661    }
1662
1663
1664    /**
1665     * GenPolynomial inverse. Required by RingElem. Throws not invertible
1666     * exception.
1667     */
1668    public GenPolynomial<C> inverse() {
1669        if (isUnit()) { // only possible if ldbcf is unit
1670            C c = leadingBaseCoefficient().inverse();
1671            return ring.getONE().multiply(c);
1672        }
1673        throw new NotInvertibleException("element not invertible " + this + " :: " + ring);
1674    }
1675
1676
1677    /**
1678     * GenPolynomial modular inverse. Only for univariate polynomials over
1679     * fields.
1680     * @param m GenPolynomial.
1681     * @return a with with a*this = 1 mod m.
1682     */
1683    public GenPolynomial<C> modInverse(GenPolynomial<C> m) {
1684        if (this.isZERO()) {
1685            throw new NotInvertibleException("zero is not invertible");
1686        }
1687        GenPolynomial<C>[] hegcd = this.hegcd(m);
1688        GenPolynomial<C> a = hegcd[0];
1689        if (!a.isUnit()) { // gcd != 1
1690            throw new AlgebraicNotInvertibleException("element not invertible, gcd != 1", m, a, m.divide(a));
1691        }
1692        GenPolynomial<C> b = hegcd[1];
1693        if (b.isZERO()) { // when m divides this, e.g. m.isUnit()
1694            throw new NotInvertibleException("element not invertible, divisible by modul");
1695        }
1696        return b;
1697    }
1698
1699
1700    /**
1701     * Extend variables. Used e.g. in module embedding. Extend all ExpVectors by
1702     * i elements and multiply by x_j^k.
1703     * @param pfac extended polynomial ring factory (by i variables).
1704     * @param j index of variable to be used for multiplication.
1705     * @param k exponent for x_j.
1706     * @return extended polynomial.
1707     */
1708    public GenPolynomial<C> extend(GenPolynomialRing<C> pfac, int j, long k) {
1709        if (ring.equals(pfac)) { // nothing to do
1710            return this;
1711        }
1712        GenPolynomial<C> Cp = pfac.getZERO().copy();
1713        if (this.isZERO()) {
1714            return Cp;
1715        }
1716        int i = pfac.nvar - ring.nvar;
1717        Map<ExpVector, C> C = Cp.val; //getMap();
1718        Map<ExpVector, C> A = val;
1719        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1720            ExpVector e = y.getKey();
1721            C a = y.getValue();
1722            ExpVector f = e.extend(i, j, k);
1723            C.put(f, a);
1724        }
1725        return Cp;
1726    }
1727
1728
1729    /**
1730     * Extend lower variables. Used e.g. in module embedding. Extend all
1731     * ExpVectors by i lower elements and multiply by x_j^k.
1732     * @param pfac extended polynomial ring factory (by i variables).
1733     * @param j index of variable to be used for multiplication.
1734     * @param k exponent for x_j.
1735     * @return extended polynomial.
1736     */
1737    public GenPolynomial<C> extendLower(GenPolynomialRing<C> pfac, int j, long k) {
1738        if (ring.equals(pfac)) { // nothing to do
1739            return this;
1740        }
1741        GenPolynomial<C> Cp = pfac.getZERO().copy();
1742        if (this.isZERO()) {
1743            return Cp;
1744        }
1745        int i = pfac.nvar - ring.nvar;
1746        Map<ExpVector, C> C = Cp.val; //getMap();
1747        Map<ExpVector, C> A = val;
1748        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1749            ExpVector e = y.getKey();
1750            C a = y.getValue();
1751            ExpVector f = e.extendLower(i, j, k);
1752            C.put(f, a);
1753        }
1754        return Cp;
1755    }
1756
1757
1758    /**
1759     * Contract variables. Used e.g. in module embedding. Remove i elements of
1760     * each ExpVector.
1761     * @param pfac contracted polynomial ring factory (by i variables).
1762     * @return Map of exponents and contracted polynomials. <b>Note:</b> could
1763     *         return SortedMap
1764     */
1765    public Map<ExpVector, GenPolynomial<C>> contract(GenPolynomialRing<C> pfac) {
1766        GenPolynomial<C> zero = pfac.getZERO();
1767        TermOrder t = new TermOrder(TermOrder.INVLEX);
1768        Map<ExpVector, GenPolynomial<C>> B = new TreeMap<ExpVector, GenPolynomial<C>>(t.getAscendComparator());
1769        if (this.isZERO()) {
1770            return B;
1771        }
1772        int i = ring.nvar - pfac.nvar;
1773        Map<ExpVector, C> A = val;
1774        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1775            ExpVector e = y.getKey();
1776            C a = y.getValue();
1777            ExpVector f = e.contract(0, i);
1778            ExpVector g = e.contract(i, e.length() - i);
1779            GenPolynomial<C> p = B.get(f);
1780            if (p == null) {
1781                p = zero;
1782            }
1783            p = p.sum(a, g);
1784            B.put(f, p);
1785        }
1786        return B;
1787    }
1788
1789
1790    /**
1791     * Contract variables to coefficient polynomial. Remove i elements of each
1792     * ExpVector, removed elements must be zero.
1793     * @param pfac contracted polynomial ring factory (by i variables).
1794     * @return contracted coefficient polynomial.
1795     */
1796    public GenPolynomial<C> contractCoeff(GenPolynomialRing<C> pfac) {
1797        Map<ExpVector, GenPolynomial<C>> ms = contract(pfac);
1798        GenPolynomial<C> c = pfac.getZERO();
1799        for (Map.Entry<ExpVector, GenPolynomial<C>> m : ms.entrySet()) {
1800            if (m.getKey().isZERO()) {
1801                c = m.getValue();
1802            } else {
1803                throw new RuntimeException("wrong coefficient contraction " + m + ", pol =  " + c);
1804            }
1805        }
1806        return c;
1807    }
1808
1809
1810    /**
1811     * Extend univariate to multivariate polynomial. This is an univariate
1812     * polynomial in variable i of the polynomial ring, it is extended to the
1813     * given polynomial ring.
1814     * @param pfac extended polynomial ring factory.
1815     * @param i index of the variable of this polynomial in pfac.
1816     * @return extended multivariate polynomial.
1817     */
1818    public GenPolynomial<C> extendUnivariate(GenPolynomialRing<C> pfac, int i) {
1819        if (i < 0 || pfac.nvar < i) {
1820            throw new IllegalArgumentException("index " + i + "out of range " + pfac.nvar);
1821        }
1822        if (ring.nvar != 1) {
1823            throw new IllegalArgumentException("polynomial not univariate " + ring.nvar);
1824        }
1825        if (this.isONE()) {
1826            return pfac.getONE();
1827        }
1828        int j = pfac.nvar - 1 - i;
1829        GenPolynomial<C> Cp = pfac.getZERO().copy();
1830        if (this.isZERO()) {
1831            return Cp;
1832        }
1833        Map<ExpVector, C> C = Cp.val; //getMap();
1834        Map<ExpVector, C> A = val;
1835        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1836            ExpVector e = y.getKey();
1837            long n = e.getVal(0);
1838            C a = y.getValue();
1839            ExpVector f = ExpVector.create(pfac.nvar, j, n);
1840            C.put(f, a); // assert not contained
1841        }
1842        return Cp;
1843    }
1844
1845
1846    /**
1847     * Make homogeneous.
1848     * @param pfac extended polynomial ring factory (by 1 variable).
1849     * @return homogeneous polynomial.
1850     */
1851    public GenPolynomial<C> homogenize(GenPolynomialRing<C> pfac) {
1852        if (ring.equals(pfac)) { // not implemented
1853            throw new UnsupportedOperationException("case with same ring not implemented");
1854        }
1855        GenPolynomial<C> Cp = pfac.getZERO().copy();
1856        if (this.isZERO()) {
1857            return Cp;
1858        }
1859        long deg = totalDegree();
1860        //int i = pfac.nvar - ring.nvar;
1861        Map<ExpVector, C> C = Cp.val; //getMap();
1862        Map<ExpVector, C> A = val;
1863        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1864            ExpVector e = y.getKey();
1865            C a = y.getValue();
1866            long d = deg - e.totalDeg();
1867            ExpVector f = e.extend(1, 0, d);
1868            C.put(f, a);
1869        }
1870        return Cp;
1871    }
1872
1873
1874    /**
1875     * Dehomogenize.
1876     * @param pfac contracted polynomial ring factory (by 1 variable).
1877     * @return in homogeneous polynomial.
1878     */
1879    public GenPolynomial<C> deHomogenize(GenPolynomialRing<C> pfac) {
1880        if (ring.equals(pfac)) { // not implemented
1881            throw new UnsupportedOperationException("case with same ring not implemented");
1882        }
1883        GenPolynomial<C> Cp = pfac.getZERO().copy();
1884        if (this.isZERO()) {
1885            return Cp;
1886        }
1887        Map<ExpVector, C> C = Cp.val; //getMap();
1888        Map<ExpVector, C> A = val;
1889        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1890            ExpVector e = y.getKey();
1891            C a = y.getValue();
1892            ExpVector f = e.contract(1, pfac.nvar);
1893            C.put(f, a);
1894        }
1895        return Cp;
1896    }
1897
1898
1899    /**
1900     * Reverse variables. Used e.g. in opposite rings.
1901     * @return polynomial with reversed variables.
1902     */
1903    public GenPolynomial<C> reverse(GenPolynomialRing<C> oring) {
1904        GenPolynomial<C> Cp = oring.getZERO().copy();
1905        if (this.isZERO()) {
1906            return Cp;
1907        }
1908        int k = -1;
1909        if (oring.tord.getEvord2() != 0 && oring.partial) {
1910            k = oring.tord.getSplit();
1911        }
1912
1913        Map<ExpVector, C> C = Cp.val; //getMap();
1914        Map<ExpVector, C> A = val;
1915        ExpVector f;
1916        for (Map.Entry<ExpVector, C> y : A.entrySet()) {
1917            ExpVector e = y.getKey();
1918            if (k >= 0) {
1919                f = e.reverse(k);
1920            } else {
1921                f = e.reverse();
1922            }
1923            C a = y.getValue();
1924            C.put(f, a);
1925        }
1926        return Cp;
1927    }
1928
1929
1930    /**
1931     * Iterator over coefficients.
1932     * @return val.values().iterator().
1933     */
1934    public Iterator<C> coefficientIterator() {
1935        return val.values().iterator();
1936    }
1937
1938
1939    /**
1940     * Iterator over exponents.
1941     * @return val.keySet().iterator().
1942     */
1943    public Iterator<ExpVector> exponentIterator() {
1944        return val.keySet().iterator();
1945    }
1946
1947
1948    /**
1949     * Iterator over monomials.
1950     * @return a PolyIterator.
1951     */
1952    public Iterator<Monomial<C>> iterator() {
1953        return new PolyIterator<C>(val);
1954    }
1955
1956
1957    /**
1958     * Map a unary function to the coefficients.
1959     * @param f evaluation functor.
1960     * @return new polynomial with coefficients f(this(e)).
1961     */
1962    public GenPolynomial<C> map(final UnaryFunctor<? super C, C> f) {
1963        GenPolynomial<C> n = ring.getZERO().copy();
1964        SortedMap<ExpVector, C> nv = n.val;
1965        for (Monomial<C> m : this) {
1966            //logger.info("m = " + m);
1967            C c = f.eval(m.c);
1968            if (c != null && !c.isZERO()) {
1969                nv.put(m.e, c);
1970            }
1971        }
1972        return n;
1973    }
1974
1975}