# Ticket #9240: trac_9240-factorial_evaluation.patch

File trac_9240-factorial_evaluation.patch, 8.4 KB (added by burcin, 2 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 from sage.symbolic.function import GinacFunction, BuiltinFunction from sage.symbolic.expression import Expression from sage.symbolic.pynac import register_symbol, symbol_table from sage.symbolic.pynac import py_factorial_py from sage.libs.pari.gen import pari from sage.symbolic.all import SR from sage.rings.all import Integer, Rational, RealField, CC, RR, \ 120 sage: gamma(1/2) sqrt(pi) sage: gamma(-4/3) gamma(-4/3) sage: gamma(-1) Infinity sage: gamma(0) Infinity :: """ GinacFunction.__init__(self, "factorial", latex_name='{\\rm factorial}', conversions=dict(maxima='factorial', mathematica='Factorial')) def _eval_(self, x): """ Evaluate the factorial function. Note that this method overrides the eval method defined in GiNaC which calls numeric evaluation on all numeric input. We preserve exact results if the input is a rational number. EXAMPLES:: sage: k = var('k') sage: k.factorial() factorial(k) sage: SR(1/2).factorial() 1/2*sqrt(pi) sage: SR(3/4).factorial() gamma(7/4) sage: SR(5).factorial() 120 sage: SR(3245908723049857203948572398475r).factorial() factorial(3245908723049857203948572398475L) sage: SR(3245908723049857203948572398475).factorial() factorial(3245908723049857203948572398475) """ if isinstance(x, Rational): return gamma(x+1) elif isinstance(x, (Integer, int)) or \ (not isinstance(x, Expression) and is_inexact(x)): return py_factorial_py(x) return None factorial = Function_factorial() class Function_binomial(GinacFunction):
• ## sage/rings/integer.pyx

diff --git a/sage/rings/integer.pyx b/sage/rings/integer.pyx
 a def factorial(self): r""" Return the factorial n! = 1 \cdot 2 \cdot 3 \cdots n. Self must fit in an unsigned long int. If the input does not fit in an unsigned long int a symbolic expression is returned. EXAMPLES:: 4 24 5 120 6 720 sage: 234234209384023842034.factorial() factorial(234234209384023842034) """ if mpz_sgn(self.value) < 0: raise ValueError, "factorial -- self = (%s) must be nonnegative"%self if not mpz_fits_uint_p(self.value): raise ValueError, "factorial not implemented for n >= 2^32.\nThis is probably OK, since the answer would have billions of digits." from sage.functions.all import factorial return factorial(self, hold=True) cdef Integer z = PY_NEW(Integer)
• ## sage/symbolic/expression.pyx

diff --git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
 a Check that the problem with applying full_simplify() to gamma functions (Trac 9240) has been fixed:: sage: gamma(37) 371993326789901217467999448150835200000000 sage: gamma(0) Infinity sage: gamma(-4) Infinity sage: gamma(-4/3) gamma(-4/3) sage: gamma(1/3) gamma(1/3) sage: gamma(1/3).full_simplify() gamma(1/3) sage: gamma(4/3) gamma(4/3) sage: gamma(4/3).full_simplify() 1/3*gamma(1/3)
• ## sage/symbolic/function.pyx

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

diff --git a/sage/symbolic/pynac.pyx b/sage/symbolic/pynac.pyx
 a # SPECIAL FUNCTIONS ################################################################# cdef public object py_tgamma(object x) except +: """ The gamma function exported to pynac. TESTS:: sage: from sage.symbolic.pynac import py_tgamma_for_doctests as py_tgamma sage: py_tgamma(4) 6 sage: py_tgamma(1/2) 1.77245385090552 """ if PY_TYPE_CHECK_EXACT(x, int) or PY_TYPE_CHECK_EXACT(x, long): x = float(x) if PY_TYPE_CHECK_EXACT(x, float): return sage_tgammal(x) #FIXME: complex # try / except blocks are faster than # if hasattr(x, 'gamma') try: return x.gamma() res = x.gamma() except AttributeError: return CC(x).gamma() # the result should be numeric, however the gamma method of rationals may # return symbolic expressions. for example (1/2).gamma() -> sqrt(pi). if isinstance(res, Expression): try: return RR(res) except ValueError: return CC(res) return res def py_tgamma_for_doctests(x): """ This function is for testing py_tgamma(). TESTS:: sage: from sage.symbolic.pynac import py_tgamma_for_doctests sage: py_tgamma_for_doctests(3) 2 """ return py_tgamma(x) from sage.rings.arith import factorial cdef public object py_factorial(object x) except +: """ The factorial function exported to pynac. TESTS:: sage: from sage.symbolic.pynac import py_factorial_py as py_factorial sage: py_factorial(4) 24 sage: py_factorial(-2/3) 2.67893853470775 """ # factorial(x) is only defined for non-negative integers x # so we first test if x can be coerced into ZZ and is non-negative. # If this is not the case then we return the symbolic expression gamma(x+1) else: return py_tgamma(x+1) def py_factorial_py(x): """ This function is a python wrapper around py_factorial(). This wrapper is needed when we override the eval() method for GiNaC's factorial function in sage.functions.other.Function_factorial. TESTS:: sage: from sage.symbolic.pynac import py_factorial_py sage: py_factorial_py(3) 6 """ return py_factorial(x) cdef public object py_doublefactorial(object x) except +: n = Integer(x) if n < -1: gsl_sf_lngamma_complex_e(PyComplex_RealAsDouble(x),PyComplex_ImagAsDouble(x), &lnr, &arg) res = gsl_complex_polar(lnr.val, arg.val) return PyComplex_FromDoubles(res.dat[0], res.dat[1]) elif hasattr(x, 'log_gamma'): return x.log_gamma() elif isinstance(x, Integer): return x.gamma().log().n() elif hasattr(x, 'gamma'): # try / except blocks are faster than # if hasattr(x, 'log_gamma') try: return x.log_gamma() except AttributeError: pass try: return x.gamma().log() except AttributeError: pass return CC(x).gamma().log() def py_lgamma_for_doctests(x):