001/*
002 * $Id$
003 */
004
005package edu.jas.ps;
006
007import java.util.List;
008import java.util.ArrayList;
009
010import edu.jas.poly.AlgebraicNumber;
011import edu.jas.poly.ExpVector;
012import edu.jas.poly.GenPolynomial;
013import edu.jas.vector.GenVector;
014import edu.jas.vector.GenVectorModul;
015import edu.jas.structure.BinaryFunctor;
016import edu.jas.structure.RingElem;
017import edu.jas.structure.Selector;
018import edu.jas.structure.UnaryFunctor;
019
020
021/**
022 * Univariate power series implementation. Uses inner classes and lazy evaluated
023 * generating function for coefficients. All ring element methods use lazy
024 * evaluation except where noted otherwise. Eager evaluated methods are
025 * <code>toString()</code>, <code>compareTo()</code>, <code>equals()</code>,
026 * <code>evaluate()</code>, or they use the <code>order()</code> method, like
027 * <code>signum()</code>, <code>abs()</code>, <code>divide()</code>,
028 * <code>remainder()</code> and <code>gcd()</code>.
029 * @param <C> ring element type
030 * @author Heinz Kredel
031 */
032
033public class UnivPowerSeries<C extends RingElem<C>> implements RingElem<UnivPowerSeries<C>> {
034
035
036    /**
037     * Power series ring factory.
038     */
039    public final UnivPowerSeriesRing<C> ring;
040
041
042    /**
043     * Data structure / generating function for coefficients. Cannot be final
044     * because of fixPoint, must be accessible in factory.
045     */
046    /*package*/Coefficients<C> lazyCoeffs;
047
048
049    /**
050     * Truncation of computations.
051     */
052    private int truncate = 11;
053
054
055    /**
056     * Order of power series.
057     */
058    private int order = -1; // == unknown
059
060
061    /**
062     * Private constructor.
063     */
064    @SuppressWarnings("unused")
065    private UnivPowerSeries() {
066        throw new IllegalArgumentException("do not use no-argument constructor");
067    }
068
069
070    /**
071     * Package constructor. Use in fixPoint only, must be accessible in factory.
072     * @param ring power series ring.
073     */
074    /*package*/ UnivPowerSeries(UnivPowerSeriesRing<C> ring) {
075        this.ring = ring;
076        this.lazyCoeffs = null;
077    }
078
079
080    /**
081     * Constructor.
082     * @param ring power series ring.
083     * @param lazyCoeffs generating function for coefficients.
084     */
085    public UnivPowerSeries(UnivPowerSeriesRing<C> ring, Coefficients<C> lazyCoeffs) {
086        if (lazyCoeffs == null || ring == null) {
087            throw new IllegalArgumentException(
088                            "null not allowed: ring = " + ring + ", lazyCoeffs = " + lazyCoeffs);
089        }
090        this.ring = ring;
091        this.lazyCoeffs = lazyCoeffs;
092        this.truncate = ring.truncate;
093    }
094
095
096    /**
097     * Get the corresponding element factory.
098     * @return factory for this Element.
099     * @see edu.jas.structure.Element#factory()
100     */
101    public UnivPowerSeriesRing<C> factory() {
102        return ring;
103    }
104
105
106    /**
107     * Clone this power series.
108     * @see java.lang.Object#clone()
109     */
110    @Override
111    public UnivPowerSeries<C> copy() {
112        return new UnivPowerSeries<C>(ring, lazyCoeffs);
113    }
114
115
116    /**
117     * String representation of power series.
118     * @see java.lang.Object#toString()
119     */
120    @Override
121    public String toString() {
122        return toString(truncate);
123    }
124
125
126    /**
127     * To String with given truncate.
128     * @return string representation of this to given truncate.
129     */
130    public String toString(int truncate) {
131        StringBuffer sb = new StringBuffer();
132        UnivPowerSeries<C> s = this;
133        String var = ring.var;
134        //System.out.println("cache = " + s.coeffCache);
135        for (int i = 0; i < truncate; i++) {
136            C c = s.coefficient(i);
137            int si = c.signum();
138            if (si != 0) {
139                if (si > 0) {
140                    if (sb.length() > 0) {
141                        sb.append(" + ");
142                    }
143                } else {
144                    c = c.negate();
145                    sb.append(" - ");
146                }
147                if (!c.isONE() || i == 0) {
148                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
149                        sb.append("{ ");
150                    }
151                    sb.append(c.toString());
152                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
153                        sb.append(" }");
154                    }
155                    if (i > 0) {
156                        sb.append(" * ");
157                    }
158                }
159                if (i == 0) {
160                    //skip; sb.append(" ");
161                } else if (i == 1) {
162                    sb.append(var);
163                } else {
164                    sb.append(var + "^" + i);
165                }
166                //sb.append(c.toString() + ", ");
167            }
168            //System.out.println("cache = " + s.coeffCache);
169        }
170        if (sb.length() == 0) {
171            sb.append("0");
172        }
173        sb.append(" + BigO(" + var + "^" + truncate + ")");
174        //sb.append("...");
175        return sb.toString();
176    }
177
178
179    /**
180     * Get a scripting compatible string representation.
181     * @return script compatible representation for this Element.
182     * @see edu.jas.structure.Element#toScript()
183     */
184    @Override
185    public String toScript() {
186        // Python case
187        StringBuffer sb = new StringBuffer("");
188        UnivPowerSeries<C> s = this;
189        String var = ring.var;
190        //System.out.println("cache = " + s.coeffCache);
191        for (int i = 0; i < truncate; i++) {
192            C c = s.coefficient(i);
193            int si = c.signum();
194            if (si != 0) {
195                if (si > 0) {
196                    if (sb.length() > 0) {
197                        sb.append(" + ");
198                    }
199                } else {
200                    c = c.negate();
201                    sb.append(" - ");
202                }
203                if (!c.isONE() || i == 0) {
204                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
205                        sb.append("{ ");
206                    }
207                    sb.append(c.toScript());
208                    if (c instanceof GenPolynomial || c instanceof AlgebraicNumber) {
209                        sb.append(" }");
210                    }
211                    if (i > 0) {
212                        sb.append(" * ");
213                    }
214                }
215                if (i == 0) {
216                    //skip; sb.append(" ");
217                } else if (i == 1) {
218                    sb.append(var);
219                } else {
220                    sb.append(var + "**" + i);
221                }
222                //sb.append(c.toString() + ", ");
223            }
224            //System.out.println("cache = " + s.coeffCache);
225        }
226        if (sb.length() == 0) {
227            sb.append("0");
228        }
229        // sb.append("," + truncate + "");
230        return sb.toString();
231    }
232
233
234    /**
235     * Get a scripting compatible string representation of the factory.
236     * @return script compatible representation for this ElemFactory.
237     * @see edu.jas.structure.Element#toScriptFactory()
238     */
239    @Override
240    public String toScriptFactory() {
241        // Python case
242        return factory().toScript();
243    }
244
245
246    /**
247     * Get coefficient.
248     * @param index number of requested coefficient.
249     * @return coefficient at index.
250     */
251    public C coefficient(int index) {
252        if (index < 0) {
253            throw new IndexOutOfBoundsException("negative index not allowed");
254        }
255        //System.out.println("cache = " + coeffCache);
256        return lazyCoeffs.get(index);
257    }
258
259
260    /**
261     * Get a GenPolynomial&lt;C&gt; from this.
262     * @return a GenPolynomial&lt;C&gt; from this up to truncate parts.
263     */
264    public GenPolynomial<C> asPolynomial() {
265        GenPolynomial<C> p = ring.polyRing().getZERO();
266        for (int i = 0; i < truncate; i++) {
267            C c = coefficient(i);
268            ExpVector e = ExpVector.create(1, 0, i);
269            p = p.sum(c, e);
270        }
271        return p;
272    }
273
274
275    /**
276     * Get a GenVector&lt;C&gt; from this.
277     * @return a GenVector&lt;C&gt; from this up to truncate parts.
278     */
279    public GenVector<C> asVector() {
280        GenVectorModul<C> vr = new GenVectorModul<C>(ring.coFac, truncate);
281        List<C> v = new ArrayList<>(truncate);
282        for (int i = 0; i < truncate; i++) {
283            C c = coefficient(i);
284            v.add(c);
285        }
286        GenVector<C> vv = new GenVector<C>(vr,v);
287        return vv;
288    }
289
290
291    /**
292     * Leading base coefficient.
293     * @return first coefficient.
294     */
295    public C leadingCoefficient() {
296        return coefficient(0);
297    }
298
299
300    /**
301     * Reductum.
302     * @return this - leading monomial.
303     */
304    public UnivPowerSeries<C> reductum() {
305        final int m = order();
306        //System.out.println("order = " + m);
307
308        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
309
310
311            @Override
312            public C generate(int i) {
313                if (m == i) {
314                    return ring.coFac.getZERO();
315                } else {
316                    return coefficient(i);
317                }
318            }
319        });
320    }
321
322
323    /**
324     * Prepend a new leading coefficient.
325     * @param h new coefficient.
326     * @return new power series.
327     */
328    public UnivPowerSeries<C> prepend(final C h) {
329        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
330
331
332            @Override
333            public C generate(int i) {
334                if (i == 0) {
335                    return h;
336                }
337                return coefficient(i - 1);
338            }
339        });
340    }
341
342
343    /**
344     * Shift coefficients.
345     * @param k shift index.
346     * @return new power series with coefficient(i) = old.coefficient(i-k).
347     */
348    public UnivPowerSeries<C> shift(final int k) {
349        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
350
351
352            @Override
353            public C generate(int i) {
354                if (i - k < 0) {
355                    return ring.coFac.getZERO();
356                }
357                return coefficient(i - k);
358            }
359        });
360    }
361
362
363    /**
364     * Select coefficients.
365     * @param sel selector functor.
366     * @return new power series with selected coefficients.
367     */
368    public UnivPowerSeries<C> select(final Selector<? super C> sel) {
369        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
370
371
372            @Override
373            public C generate(int i) {
374                C c = coefficient(i);
375                if (sel.select(c)) {
376                    return c;
377                }
378                return ring.coFac.getZERO();
379            }
380        });
381    }
382
383
384    /**
385     * Shift select coefficients. Not selected coefficients are removed from the
386     * result series.
387     * @param sel selector functor.
388     * @return new power series with shifted selected coefficients.
389     */
390    public UnivPowerSeries<C> shiftSelect(final Selector<? super C> sel) {
391        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
392
393
394            int pos = 0;
395
396
397            @Override
398            public C generate(int i) {
399                C c;
400                if (i > 0) {
401                    c = get(i - 1); // ensure coeffs are all generated
402                }
403                do {
404                    c = coefficient(pos++);
405                } while (!sel.select(c));
406                return c;
407            }
408        });
409    }
410
411
412    /**
413     * Map a unary function to this power series.
414     * @param f evaluation functor.
415     * @return new power series with coefficients f(this(i)).
416     */
417    public UnivPowerSeries<C> map(final UnaryFunctor<? super C, C> f) {
418        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
419
420
421            @Override
422            public C generate(int i) {
423                return f.eval(coefficient(i));
424            }
425        });
426    }
427
428
429    /**
430     * Map a binary function to this and another power series.
431     * @param f evaluation functor with coefficients f(this(i),other(i)).
432     * @param ps other power series.
433     * @return new power series.
434     */
435    public <C2 extends RingElem<C2>> UnivPowerSeries<C> zip(final BinaryFunctor<? super C, ? super C2, C> f,
436                    final UnivPowerSeries<C2> ps) {
437        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
438
439
440            @Override
441            public C generate(int i) {
442                return f.eval(coefficient(i), ps.coefficient(i));
443            }
444        });
445    }
446
447
448    /**
449     * Sum of two power series.
450     * @param ps other power series.
451     * @return this + ps.
452     */
453    public UnivPowerSeries<C> sum(UnivPowerSeries<C> ps) {
454        return zip(new Sum<C>(), ps);
455    }
456
457
458    /**
459     * Subtraction of two power series.
460     * @param ps other power series.
461     * @return this - ps.
462     */
463    public UnivPowerSeries<C> subtract(UnivPowerSeries<C> ps) {
464        return zip(new Subtract<C>(), ps);
465    }
466
467
468    /**
469     * Multiply by coefficient.
470     * @param c coefficient.
471     * @return this * c.
472     */
473    public UnivPowerSeries<C> multiply(C c) {
474        return map(new Multiply<C>(c));
475    }
476
477
478    /**
479     * Monic.
480     * @return 1/orderCoeff() * this.
481     */
482    public UnivPowerSeries<C> monic() {
483        int i = order();
484        C a = coefficient(i);
485        if (a.isONE()) {
486            return this;
487        }
488        if (a.isZERO()) { // sic
489            return this;
490        }
491        final C b = a.inverse();
492        return map(new UnaryFunctor<C, C>() {
493
494
495            @Override
496            public C eval(C c) {
497                return b.multiply(c);
498            }
499        });
500    }
501
502
503    /**
504     * Negate.
505     * @return - this.
506     */
507    public UnivPowerSeries<C> negate() {
508        return map(new Negate<C>());
509    }
510
511
512    /**
513     * Absolute value.
514     * @return abs(this).
515     */
516    public UnivPowerSeries<C> abs() {
517        if (signum() < 0) {
518            return negate();
519        }
520        return this;
521    }
522
523
524    /**
525     * Evaluate at given point.
526     * @return ps(c).
527     */
528    public C evaluate(C e) {
529        C v = coefficient(0);
530        C p = e;
531        for (int i = 1; i < truncate; i++) {
532            C c = coefficient(i).multiply(p);
533            v = v.sum(c);
534            p = p.multiply(e);
535        }
536        return v;
537    }
538
539
540    /**
541     * Order.
542     * @return index of first non zero coefficient.
543     */
544    public int order() {
545        if (order < 0) { // compute it
546            for (int i = 0; i < truncate; i++) {
547                if (!coefficient(i).isZERO()) {
548                    order = i;
549                    return order;
550                }
551            }
552            order = truncate/* + 1*/;
553        }
554        return order;
555    }
556
557
558    /**
559     * Truncate.
560     * @return truncate index of power series.
561     */
562    public int truncate() {
563        return truncate;
564    }
565
566
567    /**
568     * Set truncate.
569     * @param t new truncate index.
570     * @return old truncate index of power series.
571     */
572    public int setTruncate(int t) {
573        if (t < 0) {
574            throw new IllegalArgumentException("negative truncate not allowed");
575        }
576        int ot = truncate;
577        truncate = t;
578        return ot;
579    }
580
581
582    /**
583     * Signum.
584     * @return sign of first non zero coefficient.
585     */
586    public int signum() {
587        return coefficient(order()).signum();
588    }
589
590
591    /**
592     * Compare to. <b>Note: </b> compare only up to truncate.
593     * @return sign of first non zero coefficient of this-ps.
594     */
595    @Override
596    public int compareTo(UnivPowerSeries<C> ps) {
597        int m = order();
598        int n = ps.order();
599        int pos = (m <= n) ? m : n;
600        int s = 0;
601        do {
602            s = coefficient(pos).compareTo(ps.coefficient(pos));
603            pos++;
604        } while (s == 0 && pos < truncate);
605        return s;
606    }
607
608
609    /**
610     * Is power series zero. <b>Note: </b> compare only up to truncate.
611     * @return If this is 0 then true is returned, else false.
612     * @see edu.jas.structure.RingElem#isZERO()
613     */
614    public boolean isZERO() {
615        return (compareTo(ring.ZERO) == 0);
616    }
617
618
619    /**
620     * Is power series one. <b>Note: </b> compare only up to truncate.
621     * @return If this is 1 then true is returned, else false.
622     * @see edu.jas.structure.RingElem#isONE()
623     */
624    public boolean isONE() {
625        return (compareTo(ring.ONE) == 0);
626    }
627
628
629    /**
630     * Comparison with any other object. <b>Note: </b> compare only up to
631     * truncate.
632     * @see java.lang.Object#equals(java.lang.Object)
633     */
634    @Override
635    @SuppressWarnings("unchecked")
636    public boolean equals(Object B) {
637        UnivPowerSeries<C> a = null;
638        try {
639            a = (UnivPowerSeries<C>) B;
640        } catch (ClassCastException ignored) {
641        }
642        if (a == null) {
643            return false;
644        }
645        return compareTo(a) == 0;
646    }
647
648
649    /**
650     * Hash code for this polynomial. <b>Note: </b> only up to truncate.
651     * @see java.lang.Object#hashCode()
652     */
653    @Override
654    public int hashCode() {
655        int h = 0;
656        //h = ( ring.hashCode() << 23 );
657        //h += val.hashCode();
658        for (int i = 0; i < truncate; i++) {
659            C c = coefficient(i);
660            if (!c.isZERO()) {
661                h += i;
662                h = Math.abs(h << 1);
663            }
664            h += c.hashCode();
665            h = Math.abs(h); // << 3);
666        }
667        return h;
668    }
669
670
671    /**
672     * Is unit.
673     * @return true, if this power series is invertible, else false.
674     */
675    public boolean isUnit() {
676        return leadingCoefficient().isUnit();
677    }
678
679
680    /**
681     * Multiply by another power series.
682     * @return this * ps.
683     */
684    public UnivPowerSeries<C> multiply(final UnivPowerSeries<C> ps) {
685        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
686
687
688            @Override
689            public C generate(int i) {
690                C c = null; //fac.getZERO();
691                for (int k = 0; k <= i; k++) {
692                    C m = coefficient(k).multiply(ps.coefficient(i - k));
693                    if (c == null) {
694                        c = m;
695                    } else {
696                        c = c.sum(m);
697                    }
698                }
699                return c;
700            }
701        });
702    }
703
704
705    /**
706     * Inverse power series.
707     * @return ps with this * ps = 1.
708     */
709    public UnivPowerSeries<C> inverse() {
710        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
711
712
713            @Override
714            public C generate(int i) {
715                C d = leadingCoefficient().inverse(); // may fail
716                if (i == 0) {
717                    return d;
718                }
719                C c = null; //fac.getZERO();
720                for (int k = 0; k < i; k++) {
721                    C m = get(k).multiply(coefficient(i - k));
722                    if (c == null) {
723                        c = m;
724                    } else {
725                        c = c.sum(m);
726                    }
727                }
728                c = c.multiply(d.negate());
729                return c;
730            }
731        });
732    }
733
734
735    /**
736     * Divide by another power series.
737     * @return this / ps.
738     */
739    public UnivPowerSeries<C> divide(UnivPowerSeries<C> ps) {
740        if (ps.isUnit()) {
741            return multiply(ps.inverse());
742        }
743        int m = order();
744        int n = ps.order();
745        if (m < n) {
746            return ring.getZERO();
747        }
748        if (!ps.coefficient(n).isUnit()) {
749            throw new ArithmeticException(
750                            "division by non unit coefficient " + ps.coefficient(n) + ", n = " + n);
751        }
752        // now m >= n
753        UnivPowerSeries<C> st, sps, q, sq;
754        if (m == 0) {
755            st = this;
756        } else {
757            st = this.shift(-m);
758        }
759        if (n == 0) {
760            sps = ps;
761        } else {
762            sps = ps.shift(-n);
763        }
764        q = st.multiply(sps.inverse());
765        if (m == n) {
766            sq = q;
767        } else {
768            sq = q.shift(m - n);
769        }
770        return sq;
771    }
772
773
774    /**
775     * Power series remainder.
776     * @param ps nonzero power series with invertible leading coefficient.
777     * @return remainder with this = quotient * ps + remainder.
778     */
779    public UnivPowerSeries<C> remainder(UnivPowerSeries<C> ps) {
780        int m = order();
781        int n = ps.order();
782        if (m >= n) {
783            return ring.getZERO();
784        }
785        return this;
786    }
787
788
789    /**
790     * Quotient and remainder by division of this by S.
791     * @param S a UnivPowerSeries
792     * @return [this/S, this - (this/S)*S].
793     */
794    @SuppressWarnings("unchecked")
795    public UnivPowerSeries<C>[] quotientRemainder(UnivPowerSeries<C> S) {
796        return new UnivPowerSeries[] { divide(S), remainder(S) };
797    }
798
799
800    /**
801     * Differentiate.
802     * @return differentiate(this).
803     */
804    public UnivPowerSeries<C> differentiate() {
805        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
806
807
808            @Override
809            public C generate(int i) {
810                C v = coefficient(i + 1);
811                v = v.multiply(ring.coFac.fromInteger(i + 1));
812                return v;
813            }
814        });
815    }
816
817
818    /**
819     * Integrate with given constant.
820     * @param c integration constant.
821     * @return integrate(this).
822     */
823    public UnivPowerSeries<C> integrate(final C c) {
824        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
825
826
827            @Override
828            public C generate(int i) {
829                if (i == 0) {
830                    return c;
831                }
832                C v = coefficient(i - 1);
833                v = v.divide(ring.coFac.fromInteger(i));
834                return v;
835            }
836        });
837    }
838
839
840    /**
841     * Power series greatest common divisor.
842     * @param ps power series.
843     * @return gcd(this,ps).
844     */
845    public UnivPowerSeries<C> gcd(UnivPowerSeries<C> ps) {
846        if (ps.isZERO()) {
847            return this;
848        }
849        if (this.isZERO()) {
850            return ps;
851        }
852        int m = order();
853        int n = ps.order();
854        int ll = (m < n) ? m : n;
855        return ring.getONE().shift(ll);
856    }
857
858
859    /**
860     * Power series extended greatest common divisor. <b>Note:</b> not
861     * implemented.
862     * @param S power series.
863     * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
864     */
865    //SuppressWarnings("unchecked")
866    public UnivPowerSeries<C>[] egcd(UnivPowerSeries<C> S) {
867        throw new UnsupportedOperationException("egcd for power series not implemented");
868    }
869
870}
871
872
873/* arithmetic method functors */
874
875
876/**
877 * Internal summation functor.
878 */
879class Sum<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
880
881
882    public C eval(C c1, C c2) {
883        return c1.sum(c2);
884    }
885}
886
887
888/**
889 * Internal subtraction functor.
890 */
891class Subtract<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
892
893
894    public C eval(C c1, C c2) {
895        return c1.subtract(c2);
896    }
897}
898
899
900/**
901 * Internal scalar multiplication functor.
902 */
903class Multiply<C extends RingElem<C>> implements UnaryFunctor<C, C> {
904
905
906    C x;
907
908
909    public Multiply(C x) {
910        this.x = x;
911    }
912
913
914    public C eval(C c) {
915        return c.multiply(x);
916    }
917}
918
919
920/**
921 * Internal negation functor.
922 */
923class Negate<C extends RingElem<C>> implements UnaryFunctor<C, C> {
924
925
926    public C eval(C c) {
927        return c.negate();
928    }
929}
930
931
932/* only for sequential access:
933class Abs<C extends RingElem<C>> implements UnaryFunctor<C,C> {
934        int sign = 0;
935        public C eval(C c) {
936            int s = c.signum();
937            if ( s == 0 ) {
938               return c;
939            }
940            if ( sign > 0 ) {
941               return c;
942            } else if ( sign < 0 ) {
943               return c.negate();
944            }
945            // first non zero coefficient:
946            sign = s;
947            if ( s > 0 ) {
948               return c;
949            }
950            return c.negate();
951        }
952}
953*/