001/*
002 * $Id$
003 */
004
005package edu.jas.ps;
006
007
008import java.util.BitSet;
009import java.util.HashMap;
010import java.util.HashSet;
011import java.util.Iterator;
012import java.util.List;
013import java.util.Map;
014import java.util.Set;
015import java.util.TreeMap;
016
017import edu.jas.poly.AlgebraicNumber;
018import edu.jas.poly.ExpVector;
019import edu.jas.poly.GenPolynomial;
020import edu.jas.structure.BinaryFunctor;
021import edu.jas.structure.RingElem;
022import edu.jas.structure.Selector;
023import edu.jas.structure.UnaryFunctor;
024import edu.jas.util.MapEntry;
025
026
027/**
028 * Multivariate power series implementation. Uses inner classes and lazy
029 * evaluated generating function for coefficients. All ring element methods use
030 * lazy evaluation except where noted otherwise. Eager evaluated methods are
031 * <code>toString()</code>, <code>compareTo()</code>, <code>equals()</code>,
032 * <code>evaluate()</code>, or methods which use the <code>order()</code> or
033 * <code>orderExpVector()</code> methods, like <code>signum()</code>,
034 * <code>abs()</code>, <code>divide()</code>, <code>remainder()</code> and
035 * <code>gcd()</code>. <b>Note: </b> Currently the term order is fixed to the
036 * order defined by the iterator over exponent vectors in class
037 * <code>ExpVectorIterator</code>.
038 * @param <C> ring element type
039 * @author Heinz Kredel
040 */
041
042public class MultiVarPowerSeries<C extends RingElem<C>> implements RingElem<MultiVarPowerSeries<C>> {
043
044
045    /**
046     * Power series ring factory.
047     */
048    public final MultiVarPowerSeriesRing<C> ring;
049
050
051    /**
052     * Data structure / generating function for coefficients. Cannot be final
053     * because of fixPoint, must be accessible in factory.
054     */
055    /*package*/MultiVarCoefficients<C> lazyCoeffs;
056
057
058    /**
059     * Truncation of computations.
060     */
061    private int truncate;
062
063
064    /**
065     * Order of power series.
066     */
067    private int order = -1; // == unknown
068
069
070    /**
071     * ExpVector of order of power series.
072     */
073    private ExpVector evorder = null; // == unknown
074
075
076    /**
077     * Private constructor.
078     */
079    @SuppressWarnings("unused")
080    private MultiVarPowerSeries() {
081        throw new IllegalArgumentException("do not use no-argument constructor");
082    }
083
084
085    /**
086     * Package constructor. Use in fixPoint only, must be accessible in factory.
087     * @param ring power series ring.
088     */
089    /*package*/ MultiVarPowerSeries(MultiVarPowerSeriesRing<C> ring) {
090        this.ring = ring;
091        this.lazyCoeffs = null;
092        this.truncate = ring.truncate;
093    }
094
095
096    /**
097     * Constructor.
098     * @param ring power series ring.
099     * @param lazyCoeffs generating function for coefficients.
100     */
101    public MultiVarPowerSeries(MultiVarPowerSeriesRing<C> ring, MultiVarCoefficients<C> lazyCoeffs) {
102        this(ring, lazyCoeffs, ring.truncate);
103    }
104
105
106    /**
107     * Constructor.
108     * @param ring power series ring.
109     * @param lazyCoeffs generating function for coefficients.
110     * @param trunc truncate parameter for this power series.
111     */
112    public MultiVarPowerSeries(MultiVarPowerSeriesRing<C> ring, MultiVarCoefficients<C> lazyCoeffs,
113                    int trunc) {
114        if (lazyCoeffs == null || ring == null) {
115            throw new IllegalArgumentException(
116                            "null not allowed: ring = " + ring + ", lazyCoeffs = " + lazyCoeffs);
117        }
118        this.ring = ring;
119        this.lazyCoeffs = lazyCoeffs;
120        this.truncate = Math.min(trunc, ring.truncate);
121    }
122
123
124    /**
125     * Get the corresponding element factory.
126     * @return factory for this Element.
127     * @see edu.jas.structure.Element#factory()
128     */
129    public MultiVarPowerSeriesRing<C> factory() {
130        return ring;
131    }
132
133
134    /**
135     * Clone this power series.
136     * @see java.lang.Object#clone()
137     */
138    @Override
139    public MultiVarPowerSeries<C> copy() {
140        return new MultiVarPowerSeries<C>(ring, lazyCoeffs);
141    }
142
143
144    /**
145     * String representation of power series.
146     * @see java.lang.Object#toString()
147     */
148    @Override
149    public String toString() {
150        return toString(truncate);
151    }
152
153
154    /**
155     * To String with given truncate.
156     * @param trunc truncate parameter for this power series.
157     * @return string representation of this to given truncate.
158     */
159    public String toString(int trunc) {
160        StringBuffer sb = new StringBuffer();
161        MultiVarPowerSeries<C> s = this;
162        String[] vars = ring.vars;
163        //System.out.println("cache1 = " + s.lazyCoeffs.coeffCache);
164        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, trunc)) {
165            C c = s.coefficient(i);
166            //System.out.println("i = " + i + ", c = " +c);
167            int si = c.signum();
168            if (si != 0) {
169                if (si > 0) {
170                    if (sb.length() > 0) {
171                        sb.append(" + ");
172                    }
173                } else {
174                    c = c.negate();
175                    sb.append(" - ");
176                }
177                if (!c.isONE() || i.isZERO()) {
178                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
179                        sb.append("{ ");
180                    }
181                    sb.append(c.toString());
182                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
183                        sb.append(" }");
184                    }
185                    if (!i.isZERO()) {
186                        sb.append(" * ");
187                    }
188                }
189                if (i.isZERO()) {
190                    //skip; sb.append(" ");
191                } else {
192                    sb.append(i.toString(vars));
193                }
194                //sb.append(c.toString() + ", ");
195            }
196            //System.out.println("cache = " + s.coeffCache);
197        }
198        if (sb.length() == 0) {
199            sb.append("0");
200        }
201        sb.append(" + BigO( (" + ring.varsToString() + ")^" + (trunc + 1) + "(" + (ring.truncate + 1)
202                        + ") )");
203        //sb.append("...");
204        //System.out.println("cache2 = " + s.lazyCoeffs.coeffCache);
205        return sb.toString();
206    }
207
208
209    /**
210     * Get a scripting compatible string representation.
211     * @return script compatible representation for this Element.
212     * @see edu.jas.structure.Element#toScript()
213     */
214    @Override
215    public String toScript() {
216        // Python case
217        StringBuffer sb = new StringBuffer("");
218        MultiVarPowerSeries<C> s = this;
219        String[] vars = ring.vars;
220        //System.out.println("cache = " + s.coeffCache);
221        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, truncate)) {
222            C c = s.coefficient(i);
223            int si = c.signum();
224            if (si != 0) {
225                if (si > 0) {
226                    if (sb.length() > 0) {
227                        sb.append(" + ");
228                    }
229                } else {
230                    c = c.negate();
231                    sb.append(" - ");
232                }
233                if (!c.isONE() || i.isZERO()) {
234                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
235                        sb.append("( ");
236                    }
237                    sb.append(c.toScript());
238                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
239                        sb.append(" )");
240                    }
241                    if (!i.isZERO()) {
242                        sb.append(" * ");
243                    }
244                }
245                if (i.isZERO()) {
246                    //skip; sb.append(" ");
247                } else {
248                    sb.append(i.toScript(vars));
249                }
250                //sb.append(c.toString() + ", ");
251            }
252            //System.out.println("cache = " + s.coeffCache);
253        }
254        if (sb.length() == 0) {
255            sb.append("0");
256        }
257        sb.append(" + BigO( (" + ring.varsToString() + ")**" + (truncate + 1) + " )");
258        // sb.append("," + truncate + "");
259        return sb.toString();
260    }
261
262
263    /**
264     * Get a scripting compatible string representation of the factory.
265     * @return script compatible representation for this ElemFactory.
266     * @see edu.jas.structure.Element#toScriptFactory()
267     */
268    @Override
269    public String toScriptFactory() {
270        // Python case
271        return factory().toScript();
272    }
273
274
275    /**
276     * Get coefficient.
277     * @param index number of requested coefficient.
278     * @return coefficient at index.
279     */
280    public C coefficient(ExpVector index) {
281        if (index == null) {
282            throw new IndexOutOfBoundsException("null index not allowed");
283        }
284        return lazyCoeffs.get(index);
285    }
286
287
288    /**
289     * Homogeneous part.
290     * @param tdeg requested degree.
291     * @return polynomial part of given degree.
292     */
293    public GenPolynomial<C> homogeneousPart(long tdeg) {
294        return lazyCoeffs.getHomPart(tdeg);
295    }
296
297
298    /**
299     * Get a GenPolynomial&lt;C&gt; from this.
300     * @return a GenPolynomial&lt;C&gt; from this up to truncate homogeneous
301     *         parts.
302     */
303    public GenPolynomial<C> asPolynomial() {
304        GenPolynomial<C> p = homogeneousPart(0L);
305        for (int i = 1; i <= truncate; i++) {
306            p = p.sum(homogeneousPart(i));
307        }
308        return p;
309    }
310
311
312    /**
313     * Leading base coefficient.
314     * @return first coefficient.
315     */
316    public C leadingCoefficient() {
317        return coefficient(ring.EVZERO);
318    }
319
320
321    /**
322     * Prepend a new leading coefficient.
323     * @param r variable for the direction.
324     * @param h new coefficient.
325     * @return new power series.
326     */
327    public MultiVarPowerSeries<C> prepend(final C h, final int r) {
328        if (r < 0 || ring.nvar < r) {
329            throw new IllegalArgumentException("variable index out of bound");
330        }
331        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
332
333
334            @Override
335            public C generate(ExpVector i) {
336                if (i.isZERO()) {
337                    return h;
338                }
339                ExpVector e = i.subst(r, i.getVal(r) - 1);
340                if (e.signum() < 0) {
341                    return pfac.coFac.getZERO();
342                }
343                return coefficient(e);
344            }
345        });
346    }
347
348
349    /**
350     * Shift coefficients.
351     * @param k shift index.
352     * @param r variable for the direction.
353     * @return new power series with coefficient(i) = old.coefficient(i-k).
354     */
355    public MultiVarPowerSeries<C> shift(final int k, final int r) {
356        if (r < 0 || ring.nvar < r) {
357            throw new IllegalArgumentException("variable index out of bound");
358        }
359        int nt = Math.min(truncate + k, ring.truncate);
360        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
361
362
363            @Override
364            public C generate(ExpVector i) {
365                long d = i.getVal(r);
366                if (d - k < 0) {
367                    return ring.coFac.getZERO();
368                }
369                ExpVector e = i.subst(r, i.getVal(r) - k);
370                return coefficient(e);
371            }
372        }, nt);
373    }
374
375
376    /**
377     * Reductum.
378     * @param r variable for taking the reductum.
379     * @return this - leading monomial in the direction of r.
380     */
381    public MultiVarPowerSeries<C> reductum(final int r) {
382        if (r < 0 || ring.nvar < r) {
383            throw new IllegalArgumentException("variable index out of bound");
384        }
385        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
386
387
388            @Override
389            public C generate(ExpVector i) {
390                ExpVector e = i.subst(r, i.getVal(r) + 1);
391                return coefficient(e);
392            }
393        });
394    }
395
396
397    /**
398     * Reductum.
399     * @return this - leading monomial.
400     */
401    public MultiVarPowerSeries<C> reductum() {
402        Map.Entry<ExpVector, C> m = orderMonomial();
403        if (m == null) {
404            return ring.getZERO();
405        }
406        ExpVector e = m.getKey();
407        long d = e.totalDeg();
408        MultiVarCoefficients<C> mc = lazyCoeffs;
409        HashMap<Long, GenPolynomial<C>> cc = new HashMap<Long, GenPolynomial<C>>(mc.coeffCache);
410        GenPolynomial<C> p = cc.get(d);
411        if (p != null && !p.isZERO()) {
412            p = p.subtract(m.getValue(), e); // p contains this term after orderMonomial()
413            cc.put(d, p);
414        }
415        HashSet<ExpVector> z = new HashSet<ExpVector>(mc.zeroCache);
416        if (!mc.homCheck.get((int) d)) {
417            z.add(e);
418            //System.out.println("e = " + e);
419        }
420
421        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(mc.pfac, cc, z, mc.homCheck) {
422
423
424            @Override
425            public C generate(ExpVector i) {
426                return coefficient(i);
427            }
428        });
429    }
430
431
432    /**
433     * Shift coefficients. Multiply by exponent vector.
434     * @param k shift ExpVector.
435     * @return new power series with coefficient(i) = old.coefficient(i-k).
436     */
437    public MultiVarPowerSeries<C> shift(final ExpVector k) {
438        if (k == null) {
439            throw new IllegalArgumentException("null ExpVector not allowed");
440        }
441        if (k.signum() == 0) {
442            return this;
443        }
444        int nt = Math.min(truncate + (int) k.totalDeg(), ring.truncate);
445        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
446
447
448            @Override
449            public C generate(ExpVector i) {
450                ExpVector d = i.subtract(k);
451                if (d.signum() < 0) {
452                    return ring.coFac.getZERO();
453                }
454                return coefficient(d);
455            }
456        }, nt);
457    }
458
459
460    /**
461     * Multiply by exponent vector and coefficient.
462     * @param k shift ExpVector.
463     * @param c coefficient multiplier.
464     * @return new power series with coefficient(i) = old.coefficient(i-k)*c.
465     */
466    public MultiVarPowerSeries<C> multiply(final C c, final ExpVector k) {
467        if (k == null) {
468            throw new IllegalArgumentException("null ExpVector not allowed");
469        }
470        if (k.signum() == 0) {
471            return this.multiply(c);
472        }
473        if (c.signum() == 0) {
474            return ring.getZERO();
475        }
476        int nt = Math.min(ring.truncate, truncate + (int) k.totalDeg());
477
478        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
479
480
481            @Override
482            public C generate(ExpVector i) {
483                ExpVector d = i.subtract(k);
484                if (d.signum() < 0) {
485                    return ring.coFac.getZERO();
486                }
487                long tdegd = d.totalDeg();
488                if (lazyCoeffs.homCheck.get((int) tdegd)) {
489                    GenPolynomial<C> p = homogeneousPart(tdegd).multiply(c, k);
490                    long tdegi = i.totalDeg();
491                    coeffCache.put(tdegi, p); // overwrite
492                    homCheck.set((int) tdegi);
493                    C b = p.coefficient(i);
494                    //System.out.println("b = " + b + ", i = " + i + ", tdegi = " + tdegi+ ", tdegd = " + tdegd);
495                    //System.out.println("p = " + p + ", i = " + i);
496                    return b;
497                }
498                return coefficient(d).multiply(c);
499            }
500        }, nt);
501    }
502
503
504    /**
505     * Sum monomial.
506     * @param m ExpVector , coefficient pair
507     * @return this + ONE.multiply(m.coefficient,m.exponent).
508     */
509    public MultiVarPowerSeries<C> sum(Map.Entry<ExpVector, C> m) {
510        if (m == null) {
511            throw new IllegalArgumentException("null Map.Entry not allowed");
512        }
513        return sum(m.getValue(), m.getKey());
514    }
515
516
517    /**
518     * Sum exponent vector and coefficient.
519     * @param k ExpVector.
520     * @param c coefficient.
521     * @return this + ONE.multiply(c,k).
522     */
523    public MultiVarPowerSeries<C> sum(final C c, final ExpVector k) {
524        if (k == null) {
525            throw new IllegalArgumentException("null ExpVector not allowed");
526        }
527        if (c.signum() == 0) {
528            return this;
529        }
530        long d = k.totalDeg();
531        MultiVarCoefficients<C> mc = lazyCoeffs;
532        HashMap<Long, GenPolynomial<C>> cc = new HashMap<Long, GenPolynomial<C>>(mc.coeffCache);
533        GenPolynomial<C> p = cc.get(d);
534        if (p == null) {
535            p = mc.pfac.getZERO();
536        }
537        p = p.sum(c, k);
538        //System.out.println("p = " + p);
539        cc.put(d, p);
540        HashSet<ExpVector> z = new HashSet<ExpVector>(mc.zeroCache);
541        //System.out.println("z = " + z);
542        if (p.coefficient(k).isZERO() && !mc.homCheck.get((int) d)) {
543            z.add(k);
544        }
545
546        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(mc.pfac, cc, z, mc.homCheck) {
547
548
549            @Override
550            public C generate(ExpVector i) {
551                return coefficient(i);
552            }
553        });
554    }
555
556
557    /**
558     * Subtract exponent vector and coefficient.
559     * @param k ExpVector.
560     * @param c coefficient.
561     * @return this - ONE.multiply(c,k).
562     */
563    public MultiVarPowerSeries<C> subtract(final C c, final ExpVector k) {
564        if (k == null) {
565            throw new IllegalArgumentException("null ExpVector not allowed");
566        }
567        if (c.signum() == 0) {
568            return this;
569        }
570        long d = k.totalDeg();
571        MultiVarCoefficients<C> mc = lazyCoeffs;
572        HashMap<Long, GenPolynomial<C>> cc = new HashMap<Long, GenPolynomial<C>>(mc.coeffCache);
573        GenPolynomial<C> p = cc.get(d);
574        if (p == null) {
575            p = mc.pfac.getZERO();
576        }
577        p = p.subtract(c, k);
578        cc.put(d, p);
579        HashSet<ExpVector> z = new HashSet<ExpVector>(mc.zeroCache);
580        //System.out.println("z = " + z);
581        if (p.coefficient(k).isZERO() && !mc.homCheck.get((int) d)) {
582            z.add(k);
583        }
584        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(mc.pfac, cc, z, mc.homCheck) {
585
586
587            @Override
588            public C generate(ExpVector i) {
589                return coefficient(i);
590            }
591        });
592    }
593
594
595    /**
596     * Sum exponent vector and coefficient.
597     * @param mvc cached coefficients.
598     * @return this + mvc.
599     */
600    public MultiVarPowerSeries<C> sum(MultiVarCoefficients<C> mvc) {
601        MultiVarCoefficients<C> mc = lazyCoeffs;
602        TreeMap<Long, GenPolynomial<C>> cc = new TreeMap<Long, GenPolynomial<C>>(mc.coeffCache);
603        TreeMap<Long, GenPolynomial<C>> ccv = new TreeMap<Long, GenPolynomial<C>>(mvc.coeffCache);
604        long d1 = (cc.size() > 0 ? cc.lastKey() : 0);
605        long d2 = (ccv.size() > 0 ? ccv.lastKey() : 0);
606        HashSet<ExpVector> z = new HashSet<ExpVector>(mc.zeroCache);
607        z.addAll(mvc.zeroCache);
608        long d = Math.max(d1, d2);
609        BitSet hc = new BitSet((int) d);
610        for (long i = 0; i <= d; i++) {
611            GenPolynomial<C> p1 = cc.get(i);
612            GenPolynomial<C> p2 = mvc.coeffCache.get(i);
613            if (p1 == null) {
614                p1 = mc.pfac.getZERO();
615            }
616            if (p2 == null) {
617                p2 = mc.pfac.getZERO();
618            }
619            GenPolynomial<C> p = p1.sum(p2);
620            //System.out.println("p = " + p);
621            cc.put(i, p);
622            if (mc.homCheck.get((int) i) && mvc.homCheck.get((int) i)) {
623                hc.set((int) i);
624            } else {
625                Set<ExpVector> ev = new HashSet<ExpVector>(p1.getMap().keySet());
626                ev.addAll(p2.getMap().keySet());
627                ev.removeAll(p.getMap().keySet());
628                z.addAll(ev);
629            }
630        }
631        //System.out.println("cc = " + cc);
632
633        return new MultiVarPowerSeries<C>(ring,
634                        new MultiVarCoefficients<C>(mc.pfac, new HashMap<Long, GenPolynomial<C>>(cc), z, hc) {
635
636
637                            @Override
638                            public C generate(ExpVector i) {
639                                return coefficient(i);
640                            }
641                        });
642    }
643
644
645    /**
646     * Select coefficients.
647     * @param sel selector functor.
648     * @return new power series with selected coefficients.
649     */
650    public MultiVarPowerSeries<C> select(final Selector<? super C> sel) {
651        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
652
653
654            @Override
655            public C generate(ExpVector i) {
656                C c = coefficient(i);
657                if (sel.select(c)) {
658                    return c;
659                }
660                return ring.coFac.getZERO();
661            }
662        });
663    }
664
665
666    /**
667     * Shift select coefficients. Not selected coefficients are removed from the
668     * result series.
669     * @param sel selector functor.
670     * @return new power series with shifted selected coefficients.
671     */
672    public MultiVarPowerSeries<C> shiftSelect(final Selector<? super C> sel) {
673        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
674
675
676            ExpVectorIterable ib = new ExpVectorIterable(ring.nvar, true, truncate);
677
678
679            Iterator<ExpVector> pos = ib.iterator();
680
681
682            @Override
683            public C generate(ExpVector i) {
684                C c;
685                if (i.signum() > 0) {
686                    int[] deps = i.dependencyOnVariables();
687                    ExpVector x = i.subst(deps[0], i.getVal(deps[0]) - 1L);
688                    c = get(x); // ensure all coefficients are generated
689                }
690                do {
691                    c = null;
692                    if (pos.hasNext()) {
693                        ExpVector e = pos.next();
694                        c = coefficient(e);
695                    } else {
696                        break;
697                    }
698                } while (!sel.select(c));
699                if (c == null) { // not correct because not known
700                    c = ring.coFac.getZERO();
701                }
702                return c;
703            }
704        });
705    }
706
707
708    /**
709     * Map a unary function to this power series.
710     * @param f evaluation functor.
711     * @return new power series with coefficients f(this(i)).
712     */
713    public MultiVarPowerSeries<C> map(final UnaryFunctor<? super C, C> f) {
714        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
715
716
717            @Override
718            public C generate(ExpVector i) {
719                return f.eval(coefficient(i));
720            }
721        });
722    }
723
724
725    /**
726     * Map a binary function to this and another power series.
727     * @param f evaluation functor with coefficients f(this(i),other(i)).
728     * @param ps other power series.
729     * @return new power series.
730     */
731    public MultiVarPowerSeries<C> zip(final BinaryFunctor<? super C, ? super C, C> f,
732                    final MultiVarPowerSeries<C> ps) {
733        int m = Math.min(ring.truncate, Math.max(truncate, ps.truncate()));
734
735        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
736
737
738            @Override
739            public C generate(ExpVector i) {
740                return f.eval(coefficient(i), ps.coefficient(i));
741            }
742        }, m);
743    }
744
745
746    /**
747     * Sum of two power series, using zip().
748     * @param ps other power series.
749     * @return this + ps.
750     */
751    public MultiVarPowerSeries<C> sumZip(MultiVarPowerSeries<C> ps) {
752        return zip(new BinaryFunctor<C, C, C>() {
753
754
755            @Override
756            public C eval(C c1, C c2) {
757                return c1.sum(c2);
758            }
759        }, ps);
760    }
761
762
763    /**
764     * Subtraction of two power series, using zip().
765     * @param ps other power series.
766     * @return this - ps.
767     */
768    public MultiVarPowerSeries<C> subtractZip(MultiVarPowerSeries<C> ps) {
769        return zip(new BinaryFunctor<C, C, C>() {
770
771
772            @Override
773            public C eval(C c1, C c2) {
774                return c1.subtract(c2);
775            }
776        }, ps);
777    }
778
779
780    /**
781     * Multiply by coefficient.
782     * @param a coefficient.
783     * @return this * a.
784     */
785    public MultiVarPowerSeries<C> multiply(final C a) {
786        if (a.isZERO()) {
787            return ring.getZERO();
788        }
789        if (a.isONE()) {
790            return this;
791        }
792        return map(new UnaryFunctor<C, C>() {
793
794
795            @Override
796            public C eval(C c) {
797                return c.multiply(a);
798            }
799        });
800    }
801
802
803    /**
804     * Monic.
805     * @return 1/orderCoeff() * this.
806     */
807    public MultiVarPowerSeries<C> monic() {
808        ExpVector e = orderExpVector();
809        if (e == null) {
810            return this;
811        }
812        C a = coefficient(e);
813        if (a.isONE()) {
814            return this;
815        }
816        if (a.isZERO()) { // sic
817            return this;
818        }
819        final C b = a.inverse();
820        return map(new UnaryFunctor<C, C>() {
821
822
823            @Override
824            public C eval(C c) {
825                return b.multiply(c);
826            }
827        });
828    }
829
830
831    /**
832     * Negate.
833     * @return - this.
834     */
835    public MultiVarPowerSeries<C> negate() {
836        return map(new UnaryFunctor<C, C>() {
837
838
839            @Override
840            public C eval(C c) {
841                return c.negate();
842            }
843        });
844    }
845
846
847    /**
848     * Absolute value.
849     * @return abs(this).
850     */
851    public MultiVarPowerSeries<C> abs() {
852        if (signum() < 0) {
853            return negate();
854        }
855        return this;
856    }
857
858
859    /**
860     * Evaluate at given point.
861     * @return ps(a).
862     */
863    public C evaluate(List<C> a) {
864        C v = ring.coFac.getZERO();
865        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, truncate)) {
866            C c = coefficient(i);
867            if (c.isZERO()) {
868                continue;
869            }
870            c = c.multiply(i.evaluate(ring.coFac, a));
871            v = v.sum(c);
872        }
873        return v;
874    }
875
876
877    /**
878     * Order.
879     * @return index of first non zero coefficient.
880     */
881    public int order() {
882        if (order >= 0) {
883            return order;
884        }
885        // must compute it
886        GenPolynomial<C> p = null;
887        int t = 0;
888        while (lazyCoeffs.homCheck.get(t)) {
889            p = lazyCoeffs.coeffCache.get((long) t);
890            if (p == null || p.isZERO()) { // ??
891                t++;
892                continue;
893            }
894            order = t;
895            evorder = p.trailingExpVector();
896            //System.out.println("order = " + t);
897            return order;
898        }
899        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, truncate)) {
900            if (!coefficient(i).isZERO()) {
901                order = (int) i.totalDeg(); //ord;
902                evorder = i;
903                //System.out.println("order = " + order + ", evorder = " + evorder);
904                return order;
905            }
906        }
907        order = truncate + 1;
908        // evorder is null
909        return order;
910    }
911
912
913    /**
914     * Order ExpVector.
915     * @return ExpVector of first non zero coefficient.
916     */
917    public ExpVector orderExpVector() {
918        //int x = 
919        order(); // ensure evorder is set
920        return evorder;
921    }
922
923
924    /**
925     * Order monomial.
926     * @return first map entry or null.
927     */
928    public Map.Entry<ExpVector, C> orderMonomial() {
929        ExpVector e = orderExpVector();
930        if (e == null) {
931            return null;
932        }
933        //JAVA6only:
934        //return new AbstractMap.SimpleImmutableEntry<ExpVector, C>(e, coefficient(e));
935        return new MapEntry<ExpVector, C>(e, coefficient(e));
936    }
937
938
939    /**
940     * Truncate.
941     * @return truncate index of power series.
942     */
943    public int truncate() {
944        return truncate;
945    }
946
947
948    /**
949     * Set truncate.
950     * @param t new truncate index.
951     * @return old truncate index of power series.
952     */
953    public int setTruncate(int t) {
954        if (t < 0) {
955            throw new IllegalArgumentException("negative truncate not allowed");
956        }
957        int ot = truncate;
958        if (order >= 0) {
959            if (order > truncate) {
960                order = -1; // reset
961                evorder = null;
962            }
963        }
964        truncate = t;
965        return ot;
966    }
967
968
969    /**
970     * Ecart.
971     * @return ecart.
972     */
973    public long ecart() {
974        ExpVector e = orderExpVector();
975        if (e == null) {
976            return 0L;
977        }
978        long d = e.totalDeg();
979        long hd = d;
980        for (long i = d + 1L; i <= truncate; i++) {
981            if (!homogeneousPart(i).isZERO()) {
982                hd = i;
983            }
984        }
985        //System.out.println("d = " + d + ", hd = " + hd + ", coeffCache = " + lazyCoeffs.coeffCache + ", this = " + this);
986        return hd - d;
987    }
988
989
990    /**
991     * Signum.
992     * @return sign of first non zero coefficient.
993     */
994    public int signum() {
995        ExpVector ev = orderExpVector();
996        if (ev != null) {
997            return coefficient(ev).signum();
998        }
999        return 0;
1000    }
1001
1002
1003    /**
1004     * Compare to. <b>Note: </b> compare only up to max(truncates).
1005     * @return sign of first non zero coefficient of this-ps.
1006     */
1007    @Override
1008    public int compareTo(MultiVarPowerSeries<C> ps) {
1009        final int m = truncate();
1010        final int n = ps.truncate();
1011        final int pos = Math.min(ring.truncate, Math.min(m, n));
1012        int s = 0;
1013        //System.out.println("coeffCache_c1 = " + lazyCoeffs.coeffCache);
1014        //System.out.println("coeffCache_c2 = " + ps.lazyCoeffs.coeffCache);
1015        // test homogeneous parts first is slower
1016        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, pos)) {
1017            s = coefficient(i).compareTo(ps.coefficient(i));
1018            if (s != 0) {
1019                //System.out.println("i = " + i + ", coeff = " + coefficient(i) + ", ps.coeff = " + ps.coefficient(i));
1020                return s;
1021            }
1022        }
1023        for (int j = pos + 1; j <= Math.min(ring.truncate, Math.max(m, n)); j++) {
1024            for (ExpVector i : new ExpVectorIterable(ring.nvar, j)) {
1025                s = coefficient(i).compareTo(ps.coefficient(i));
1026                //System.out.println("i = " + i + ", coeff = " + coefficient(i) + ", ps.coeff = " + ps.coefficient(i));
1027                if (s != 0) {
1028                    //System.out.println("i = " + i + ", coeff = " + coefficient(i) + ", ps.coeff = " + ps.coefficient(i));
1029                    return s;
1030                }
1031            }
1032        }
1033        return s;
1034    }
1035
1036
1037    /**
1038     * Is power series zero. <b>Note: </b> compare only up to truncate.
1039     * @return If this is 0 then true is returned, else false.
1040     * @see edu.jas.structure.RingElem#isZERO()
1041     */
1042    public boolean isZERO() {
1043        return (signum() == 0);
1044    }
1045
1046
1047    /**
1048     * Is power series one. <b>Note: </b> compare only up to truncate.
1049     * @return If this is 1 then true is returned, else false.
1050     * @see edu.jas.structure.RingElem#isONE()
1051     */
1052    public boolean isONE() {
1053        if (!leadingCoefficient().isONE()) {
1054            return false;
1055        }
1056        return (compareTo(ring.ONE) == 0);
1057        //return reductum().isZERO();
1058    }
1059
1060
1061    /**
1062     * Comparison with any other object. <b>Note: </b> compare only up to
1063     * truncate.
1064     * @see java.lang.Object#equals(java.lang.Object)
1065     */
1066    @Override
1067    @SuppressWarnings("unchecked")
1068    public boolean equals(Object B) {
1069        if (B == null) {
1070            return false;
1071        }
1072        if (!(B instanceof MultiVarPowerSeries)) {
1073            return false;
1074        }
1075        MultiVarPowerSeries<C> a = (MultiVarPowerSeries<C>) B;
1076        return compareTo(a) == 0;
1077    }
1078
1079
1080    /**
1081     * Hash code for this polynomial. <b>Note: </b> only up to truncate.
1082     * @see java.lang.Object#hashCode()
1083     */
1084    @Override
1085    public int hashCode() {
1086        int h = 0;
1087        //h = ( ring.hashCode() << 23 );
1088        //h += val.hashCode();
1089        for (ExpVector i : new ExpVectorIterable(ring.nvar, true, truncate)) {
1090            C c = coefficient(i);
1091            //System.out.print(" i = " + i + " c = " + c);
1092            //System.out.println(" #i = " + i.hashCode() + " #c = " + c.hashCode() + "# = " + h);
1093            if (!c.isZERO()) {
1094                h += i.hashCode();
1095                h = Math.abs(h << 1);
1096            }
1097            h += c.hashCode();
1098            h = Math.abs(h); // << 3);
1099        }
1100        return h;
1101    }
1102
1103
1104    /**
1105     * Is unit.
1106     * @return true, if this power series is invertible, else false.
1107     */
1108    public boolean isUnit() {
1109        return leadingCoefficient().isUnit();
1110    }
1111
1112
1113    /**
1114     * Sum a another power series.
1115     * @param ps other power series.
1116     * @return this + ps.
1117     */
1118    public MultiVarPowerSeries<C> sum(final MultiVarPowerSeries<C> ps) {
1119        //final MultiVarPowerSeries<C> ps1 = this; // method name was ambiguous in generate
1120        int nt = Math.min(ring.truncate, Math.max(truncate(), ps.truncate()));
1121
1122        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1123
1124
1125            @Override
1126            public C generate(ExpVector e) {
1127                long tdeg = e.totalDeg();
1128                if (lazyCoeffs.homCheck.get((int) tdeg)) {
1129                    // generate respective homogeneous polynomial
1130                    GenPolynomial<C> p = homogeneousPart(tdeg).sum(ps.homogeneousPart(tdeg));
1131                    coeffCache.put(tdeg, p); // overwrite
1132                    homCheck.set((int) tdeg);
1133                    C c = p.coefficient(e);
1134                    //System.out.println("c = " + c + ", e = " + e + ", tdeg = " + tdeg);
1135                    return c;
1136                }
1137                return coefficient(e).sum(ps.coefficient(e));
1138            }
1139        }, nt);
1140    }
1141
1142
1143    /**
1144     * Subtract a another power series.
1145     * @param ps other power series.
1146     * @return this - ps.
1147     */
1148    public MultiVarPowerSeries<C> subtract(final MultiVarPowerSeries<C> ps) {
1149        //final MultiVarPowerSeries<C> ps1 = this; // method name was ambiguous in generate
1150        int nt = Math.min(ring.truncate, Math.max(truncate(), ps.truncate()));
1151
1152        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1153
1154
1155            @Override
1156            public C generate(ExpVector e) {
1157                long tdeg = e.totalDeg();
1158                if (lazyCoeffs.homCheck.get((int) tdeg)) {
1159                    // generate respective homogeneous polynomial
1160                    GenPolynomial<C> p = homogeneousPart(tdeg).subtract(ps.homogeneousPart(tdeg));
1161                    coeffCache.put(tdeg, p); // overwrite
1162                    homCheck.set((int) tdeg);
1163                    C c = p.coefficient(e);
1164                    //System.out.println("p = " + p + ", e = " + e + ", tdeg = " + tdeg);
1165                    return c;
1166                }
1167                return coefficient(e).subtract(ps.coefficient(e));
1168            }
1169        }, nt);
1170    }
1171
1172
1173    /**
1174     * Multiply by another power series.
1175     * @param ps other power series.
1176     * @return this * ps.
1177     */
1178    public MultiVarPowerSeries<C> multiply(final MultiVarPowerSeries<C> ps) {
1179        //final MultiVarPowerSeries<C> ps1 = this; // method name was ambiguous in generate
1180        int nt = Math.min(ring.truncate, truncate() + ps.truncate());
1181
1182        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1183
1184
1185            @Override
1186            public C generate(ExpVector e) {
1187                long tdeg = e.totalDeg();
1188                // generate respective homogeneous polynomial
1189                GenPolynomial<C> p = null; //fac.getZERO();
1190                for (int k = 0; k <= tdeg; k++) {
1191                    GenPolynomial<C> m = homogeneousPart(k).multiply(ps.homogeneousPart(tdeg - k));
1192                    if (p == null) {
1193                        p = m;
1194                    } else {
1195                        p = p.sum(m);
1196                    }
1197                }
1198                coeffCache.put(tdeg, p); // overwrite
1199                homCheck.set((int) tdeg);
1200                C c = p.coefficient(e);
1201                return c;
1202            }
1203        }, nt);
1204    }
1205
1206
1207    /**
1208     * Inverse power series.
1209     * @return ps with this * ps = 1.
1210     */
1211    public MultiVarPowerSeries<C> inverse() {
1212        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1213
1214
1215            @Override
1216            public C generate(ExpVector e) {
1217                long tdeg = e.totalDeg();
1218                C d = leadingCoefficient().inverse(); // may fail
1219                if (tdeg == 0) {
1220                    return d;
1221                }
1222                GenPolynomial<C> p = null; //fac.getZERO();
1223                for (int k = 0; k < tdeg; k++) {
1224                    GenPolynomial<C> m = getHomPart(k).multiply(homogeneousPart(tdeg - k));
1225                    if (p == null) {
1226                        p = m;
1227                    } else {
1228                        p = p.sum(m);
1229                    }
1230                }
1231                p = p.multiply(d.negate());
1232                //System.out.println("tdeg = " + tdeg + ", p = " + p);
1233                coeffCache.put(tdeg, p); // overwrite
1234                homCheck.set((int) tdeg);
1235                C c = p.coefficient(e);
1236                return c;
1237            }
1238        });
1239    }
1240
1241
1242    /**
1243     * Divide by another power series.
1244     * @param ps nonzero power series with invertible coefficient.
1245     * @return this / ps.
1246     */
1247    public MultiVarPowerSeries<C> divide(MultiVarPowerSeries<C> ps) {
1248        if (ps.isUnit()) {
1249            return multiply(ps.inverse());
1250        }
1251        int m = order();
1252        int n = ps.order();
1253        if (m < n) {
1254            return ring.getZERO();
1255        }
1256        ExpVector em = orderExpVector();
1257        ExpVector en = ps.orderExpVector();
1258        if (!ps.coefficient(en).isUnit()) {
1259            throw new ArithmeticException("division by non unit coefficient " + ps.coefficient(ps.evorder)
1260                            + ", evorder = " + ps.evorder);
1261        }
1262        // now m >= n
1263        MultiVarPowerSeries<C> st, sps, q, sq;
1264        if (m == 0) {
1265            st = this;
1266        } else {
1267            st = this.shift(em.negate());
1268        }
1269        if (n == 0) {
1270            sps = ps;
1271        } else {
1272            sps = ps.shift(en.negate());
1273        }
1274        q = st.multiply(sps.inverse());
1275        sq = q.shift(em.subtract(en));
1276        return sq;
1277    }
1278
1279
1280    /**
1281     * Power series remainder.
1282     * @param ps nonzero power series with invertible leading coefficient.
1283     * @return remainder with this = quotient * ps + remainder.
1284     */
1285    public MultiVarPowerSeries<C> remainder(MultiVarPowerSeries<C> ps) {
1286        int m = order();
1287        int n = ps.order();
1288        if (m >= n) {
1289            return ring.getZERO();
1290        }
1291        return this;
1292    }
1293
1294
1295    /**
1296     * Quotient and remainder by division of this by S.
1297     * @param S a MultiVarPowerSeries
1298     * @return [this/S, this - (this/S)*S].
1299     */
1300    @SuppressWarnings("unchecked")
1301    public MultiVarPowerSeries<C>[] quotientRemainder(MultiVarPowerSeries<C> S) {
1302        return new MultiVarPowerSeries[] { divide(S), remainder(S) };
1303    }
1304
1305
1306    /**
1307     * Differentiate with respect to variable r.
1308     * @param r variable for the direction.
1309     * @return differentiate(this).
1310     */
1311    public MultiVarPowerSeries<C> differentiate(final int r) {
1312        if (r < 0 || ring.nvar < r) {
1313            throw new IllegalArgumentException("variable index out of bound");
1314        }
1315        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1316
1317
1318            @Override
1319            public C generate(ExpVector i) {
1320                long d = i.getVal(r);
1321                ExpVector e = i.subst(r, d + 1);
1322                C v = coefficient(e);
1323                v = v.multiply(ring.coFac.fromInteger(d + 1));
1324                return v;
1325            }
1326        });
1327    }
1328
1329
1330    /**
1331     * Integrate with respect to variable r and with given constant.
1332     * @param c integration constant.
1333     * @param r variable for the direction.
1334     * @return integrate(this).
1335     */
1336    public MultiVarPowerSeries<C> integrate(final C c, final int r) {
1337        if (r < 0 || ring.nvar < r) {
1338            throw new IllegalArgumentException("variable index out of bound");
1339        }
1340        int nt = Math.min(ring.truncate, truncate + 1);
1341        return new MultiVarPowerSeries<C>(ring, new MultiVarCoefficients<C>(ring) {
1342
1343
1344            @Override
1345            public C generate(ExpVector i) {
1346                if (i.isZERO()) {
1347                    return c;
1348                }
1349                long d = i.getVal(r);
1350                if (d > 0) {
1351                    ExpVector e = i.subst(r, d - 1);
1352                    C v = coefficient(e);
1353                    v = v.divide(ring.coFac.fromInteger(d));
1354                    return v;
1355                }
1356                return ring.coFac.getZERO();
1357            }
1358        }, nt);
1359    }
1360
1361
1362    /**
1363     * Power series greatest common divisor.
1364     * @param ps power series.
1365     * @return gcd(this,ps).
1366     */
1367    public MultiVarPowerSeries<C> gcd(MultiVarPowerSeries<C> ps) {
1368        if (ps.isZERO()) {
1369            return this;
1370        }
1371        if (this.isZERO()) {
1372            return ps;
1373        }
1374        ExpVector em = orderExpVector();
1375        ExpVector en = ps.orderExpVector();
1376        return ring.getONE().shift(em.gcd(en));
1377    }
1378
1379
1380    /**
1381     * Power series extended greatest common divisor. <b>Note:</b> not
1382     * implemented.
1383     * @param S power series.
1384     * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
1385     */
1386    //SuppressWarnings("unchecked")
1387    public MultiVarPowerSeries<C>[] egcd(MultiVarPowerSeries<C> S) {
1388        throw new UnsupportedOperationException("egcd for power series not implemented");
1389    }
1390
1391}