001/* 002 * $Id: UnivPowerSeries.java 5934 2018-09-30 11:23:44Z 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( 084 "null not allowed: ring = " + ring + ", lazyCoeffs = " + 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<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 } 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( 718 "division by non unit coefficient " + ps.coefficient(n) + ", n = " + 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 @SuppressWarnings("unchecked") 763 public UnivPowerSeries<C>[] quotientRemainder(UnivPowerSeries<C> S) { 764 return new UnivPowerSeries[] { divide(S), remainder(S) }; 765 } 766 767 768 /** 769 * Differentiate. 770 * @return differentiate(this). 771 */ 772 public UnivPowerSeries<C> differentiate() { 773 return new UnivPowerSeries<C>(ring, new Coefficients<C>() { 774 775 776 @Override 777 public C generate(int i) { 778 C v = coefficient(i + 1); 779 v = v.multiply(ring.coFac.fromInteger(i + 1)); 780 return v; 781 } 782 }); 783 } 784 785 786 /** 787 * Integrate with given constant. 788 * @param c integration constant. 789 * @return integrate(this). 790 */ 791 public UnivPowerSeries<C> integrate(final C c) { 792 return new UnivPowerSeries<C>(ring, new Coefficients<C>() { 793 794 795 @Override 796 public C generate(int i) { 797 if (i == 0) { 798 return c; 799 } 800 C v = coefficient(i - 1); 801 v = v.divide(ring.coFac.fromInteger(i)); 802 return v; 803 } 804 }); 805 } 806 807 808 /** 809 * Power series greatest common divisor. 810 * @param ps power series. 811 * @return gcd(this,ps). 812 */ 813 public UnivPowerSeries<C> gcd(UnivPowerSeries<C> ps) { 814 if (ps.isZERO()) { 815 return this; 816 } 817 if (this.isZERO()) { 818 return ps; 819 } 820 int m = order(); 821 int n = ps.order(); 822 int ll = (m < n) ? m : n; 823 return ring.getONE().shift(ll); 824 } 825 826 827 /** 828 * Power series extended greatest common divisor. <b>Note:</b> not 829 * implemented. 830 * @param S power series. 831 * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S). 832 */ 833 //SuppressWarnings("unchecked") 834 public UnivPowerSeries<C>[] egcd(UnivPowerSeries<C> S) { 835 throw new UnsupportedOperationException("egcd for power series not implemented"); 836 } 837 838} 839 840 841/* arithmetic method functors */ 842 843 844/** 845 * Internal summation functor. 846 */ 847class Sum<C extends RingElem<C>> implements BinaryFunctor<C, C, C> { 848 849 850 public C eval(C c1, C c2) { 851 return c1.sum(c2); 852 } 853} 854 855 856/** 857 * Internal subtraction functor. 858 */ 859class Subtract<C extends RingElem<C>> implements BinaryFunctor<C, C, C> { 860 861 862 public C eval(C c1, C c2) { 863 return c1.subtract(c2); 864 } 865} 866 867 868/** 869 * Internal scalar multiplication functor. 870 */ 871class Multiply<C extends RingElem<C>> implements UnaryFunctor<C, C> { 872 873 874 C x; 875 876 877 public Multiply(C x) { 878 this.x = x; 879 } 880 881 882 public C eval(C c) { 883 return c.multiply(x); 884 } 885} 886 887 888/** 889 * Internal negation functor. 890 */ 891class Negate<C extends RingElem<C>> implements UnaryFunctor<C, C> { 892 893 894 public C eval(C c) { 895 return c.negate(); 896 } 897} 898 899 900/* only for sequential access: 901class Abs<C extends RingElem<C>> implements UnaryFunctor<C,C> { 902 int sign = 0; 903 public C eval(C c) { 904 int s = c.signum(); 905 if ( s == 0 ) { 906 return c; 907 } 908 if ( sign > 0 ) { 909 return c; 910 } else if ( sign < 0 ) { 911 return c.negate(); 912 } 913 // first non zero coefficient: 914 sign = s; 915 if ( s > 0 ) { 916 return c; 917 } 918 return c.negate(); 919 } 920} 921*/