Ticket #12068: trac_12068-numer_denom_ginac-folded-fh.patch

File trac_12068-numer_denom_ginac-folded-fh.patch, 13.7 KB (added by hivert, 10 years ago)
  • sage/libs/ginac/decl.pxi

    # HG changeset patch
    # User Florent Hivert <Florent.Hivert@univ-rouen.fr>
    # Date 1322130837 -3600
    # Node ID 1911f8e3cefb2ef0e7f99dd9c21ab43db4e5ca28
    # Parent  fcbb634bae62c51df1fa4fc0c2af133cbe47e9f1
    #12068: Numerator for symbolic expression should'nt use maxima
    * * *
    trac 12068: denominator() for symbolic expressions recognizes more cases
    
    diff --git a/sage/libs/ginac/decl.pxi b/sage/libs/ginac/decl.pxi
    a b cdef extern from "ginac_wrap.h": 
    3636
    3737    object GSymbol_to_str "_to_PyString<symbol>"(GSymbol *s)
    3838
    39     ctypedef struct GExPair "std::pair<ex, ex>"
     39    ctypedef struct GExPair "std::pair<ex, ex>":
     40        pass
    4041    ctypedef struct GExMap "exmap":
    4142        void insert(GExPair e)
    4243
    cdef extern from "ginac_wrap.h": 
    6970        GEx subs_map "subs" (GExMap map) except +
    7071        GEx coeff(GEx expr, int n)    except +
    7172        GEx lcoeff(GEx expr)          except +
    72         GEx tcoeff(GEx expr)          except +       
     73        GEx tcoeff(GEx expr)          except +
     74        GEx numer()                   except +
     75        GEx denom()                   except +
     76        GEx numer_denom()             except +
    7377        int degree(GEx expr)          except +
    7478        int ldegree(GEx expr)         except +
    7579        GEx rhs()                     except +
    cdef extern from "ginac_wrap.h": 
    8589
    8690    GExPair make_pair "std::make_pair" (GEx, GEx)
    8791
     92    ctypedef struct GNumeric "numeric":
     93        bint is_positive() except +
     94        bint is_negative() except +
     95
    8896    # Numericals
    8997    bint is_a_numeric "is_a<numeric>" (GEx e)
     98    GNumeric ex_to_numeric "ex_to<numeric>" (GEx e)
    9099    # given a GEx that is known to be a numeric, return reference to
    91100    # the underlying PyObject*.
    92101    object py_object_from_numeric(GEx e)     except +
    cdef extern from "ginac_wrap.h": 
    143152    unsigned info_nonnegint     "GiNaC::info_flags::nonnegint"
    144153    unsigned info_even          "GiNaC::info_flags::even"
    145154    unsigned info_odd           "GiNaC::info_flags::odd"
     155    unsigned info_rational_function "GiNaC::info_flags::rational_function"
    146156
    147157    # Constants
    148158    GEx g_Pi "Pi"
    cdef extern from "ginac_wrap.h": 
    216226    bint is_a_function "is_a<function>" (GEx e)
    217227    bint is_a_ncmul "is_a<ncmul>" (GEx e)
    218228
    219 
    220229    # Arithmetic
    221230    int ginac_error()
    222231    GEx gadd "ADD_WRAP" (GEx left, GEx right) except +
  • sage/symbolic/expression.pyx

    diff --git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
    a b cdef class Expression(CommutativeRingEle 
    18671867        assured that if True or False is returned (and proof is False) then
    18681868        the answer is correct.
    18691869
    1870         INPUT::
     1870        INPUT:
    18711871       
    18721872           ntests -- (default 20) the number of iterations to run
    18731873           domain -- (optional) the domain from which to draw the random values
    cdef class Expression(CommutativeRingEle 
    64156415            ((x - 1)*x + y^2)/(x^2 - 7) + (b + c)/a + 1/(x + 1)
    64166416        """
    64176417        return self.parent()(self._maxima_().combine())
    6418    
    6419     def numerator(self):
    6420         """
    6421         Returns the numerator of this symbolic expression.  If the
    6422         expression is not a quotient, then this will return the
    6423         expression itself.
    6424        
    6425         EXAMPLES::
    6426        
     6418
     6419    def numerator(self, bint normalize = True):
     6420        """
     6421        Returns the numerator of this symbolic expression
     6422
     6423        INPUT:
     6424
     6425        - ``normalize`` -- (default: ``True``) a boolean.
     6426
     6427        If ``normalize`` is ``True``, the expression is first normalized to
     6428        have it as a fraction before getting the numerator.
     6429
     6430        If ``normalize`` is ``False``, the expression is kept and if it is not
     6431        a quotient, then this will return the expression itself.
     6432
     6433        EXAMPLES::
     6434
    64276435            sage: a, x, y = var('a,x,y')
    64286436            sage: f = x*(x-a)/((x^2 - y)*(x-a)); f
    64296437            x/(x^2 - y)
    cdef class Expression(CommutativeRingEle 
    64316439            x
    64326440            sage: f.denominator()
    64336441            x^2 - y
     6442            sage: f.numerator(normalize=False)
     6443            x
     6444            sage: f.denominator(normalize=False)
     6445            x^2 - y
    64346446
    64356447            sage: y = var('y')
    64366448            sage: g = x + y/(x + 2); g
    64376449            x + y/(x + 2)
    64386450            sage: g.numerator()
     6451            x^2 + 2*x + y
     6452            sage: g.denominator()
     6453            x + 2
     6454            sage: g.numerator(normalize=False)
    64396455            x + y/(x + 2)
    6440             sage: g.denominator()
     6456            sage: g.denominator(normalize=False)
    64416457            1
    64426458
    6443         """
    6444         return self.parent()(self._maxima_().num())
    6445 
    6446     def denominator(self):
    6447         """
    6448         Returns the denominator of this symbolic expression.  If the
    6449         expression is not a quotient, then this will just return 1.
    6450        
    6451         EXAMPLES::
    6452        
     6459        TESTS::
     6460
     6461            sage: ((x+y)^2/(x-y)^3*x^3).numerator(normalize=False)
     6462            (x + y)^2*x^3
     6463            sage: ((x+y)^2*x^3).numerator(normalize=False)
     6464            (x + y)^2*x^3
     6465            sage: (y/x^3).numerator(normalize=False)
     6466            y
     6467            sage: t = y/x^3/(x+y)^(1/2); t
     6468            y/(sqrt(x + y)*x^3)
     6469            sage: t.numerator(normalize=False)
     6470            y
     6471            sage: (1/x^3).numerator(normalize=False)
     6472            1
     6473            sage: (x^3).numerator(normalize=False)
     6474            x^3
     6475            sage: (y*x^sin(x)).numerator(normalize=False)
     6476            Traceback (most recent call last):
     6477            ...
     6478            TypeError: self is not a rational expression
     6479        """
     6480        cdef GExVector vec
     6481        cdef GEx oper, power
     6482        if normalize:
     6483            return new_Expression_from_GEx(self._parent, self._gobj.numer())
     6484        elif is_a_mul(self._gobj):
     6485            for i from 0 <= i < self._gobj.nops():
     6486                oper = self._gobj.op(i)
     6487                if not is_a_power(oper):
     6488                    vec.push_back(oper)
     6489                else:
     6490                    power = oper.op(1)
     6491                    if not is_a_numeric(power):
     6492                        raise TypeError, "self is not a rational expression"
     6493                    elif ex_to_numeric(power).is_positive():
     6494                        vec.push_back(oper)
     6495            return new_Expression_from_GEx(self._parent,
     6496                                           g_mul_construct(vec, True))
     6497        elif is_a_power(self._gobj):
     6498            power = self._gobj.op(1)
     6499            if is_a_numeric(power) and ex_to_numeric(power).is_negative():
     6500                return self._parent.one()
     6501        return self
     6502
     6503    def denominator(self, bint normalize=True):
     6504        """
     6505        Returns the denominator of this symbolic expression
     6506
     6507        INPUT:
     6508
     6509        - ``normalize`` -- (default: ``True``) a boolean.
     6510
     6511        If ``normalize`` is ``True``, the expression is first normalized to
     6512        have it as a fraction before getting the denominator.
     6513
     6514        If ``normalize`` is ``False``, the expression is kept and if it is not
     6515        a quotient, then this will just return 1.
     6516
     6517        EXAMPLES::
     6518
    64536519            sage: x, y, z, theta = var('x, y, z, theta')
    64546520            sage: f = (sqrt(x) + sqrt(y) + sqrt(z))/(x^10 - y^10 - sqrt(theta))
     6521            sage: f.numerator()
     6522            sqrt(x) + sqrt(y) + sqrt(z)
    64556523            sage: f.denominator()
     6524            -sqrt(theta) + x^10 - y^10
     6525
     6526            sage: f.numerator(normalize=False)
     6527            -(sqrt(x) + sqrt(y) + sqrt(z))
     6528            sage: f.denominator(normalize=False)
    64566529            sqrt(theta) - x^10 + y^10
    64576530
    64586531            sage: y = var('y')
    64596532            sage: g = x + y/(x + 2); g
    64606533            x + y/(x + 2)
    6461             sage: g.numerator()
     6534            sage: g.numerator(normalize=False)
    64626535            x + y/(x + 2)
    6463             sage: g.denominator()
     6536            sage: g.denominator(normalize=False)
    64646537            1
    6465         """
    6466         return self.parent()(self._maxima_().denom())
     6538
     6539        TESTS::
     6540
     6541            sage: ((x+y)^2/(x-y)^3*x^3).denominator(normalize=False)
     6542            (x - y)^3
     6543            sage: ((x+y)^2*x^3).denominator(normalize=False)
     6544            1
     6545            sage: (y/x^3).denominator(normalize=False)
     6546            x^3
     6547            sage: t = y/x^3/(x+y)^(1/2); t
     6548            y/(sqrt(x + y)*x^3)
     6549            sage: t.denominator(normalize=False)
     6550            sqrt(x + y)*x^3
     6551            sage: (1/x^3).denominator(normalize=False)
     6552            x^3
     6553            sage: (x^3).denominator(normalize=False)
     6554            1
     6555            sage: (y*x^sin(x)).denominator(normalize=False)
     6556            Traceback (most recent call last):
     6557            ...
     6558            TypeError: self is not a rational expression
     6559        """
     6560        cdef GExVector vec
     6561        cdef GEx oper, ex, power
     6562        if normalize:
     6563            return new_Expression_from_GEx(self._parent, self._gobj.denom())
     6564        elif is_a_mul(self._gobj):
     6565            for i from 0 <= i < self._gobj.nops():
     6566                oper = self._gobj.op(i)
     6567                if is_a_power(oper):
     6568                    ex = oper.op(0)
     6569                    power = oper.op(1)
     6570                    if not is_a_numeric(power):
     6571                        raise TypeError, "self is not a rational expression"
     6572                    elif ex_to_numeric(power).is_negative():
     6573                        vec.push_back(g_pow(ex, g_abs(power)))
     6574            return new_Expression_from_GEx(self._parent,
     6575                                           g_mul_construct(vec, False))
     6576        elif is_a_power(self._gobj):
     6577            power = self._gobj.op(1)
     6578            if is_a_numeric(power) and ex_to_numeric(power).is_negative():
     6579                return new_Expression_from_GEx(self._parent,
     6580                        g_pow(self._gobj.op(0), g_abs(power)))
     6581
     6582        return self._parent.one()
     6583
     6584    def numerator_denominator(self, bint normalize=True):
     6585        """
     6586        Returns the numerator and the denominator of this symbolic expression
     6587
     6588        INPUT:
     6589
     6590        - ``normalize`` -- (default: ``True``) a boolean.
     6591
     6592        If ``normalize`` is ``True``, the expression is first normalized to
     6593        have it as a fraction before getting the numerator and denominator.
     6594
     6595        If ``normalize`` is ``False``, the expression is kept and if it is not
     6596        a quotient, then this will return the expression itself together with
     6597        1.
     6598
     6599        EXAMPLE::
     6600
     6601            sage: x, y, a = var("x y a")
     6602            sage: ((x+y)^2/(x-y)^3*x^3).numerator_denominator()
     6603            ((x + y)^2*x^3, (x - y)^3)
     6604
     6605            sage: ((x+y)^2/(x-y)^3*x^3).numerator_denominator(False)
     6606            ((x + y)^2*x^3, (x - y)^3)
     6607
     6608            sage: g = x + y/(x + 2)
     6609            sage: g.numerator_denominator()
     6610            (x^2 + 2*x + y, x + 2)
     6611            sage: g.numerator_denominator(normalize=False)
     6612            (x + y/(x + 2), 1)
     6613
     6614            sage: g = x^2*(x + 2)
     6615            sage: g.numerator_denominator()
     6616            ((x + 2)*x^2, 1)
     6617            sage: g.numerator_denominator(normalize=False)
     6618            ((x + 2)*x^2, 1)
     6619
     6620        TESTS::
     6621
     6622            sage: ((x+y)^2/(x-y)^3*x^3).numerator_denominator(normalize=False)
     6623            ((x + y)^2*x^3, (x - y)^3)
     6624            sage: ((x+y)^2*x^3).numerator_denominator(normalize=False)
     6625            ((x + y)^2*x^3, 1)
     6626            sage: (y/x^3).numerator_denominator(normalize=False)
     6627            (y, x^3)
     6628            sage: t = y/x^3/(x+y)^(1/2); t
     6629            y/(sqrt(x + y)*x^3)
     6630            sage: t.numerator_denominator(normalize=False)
     6631            (y, sqrt(x + y)*x^3)
     6632            sage: (1/x^3).numerator_denominator(normalize=False)
     6633            (1, x^3)
     6634            sage: (x^3).numerator_denominator(normalize=False)
     6635            (x^3, 1)
     6636            sage: (y*x^sin(x)).numerator_denominator(normalize=False)
     6637            Traceback (most recent call last):
     6638            ...
     6639            TypeError: self is not a rational expression
     6640        """
     6641        cdef GExVector vecnumer, vecdenom
     6642        cdef GEx oper, ex, power
     6643        cdef int py_pow
     6644        cdef GNumeric power_num
     6645        if normalize:
     6646            ex = self._gobj.numer_denom()
     6647            return (new_Expression_from_GEx(self._parent, ex.op(0)),
     6648                    new_Expression_from_GEx(self._parent, ex.op(1)))
     6649        elif is_a_mul(self._gobj):
     6650            for i from 0 <= i < self._gobj.nops():
     6651                oper = self._gobj.op(i)
     6652                if is_a_power(oper):   # oper = ex^power
     6653                    ex = oper.op(0)
     6654                    power = oper.op(1)
     6655                    if not is_a_numeric(power):
     6656                        raise TypeError, "self is not a rational expression"
     6657                    elif is_a_numeric(power):
     6658                        power_num = ex_to_numeric(power)
     6659                        if power_num.is_positive():
     6660                            vecnumer.push_back(oper)
     6661                        else:
     6662                            vecdenom.push_back(g_pow(ex, g_abs(power)))
     6663                else:
     6664                    vecnumer.push_back(oper)
     6665            return (new_Expression_from_GEx(self._parent,
     6666                                            g_mul_construct(vecnumer, False)),
     6667                    new_Expression_from_GEx(self._parent,
     6668                                            g_mul_construct(vecdenom, False)))
     6669        elif is_a_power(self._gobj):
     6670            power = self._gobj.op(1)
     6671            if is_a_numeric(power) and ex_to_numeric(power).is_positive():
     6672                return (self, self._parent.one())
     6673            else:
     6674                return (self._parent.one(),
     6675                        new_Expression_from_GEx(self._parent,
     6676                               g_pow(self._gobj.op(0), g_abs(power))))
     6677        else:
     6678            return (self, self._parent.one())
    64676679
    64686680    def partial_fraction(self, var=None):
    64696681        r"""