Ticket #7748: trac_7748-incomplete_gamma.patch

File trac_7748-incomplete_gamma.patch, 6.5 KB (added by fstan, 12 years ago)

Make incomplete gamma function symbolic.

  • sage/functions/all.py

    # HG changeset patch
    # User Flavia Stan <flavia.stan@gmail.com>
    # Date 1265542370 -3600
    # Node ID 4c0e281ec210347ab2b451e54e294602b7fa75d1
    # Parent  76deac46abf330fa30d2ec1eb6734acaa93f079b
    Trac 7748: Make incomplete gamma function symbolic.
    
    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    1717
    1818from other import ( ceil, floor, gamma, factorial,
    1919                    abs_symbolic, erf, sqrt,
     20                    gamma_inc, incomplete_gamma,
    2021                    real_part, real,
    2122                    imag_part, imag, imaginary, conjugate)
    2223
     
    2425
    2526
    2627from transcendental import (exponential_integral_1,
    27                             gamma_inc, incomplete_gamma,
    2828                            zeta, zeta_symmetric,
    2929                            Li, Ei,
    3030                            dickman_rho)
  • sage/functions/other.py

    diff --git a/sage/functions/other.py b/sage/functions/other.py
    a b  
    55from sage.symbolic.expression import Expression
    66from sage.libs.pari.gen import pari
    77from sage.symbolic.all import SR
    8 from sage.rings.all import Integer, Rational, RealField, CC, RR
     8from sage.rings.all import Integer, Rational, RealField, CC, RR, \
     9     is_ComplexNumber, ComplexField
    910from sage.misc.latex import latex
    1011import math
    1112
     13import sage.structure.element
     14coercion_model = sage.structure.element.get_coercion_model()
     15
     16from sage.structure.coerce import parent
     17from sage.symbolic.constants import pi
     18from sage.symbolic.function import is_inexact
     19from sage.functions.log import exp
     20from sage.functions.transcendental import Ei
     21
    1222one_half = ~SR(2)
    1323
    1424class Function_erf(BuiltinFunction):
     
    538548
    539549gamma = Function_gamma()
    540550
     551class Function_gamma_inc(BuiltinFunction):
     552    def __init__(self):
     553        """
     554        The incomplete gamma function.
     555
     556        EXAMPLES::
     557
     558            sage: gamma_inc(CDF(0,1), 3)
     559            0.00320857499337 + 0.0124061858119*I
     560            sage: gamma_inc(RDF(1), 3)
     561            0.0497870683678639
     562            sage: gamma_inc(3,2)
     563            gamma(3, 2)
     564            sage: gamma_inc(x,0)
     565            gamma(x)
     566            sage: latex(gamma_inc(3,2))
     567            \Gamma\left(3, 2\right)
     568            sage: loads(dumps((gamma_inc(3,2))))
     569            gamma(3, 2)
     570        """
     571        BuiltinFunction.__init__(self, "gamma", nargs=2, latex_name=r"\Gamma")
     572
     573    def _eval_(self, x, y):
     574        """
     575        EXAMPLES::
     576
     577            sage: gamma_inc(2.,0)
     578            1.00000000000000
     579            sage: gamma_inc(2,0)
     580            1
     581            sage: gamma_inc(1/2,2)
     582            -(erf(sqrt(2)) - 1)*sqrt(pi)
     583            sage: gamma_inc(1/2,1)
     584            -(erf(1) - 1)*sqrt(pi)
     585            sage: gamma_inc(1/2,0)
     586            sqrt(pi)
     587            sage: gamma_inc(x,0)
     588            gamma(x)
     589            sage: gamma_inc(1,2)
     590            e^(-2)
     591            sage: gamma_inc(0,2)
     592            -Ei(-2)
     593        """
     594        if not isinstance(x, Expression) and not isinstance(y, Expression) and \
     595               (is_inexact(x) or is_inexact(y)):
     596            x, y = coercion_model.canonical_coercion(x, y)
     597            return self._evalf_(x, y, parent(x))
     598
     599        if y == 0:
     600            return gamma(x)
     601        if x == 1:
     602            return exp(-y)
     603        if x == 0:
     604            return -Ei(-y)
     605        if x == Rational(1)/2: #only for x>0
     606            return sqrt(pi)*(1-erf(sqrt(y)))
     607        return None
     608
     609    def _evalf_(self, x, y, parent=None):
     610        """
     611        EXAMPLES::
     612
     613            sage: gamma_inc(0,2)
     614            -Ei(-2)
     615            sage: gamma_inc(0,2.)
     616            0.0489005107080611
     617            sage: gamma_inc(3,2).n()
     618            1.35335283236613
     619        """
     620        try:
     621            return x.gamma_inc(y)
     622        except AttributeError:
     623            if not (is_ComplexNumber(x)):
     624                if is_ComplexNumber(y):
     625                    C = y.parent()
     626                else:
     627                    C = ComplexField()
     628                    x = C(x)
     629            return x.gamma_inc(y)
     630           
     631# synonym.
     632incomplete_gamma = gamma_inc=Function_gamma_inc()
     633
     634
    541635class Function_factorial(GinacFunction):
    542636    def __init__(self):
    543637        r"""
  • sage/functions/transcendental.py

    diff --git a/sage/functions/transcendental.py b/sage/functions/transcendental.py
    a b  
    100100    else:
    101101        return [float(z) for z in pari(x).eint1(n)]
    102102
    103 def gamma_inc(s, t):
    104     """
    105     Incomplete Gamma function Gamma(s,t).
    106    
    107     EXAMPLES::
    108    
    109         sage: gamma_inc(CDF(0,1), 3)
    110         0.00320857499337 + 0.0124061858119*I
    111         sage: gamma_inc(3, 3)
    112         0.846380162253687
    113         sage: gamma_inc(RDF(1), 3)
    114         0.0497870683678639
    115     """
    116     try:
    117         return s.gamma_inc(t)
    118     except AttributeError:
    119         if not (is_ComplexNumber(s)):
    120             if is_ComplexNumber(t):
    121                 C = t.parent()
    122             else:
    123                 C = ComplexField()
    124             s = C(s)
    125         return s.gamma_inc(t)
    126 
    127 
    128 # synonym.
    129 incomplete_gamma = gamma_inc
    130103
    131104def zeta(s):
    132105    """
  • sage/symbolic/random_tests.py

    diff --git a/sage/symbolic/random_tests.py b/sage/symbolic/random_tests.py
    a b  
    202202
    203203        sage: from sage.symbolic.random_tests import *
    204204        sage: random_expr(50, nvars=3, coeff_generator=CDF.random_element)
    205         sinh(arcsech(-heaviside(v2)/erf(-(0.615863165633 + 0.879368031485*I)*v1^2*v3) - gamma(pi) + erf(-(0.708874026302 - 0.954135400334*I)*v3)))^coth(-cosh(-polylog((v2 + 0.913564344312 + 0.0898040160336*I)^(-(0.723896589334 - 0.799038508886*I)*v2), -v1 - v3))/dilog(-(0.0263902659909 + 0.153261789843*I)*arctan2(pi, arccot(pi))))
     205        sinh(arcsech(-(0.314177274493 + 0.144437996366*I)/erf(-v1^2*e/v3) + erf(-(0.708874026302 - 0.954135400334*I)*v3) + 0.0275857401668 - 0.479027260657*I))^(-cosh(-polylog((v2^2 + (0.067987275089 + 1.08529153495*I)*v3)^(-v1 - v3), (5.29385548262 + 2.57440711353*I)*e/arctan2(arccsch(catalan), unit_step(-0.436810529675 + 0.736945423566*I)))))
    206206        sage: random_expr(5, verbose=True)
    207207        About to apply <built-in function add> to [v1, v1]
    208         About to apply <built-in function div> to [-1/3, 2*v1]
    209         -1/6/v1
     208        About to apply <built-in function mul> to [v1, -1/2]
     209        About to apply <built-in function mul> to [2*v1, -1/2*v1]
     210        -v1^2
    210211    """
    211212    vars = [(1.0, sage.calculus.calculus.var('v%d' % (n+1))) for n in range(nvars)]
    212213    if ncoeffs is None: