sympy/functions/combinatorial/factorials.py

   1
   2 from sympy.core.basic import Basic, S
   3 from sympy.core.function import SingleValuedFunction
   4 from sympy.ntheory import sieve
   5 from math import sqrt
   6
   7 ###############################################################################
   8 ######################## FACTORIAL and MULTI-FACTORIAL ########################
   9 ###############################################################################
  10
  11 class Factorial(SingleValuedFunction):
  12     """Implementation of factorial function over nonnegative integers.
  13        For the sake of convenience and simplicity of procedures using
  14        this function it is defined for negative integers and returns
  15        zero in this case.
  16
  17        The factorial is very important in combinatorics where it gives
  18        the number of ways in which 'n' objects can be permuted. It also
  19        arises in calculus, probability, number theory etc.
  20
  21        There is strict relation of factorial with gamma function. In
  22        fact n! = gamma(n+1) for nonnegarive integers. Rewrite of this
  23        kind is very useful in case of combinatorial simplification.
  24
  25        Computation of the factorial is done using two algorithms. For
  26        small arguments naive product is evaluated. However for bigger
  27        input algorithm Prime-Swing is used. It is the fastest algorithm
  28        known and computes n! via prime factorization of special class
  29        of numbers, called here the 'Swing Numbers'.
  30
  31        >>> from sympy import *
  32        >>> n = Symbol('n', integer=True)
  33
  34        >>> factorial(-2)
  35        0
  36
  37        >>> factorial(0)
  38        1
  39
  40        >>> factorial(7)
  41        5040
  42
  43        >>> factorial(n)
  44        n!
  45
  46        >>> factorial(2*n)
  47        (2*n)!
  48
  49     """
  50
  51     nofargs = 1
  52
  53     _small_swing = [
  54         1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
  55         12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
  56         35102025,5014575,145422675,9694845,300540195,300540195
  57     ]
  58
  59     @classmethod
  60     def _swing(cls, n):
  61         if n < 33:
  62             return cls._small_swing[n]
  63         else:
  64             N, primes = int(sqrt(n)), []
  65
  66             for prime in sieve.primerange(3, N+1):
  67                 p, q = 1, n
  68
  69                 while True:
  70                     q /= prime
  71
  72                     if q > 0:
  73                         if q & 1 == 1:
  74                             p *= prime
  75                     else:
  76                         break
  77
  78                 if p > 1:
  79                     primes.append(p)
  80
  81             for prime in sieve.primerange(N+1, n/3 + 1):
  82                 if (n / prime) & 1 == 1:
  83                     primes.append(prime)
  84
  85             L_product = R_product = 1
  86
  87             for prime in sieve.primerange(n/2 + 1, n+1):
  88                 L_product *= prime
  89
  90             for prime in primes:
  91                 R_product *= prime
  92
  93             return L_product*R_product
  94
  95     @classmethod
  96     def _recursive(cls, n):
  97         if n < 2:
  98             return 1
  99         else:
 100             return (cls._recursive(n/2)**2)*cls._swing(n)
 101
 102     @classmethod
 103     def canonize(cls, n):
 104         n = Basic.sympify(n)
 105
 106         if isinstance(n, Basic.Number):
 107             if isinstance(n, Basic.Zero):
 108                 return S.One
 109             elif isinstance(n, Basic.Integer):
 110                 if n.is_negative:
 111                     return S.Zero
 112                 else:
 113                     n, result = n.p, 1
 114
 115                     if n < 20:
 116                         for i in range(2, n+1):
 117                             result *= i
 118                     else:
 119                         N, bits = n, 0
 120
 121                         while N != 0:
 122                             if N & 1 == 1:
 123                                 bits += 1
 124
 125                             N = N >> 1
 126
 127                         result = cls._recursive(n)*2**(n-bits)
 128
 129                     return Basic.Integer(result)
 130
 131         if n.is_integer:
 132             if n.is_negative:
 133                 return S.Zero
 134         else:
 135             return Basic.gamma(n+1)
 136
 137
 138     @classmethod    # ?
 139     def _eval_rewrite_as_gamma(self, arg):
 140         return Basic.gamma(1 + arg)
 141
 142     def tostr(self, level=0):
 143         return '%s!' % self[0].tostr(self.precedence)
 144
 145     def _eval_is_integer(self):
 146         return self[0].is_integer
 147
 148 class MultiFactorial(SingleValuedFunction):
 149     pass
 150
 151
 152 factorial   = Factorial
 153
 154 ###############################################################################
 155 ######################## RISING and FALLING FACTORIALS ########################
 156 ###############################################################################
 157
 158 class RisingFactorial(SingleValuedFunction):
 159     """Rising factorial (also called Pochhammer symbol) is a double valued
 160        function arising in concrete mathematics, hypergeometric functions
 161        and series expanansions. It is defined by
 162
 163                    rf(x, k) = x * (x+1) * ... * (x + k-1)
 164
 165        where 'x' can be arbitrary expression and 'k' is an integer. For
 166        more information check "Concrete mathematics" by Graham, pp. 66
 167        or visit http://mathworld.wolfram.com/RisingFactorial.html page.
 168
 169        >>> from sympy import *
 170        >>> x = Symbol('x')
 171
 172        >>> rf(x, 0)
 173        1
 174
 175        >>> rf(1, 5)
 176        120
 177
 178        >>> rf(x, 5)
 179        x*(1 + x)*(2 + x)*(3 + x)*(4 + x)
 180
 181     """
 182
 183     nofargs = 2
 184
 185     @classmethod
 186     def canonize(cls, x, k):
 187         x = Basic.sympify(x)
 188         k = Basic.sympify(k)
 189
 190         if isinstance(x, Basic.NaN):
 191             return S.NaN
 192         elif isinstance(x, Basic.One):
 193             return factorial(k)
 194         elif isinstance(k, Basic.Integer):
 195             if isinstance(k, Basic.NaN):
 196                 return S.NaN
 197             elif isinstance(k, Basic.Zero):
 198                 return S.One
 199             else:
 200                 if k.is_positive:
 201                     if isinstance(x, Basic.Infinity):
 202                         return S.Infinity
 203                     elif isinstance(x, Basic.NegativeInfinity):
 204                         if k.is_odd:
 205                             return S.NegativeInfinity
 206                         else:
 207                             return S.Infinity
 208                     else:
 209                         return reduce(lambda r, i: r*(x+i), xrange(0, int(k)), 1)
 210                 else:
 211                     if isinstance(x, Basic.Infinity):
 212                         return S.Infinity
 213                     elif isinstance(x, Basic.NegativeInfinity):
 214                         return S.Infinity
 215                     else:
 216                         return 1/reduce(lambda r, i: r*(x-i), xrange(1, abs(int(k))+1), 1)
 217
 218     def _eval_rewrite_as_gamma(self, x, k):
 219         return Basic.gamma(x + k) / Basic.gamma(x)
 220
 221 class FallingFactorial(SingleValuedFunction):
 222     """Falling factorial (related to rising factorial) is a double valued
 223        function arising in concrete mathematics, hypergeometric functions
 224        and series expanansions. It is defined by
 225
 226                    ff(x, k) = x * (x-1) * ... * (x - k+1)
 227
 228        where 'x' can be arbitrary expression and 'k' is an integer. For
 229        more information check "Concrete mathematics" by Graham, pp. 66
 230        or visit http://mathworld.wolfram.com/FallingFactorial.html page.
 231
 232        >>> from sympy import *
 233        >>> x = Symbol('x')
 234
 235        >>> ff(x, 0)
 236        1
 237
 238        >>> ff(5, 5)
 239        120
 240
 241        >>> ff(x, 5)
 242        x*(1 - x)*(2 - x)*(3 - x)*(4 - x)
 243
 244     """
 245
 246     nofargs = 2
 247
 248     @classmethod
 249     def canonize(cls, x, k):
 250         x = Basic.sympify(x)
 251         k = Basic.sympify(k)
 252
 253         if isinstance(x, Basic.NaN):
 254             return S.NaN
 255         elif isinstance(k, Basic.Integer):
 256             if isinstance(k, Basic.NaN):
 257                 return S.NaN
 258             elif isinstance(k, Basic.Zero):
 259                 return S.One
 260             else:
 261                 result = S.One
 262
 263                 if k.is_positive:
 264                     if isinstance(x, Basic.Infinity):
 265                         return S.Infinity
 266                     elif isinstance(x, Basic.NegativeInfinity):
 267                         if k.is_odd:
 268                             return S.NegativeInfinity
 269                         else:
 270                             return S.Infinity
 271                     else:
 272                         return reduce(lambda r, i: r*(x-i), xrange(0, int(k)), 1)
 273                 else:
 274                     if isinstance(x, Basic.Infinity):
 275                         return S.Infinity
 276                     elif isinstance(x, Basic.NegativeInfinity):
 277                         return S.Infinity
 278                     else:
 279                         return 1/reduce(lambda r, i: r*(x+i), xrange(1, abs(int(k))+1), 1)
 280
 281
 282     def _eval_rewrite_as_gamma(self, x, k):
 283         return (-1)**k * Basic.gamma(-x + k) / Basic.gamma(-x)
 284
 285 rf = RisingFactorial
 286 ff = FallingFactorial
 287
 288 ###############################################################################
 289 ########################### BINOMIAL COEFFICIENTS #############################
 290 ###############################################################################
 291
 292 class Binomial(SingleValuedFunction):
 293     """Implementation of the binomial coefficient. It can be defined
 294        in two ways depending on its desired interpretation:
 295
 296            C(n,k) = n!/(k!(n-k)!)   or   C(n, k) = ff(n, k)/k!
 297
 298        First formula has strict combinatorial meaning, definig the
 299        number of ways we can choose 'k' elements from 'n' element
 300        set. In this case both arguments are nonnegative integers
 301        and binomial is computed using efficient algorithm based
 302        on prime factorisation.
 303
 304        The other definition is generalisation for arbitaty 'n',
 305        however 'k' must be also nonnegative. This case is very
 306        useful in case for evaluating summations.
 307
 308        For the sake of convenience for negative 'k' this function
 309        will return zero no matter what valued is the other argument.
 310
 311        >>> from sympy import *
 312        >>> n = symbols('n', integer=True)
 313
 314        >>> binomial(15, 8)
 315        6435
 316
 317        >>> binomial(n, -1)
 318        0
 319
 320        >>> [ binomial(0, i) for i in range(1)]
 321        [1]
 322        >>> [ binomial(1, i) for i in range(2)]
 323        [1, 1]
 324        >>> [ binomial(2, i) for i in range(3)]
 325        [1, 2, 1]
 326        >>> [ binomial(3, i) for i in range(4)]
 327        [1, 3, 3, 1]
 328        >>> [ binomial(4, i) for i in range(5)]
 329        [1, 4, 6, 4, 1]
 330
 331        >>> binomial(Rational(5,4), 3)
 332        -5/128
 333
 334        >>> binomial(n, 3)
 335        (1/6)*n*(1 - n)*(2 - n)
 336
 337     """
 338
 339     nofargs = 2
 340
 341     @classmethod
 342     def canonize(cls, r, k):
 343         r, k = map(Basic.sympify, (r, k))
 344
 345         if isinstance(k, Basic.Number):
 346             if isinstance(k, Basic.Zero):
 347                 return S.One
 348             elif isinstance(k, Basic.Integer):
 349                 if k.is_negative:
 350                     return S.Zero
 351                 else:
 352                     if isinstance(r, Basic.Integer) and r.is_nonnegative:
 353                         r, k = int(r), int(k)
 354
 355                         if k > r:
 356                             return S.Zero
 357                         elif k > r / 2:
 358                             k = r - k
 359
 360                         M, result = int(sqrt(r)), 1
 361
 362                         for prime in sieve.primerange(2, r+1):
 363                             if prime > r - k:
 364                                 result *= prime
 365                             elif prime > r / 2:
 366                                 continue
 367                             elif prime > M:
 368                                 if r % prime < k % prime:
 369                                     result *= prime
 370                             else:
 371                                 R, K = r, k
 372                                 exp = a = 0
 373
 374                                 while R > 0:
 375                                     a = int((R % prime) < (K % prime + a))
 376                                     R, K = R / prime, K / prime
 377                                     exp = a + exp
 378
 379                                 if exp > 0:
 380                                     result *= prime**exp
 381
 382                         return Basic.Integer(result)
 383                     else:
 384                         result = r - k + 1
 385
 386                         for i in xrange(2, k+1):
 387                             result *= r-k+i
 388                             result /= i
 389
 390                         return result
 391
 392         if k.is_integer:
 393             if k.is_negative:
 394                 return S.Zero
 395         else:
 396             return Basic.gamma(r+1)/(Basic.gamma(r-k+1)*Basic.gamma(k+1))
 397
 398
 399     def _eval_rewrite_as_gamma(self, r, k):
 400         return Basic.gamma(r+1) / (Basic.gamma(r-k+1)*Basic.gamma(k+1))
 401
 402     def _eval_is_integer(self):
 403         return self[0].is_integer and self[1].is_integer
 404
 405
 406 binomial = Binomial