Ticket #9240: trac_9240-factorial_evaluation.patch

File trac_9240-factorial_evaluation.patch, 8.4 KB (added by burcin, 10 years ago)

apply after trac_9240_full_simplify_gamma.patch

  • sage/functions/other.py

    # HG changeset patch
    # User Burcin Erocal <burcin@erocal.org>
    # Date 1306263152 -7200
    # Node ID ceef058f2d756de71375e35a8464975eba83bb7d
    # Parent  3de52b3e5a84b45ac98d6cdf48eb3d3ff9859e3e
    trac 9240: fix evaluation of factorial
    
    The evaluation of factorial() is defined by GiNaC, to get a numeric result for
    any numeric input. This definition makes sense since GiNaC limits the numeric
    evaluation function of factorial to accept positive integer arguments only.
    Pynac extends this function to other values using the gamma function, without
    changing the evaluation behavior. This led to several inconsistencies.
    
    We now override the evaluation for factorial() in Sage to return exact values
    when possible. For example factorial(1/2) -> 1/2*sqrt(pi).
    
    This patch also adds doctests to py_tgamma() and py_factorial() and makes
    py_tgamma() return numeric results in some corner cases.
    
    diff --git a/sage/functions/other.py b/sage/functions/other.py
    a b  
    44from sage.symbolic.function import GinacFunction, BuiltinFunction
    55from sage.symbolic.expression import Expression
    66from sage.symbolic.pynac import register_symbol, symbol_table
     7from sage.symbolic.pynac import py_factorial_py
    78from sage.libs.pari.gen import pari
    89from sage.symbolic.all import SR
    910from sage.rings.all import Integer, Rational, RealField, CC, RR, \
     
    733734            120
    734735            sage: gamma(1/2)
    735736            sqrt(pi)
     737            sage: gamma(-4/3)
     738            gamma(-4/3)
    736739            sage: gamma(-1)
    737740            Infinity
     741            sage: gamma(0)
     742            Infinity
    738743
    739744        ::
    740745
     
    10451050        """
    10461051        GinacFunction.__init__(self, "factorial", latex_name='{\\rm factorial}',
    10471052                conversions=dict(maxima='factorial', mathematica='Factorial'))
    1048  
     1053
     1054    def _eval_(self, x):
     1055        """
     1056        Evaluate the factorial function.
     1057
     1058        Note that this method overrides the eval method defined in GiNaC
     1059        which calls numeric evaluation on all numeric input. We preserve
     1060        exact results if the input is a rational number.
     1061
     1062        EXAMPLES::
     1063
     1064            sage: k = var('k')
     1065            sage: k.factorial()
     1066            factorial(k)
     1067            sage: SR(1/2).factorial()
     1068            1/2*sqrt(pi)
     1069            sage: SR(3/4).factorial()
     1070            gamma(7/4)
     1071            sage: SR(5).factorial()
     1072            120
     1073            sage: SR(3245908723049857203948572398475r).factorial()
     1074            factorial(3245908723049857203948572398475L)
     1075            sage: SR(3245908723049857203948572398475).factorial()
     1076            factorial(3245908723049857203948572398475)
     1077        """
     1078        if isinstance(x, Rational):
     1079            return gamma(x+1)
     1080        elif isinstance(x, (Integer, int)) or \
     1081                (not isinstance(x, Expression) and is_inexact(x)):
     1082            return py_factorial_py(x)
     1083
     1084        return None
     1085
    10491086factorial = Function_factorial()
    10501087
    10511088class Function_binomial(GinacFunction):
  • sage/rings/integer.pyx

    diff --git a/sage/rings/integer.pyx b/sage/rings/integer.pyx
    a b  
    36963696    def factorial(self):
    36973697        r"""
    36983698        Return the factorial `n! = 1 \cdot 2 \cdot 3 \cdots n`.
    3699         Self must fit in an ``unsigned long int``.
     3699
     3700        If the input does not fit in an ``unsigned long int`` a symbolic
     3701        expression is returned.
    37003702       
    37013703        EXAMPLES::
    37023704       
     
    37093711            4 24
    37103712            5 120
    37113713            6 720
     3714            sage: 234234209384023842034.factorial()
     3715            factorial(234234209384023842034)
    37123716        """
    37133717        if mpz_sgn(self.value) < 0:
    37143718            raise ValueError, "factorial -- self = (%s) must be nonnegative"%self
    37153719
    37163720        if not mpz_fits_uint_p(self.value):
    3717             raise ValueError, "factorial not implemented for n >= 2^32.\nThis is probably OK, since the answer would have billions of digits."
     3721            from sage.functions.all import factorial
     3722            return factorial(self, hold=True)
    37183723
    37193724        cdef Integer z = PY_NEW(Integer)
    37203725
  • sage/symbolic/expression.pyx

    diff --git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
    a b  
    67066706        Check that the problem with applying full_simplify() to gamma functions (Trac 9240)
    67076707        has been fixed::
    67086708
    6709             sage: gamma(37)
    6710             371993326789901217467999448150835200000000
    6711             sage: gamma(0)
    6712             Infinity
    6713             sage: gamma(-4)
    6714             Infinity
    6715             sage: gamma(-4/3)
    6716             gamma(-4/3)
     6709            sage: gamma(1/3)
     6710            gamma(1/3)
    67176711            sage: gamma(1/3).full_simplify()
    67186712            gamma(1/3)
     6713            sage: gamma(4/3)
     6714            gamma(4/3)
    67196715            sage: gamma(4/3).full_simplify()
    67206716            1/3*gamma(1/3)
    67216717
  • sage/symbolic/function.pyx

    diff --git a/sage/symbolic/function.pyx b/sage/symbolic/function.pyx
    a b  
    127127    cdef _register_function(self):
    128128        cdef GFunctionOpt opt
    129129        opt = g_function_options_args(self._name, self._nargs)
    130         opt.set_python_func()
    131130
    132131        if hasattr(self, '_eval_'):
    133132            opt.eval_func(self)
  • sage/symbolic/pynac.pyx

    diff --git a/sage/symbolic/pynac.pyx b/sage/symbolic/pynac.pyx
    a b  
    10301030# SPECIAL FUNCTIONS
    10311031#################################################################
    10321032cdef public object py_tgamma(object x) except +:
     1033    """
     1034    The gamma function exported to pynac.
     1035
     1036    TESTS::
     1037
     1038        sage: from sage.symbolic.pynac import py_tgamma_for_doctests as py_tgamma
     1039        sage: py_tgamma(4)
     1040        6
     1041        sage: py_tgamma(1/2)
     1042        1.77245385090552
     1043    """
     1044    if PY_TYPE_CHECK_EXACT(x, int) or PY_TYPE_CHECK_EXACT(x, long):
     1045        x = float(x)
    10331046    if PY_TYPE_CHECK_EXACT(x, float):
    10341047        return sage_tgammal(x)
    1035     #FIXME: complex
     1048
     1049    # try / except blocks are faster than
     1050    # if hasattr(x, 'gamma')
    10361051    try:
    1037         return x.gamma()
     1052        res = x.gamma()
    10381053    except AttributeError:
    10391054        return CC(x).gamma()
    10401055
     1056    # the result should be numeric, however the gamma method of rationals may
     1057    # return symbolic expressions. for example (1/2).gamma() -> sqrt(pi).
     1058    if isinstance(res, Expression):
     1059        try:
     1060            return RR(res)
     1061        except ValueError:
     1062            return CC(res)
     1063    return res
     1064
     1065def py_tgamma_for_doctests(x):
     1066    """
     1067    This function is for testing py_tgamma().
     1068
     1069    TESTS::
     1070
     1071        sage: from sage.symbolic.pynac import py_tgamma_for_doctests
     1072        sage: py_tgamma_for_doctests(3)
     1073        2
     1074    """
     1075    return py_tgamma(x)
     1076
    10411077from sage.rings.arith import factorial
    10421078cdef public object py_factorial(object x) except +:
     1079    """
     1080    The factorial function exported to pynac.
     1081
     1082    TESTS::
     1083
     1084        sage: from sage.symbolic.pynac import py_factorial_py as py_factorial
     1085        sage: py_factorial(4)
     1086        24
     1087        sage: py_factorial(-2/3)
     1088        2.67893853470775
     1089    """
    10431090    # factorial(x) is only defined for non-negative integers x
    10441091    # so we first test if x can be coerced into ZZ and is non-negative.
    10451092    # If this is not the case then we return the symbolic expression gamma(x+1)
     
    10551102    else:
    10561103        return py_tgamma(x+1)
    10571104
     1105def py_factorial_py(x):
     1106    """
     1107    This function is a python wrapper around py_factorial(). This wrapper
     1108    is needed when we override the eval() method for GiNaC's factorial
     1109    function in sage.functions.other.Function_factorial.
     1110
     1111    TESTS::
     1112
     1113        sage: from sage.symbolic.pynac import py_factorial_py
     1114        sage: py_factorial_py(3)
     1115        6
     1116    """
     1117    return py_factorial(x)
     1118
    10581119cdef public object py_doublefactorial(object x) except +:
    10591120    n = Integer(x)
    10601121    if n < -1:
     
    14671528        gsl_sf_lngamma_complex_e(PyComplex_RealAsDouble(x),PyComplex_ImagAsDouble(x), &lnr, &arg)
    14681529        res = gsl_complex_polar(lnr.val, arg.val)
    14691530        return PyComplex_FromDoubles(res.dat[0], res.dat[1])
    1470     elif hasattr(x, 'log_gamma'):
    1471          return x.log_gamma()
    14721531    elif isinstance(x, Integer):
    14731532        return x.gamma().log().n()
    1474     elif hasattr(x, 'gamma'):
     1533
     1534    # try / except blocks are faster than
     1535    # if hasattr(x, 'log_gamma')
     1536    try:
     1537        return x.log_gamma()
     1538    except AttributeError:
     1539        pass
     1540
     1541    try:
    14751542        return x.gamma().log()
     1543    except AttributeError:
     1544        pass
     1545
    14761546    return CC(x).gamma().log()
    14771547
    14781548def py_lgamma_for_doctests(x):