Ticket #9240: trac_9240-factorial_evaluation.patch

sage/functions/other.py

@@ -4 +4 @@
 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, \
@@ -733 +734 @@
             120
             sage: gamma(1/2)
             sqrt(pi)
+            sage: gamma(-4/3)
+            gamma(-4/3)
             sage: gamma(-1)
             Infinity
+            sage: gamma(0)
+            Infinity
 
         ::
 
@@ -1045 +1050 @@
         """
         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

@@ -3696 +3696 @@
     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::
         
@@ -3709 +3711 @@
             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

@@ -6706 +6706 @@
         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

@@ -127 +127 @@
     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

@@ -1030 +1030 @@
 # 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)
@@ -1055 +1102 @@
     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:
@@ -1467 +1528 @@
         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):