001/*
002 * $Id: UnivPowerSeries.java 4125 2012-08-19 19:05:22Z 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    //JAVA6only: @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    //JAVA6only: @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     * Negate.
453     * @return - this.
454     */
455    public UnivPowerSeries<C> negate() {
456        return map(new Negate<C>());
457    }
458
459
460    /**
461     * Absolute value.
462     * @return abs(this).
463     */
464    public UnivPowerSeries<C> abs() {
465        if (signum() < 0) {
466            return negate();
467        }
468        return this;
469    }
470
471
472    /**
473     * Evaluate at given point.
474     * @return ps(c).
475     */
476    public C evaluate(C e) {
477        C v = coefficient(0);
478        C p = e;
479        for (int i = 1; i < truncate; i++) {
480            C c = coefficient(i).multiply(p);
481            v = v.sum(c);
482            p = p.multiply(e);
483        }
484        return v;
485    }
486
487
488    /**
489     * Order.
490     * @return index of first non zero coefficient.
491     */
492    public int order() {
493        if (order < 0) { // compute it
494            for (int i = 0; i <= truncate; i++) {
495                if (!coefficient(i).isZERO()) {
496                    order = i;
497                    return order;
498                }
499            }
500            order = truncate + 1;
501        }
502        return order;
503    }
504
505
506    /**
507     * Truncate.
508     * @return truncate index of power series.
509     */
510    public int truncate() {
511        return truncate;
512    }
513
514
515    /**
516     * Set truncate.
517     * @param t new truncate index.
518     * @return old truncate index of power series.
519     */
520    public int setTruncate(int t) {
521        if (t < 0) {
522            throw new IllegalArgumentException("negative truncate not allowed");
523        }
524        int ot = truncate;
525        truncate = t;
526        return ot;
527    }
528
529
530    /**
531     * Signum.
532     * @return sign of first non zero coefficient.
533     */
534    public int signum() {
535        return coefficient(order()).signum();
536    }
537
538
539    /**
540     * Compare to. <b>Note: </b> compare only up to truncate.
541     * @return sign of first non zero coefficient of this-ps.
542     */
543    //JAVA6only: @Override
544    public int compareTo(UnivPowerSeries<C> ps) {
545        int m = order();
546        int n = ps.order();
547        int pos = (m <= n) ? m : n;
548        int s = 0;
549        do {
550            s = coefficient(pos).compareTo(ps.coefficient(pos));
551            pos++;
552        } while (s == 0 && pos <= truncate);
553        return s;
554    }
555
556
557    /**
558     * Is power series zero. <b>Note: </b> compare only up to truncate.
559     * @return If this is 0 then true is returned, else false.
560     * @see edu.jas.structure.RingElem#isZERO()
561     */
562    public boolean isZERO() {
563        return (compareTo(ring.ZERO) == 0);
564    }
565
566
567    /**
568     * Is power series one. <b>Note: </b> compare only up to truncate.
569     * @return If this is 1 then true is returned, else false.
570     * @see edu.jas.structure.RingElem#isONE()
571     */
572    public boolean isONE() {
573        return (compareTo(ring.ONE) == 0);
574    }
575
576
577    /**
578     * Comparison with any other object. <b>Note: </b> compare only up to
579     * truncate.
580     * @see java.lang.Object#equals(java.lang.Object)
581     */
582    @Override
583    @SuppressWarnings("unchecked")
584    public boolean equals(Object B) {
585        UnivPowerSeries<C> a = null;
586        try {
587            a = (UnivPowerSeries<C>) B;
588        } catch (ClassCastException ignored) {
589        }
590        if (a == null) {
591            return false;
592        }
593        return compareTo(a) == 0;
594    }
595
596
597    /**
598     * Hash code for this polynomial. <b>Note: </b> only up to truncate.
599     * @see java.lang.Object#hashCode()
600     */
601    @Override
602    public int hashCode() {
603        int h = 0;
604        //h = ( ring.hashCode() << 23 );
605        //h += val.hashCode();
606        for (int i = 0; i <= truncate; i++) {
607            h += coefficient(i).hashCode();
608            h = (h << 23);
609        }
610        return h;
611    }
612
613
614    /**
615     * Is unit.
616     * @return true, if this power series is invertible, else false.
617     */
618    public boolean isUnit() {
619        return leadingCoefficient().isUnit();
620    }
621
622
623    /**
624     * Multiply by another power series.
625     * @return this * ps.
626     */
627    public UnivPowerSeries<C> multiply(final UnivPowerSeries<C> ps) {
628        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
629
630
631            @Override
632            public C generate(int i) {
633                C c = null; //fac.getZERO();
634                for (int k = 0; k <= i; k++) {
635                    C m = coefficient(k).multiply(ps.coefficient(i - k));
636                    if (c == null) {
637                        c = m;
638                    } else {
639                        c = c.sum(m);
640                    }
641                }
642                return c;
643            }
644        });
645    }
646
647
648    /**
649     * Inverse power series.
650     * @return ps with this * ps = 1.
651     */
652    public UnivPowerSeries<C> inverse() {
653        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
654
655
656            @Override
657            public C generate(int i) {
658                C d = leadingCoefficient().inverse(); // may fail
659                if (i == 0) {
660                    return d;
661                }
662                C c = null; //fac.getZERO();
663                for (int k = 0; k < i; k++) {
664                    C m = get(k).multiply(coefficient(i - k));
665                    if (c == null) {
666                        c = m;
667                    } else {
668                        c = c.sum(m);
669                    }
670                }
671                c = c.multiply(d.negate());
672                return c;
673            }
674        });
675    }
676
677
678    /**
679     * Divide by another power series.
680     * @return this / ps.
681     */
682    public UnivPowerSeries<C> divide(UnivPowerSeries<C> ps) {
683        if (ps.isUnit()) {
684            return multiply(ps.inverse());
685        }
686        int m = order();
687        int n = ps.order();
688        if (m < n) {
689            return ring.getZERO();
690        }
691        if (!ps.coefficient(n).isUnit()) {
692            throw new ArithmeticException("division by non unit coefficient " + ps.coefficient(n) + ", n = "
693                            + n);
694        }
695        // now m >= n
696        UnivPowerSeries<C> st, sps, q, sq;
697        if (m == 0) {
698            st = this;
699        } else {
700            st = this.shift(-m);
701        }
702        if (n == 0) {
703            sps = ps;
704        } else {
705            sps = ps.shift(-n);
706        }
707        q = st.multiply(sps.inverse());
708        if (m == n) {
709            sq = q;
710        } else {
711            sq = q.shift(m - n);
712        }
713        return sq;
714    }
715
716
717    /**
718     * Power series remainder.
719     * @param ps nonzero power series with invertible leading coefficient.
720     * @return remainder with this = quotient * ps + remainder.
721     */
722    public UnivPowerSeries<C> remainder(UnivPowerSeries<C> ps) {
723        int m = order();
724        int n = ps.order();
725        if (m >= n) {
726            return ring.getZERO();
727        }
728        return this;
729    }
730
731
732    /**
733     * Differentiate.
734     * @return differentiate(this).
735     */
736    public UnivPowerSeries<C> differentiate() {
737        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
738
739
740            @Override
741            public C generate(int i) {
742                C v = coefficient(i + 1);
743                v = v.multiply(ring.coFac.fromInteger(i + 1));
744                return v;
745            }
746        });
747    }
748
749
750    /**
751     * Integrate with given constant.
752     * @param c integration constant.
753     * @return integrate(this).
754     */
755    public UnivPowerSeries<C> integrate(final C c) {
756        return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
757
758
759            @Override
760            public C generate(int i) {
761                if (i == 0) {
762                    return c;
763                }
764                C v = coefficient(i - 1);
765                v = v.divide(ring.coFac.fromInteger(i));
766                return v;
767            }
768        });
769    }
770
771
772    /**
773     * Power series greatest common divisor.
774     * @param ps power series.
775     * @return gcd(this,ps).
776     */
777    public UnivPowerSeries<C> gcd(UnivPowerSeries<C> ps) {
778        if (ps.isZERO()) {
779            return this;
780        }
781        if (this.isZERO()) {
782            return ps;
783        }
784        int m = order();
785        int n = ps.order();
786        int ll = (m < n) ? m : n;
787        return ring.getONE().shift(ll);
788    }
789
790
791    /**
792     * Power series extended greatest common divisor. <b>Note:</b> not
793     * implemented.
794     * @param S power series.
795     * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
796     */
797    //SuppressWarnings("unchecked")
798    public UnivPowerSeries<C>[] egcd(UnivPowerSeries<C> S) {
799        throw new UnsupportedOperationException("egcd for power series not implemented");
800    }
801
802}
803
804
805/* arithmetic method functors */
806
807
808/**
809 * Internal summation functor.
810 */
811class Sum<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
812
813
814    public C eval(C c1, C c2) {
815        return c1.sum(c2);
816    }
817}
818
819
820/**
821 * Internal subtraction functor.
822 */
823class Subtract<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
824
825
826    public C eval(C c1, C c2) {
827        return c1.subtract(c2);
828    }
829}
830
831
832/**
833 * Internal scalar multiplication functor.
834 */
835class Multiply<C extends RingElem<C>> implements UnaryFunctor<C, C> {
836
837
838    C x;
839
840
841    public Multiply(C x) {
842        this.x = x;
843    }
844
845
846    public C eval(C c) {
847        return c.multiply(x);
848    }
849}
850
851
852/**
853 * Internal negation functor.
854 */
855class Negate<C extends RingElem<C>> implements UnaryFunctor<C, C> {
856
857
858    public C eval(C c) {
859        return c.negate();
860    }
861}
862
863
864/* only for sequential access:
865class Abs<C extends RingElem<C>> implements UnaryFunctor<C,C> {
866        int sign = 0;
867        public C eval(C c) {
868            int s = c.signum();
869            if ( s == 0 ) {
870               return c;
871            }
872            if ( sign > 0 ) {
873               return c;
874            } else if ( sign < 0 ) {
875               return c.negate();
876            }
877            // first non zero coefficient:
878            sign = s;
879            if ( s > 0 ) {
880               return c;
881            }
882            return c.negate();
883        }
884}
885*/