001 /*
002 * $Id: UnivPowerSeries.java 3571 2011-03-18 22:02:51Z kredel $
003 */
004
005 package edu.jas.ps;
006
007
008 import edu.jas.poly.AlgebraicNumber;
009 import edu.jas.poly.ExpVector;
010 import edu.jas.poly.GenPolynomial;
011 import edu.jas.structure.BinaryFunctor;
012 import edu.jas.structure.RingElem;
013 import edu.jas.structure.Selector;
014 import 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
029 public 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> clone() {
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<C> from this.
258 * @return a GenPolynomial<C> 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 } else {
310 return coefficient(i - 1);
311 }
312 }
313 });
314 }
315
316
317 /**
318 * Shift coefficients.
319 * @param k shift index.
320 * @return new power series with coefficient(i) = old.coefficient(i-k).
321 */
322 public UnivPowerSeries<C> shift(final int k) {
323 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
324
325
326 @Override
327 public C generate(int i) {
328 if (i - k < 0) {
329 return ring.coFac.getZERO();
330 }
331 return coefficient(i - k);
332 }
333 });
334 }
335
336
337 /**
338 * Select coefficients.
339 * @param sel selector functor.
340 * @return new power series with selected coefficients.
341 */
342 public UnivPowerSeries<C> select(final Selector<? super C> sel) {
343 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
344
345
346 @Override
347 public C generate(int i) {
348 C c = coefficient(i);
349 if (sel.select(c)) {
350 return c;
351 } else {
352 return ring.coFac.getZERO();
353 }
354 }
355 });
356 }
357
358
359 /**
360 * Shift select coefficients. Not selected coefficients are removed from the
361 * result series.
362 * @param sel selector functor.
363 * @return new power series with shifted selected coefficients.
364 */
365 public UnivPowerSeries<C> shiftSelect(final Selector<? super C> sel) {
366 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
367
368
369 int pos = 0;
370
371
372 @Override
373 public C generate(int i) {
374 C c;
375 if (i > 0) {
376 c = get(i - 1);
377 }
378 do {
379 c = coefficient(pos++);
380 } while (!sel.select(c));
381 return c;
382 }
383 });
384 }
385
386
387 /**
388 * Map a unary function to this power series.
389 * @param f evaluation functor.
390 * @return new power series with coefficients f(this(i)).
391 */
392 public UnivPowerSeries<C> map(final UnaryFunctor<? super C, C> f) {
393 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
394
395
396 @Override
397 public C generate(int i) {
398 return f.eval(coefficient(i));
399 }
400 });
401 }
402
403
404 /**
405 * Map a binary function to this and another power series.
406 * @param f evaluation functor with coefficients f(this(i),other(i)).
407 * @param ps other power series.
408 * @return new power series.
409 */
410 public <C2 extends RingElem<C2>> UnivPowerSeries<C> zip(final BinaryFunctor<? super C, ? super C2, C> f,
411 final UnivPowerSeries<C2> ps) {
412 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
413
414
415 @Override
416 public C generate(int i) {
417 return f.eval(coefficient(i), ps.coefficient(i));
418 }
419 });
420 }
421
422
423 /**
424 * Sum of two power series.
425 * @param ps other power series.
426 * @return this + ps.
427 */
428 public UnivPowerSeries<C> sum(UnivPowerSeries<C> ps) {
429 return zip(new Sum<C>(), ps);
430 }
431
432
433 /**
434 * Subtraction of two power series.
435 * @param ps other power series.
436 * @return this - ps.
437 */
438 public UnivPowerSeries<C> subtract(UnivPowerSeries<C> ps) {
439 return zip(new Subtract<C>(), ps);
440 }
441
442
443 /**
444 * Multiply by coefficient.
445 * @param c coefficient.
446 * @return this * c.
447 */
448 public UnivPowerSeries<C> multiply(C c) {
449 return map(new Multiply<C>(c));
450 }
451
452
453 /**
454 * Negate.
455 * @return - this.
456 */
457 public UnivPowerSeries<C> negate() {
458 return map(new Negate<C>());
459 }
460
461
462 /**
463 * Absolute value.
464 * @return abs(this).
465 */
466 public UnivPowerSeries<C> abs() {
467 if (signum() < 0) {
468 return negate();
469 } else {
470 return this;
471 }
472 }
473
474
475 /**
476 * Evaluate at given point.
477 * @return ps(c).
478 */
479 public C evaluate(C e) {
480 C v = coefficient(0);
481 C p = e;
482 for (int i = 1; i < truncate; i++) {
483 C c = coefficient(i).multiply(p);
484 v = v.sum(c);
485 p = p.multiply(e);
486 }
487 return v;
488 }
489
490
491 /**
492 * Order.
493 * @return index of first non zero coefficient.
494 */
495 public int order() {
496 if (order < 0) { // compute it
497 for (int i = 0; i <= truncate; i++) {
498 if (!coefficient(i).isZERO()) {
499 order = i;
500 return order;
501 }
502 }
503 order = truncate + 1;
504 }
505 return order;
506 }
507
508
509 /**
510 * Truncate.
511 * @return truncate index of power series.
512 */
513 public int truncate() {
514 return truncate;
515 }
516
517
518 /**
519 * Set truncate.
520 * @param t new truncate index.
521 * @return old truncate index of power series.
522 */
523 public int setTruncate(int t) {
524 if (t < 0) {
525 throw new IllegalArgumentException("negative truncate not allowed");
526 }
527 int ot = truncate;
528 truncate = t;
529 return ot;
530 }
531
532
533 /**
534 * Signum.
535 * @return sign of first non zero coefficient.
536 */
537 public int signum() {
538 return coefficient(order()).signum();
539 }
540
541
542 /**
543 * Compare to. <b>Note: </b> compare only up to truncate.
544 * @return sign of first non zero coefficient of this-ps.
545 */
546 //JAVA6only: @Override
547 public int compareTo(UnivPowerSeries<C> ps) {
548 int m = order();
549 int n = ps.order();
550 int pos = (m <= n) ? m : n;
551 int s = 0;
552 do {
553 s = coefficient(pos).compareTo(ps.coefficient(pos));
554 pos++;
555 } while (s == 0 && pos <= truncate);
556 return s;
557 }
558
559
560 /**
561 * Is power series zero. <b>Note: </b> compare only up to truncate.
562 * @return If this is 0 then true is returned, else false.
563 * @see edu.jas.structure.RingElem#isZERO()
564 */
565 public boolean isZERO() {
566 return (compareTo(ring.ZERO) == 0);
567 }
568
569
570 /**
571 * Is power series one. <b>Note: </b> compare only up to truncate.
572 * @return If this is 1 then true is returned, else false.
573 * @see edu.jas.structure.RingElem#isONE()
574 */
575 public boolean isONE() {
576 return (compareTo(ring.ONE) == 0);
577 }
578
579
580 /**
581 * Comparison with any other object. <b>Note: </b> compare only up to
582 * truncate.
583 * @see java.lang.Object#equals(java.lang.Object)
584 */
585 @Override
586 @SuppressWarnings("unchecked")
587 public boolean equals(Object B) {
588 if (!(B instanceof UnivPowerSeries)) {
589 return false;
590 }
591 UnivPowerSeries<C> a = null;
592 try {
593 a = (UnivPowerSeries<C>) B;
594 } catch (ClassCastException ignored) {
595 }
596 if (a == null) {
597 return false;
598 }
599 return compareTo(a) == 0;
600 }
601
602
603 /**
604 * Hash code for this polynomial. <b>Note: </b> only up to truncate.
605 * @see java.lang.Object#hashCode()
606 */
607 @Override
608 public int hashCode() {
609 int h = 0;
610 //h = ( ring.hashCode() << 23 );
611 //h += val.hashCode();
612 for (int i = 0; i <= truncate; i++) {
613 h += coefficient(i).hashCode();
614 h = (h << 23);
615 };
616 return h;
617 }
618
619
620 /**
621 * Is unit.
622 * @return true, if this power series is invertible, else false.
623 */
624 public boolean isUnit() {
625 return leadingCoefficient().isUnit();
626 }
627
628
629 /**
630 * Multiply by another power series.
631 * @return this * ps.
632 */
633 public UnivPowerSeries<C> multiply(final UnivPowerSeries<C> ps) {
634 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
635
636
637 @Override
638 public C generate(int i) {
639 C c = null; //fac.getZERO();
640 for (int k = 0; k <= i; k++) {
641 C m = coefficient(k).multiply(ps.coefficient(i - k));
642 if (c == null) {
643 c = m;
644 } else {
645 c = c.sum(m);
646 }
647 }
648 return c;
649 }
650 });
651 }
652
653
654 /**
655 * Inverse power series.
656 * @return ps with this * ps = 1.
657 */
658 public UnivPowerSeries<C> inverse() {
659 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
660
661
662 @Override
663 public C generate(int i) {
664 C d = leadingCoefficient().inverse(); // may fail
665 if (i == 0) {
666 return d;
667 }
668 C c = null; //fac.getZERO();
669 for (int k = 0; k < i; k++) {
670 C m = get(k).multiply(coefficient(i - k));
671 if (c == null) {
672 c = m;
673 } else {
674 c = c.sum(m);
675 }
676 }
677 c = c.multiply(d.negate());
678 return c;
679 }
680 });
681 }
682
683
684 /**
685 * Divide by another power series.
686 * @return this / ps.
687 */
688 public UnivPowerSeries<C> divide(UnivPowerSeries<C> ps) {
689 if (ps.isUnit()) {
690 return multiply(ps.inverse());
691 }
692 int m = order();
693 int n = ps.order();
694 if (m < n) {
695 return ring.getZERO();
696 }
697 if (!ps.coefficient(n).isUnit()) {
698 throw new ArithmeticException("division by non unit coefficient " + ps.coefficient(n) + ", n = "
699 + n);
700 }
701 // now m >= n
702 UnivPowerSeries<C> st, sps, q, sq;
703 if (m == 0) {
704 st = this;
705 } else {
706 st = this.shift(-m);
707 }
708 if (n == 0) {
709 sps = ps;
710 } else {
711 sps = ps.shift(-n);
712 }
713 q = st.multiply(sps.inverse());
714 if (m == n) {
715 sq = q;
716 } else {
717 sq = q.shift(m - n);
718 }
719 return sq;
720 }
721
722
723 /**
724 * Power series remainder.
725 * @param ps nonzero power series with invertible leading coefficient.
726 * @return remainder with this = quotient * ps + remainder.
727 */
728 public UnivPowerSeries<C> remainder(UnivPowerSeries<C> ps) {
729 int m = order();
730 int n = ps.order();
731 if (m >= n) {
732 return ring.getZERO();
733 }
734 return this;
735 }
736
737
738 /**
739 * Differentiate.
740 * @return differentiate(this).
741 */
742 public UnivPowerSeries<C> differentiate() {
743 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
744
745
746 @Override
747 public C generate(int i) {
748 C v = coefficient(i + 1);
749 v = v.multiply(ring.coFac.fromInteger(i + 1));
750 return v;
751 }
752 });
753 }
754
755
756 /**
757 * Integrate with given constant.
758 * @param c integration constant.
759 * @return integrate(this).
760 */
761 public UnivPowerSeries<C> integrate(final C c) {
762 return new UnivPowerSeries<C>(ring, new Coefficients<C>() {
763
764
765 @Override
766 public C generate(int i) {
767 if (i == 0) {
768 return c;
769 }
770 C v = coefficient(i - 1);
771 v = v.divide(ring.coFac.fromInteger(i));
772 return v;
773 }
774 });
775 }
776
777
778 /**
779 * Power series greatest common divisor.
780 * @param ps power series.
781 * @return gcd(this,ps).
782 */
783 public UnivPowerSeries<C> gcd(UnivPowerSeries<C> ps) {
784 if (ps.isZERO()) {
785 return this;
786 }
787 if (this.isZERO()) {
788 return ps;
789 }
790 int m = order();
791 int n = ps.order();
792 int ll = (m < n) ? m : n;
793 return ring.getONE().shift(ll);
794 }
795
796
797 /**
798 * Power series extended greatest common divisor. <b>Note:</b> not
799 * implemented.
800 * @param S power series.
801 * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
802 */
803 //SuppressWarnings("unchecked")
804 public UnivPowerSeries<C>[] egcd(UnivPowerSeries<C> S) {
805 throw new UnsupportedOperationException("egcd for power series not implemented");
806 }
807
808 }
809
810
811 /* arithmetic method functors */
812
813
814 /**
815 * Internal summation functor.
816 */
817 class Sum<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
818
819
820 public C eval(C c1, C c2) {
821 return c1.sum(c2);
822 }
823 }
824
825
826 /**
827 * Internal subtraction functor.
828 */
829 class Subtract<C extends RingElem<C>> implements BinaryFunctor<C, C, C> {
830
831
832 public C eval(C c1, C c2) {
833 return c1.subtract(c2);
834 }
835 }
836
837
838 /**
839 * Internal scalar multiplication functor.
840 */
841 class Multiply<C extends RingElem<C>> implements UnaryFunctor<C, C> {
842
843
844 C x;
845
846
847 public Multiply(C x) {
848 this.x = x;
849 }
850
851
852 public C eval(C c) {
853 return c.multiply(x);
854 }
855 }
856
857
858 /**
859 * Internal negation functor.
860 */
861 class Negate<C extends RingElem<C>> implements UnaryFunctor<C, C> {
862
863
864 public C eval(C c) {
865 return c.negate();
866 }
867 }
868
869
870 /* only for sequential access:
871 class Abs<C extends RingElem<C>> implements UnaryFunctor<C,C> {
872 int sign = 0;
873 public C eval(C c) {
874 int s = c.signum();
875 if ( s == 0 ) {
876 return c;
877 }
878 if ( sign > 0 ) {
879 return c;
880 } else if ( sign < 0 ) {
881 return c.negate();
882 }
883 // first non zero coefficient:
884 sign = s;
885 if ( s > 0 ) {
886 return c;
887 }
888 return c.negate();
889 }
890 }
891 */