sympy/specfun/factorials.py

   1 from sympy.core import *
   2 from sympy.core.basic import S
   3
   4 from sympy.utilities.memoization import recurrence_memo
   5
   6 # Lanczos approximation for low-precision numerical factorial
   7 _lanczos_coef = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
   8   771.32342877765313, -176.61502916214059, 12.507343278686905,
   9     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
  10
  11 def _lanczos(z):
  12     from cmath import pi, sin, log, exp
  13     if z.real < 0:
  14         return pi*z / (sin(pi*z) * _lanczos(-z))
  15     else:
  16         x = _lanczos_coef[0]
  17         for i in range(1, 9):
  18             x += _lanczos_coef[i]/(z+i)
  19         logw = 0.91893853320467267+(z+0.5)*log(z+7.5)+log(x)-z-7.5
  20         return exp(logw)
  21
  22 class _Factorial(SingleValuedFunction):
  23     """
  24     Factorials and multiple factorials
  25
  26     Usage
  27     =====
  28         factorial(x) gives the factorial of x, defined as
  29         x! = 1*2*3*...*x if x is a positive integer. If x is not a
  30         positive integer, the value of x! is defined in terms of the
  31         gamma function.
  32
  33         factorial(x, m) returns the m-th order multifactorial of x;
  34         i.e., 1 * ... * (x-2*m) * (x-m) * x. In particular,
  35         factorial(x, 2) returns the double factorial of x,
  36
  37           * x!! = 2*4*...*(x-2)*x if x is a positive even integer
  38           * x!! = 1*3*...*(x-2)*x if x is a positive odd integer
  39
  40         The argument m must be a positive integer.
  41
  42     Examples
  43     ========
  44         >>> factorial(5)
  45         120
  46         >>> factorial(0)
  47         1
  48         >>> factorial(Rational(5,2))
  49         (15/8)*pi**(1/2)
  50
  51         >>> factorial(5, 2)
  52         15
  53         >>> factorial(6, 2)
  54         48
  55
  56     """
  57
  58     @staticmethod
  59     @recurrence_memo([1])
  60     def _fac1(n, prev):
  61         return n * prev[-1]
  62
  63     # generators for multifactorials
  64     _generators = {}
  65
  66     @classmethod
  67     def canonize(cls, x, m=1):
  68
  69         # the usual case
  70         if m == 1 and x.is_integer:
  71             # handle factorial poles, also for symbols with assumptions
  72             if x.is_negative:
  73                 return oo
  74             if isinstance(x, Integer):
  75                 return Integer(cls._fac1(int(x)))
  76
  77         # half-integer case of the ordinary factorial
  78         if m == 1 and isinstance(x, Rational) and x.q == 2:
  79             n = (x.p + 1) / 2
  80             if n < 0:
  81                 return (-1)**(-n+1) * pi * x / factorial(-x)
  82             return Basic.sqrt(pi) * Rational(1, 2**n) * factorial(2*n-1, 2)
  83
  84         # multifactorials are only defined for integers
  85         if (not isinstance(m, Integer) and m > 0) or not \
  86             isinstance(x, Integer):
  87             return
  88
  89         x = int(x)
  90         m = int(m)
  91
  92         if x == 0:
  93             return Integer(1)
  94
  95         # e.g., for double factorials, the start value is 1 for odd x
  96         # and 2 for even x
  97         start_value = x % m
  98         if start_value == 0:
  99             start_value = m
 100
 101         # we can define the m-multifactorials for negative numbers
 102         # to satisfy f(x,m) = f(x+m,m) / (x+m)
 103         if x < 0:
 104             if start_value == m:
 105                 return oo
 106             return factorial(x+m, m) / (x+m)
 107
 108         # positive case
 109         if not (m, start_value) in cls._generators:
 110             @recurrence_memo([start_value])
 111             def _f(k, prev):
 112                 return (k*m + start_value) * prev[-1]
 113             cls._generators[m, start_value] = _f
 114
 115         return Integer(cls._generators[m, start_value]((x-start_value) // m))
 116
 117     # This should give a series expansion around x = oo. Needs fixing
 118     # def series(self, x, n):
 119     #    return sqrt(2*pi*x) * x**x * exp(-x) * (1 + O(1/x))
 120
 121     # Derivatives are given in terms of polygamma functions
 122     # (XXX: only for order m = 1)
 123     def fdiff(self, argindex=1):
 124         if argindex == 1:
 125             from sympy.functions import gamma, polygamma
 126             return gamma(self[0]+1)*polygamma(0,self[0]+1)
 127         else:
 128             raise ArgumentIndexError(self, argindex)
 129
 130     @classmethod
 131     def _eval_apply_evalf(cls, x, m=1):
 132         """Return a low-precision numerical approximation."""
 133         assert m == 1
 134         a, b = x.as_real_imag()
 135         y = _lanczos(complex(a, b))
 136         if b == 0:
 137             return Real(y.real)
 138         else:
 139             Real(y.real) + I*Real(y.imag)
 140
 141 factorial = _Factorial
 142
 143 class UnevaluatedFactorial(_Factorial):
 144
 145     @classmethod
 146     def canonize(cls, x, m=1):
 147         return None
 148
 149 unfac = UnevaluatedFactorial
 150
 151
 152
 153 # factorial_simplify helpers; could use refactoring
 154
 155 def _isfactorial(expr):
 156     #return isinstance(expr, Apply) and isinstance(expr[0], Factorial)
 157     return isinstance(expr, _Factorial)
 158
 159 def _collect_factors(expr):
 160     assert isinstance(expr, Mul)
 161     numer_args = []
 162     denom_args = []
 163     other = []
 164     for x in expr:
 165         if isinstance(x, Mul):
 166             n, d, o = _collect_factors(x)
 167             numer_args += n
 168             denom_args += d
 169             other += o
 170         elif isinstance(x, Pow):
 171             base, exp = x[:]
 172             if _isfactorial(base) and \
 173                 isinstance(exp, Rational) and exp.is_integer:
 174                 if exp > 0:
 175                     for i in xrange(exp.p): numer_args.append(base[0])
 176                 else:
 177                     for i in xrange(-exp.p): denom_args.append(base[0])
 178             else:
 179                 other.append(x)
 180         elif _isfactorial(x):
 181             numer_args.append(x[0])
 182         else:
 183             other.append(x)
 184     return numer_args, denom_args, other
 185
 186 # handle x!/(x+n)!
 187 def _simplify_quotient(na, da, other):
 188     while 1:
 189         candidates = []
 190         for i, y in enumerate(na):
 191             for j, x in enumerate(da):
 192                 #delta = simplify(y - x)
 193                 delta = y - x
 194                 if isinstance(delta, Rational) and delta.is_integer:
 195                     candidates.append((delta, i, j))
 196         if candidates:
 197             # There may be multiple candidates. Choose the quotient pair
 198             # that minimizes the work.
 199             candidates.sort(key=lambda x: abs(x[0]))
 200             delta, i, j = candidates[0]
 201             p = na[i]
 202             q = da[j]
 203             t = Rational(1)
 204             if delta > 0:
 205                 for k in xrange(1, int(delta)+1):
 206                     t *= (q+k)
 207             else:
 208                 for k in xrange(1, -int(delta)+1):
 209                     t /= (p+k)
 210             other.append(t)
 211             del na[i], da[j]
 212         else:
 213             return
 214
 215 # handle x!*(x+1) and x!/x
 216 def _simplify_recurrence(facs, other, reciprocal=False):
 217     # this should probably be rewritten more elegantly
 218     i = 0
 219     while i < len(facs):
 220         j = 0
 221         while j < len(other):
 222             othr = other[j]
 223             fact = facs[i]
 224             if reciprocal:
 225                 othr = 1/othr
 226             if   othr - fact == 1: facs[i] += 1; del other[j]; j -= 1
 227             elif -othr - fact == 1: facs[i] += 1; del other[j]; other.append(-1); j -= 1
 228             elif 1/othr - fact == 0: facs[i] -= 1; del other[j]; j -= 1
 229             elif -1/othr - fact == 0: facs[i] -= 1; del other[j]; other.append(-1); j -= 1
 230             j += 1
 231         i += 1
 232
 233 def factorial_simplify(expr):
 234     """
 235     This function takes an expression containing factorials
 236     and removes as many of them as possible by combining
 237     products and quotients of factorials into polynomials
 238     or other simpler expressions.
 239
 240     TODO: handle reflection formula, duplication formula
 241     double factorials
 242     """
 243
 244
 245     if isinstance(expr, Add):
 246         return Add(*(factorial_simplify(x) for x in expr))
 247
 248     if isinstance(expr, Pow):
 249         return Pow(factorial_simplify(expr[0]), expr[1])
 250
 251     if isinstance(expr, Mul):
 252         na, da, other = _collect_factors(expr)
 253         _simplify_quotient(na, da, other)
 254         _simplify_recurrence(na, other)
 255         _simplify_recurrence(da, other, reciprocal=True)
 256
 257         result = Rational(1)
 258         for n in na: result *= factorial(n)
 259         for d in da: result /= factorial(d)
 260         for o in other: result *= o
 261         return result
 262
 263     expr = expr.subs(unfac, factorial)
 264
 265     return expr
 266
 267 class Rising_factorial(SingleValuedFunction):
 268     """
 269     Usage
 270     =====
 271         Calculate the rising factorial (x)^(n) = x(x+1)...(x+n-1)
 272         as a quotient of factorials
 273
 274     Examples
 275     ========
 276         >>> rising_factorial(3, 2)
 277         12
 278
 279     """
 280     nofargs = 2
 281
 282     @classmethod
 283     def canonize(cls, x, n):
 284         return factorial_simplify(unfac(x+n-1) / unfac(x-1))
 285
 286
 287 rising_factorial = Rising_factorial
 288
 289
 290 class Falling_factorial(SingleValuedFunction):
 291     """
 292     Usage
 293     =====
 294         Calculate the falling factorial (x)_(n) = x(x-1)...(x-n+1)
 295         as a quotient of factorials
 296
 297     Examples
 298     ========
 299         >>> falling_factorial(5, 3)
 300         60
 301
 302     """
 303     nofargs = 2
 304
 305     @classmethod
 306     def canonize(cls, x, n):
 307         return factorial_simplify(unfac(x) / unfac(x-n))
 308
 309
 310 falling_factorial = Falling_factorial
 311
 312
 313 class Binomial2(SingleValuedFunction):
 314     """
 315     Usage
 316     =====
 317         Calculate the binomial coefficient C(n,k) = n!/(k!(n-k)!)
 318
 319     Notes
 320     =====
 321         When n and k are positive integers, the result is always
 322         a positive integer
 323
 324         If k is a positive integer, the result is a polynomial in n
 325         that is evaluated explicitly.
 326
 327     Examples
 328     ========
 329         >>> binomial2(15,8)
 330         6435
 331         >>> # Building Pascal's triangle
 332         >>> [binomial2(0,k) for k in range(1)]
 333         [1]
 334         >>> [binomial2(1,k) for k in range(2)]
 335         [1, 1]
 336         >>> [binomial2(2,k) for k in range(3)]
 337         [1, 2, 1]
 338         >>> [binomial2(3,k) for k in range(4)]
 339         [1, 3, 3, 1]
 340         >>> # n can be arbitrary if k is a positive integer
 341         >>> binomial2(Rational(5,4), 3)
 342         -5/128
 343         >>> x = Symbol('x')
 344         >>> binomial2(x, 3)
 345         (1/6)*x*(1 - x)*(2 - x)
 346
 347     """
 348     nofargs = 2
 349
 350     @classmethod
 351     def canonize(cls, n, k):
 352
 353         # TODO: move these two cases to factorial_simplify as well
 354         if n == 0 and k != 0:
 355             return Basic.sin(pi*k)/(pi*k)
 356
 357         return factorial_simplify(unfac(n) / unfac(k) / unfac(n-k))
 358
 359 binomial2 = Binomial2