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