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<C> from this. 262 * @return a GenPolynomial<C> 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<C> from this. 277 * @return a GenVector<C> 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*/