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