| 388 | def partial_fraction_decomposition(self, kind=None, factor=True): |
| 389 | r""" |
| 390 | Return a partial fraction decomposition of ``self`` of kind |
| 391 | ``kind``. |
| 392 | |
| 393 | INPUT: |
| 394 | |
| 395 | - ``kind`` - (Optional; default=None) The string |
| 396 | 'nullstellensatz' or 'algebraic_dependence', indicating what kind |
| 397 | of multivariate decomposition to return in case ``self`` is |
| 398 | multivariate. |
| 399 | The default multivariate decomposition is a Leinartas |
| 400 | decomposition. |
| 401 | - ``factor`` - (Optional; default=True) If True, then the |
| 402 | output will be a REFDSum instance (which features factored |
| 403 | denominators). |
| 404 | Note that factoring the denominator of ``self`` happens anyway as |
| 405 | part of the algorithm. |
| 406 | |
| 407 | OUTPUT: |
| 408 | |
| 409 | A sum of fraction field elements or, if factor=True, |
| 410 | a REFDSum instance that is a partial fraction decomposition of |
| 411 | ``self``. |
| 412 | |
| 413 | EXAMPLES:: |
| 414 | |
| 415 | No variables:: |
| 416 | |
| 417 | sage: f = 64/45 |
| 418 | sage: print f.partial_fraction_decomposition() |
| 419 | [(1, []), (2, [(3, 2)]), (1, [(5, 1)])] |
| 420 | sage: print f.partial_fraction_decomposition(factor=False) |
| 421 | [1, 2/9, 1/5] |
| 422 | sage: print sum(f.partial_fraction_decomposition(factor=False)) == f |
| 423 | True |
| 424 | |
| 425 | One variable:: |
| 426 | |
| 427 | sage: R.<x> = PolynomialRing(QQ) |
| 428 | sage: f = 1 + 2/x + 2*x/(x^2 - 4*x + 8) |
| 429 | sage: print f |
| 430 | (x^3 + 16)/(x^3 - 4*x^2 + 8*x) |
| 431 | sage: d1 = f.partial_fraction_decomposition() |
| 432 | sage: print d1 |
| 433 | [(1, []), (2, [(x, 1)]), (2*x, [(x^2 - 4*x + 8, 1)])] |
| 434 | sage: d1.sum().quotient() == f |
| 435 | True |
| 436 | sage: print f.partial_fraction_decomposition(factor=False) |
| 437 | [1, 2/x, 2*x/(x^2 - 4*x + 8)] |
| 438 | |
| 439 | Two variables:: |
| 440 | |
| 441 | sage: R.<x,y>= PolynomialRing(QQ) |
| 442 | sage: f = 1/x + 1/y + 1/(x + y) |
| 443 | sage: print f |
| 444 | (x^2 + 3*x*y + y^2)/(x^2*y + x*y^2) |
| 445 | sage: print f.partial_fraction_decomposition(kind='nullstellensatz') |
| 446 | [(0, []), (x^2 + 3*x*y + y^2, [(y, 1), (x, 1), (x + y, 1)])] |
| 447 | sage: print f.partial_fraction_decomposition(kind='algebraic_dependence') |
| 448 | [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])] |
| 449 | sage: print f.partial_fraction_decomposition() |
| 450 | [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])] |
| 451 | |
| 452 | Three variables over a finite field:: |
| 453 | |
| 454 | sage: R.<x,y,z>= PolynomialRing(GF(2)) |
| 455 | sage: f = x + 1/x + 1/y + 1/z + 1/(x*y + z) |
| 456 | sage: print f |
| 457 | (x^3*y^2*z + x^2*y*z^2 + x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2)/(x^2*y^2*z + x*y*z^2) |
| 458 | sage: d1 = f.partial_fraction_decomposition() |
| 459 | sage: print d1 |
| 460 | [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])] |
| 461 | sage: d2 = f.partial_fraction_decomposition() |
| 462 | sage: print d2 |
| 463 | [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])] |
| 464 | sage: print d2.sum().quotient() == f |
| 465 | True |
| 466 | |
| 467 | NOTE:: |
| 468 | |
| 469 | See the docstrings of the REFD class methods |
| 470 | leinartas_decomposition(), |
| 471 | nullstellensatz_decomposition(), and |
| 472 | algebraic_dependence_decomposition() |
| 473 | for definitions of the different kinds of partial fraction |
| 474 | decomposition options available. |
| 475 | |
| 476 | AUTHORS: |
| 477 | |
| 478 | - Robert Bradshaw (2007-05-31) |
| 479 | - Alexander Raichev (2012-06-25) |
| 480 | """ |
| 481 | f = REFD(quotient=self) |
| 482 | if f.ring().ngens() <= 1: |
| 483 | decomp = f.univariate_decomposition() |
| 484 | else: |
| 485 | if kind == 'nullstellensatz': |
| 486 | decomp = f.nullstellensatz_decomposition() |
| 487 | elif kind == 'algebraic_dependence': |
| 488 | decomp = f.algebraic_dependence_decomposition() |
| 489 | else: |
| 490 | decomp = f.leinartas_decomposition() |
| 491 | if factor: |
| 492 | return decomp |
| 493 | else: |
| 494 | return [r.quotient() for r in decomp] |
| 495 | |
| 496 | |
| 497 | @total_ordering |
| 498 | class REFD(object): |
| 499 | r""" |
| 500 | Represents a rational expression with factored denominator (REFD) |
| 501 | `p/(q_1^{e_1} \cdots q_m^{e_m})` by storing the parts `p` and |
| 502 | `[(q_1,e_1), \ldots, (q_m,e_m)]`. |
| 503 | Here `p, q_1, \ldots, q_m` are elements of a 0- or 1-variate factorial |
| 504 | polynomial ring `R`, `q_1, \ldots, q_m` are distinct irreducible elements |
| 505 | of `R`, and `e_1, \ldots, e_m` are positive integers. |
| 506 | An element `r` of `R` is represented as `[r, (,)]`. |
| 507 | |
| 508 | AUTHORS: |
| 509 | |
| 510 | - Alexander Raichev (2012-06-25) |
| 511 | """ |
| 512 | |
| 513 | def __init__(self, numerator=None, denominator_list=None, quotient=None, |
| 514 | reduce_=True): |
| 515 | r""" |
| 516 | Create a REFD instance. |
| 517 | |
| 518 | INPUT: |
| 519 | |
| 520 | - ``numerator`` - (Optional; default=None) An element `p` of a |
| 521 | 0- or 1-variate factorial polynomial ring `R`. |
| 522 | - ``denominator_list`` - (Optional; default=None) A list of the form |
| 523 | `[(q_1,e_1), \ldots, (q_m,e_m)]` where the `q_1, \ldots, q_m` are |
| 524 | distinct irreducible elements of `R` and the `e_i` are positive |
| 525 | integers. |
| 526 | - ``quotient`` - (Optional; default=None) An element of a field of |
| 527 | fractions of a factorial ring. |
| 528 | - ``reduce_`` - (Optional; default=True) If True, then represent |
| 529 | `p/(q_1^{e_1} \cdots q_m^{e_m})` in lowest terms. |
| 530 | If False, then won't attempt to divide `p` by any of the `q_i`. |
| 531 | |
| 532 | OUTPUT: |
| 533 | |
| 534 | A REFD instance representing the rational expression |
| 535 | `p/(q_1^{e_1} \cdots q_m^{e_m})`. |
| 536 | To get a non-None output, one of ``numerator`` or ``quotient`` must be |
| 537 | non-None. |
| 538 | |
| 539 | EXAMPLES:: |
| 540 | |
| 541 | sage: from sage.categories.quotient_fields import REFD |
| 542 | sage: f = REFD(64, [(3, 2), (5, 1)]) |
| 543 | sage: print f |
| 544 | (64, [(3, 2), (5, 1)]) |
| 545 | |
| 546 | :: |
| 547 | |
| 548 | sage: R.<x,y> = PolynomialRing(QQ) |
| 549 | sage: qs = [x, 1], [y, 1], [x*y+1, 1] |
| 550 | sage: f = REFD(x, qs) |
| 551 | sage: print f |
| 552 | (1, [(y, 1), (x*y + 1, 1)]) |
| 553 | sage: ff = REFD(x, qs, reduce_=False) |
| 554 | sage: print ff |
| 555 | (x, [(y, 1), (x, 1), (x*y + 1, 1)]) |
| 556 | |
| 557 | :: |
| 558 | |
| 559 | sage: f = REFD(x + y, [(x + y, 1)]) |
| 560 | sage: print f |
| 561 | (1, []) |
| 562 | |
| 563 | :: |
| 564 | |
| 565 | sage: f = 64/45 |
| 566 | sage: print REFD(quotient=f) |
| 567 | (64, [(3, 2), (5, 1)]) |
| 568 | |
| 569 | :: |
| 570 | |
| 571 | sage: R.<x> = PolynomialRing(QQ) |
| 572 | sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) |
| 573 | sage: print REFD(quotient=f) |
| 574 | (5*x^7 - 5*x^6 + 5/3*x^5 - 5/3*x^4 + 2*x^3 - 2/3*x^2 + 1/3*x - 1/3, [(x - 1, 1), (x, 1), (x^2 + 1/3, 1)]) |
| 575 | |
| 576 | :: |
| 577 | |
| 578 | sage: R.<x,y> = PolynomialRing(QQ) |
| 579 | sage: f = 2*y/(5*(x^3 - 1)*(y + 1)) |
| 580 | sage: print REFD(quotient=f) |
| 581 | (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)]) |
| 582 | |
| 583 | Singular throws a 'not implemented' error when trying to factor in |
| 584 | a multivariate polynomial ring over an inexact field :: |
| 585 | |
| 586 | sage: R.<x,y>= PolynomialRing(CC) |
| 587 | sage: f = (x + 1)/(x*y*(x*y + 1)^2) |
| 588 | sage: REFD(quotient=f) |
| 589 | Traceback (most recent call last): |
| 590 | ... |
| 591 | TypeError: Singular error: |
| 592 | ? not implemented |
| 593 | ? error occurred in or before STDIN line 17: `def sage9=factorize(sage8);` |
| 594 | |
| 595 | """ |
| 596 | from sage.rings.polynomial.polynomial_ring import is_PolynomialRing |
| 597 | from sage.rings.polynomial.multi_polynomial import is_MPolynomial |
| 598 | |
| 599 | if quotient is not None: |
| 600 | p = quotient.numerator() |
| 601 | q = quotient.denominator() |
| 602 | R = p.parent() |
| 603 | self._ring = R |
| 604 | if (is_PolynomialRing(R) or is_MPolynomial(R)) and \ |
| 605 | (q in R.base_ring()): |
| 606 | # R is a polynomial ring and q is in its base ring. |
| 607 | # No factoring needed. |
| 608 | self._numerator = quotient |
| 609 | self._denominator_list = [] |
| 610 | else: |
| 611 | # Otherwise factor q. |
| 612 | try: |
| 613 | qs = q.factor() |
| 614 | except NotImplementedError: |
| 615 | # Singular's factor() needs `proof=False`. |
| 616 | qs = q.factor(proof=False) |
| 617 | if qs.unit() != 1: |
| 618 | p = p/qs.unit() |
| 619 | qs = sorted([tuple(t) for t in qs]) # Sort for consitency. |
| 620 | self._numerator = p |
| 621 | self._denominator_list = qs |
| 622 | # Done. No reducing needed as Sage automatically reduces fractions. |
| 623 | return |
| 624 | |
| 625 | if numerator is not None: |
| 626 | R = numerator.parent() |
| 627 | else: |
| 628 | R = None |
| 629 | self._numerator = numerator |
| 630 | self._denominator_list = sorted([tuple(t) for t in denominator_list]) |
| 631 | self._ring = R |
| 632 | |
| 633 | if (numerator is not None) and reduce_: |
| 634 | # Reduce rational expression if possible by canceling |
| 635 | # irreducible factors in denominator_list. |
| 636 | p = numerator |
| 637 | qs = denominator_list |
| 638 | p_new = p |
| 639 | qs_new = [] |
| 640 | for (q, e) in qs: |
| 641 | a = 0 |
| 642 | while q.divides(p_new) and a < e: |
| 643 | a += 1 |
| 644 | p_new = R(p_new/q) |
| 645 | if e - a > 0: |
| 646 | qs_new.append((q,e-a)) |
| 647 | self._numerator = p_new |
| 648 | self._denominator_list = qs_new |
| 649 | |
| 650 | def numerator(self): |
| 651 | r""" |
| 652 | Return the numerator of ``self``. |
| 653 | """ |
| 654 | return self._numerator |
| 655 | |
| 656 | def denominator_list(self): |
| 657 | r""" |
| 658 | Return the denominator list of ``self``. |
| 659 | """ |
| 660 | return self._denominator_list |
| 661 | |
| 662 | def denominator(self): |
| 663 | r""" |
| 664 | Return the denominator of ``self``. |
| 665 | """ |
| 666 | from sage.misc.misc import prod |
| 667 | return prod([q**e for q,e in self.denominator_list()]) |
| 668 | |
| 669 | def ring(self): |
| 670 | r""" |
| 671 | Return the ring of ``self``. |
| 672 | """ |
| 673 | return self._ring |
| 674 | |
| 675 | def list(self): |
| 676 | r""" |
| 677 | Convert ``self`` into a list. |
| 678 | """ |
| 679 | return (self.numerator(), self.denominator_list()) |
| 680 | |
| 681 | def quotient(self): |
| 682 | r""" |
| 683 | Convert ``self`` into a field of fractions element. |
| 684 | """ |
| 685 | return self.numerator()/self.denominator() |
| 686 | |
| 687 | def __str__(self): |
| 688 | return str(self.list()) |
| 689 | |
| 690 | def __eq__(self, other): |
| 691 | r""" |
| 692 | Two REFD instances are equal iff they represent the same |
| 693 | rational expression. |
| 694 | |
| 695 | EXAMPLES:: |
| 696 | |
| 697 | sage: from sage.categories.quotient_fields import REFD |
| 698 | sage: R.<x,y>= PolynomialRing(QQ) |
| 699 | sage: qs = [x,1], [y,1], [x*y+1,1] |
| 700 | sage: f = REFD(x, qs) |
| 701 | sage: ff = REFD(x, qs, reduce_=False) |
| 702 | sage: f == ff |
| 703 | True |
| 704 | sage: g = REFD(y, qs) |
| 705 | sage: g == f |
| 706 | False |
| 707 | """ |
| 708 | return self.quotient() == other.quotient() |
| 709 | |
| 710 | def __ne__(self, other): |
| 711 | r""" |
| 712 | EXAMPLES:: |
| 713 | |
| 714 | sage: from sage.categories.quotient_fields import REFD |
| 715 | sage: R.<x,y>= PolynomialRing(QQ) |
| 716 | sage: qs = [x,1], [y,1], [x*y+1,1] |
| 717 | sage: f = REFD(x, qs) |
| 718 | sage: ff = REFD(x, qs, reduce_=False) |
| 719 | sage: f != ff |
| 720 | False |
| 721 | sage: g = REFD(y, qs) |
| 722 | sage: g != f |
| 723 | True |
| 724 | """ |
| 725 | return not (self == other) |
| 726 | |
| 727 | def __lt__(self, other): |
| 728 | r""" |
| 729 | REFD A is less than REFD B iff the denominator list of A |
| 730 | is shorter than that of B, or |
| 731 | the denominator list lengths are equal and |
| 732 | the denominator of A is less than that of B in the ambient ring, or |
| 733 | the denominator list lengths are equal and the denominators are equal |
| 734 | and the the numerator of A is less than that of B in the ambient ring. |
| 735 | |
| 736 | EXAMPLES:: |
| 737 | |
| 738 | sage: from sage.categories.quotient_fields import REFD |
| 739 | sage: R.<x,y>= PolynomialRing(QQ) |
| 740 | sage: qs = [x,1], [y,1], [x*y+1,1] |
| 741 | sage: f = REFD(x, qs) |
| 742 | sage: ff = REFD(x, qs, reduce_=False) |
| 743 | sage: g = REFD(x + 1, qs) |
| 744 | sage: h = REFD(x + 2, qs) |
| 745 | sage: print f < ff |
| 746 | False |
| 747 | sage: print f < g |
| 748 | True |
| 749 | sage: print g < f |
| 750 | False |
| 751 | sage: print g < h |
| 752 | True |
| 753 | """ |
| 754 | sn = self.numerator() |
| 755 | on = other.numerator() |
| 756 | sdl = self.denominator_list() |
| 757 | odl = other.denominator_list() |
| 758 | sd = self.denominator() |
| 759 | od = other.denominator() |
| 760 | |
| 761 | if self == other: |
| 762 | return False |
| 763 | elif len(sdl) < len(odl) or\ |
| 764 | (len(sdl) == len(odl) and sd < od) or\ |
| 765 | (len(sdl) == len(odl) and sd == od and sn < on): |
| 766 | return True |
| 767 | else: |
| 768 | return False |
| 769 | |
| 770 | def univariate_decomposition(self): |
| 771 | r""" |
| 772 | Return the usual 0- or 1-variate partial fraction decomposition |
| 773 | of ``self`` as a REFDSum instance. |
| 774 | Assume that ``self`` lies in the field of fractions |
| 775 | of a 0- or 1-variate factorial polynomial ring. |
| 776 | |
| 777 | EXAMPLES:: |
| 778 | |
| 779 | No variables:: |
| 780 | |
| 781 | sage: from sage.categories.quotient_fields import REFD |
| 782 | sage: f = 64/45 |
| 783 | sage: decomp = REFD(quotient=f).univariate_decomposition() |
| 784 | sage: print decomp |
| 785 | [(1, []), (2, [(3, 2)]), (1, [(5, 1)])] |
| 786 | sage: print decomp.sum().quotient() == f |
| 787 | True |
| 788 | |
| 789 | One variable:: |
| 790 | |
| 791 | sage: R.<x> = PolynomialRing(QQ) |
| 792 | sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) |
| 793 | sage: print f |
| 794 | (15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 - 2*x^2 + x - 1)/(3*x^4 - 3*x^3 + x^2 - x) |
| 795 | sage: decomp = REFD(quotient=f).univariate_decomposition() |
| 796 | sage: print decomp |
| 797 | [(5*x^3, []), (1, [(x - 1, 1)]), (1, [(x, 1)]), (1/3, [(x^2 + 1/3, 1)])] |
| 798 | sage: print decomp.sum().quotient() == f |
| 799 | True |
| 800 | |
| 801 | One variable over a finite field:: |
| 802 | |
| 803 | sage: R.<x> = PolynomialRing(GF(2)) |
| 804 | sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) |
| 805 | sage: print f |
| 806 | (x^6 + x^4 + 1)/(x^3 + x) |
| 807 | sage: decomp = REFD(quotient=f).univariate_decomposition() |
| 808 | sage: print decomp |
| 809 | [(x^3, []), (1, [(x, 1)]), (x, [(x + 1, 2)])] |
| 810 | sage: print decomp.sum().quotient() == f |
| 811 | True |
| 812 | |
| 813 | One variable over an inexact field:: |
| 814 | |
| 815 | sage: R.<x> = PolynomialRing(CC) |
| 816 | sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) |
| 817 | sage: print f |
| 818 | (15.0000000000000*x^7 - 15.0000000000000*x^6 + 5.00000000000000*x^5 - 5.00000000000000*x^4 + 6.00000000000000*x^3 - 2.00000000000000*x^2 + x - 1.00000000000000)/(3.00000000000000*x^4 - 3.00000000000000*x^3 + x^2 - x) |
| 819 | sage: decomp = REFD(quotient=f).univariate_decomposition() |
| 820 | sage: print decomp |
| 821 | [(5.00000000000000*x^3, []), (1.00000000000000, [(x - 1.00000000000000, 1)]), (-0.288675134594813*I, [(x - 0.577350269189626*I, 1)]), (1.00000000000000, [(x, 1)]), (0.288675134594813*I, [(x + 0.577350269189626*I, 1)])] |
| 822 | sage: print decomp.sum().quotient() == f # Rounding error coming |
| 823 | False |
| 824 | |
| 825 | NOTE:: |
| 826 | |
| 827 | Let `f = p/q` be a rational expression where `p` and `q` lie in a |
| 828 | 0- or 1-variate factorial polynomial ring `R`. |
| 829 | Let `q_1^{e_1} \cdots q_m^{e_m}` be the |
| 830 | unique factorization of `q` in `R` into irreducible factors. |
| 831 | Then `f` can be written uniquely as |
| 832 | |
| 833 | (*) `p_0 + \sum_{i=1}^{m} p_i/ q_i^{e_i}`, |
| 834 | |
| 835 | for some `p_j \in R`. |
| 836 | I call (*) the *usual partial fraction decomposition* of f. |
| 837 | |
| 838 | AUTHORS: |
| 839 | |
| 840 | - Robert Bradshaw (2007-05-31) |
| 841 | - Alexander Raichev (2012-06-25) |
| 842 | """ |
| 843 | from sage.misc.misc import prod |
| 844 | |
| 845 | p = self.numerator() |
| 846 | q = self.denominator() |
| 847 | whole, p = p.quo_rem(q) |
| 848 | qs = self.denominator_list() |
| 849 | decomp = [REFD(whole, [])] |
| 850 | for (a, m) in qs: |
| 851 | numer = p * prod([b**n for (b, n) in qs if b != a])\ |
| 852 | .inverse_mod(a**m) % (a**m) |
| 853 | # The inverse exists because the product and a**m |
| 854 | # are relatively prime. |
| 855 | decomp.append(REFD(numer, [(a, m)])) |
| 856 | return REFDSum(decomp) |
| 857 | |
| 858 | def nullstellensatz_certificate(self): |
| 859 | r""" |
| 860 | Let `[(q_1, e_1), \ldots, (q_m,e_m)]` be the denominator list of |
| 861 | ``self``. |
| 862 | Return a list of polynomials `h_1, \ldots, h_m` in ``self.ring()`` |
| 863 | that satisfies `h_1 q_1 + \cdots + h_m q_m = 1` if it exists. |
| 864 | Otherwise return None. |
| 865 | Only works for multivariate ``self``. |
| 866 | |
| 867 | EXAMPLES:: |
| 868 | |
| 869 | sage: from sage.categories.quotient_fields import REFD |
| 870 | sage: R.<x,y> = PolynomialRing(QQ) |
| 871 | sage: f = 1/(x^2 * (x*y + 1)) |
| 872 | sage: print f |
| 873 | 1/(x^3*y + x^2) |
| 874 | sage: ff = REFD(quotient=f) |
| 875 | sage: L = ff.nullstellensatz_certificate() |
| 876 | sage: print L |
| 877 | [y^2, -x*y + 1] |
| 878 | sage: qs = ff.denominator_list() |
| 879 | sage: sum([L[i]*qs[i][0]**qs[i][1] for i in range(len(qs))]) == 1 |
| 880 | True |
| 881 | |
| 882 | :: |
| 883 | |
| 884 | sage: f = 1/(x*y) |
| 885 | sage: L = REFD(quotient=f).nullstellensatz_certificate() |
| 886 | sage: L is None |
| 887 | True |
| 888 | |
| 889 | """ |
| 890 | |
| 891 | R = self.ring() |
| 892 | qs = self.denominator_list() |
| 893 | J = R.ideal([q**e for q, e in qs]) |
| 894 | if R(1) in J: |
| 895 | return R(1).lift(J) |
| 896 | else: |
| 897 | return None |
| 898 | |
| 899 | def nullstellensatz_decomposition(self): |
| 900 | r""" |
| 901 | Return a Nullstellensatz decomposition of ``self`` as a REFDSum |
| 902 | instance. |
| 903 | Recursive. |
| 904 | Only works for multivariate ``self``. |
| 905 | |
| 906 | EXAMPLES:: |
| 907 | |
| 908 | sage: from sage.categories.quotient_fields import REFD |
| 909 | sage: R.<x,y> = PolynomialRing(QQ) |
| 910 | sage: f = 1/(x*(x*y + 1)) |
| 911 | sage: decomp = REFD(quotient=f).nullstellensatz_decomposition() |
| 912 | sage: print decomp |
| 913 | [(0, []), (1, [(x, 1)]), (-y, [(x*y + 1, 1)])] |
| 914 | sage: decomp.sum().quotient() == f |
| 915 | True |
| 916 | sage: for r in decomp: |
| 917 | ... L = r.nullstellensatz_certificate() |
| 918 | ... L is None |
| 919 | True |
| 920 | True |
| 921 | True |
| 922 | |
| 923 | NOTE:: |
| 924 | |
| 925 | Let `f = p/q` be a rational expression where `p` and `q` lie in a |
| 926 | `d`-variate polynomial ring `K[X]` for some field `K` and `d \ge 1`. |
| 927 | Let `q_1^{e_1} \cdots q_m^{e_m}` be the |
| 928 | unique factorization of `q` in `K[X]` into irreducible factors and |
| 929 | let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i` |
| 930 | over the algebraic closure `L` of `K`. |
| 931 | By [Raic2012]_, `f` can be written as |
| 932 | |
| 933 | (*) `\sum p_A/\prod_{i \in A} q_i^{e_i}`, |
| 934 | |
| 935 | where the `p_A` are in `K[X]` and the sum is taken over all subsets |
| 936 | `A \subseteq \{1, \ldots, m\}` such that |
| 937 | `\cap_{i\in A} T_i \neq \emptyset`. |
| 938 | |
| 939 | I call (*) a *Nullstellensatz decomposition* of `f`. |
| 940 | Nullstellensatz decompositions are not unique. |
| 941 | |
| 942 | The algorithm used comes from [Raic2012]_. |
| 943 | """ |
| 944 | decomp = REFDSum() |
| 945 | L = self.nullstellensatz_certificate() |
| 946 | if L is None: |
| 947 | # No decomposing possible. |
| 948 | decomp = REFDSum([self]) |
| 949 | else: |
| 950 | # Decompose recursively. |
| 951 | p = self.numerator() |
| 952 | qs = self.denominator_list() |
| 953 | m = len(qs) |
| 954 | iteration1 = REFDSum([REFD(p*L[i],[(qs[j][0], qs[j][1]) |
| 955 | for j in range(m) if j != i]) |
| 956 | for i in range(m) if L[i] != 0]) |
| 957 | |
| 958 | # Now decompose each REFD of iteration1. |
| 959 | for r in iteration1: |
| 960 | decomp.extend(r.nullstellensatz_decomposition()) |
| 961 | |
| 962 | # Simplify and return result. |
| 963 | return decomp.combine_like_terms().whole_and_parts() |
| 964 | |
| 965 | def algebraic_dependence_certificate(self): |
| 966 | r""" |
| 967 | Return the ideal `J` of annihilating polynomials for the set |
| 968 | of polynomials ``[q**e for (q,e) in self.denominator_list()]``, |
| 969 | which could be the zero ideal. |
| 970 | The ideal `J` lies in a polynomial ring over the field |
| 971 | ``self.ring().base_ring()`` that has |
| 972 | ``m = len(self.denominator_list())`` indeterminates. |
| 973 | |
| 974 | EXAMPLES: |
| 975 | |
| 976 | sage: from sage.categories.quotient_fields import REFD |
| 977 | sage: R.<x,y> = PolynomialRing(QQ) |
| 978 | sage: f = 1/(x^2 * (x*y + 1) * y^3) |
| 979 | sage: ff = REFD(quotient=f) |
| 980 | sage: J = ff.algebraic_dependence_certificate() |
| 981 | sage: print J |
| 982 | Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 - 6*T2^5 + T2^6) of Multivariate Polynomial Ring in T0, T1, T2 over Rational Field |
| 983 | sage: g = J.gens()[0] |
| 984 | sage: qs = ff.denominator_list() |
| 985 | sage: g(*(q**e for q, e in qs)) == 0 |
| 986 | True |
| 987 | |
| 988 | :: |
| 989 | |
| 990 | sage: f = 1/(x^3 * y^2) |
| 991 | sage: J = REFD(quotient=f).algebraic_dependence_certificate() |
| 992 | sage: print J |
| 993 | Ideal (0) of Multivariate Polynomial Ring in T0, T1 over Rational Field |
| 994 | |
| 995 | """ |
| 996 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |
| 997 | |
| 998 | R = self.ring() |
| 999 | qs = self.denominator_list() |
| 1000 | if not qs: |
| 1001 | return R.ideal() # The zero ideal. |
| 1002 | m = len(qs) |
| 1003 | F = R.base_ring() |
| 1004 | Xs = list(R.gens()) |
| 1005 | d = len(Xs) |
| 1006 | |
| 1007 | # Expand R by 2*m new variables. |
| 1008 | S = 'S' |
| 1009 | while S in [str(x) for x in Xs]: |
| 1010 | S = S + 'S' |
| 1011 | Ss = [S + str(i) for i in range(m)] |
| 1012 | T = 'T' |
| 1013 | while T in [str(x) for x in Xs]: |
| 1014 | T = T + 'T' |
| 1015 | Ts = [T + str(i) for i in range(m)] |
| 1016 | |
| 1017 | Vs = [str(x) for x in Xs] + Ss + Ts |
| 1018 | RR = PolynomialRing(F, Vs) |
| 1019 | Xs = RR.gens()[:d] |
| 1020 | Ss = RR.gens()[d: d + m] |
| 1021 | Ts = RR.gens()[d + m: d + 2 * m] |
| 1022 | |
| 1023 | # Compute the appropriate elimination ideal. |
| 1024 | J = RR.ideal([ Ss[j] - RR(qs[j][0]) for j in range(m)] +\ |
| 1025 | [ Ss[j]**qs[j][1] - Ts[j] for j in range(m)]) |
| 1026 | J = J.elimination_ideal(Xs + Ss) |
| 1027 | |
| 1028 | # Coerce J into the polynomial ring in the indeteminates Ts[m:]. |
| 1029 | # I choose the negdeglex order because i find it useful in my work. |
| 1030 | RRR = PolynomialRing(F, [str(t) for t in Ts], order='negdeglex') |
| 1031 | return RRR.ideal(J) |
| 1032 | |
| 1033 | def algebraic_dependence_decomposition(self, whole_and_parts=True): |
| 1034 | r""" |
| 1035 | Return an algebraic dependence decomposition of ``self`` as a REFDSum |
| 1036 | instance. |
| 1037 | Recursive. |
| 1038 | |
| 1039 | EXAMPLES:: |
| 1040 | |
| 1041 | sage: from sage.categories.quotient_fields import REFD |
| 1042 | sage: R.<x,y> = PolynomialRing(QQ) |
| 1043 | sage: f = 1/(x^2 * (x*y + 1) * y^3) |
| 1044 | sage: decomp = REFD(quotient=f).algebraic_dependence_decomposition() |
| 1045 | sage: print decomp |
| 1046 | [(0, []), (-x, [(x*y + 1, 1)]), (x^2*y^2 - x*y + 1, [(y, 3), (x, 2)])] |
| 1047 | sage: print decomp.sum().quotient() == f |
| 1048 | True |
| 1049 | sage: for r in decomp: |
| 1050 | ... J = r.algebraic_dependence_certificate() |
| 1051 | ... Z = J.ring().ideal() # The zero ideal |
| 1052 | ... J == Z |
| 1053 | True |
| 1054 | True |
| 1055 | True |
| 1056 | |
| 1057 | NOTE:: |
| 1058 | |
| 1059 | Let `f = p/q` be a rational expression where `p` and `q` lie in a |
| 1060 | `d`-variate polynomial ring `K[X]` for some field `K`. |
| 1061 | Let `q_1^{e_1} \cdots q_m^{e_m}` be the |
| 1062 | unique factorization of `q` in `K[X]` into irreducible factors and |
| 1063 | let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i` |
| 1064 | over the algebraic closure `L` of `K`. |
| 1065 | By [Raic2012]_, `f` can be written as |
| 1066 | |
| 1067 | (*) `\sum p_A/\prod_{i \in A} q_i^{b_i}`, |
| 1068 | |
| 1069 | where the `b_i` are positive integers, the `p_A` are in `K[X]`, |
| 1070 | and the sum is taken over all subsets `A \subseteq \{1, \ldots, m\}` |
| 1071 | such that `|A| \le d` and `\{q_i : i\in A\}` is algebraically |
| 1072 | independent. |
| 1073 | |
| 1074 | I call (*) an *algebraic dependence decomposition* of `f`. |
| 1075 | Algebraic dependence decompositions are not unique. |
| 1076 | |
| 1077 | The algorithm used comes from [Raic2012]_. |
| 1078 | """ |
| 1079 | from sage.misc.misc import prod |
| 1080 | |
| 1081 | decomp = REFDSum() |
| 1082 | R = self.ring() |
| 1083 | J = self.algebraic_dependence_certificate() |
| 1084 | if not J: |
| 1085 | # No decomposing possible. |
| 1086 | decomp = REFDSum([self]) |
| 1087 | else: |
| 1088 | # Decompose recursively. |
| 1089 | p = self.numerator() |
| 1090 | qs = self.denominator_list() |
| 1091 | m = len(qs) |
| 1092 | g = J.gens()[0] # An annihilating polynomial for the qs. |
| 1093 | new_vars = J.ring().gens() |
| 1094 | # Note that each new_vars[j] corresponds to qs[j] such that |
| 1095 | # g([q**e for q,e in qs]) = 0. |
| 1096 | # Assuming here that g.parent() has negdeglex term order |
| 1097 | # so that g.lt() is indeed the monomial we want below. |
| 1098 | # Use g to rewrite r into a sum of REFDs, |
| 1099 | # each with < m distinct denominator factors. |
| 1100 | gg = (g.lt() - g)/(g.lc()) |
| 1101 | numers = map(prod, zip(gg.coefficients(), gg.monomials())) |
| 1102 | e = list(g.lt().exponents())[0:m] |
| 1103 | denoms = [(new_vars[j], e[0][j] + 1) for j in range(m)] |
| 1104 | # Write r in terms of new_vars, |
| 1105 | # cancel factors in the denominator, and combine like terms. |
| 1106 | iteration1_temp = REFDSum([REFD(a, denoms) for a in numers]).\ |
| 1107 | combine_like_terms() |
| 1108 | # Substitute in qs. |
| 1109 | qpowsub = dict([(new_vars[j],qs[j][0]**qs[j][1]) |
| 1110 | for j in range(m)]) |
| 1111 | iteration1 = REFDSum() |
| 1112 | for r in iteration1_temp: |
| 1113 | num1 = p*g.parent()(r.numerator()).subs(qpowsub) |
| 1114 | denoms1 = [] |
| 1115 | for q,e in r.denominator_list(): |
| 1116 | j = new_vars.index(q) |
| 1117 | denoms1.append((qs[j][0], qs[j][1]*e)) |
| 1118 | iteration1.append(REFD(num1, denoms1)) |
| 1119 | # Now decompose each REFD of iteration1. |
| 1120 | for r in iteration1: |
| 1121 | decomp.extend(r.algebraic_dependence_decomposition()) |
| 1122 | |
| 1123 | # Simplify and return result. |
| 1124 | return decomp.combine_like_terms().whole_and_parts() |
| 1125 | |
| 1126 | def leinartas_decomposition(self): |
| 1127 | r""" |
| 1128 | Return a Leinartas decomposition of ``self`` as a REFDSum instance. |
| 1129 | |
| 1130 | EXAMPLES:: |
| 1131 | |
| 1132 | sage: from sage.categories.quotient_fields import REFD |
| 1133 | sage: R.<x,y> = PolynomialRing(QQ) |
| 1134 | sage: f = 1/x + 1/y + 1/(x*y + 1) |
| 1135 | sage: decomp = REFD(quotient=f).leinartas_decomposition() |
| 1136 | sage: print decomp |
| 1137 | [(0, []), (1, [(x*y + 1, 1)]), (x + y, [(y, 1), (x, 1)])] |
| 1138 | sage: print decomp.sum().quotient() == f |
| 1139 | True |
| 1140 | sage: for r in decomp: |
| 1141 | ... L = r.nullstellensatz_certificate() |
| 1142 | ... print L is None |
| 1143 | ... J = r.algebraic_dependence_certificate() |
| 1144 | ... Z = J.ring().ideal() # The zero ideal |
| 1145 | ... print J == Z |
| 1146 | True |
| 1147 | True |
| 1148 | True |
| 1149 | True |
| 1150 | True |
| 1151 | True |
| 1152 | |
| 1153 | :: |
| 1154 | |
| 1155 | sage: R.<x,y,z>= PolynomialRing(GF(2, 'a')) |
| 1156 | sage: f = 1/(x * y * z * (x*y + z)) |
| 1157 | sage: decomp = REFD(quotient=f).leinartas_decomposition() |
| 1158 | sage: print decomp |
| 1159 | [(0, []), (1, [(z, 2), (x*y + z, 1)]), (1, [(z, 2), (y, 1), (x, 1)])] |
| 1160 | sage: print decomp.sum().quotient() == f |
| 1161 | True |
| 1162 | |
| 1163 | NOTE:: |
| 1164 | |
| 1165 | Let `f = p/q` be a rational expression where `p` and `q` lie in a |
| 1166 | `d`-variate polynomial ring `K[X]` for some field `K`. |
| 1167 | Let `q_1^{e_1} \cdots q_m^{e_m}` be the |
| 1168 | unique factorization of `q` in `K[X]` into irreducible factors and |
| 1169 | let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i` |
| 1170 | over the algebraic closure `L` of `K`. |
| 1171 | By [Raic2012]_, `f` can be written as |
| 1172 | |
| 1173 | (*) `\sum p_A/\prod_{i \in A} q_i^{b_i}`, |
| 1174 | |
| 1175 | where the `b_i` are positive integers, the `p_A` are in `K[X]`, |
| 1176 | and the sum is taken over all subsets `A \subseteq \{1, \ldots, m\}` |
| 1177 | such that |
| 1178 | (1) `|A| \le d`, |
| 1179 | (2) `\cap_{i\in A} T_i \neq \emptyset`, and |
| 1180 | (3) `\{q_i : i\in A\}` is algebraically independent. |
| 1181 | |
| 1182 | In particular, any rational expression in `d` variables |
| 1183 | can be represented as a sum of rational expressions |
| 1184 | whose denominators each contain at most `d` distinct irreducible |
| 1185 | factors. |
| 1186 | |
| 1187 | I call (*) a *Leinartas decomposition* of `f`. |
| 1188 | Leinartas decompositions are not unique. |
| 1189 | |
| 1190 | The algorithm used comes from [Raic2012]_. |
| 1191 | """ |
| 1192 | d = len(self.ring().gens()) |
| 1193 | if d == 1: |
| 1194 | # Sage's lift() function doesn't work in |
| 1195 | # univariate polynomial rings. |
| 1196 | # So nullstellensatz_decomposition() won't work. |
| 1197 | # So use algebraic_dependence_decomposition(), |
| 1198 | # which is sufficient. |
| 1199 | temp = REFDSum([self]) |
| 1200 | else: |
| 1201 | temp = self.nullstellensatz_decomposition() |
| 1202 | decomp = REFDSum() |
| 1203 | for r in temp: |
| 1204 | decomp.extend(r.algebraic_dependence_decomposition()) |
| 1205 | |
| 1206 | # Simplify and return result. |
| 1207 | return decomp.combine_like_terms().whole_and_parts() |
| 1208 | |
| 1209 | class REFDSum(list): |
| 1210 | r""" |
| 1211 | A list representing the sum of REFD objects with distinct |
| 1212 | denominator lists. |
| 1213 | |
| 1214 | AUTHORS: |
| 1215 | |
| 1216 | - Alexander Raichev (2012-06-25) |
| 1217 | """ |
| 1218 | def __str__(self): |
| 1219 | return str([r.list() for r in self]) |
| 1220 | |
| 1221 | def __eq__(self, other): |
| 1222 | return sorted(self) == sorted(other) |
| 1223 | |
| 1224 | def __ne__(self, other): |
| 1225 | return not (self == other) |
| 1226 | |
| 1227 | def ring(self): |
| 1228 | r""" |
| 1229 | Return the ring of ``self``. |
| 1230 | """ |
| 1231 | if self: |
| 1232 | return self[0].ring() |
| 1233 | else: |
| 1234 | return None |
| 1235 | |
| 1236 | def whole_and_parts(self): |
| 1237 | r""" |
| 1238 | Rewrite ``self`` as a REFDSum of a (possibly zero) polynomial |
| 1239 | REFD followed by reduced rational expression REFDs. |
| 1240 | Only useful for multivariate decompositions. |
| 1241 | |
| 1242 | EXAMPLES:: |
| 1243 | |
| 1244 | sage: from sage.categories.quotient_fields import REFD |
| 1245 | sage: from sage.categories.quotient_fields import REFDSum |
| 1246 | sage: R = PolynomialRing(QQ, 'x, y') |
| 1247 | sage: x, y = R.gens() |
| 1248 | sage: f = x**2 + 3*y + 1/x + 1/y |
| 1249 | sage: f = REFD(quotient=f) |
| 1250 | sage: print f |
| 1251 | (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)]) |
| 1252 | sage: print REFDSum([f]).whole_and_parts() |
| 1253 | [(x^2 + 3*y, []), (x + y, [(y, 1), (x, 1)])] |
| 1254 | |
| 1255 | """ |
| 1256 | whole = 0 |
| 1257 | parts = [] |
| 1258 | R = self.ring() |
| 1259 | for r in self: |
| 1260 | # Since r has already passed through REFD.__init__()'s reducing |
| 1261 | # procedure, r is already in lowest terms. |
| 1262 | # Check if can write r as a mixed fraction: whole + fraction. |
| 1263 | p = r.numerator() |
| 1264 | q = r.denominator() |
| 1265 | if q == 1: |
| 1266 | # r is already a whole polynomial |
| 1267 | whole += p |
| 1268 | else: |
| 1269 | # r is a fraction |
| 1270 | num = R(r.numerator()) # Coerce into R in case its a constant. |
| 1271 | denom = r.denominator() |
| 1272 | a, b = num.quo_rem(denom) |
| 1273 | whole += a |
| 1274 | parts.append(REFD(b, r.denominator_list(), reduce_=False)) |
| 1275 | return REFDSum([REFD(whole, ())] + parts) |
| 1276 | |
| 1277 | def combine_like_terms(self): |
| 1278 | r""" |
| 1279 | Combine terms in ``self`` with the same denominator. |
| 1280 | Only useful for multivariate decompositions. |
| 1281 | |
| 1282 | EXAMPLES:: |
| 1283 | |
| 1284 | sage: from sage.categories.quotient_fields import REFD |
| 1285 | sage: from sage.categories.quotient_fields import REFDSum |
| 1286 | sage: R.<x,y>= PolynomialRing(QQ) |
| 1287 | sage: f = REFD(quotient=1/(x * y * (x*y + 1))) |
| 1288 | sage: g = REFD(quotient=x/(x * y * (x*y + 1))) |
| 1289 | sage: s = REFDSum([f, g, f]) |
| 1290 | sage: t = s.combine_like_terms() |
| 1291 | sage: print s |
| 1292 | [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (1, [(y, 1), (x*y + 1, 1)]), (1, [(y, 1), (x, 1), (x*y + 1, 1)])] |
| 1293 | sage: print t |
| 1294 | [(1, [(y, 1), (x*y + 1, 1)]), (2, [(y, 1), (x, 1), (x*y + 1, 1)])] |
| 1295 | |
| 1296 | """ |
| 1297 | if not self: |
| 1298 | return self |
| 1299 | |
| 1300 | # Combine like terms. |
| 1301 | refds = sorted(self) |
| 1302 | new_refds = [] |
| 1303 | temp = refds[0] |
| 1304 | for f in refds[1:]: |
| 1305 | if temp.denominator_list() == f.denominator_list(): |
| 1306 | # Add f to temp. |
| 1307 | num = temp.numerator() + f.numerator() |
| 1308 | temp = REFD(num, temp.denominator_list()) |
| 1309 | else: |
| 1310 | # Append temp to new_refds and update temp. |
| 1311 | new_refds.append(temp) |
| 1312 | temp = f |
| 1313 | new_refds.append(temp) |
| 1314 | return REFDSum(new_refds) |
| 1315 | |
| 1316 | def sum(self): |
| 1317 | r""" |
| 1318 | Return the sum of the REFDs in ``self`` as a REFD. |
| 1319 | |
| 1320 | EXAMPLES:: |
| 1321 | |
| 1322 | sage: from sage.categories.quotient_fields import REFD |
| 1323 | sage: from sage.categories.quotient_fields import REFDSum |
| 1324 | sage: R.<x,y> = PolynomialRing(QQ) |
| 1325 | sage: qs = (x, 1), (y, 1), (x*y + 1, 1) |
| 1326 | sage: f = REFD(quotient=1/(x * y * (x*y + 1))) |
| 1327 | sage: g = REFD(quotient=x*y/(x * y * (x*y + 1))) |
| 1328 | sage: print REFDSum([f, g]).sum() |
| 1329 | (1, [(y, 1), (x, 1)]) |
| 1330 | |
| 1331 | """ |
| 1332 | if not self: |
| 1333 | return self |
| 1334 | |
| 1335 | # Compute the sum's numerator and denominator. |
| 1336 | summy = sum((f.quotient() for f in self)) |
| 1337 | num = summy.numerator() |
| 1338 | denom = summy.denominator() |
| 1339 | if denom == 1: |
| 1340 | # Done |
| 1341 | return REFD(num, []) |
| 1342 | |
| 1343 | # Compute the irreducible factors of denom. |
| 1344 | # Could use the factor() command here, but that might be slow. |
| 1345 | factors = [] |
| 1346 | for f in self: |
| 1347 | factors.extend([q for q,e in f.denominator_list()]) |
| 1348 | |
| 1349 | # Eliminate repeats from factors and sort |
| 1350 | factors = sorted(list(set(factors))) |
| 1351 | |
| 1352 | # The irreducible factors of denom lie in factors. |
| 1353 | # Use this fact to build a denominator list for denom. |
| 1354 | denom_list = [] |
| 1355 | for q in factors: |
| 1356 | e = 0 |
| 1357 | quo, rem = denom.quo_rem(q) |
| 1358 | while rem == 0: |
| 1359 | e += 1 |
| 1360 | denom = quo |
| 1361 | quo, rem = denom.quo_rem(q) |
| 1362 | if e > 0: |
| 1363 | denom_list.append((q, e)) |
| 1364 | return REFD(num, denom_list, reduce_=False) |
| 1365 | No newline at end of file |