1from typing import List 2from functools import reduce 3 4from sympy.core import S, sympify, Dummy, Mod 5from sympy.core.cache import cacheit 6from sympy.core.compatibility import HAS_GMPY 7from sympy.core.function import Function, ArgumentIndexError, PoleError 8from sympy.core.logic import fuzzy_and 9from sympy.core.numbers import Integer, pi 10from sympy.core.relational import Eq 11from sympy.ntheory import sieve 12from sympy.polys.polytools import Poly 13 14from math import sqrt as _sqrt 15 16 17class CombinatorialFunction(Function): 18 """Base class for combinatorial functions. """ 19 20 def _eval_simplify(self, **kwargs): 21 from sympy.simplify.combsimp import combsimp 22 # combinatorial function with non-integer arguments is 23 # automatically passed to gammasimp 24 expr = combsimp(self) 25 measure = kwargs['measure'] 26 if measure(expr) <= kwargs['ratio']*measure(self): 27 return expr 28 return self 29 30 31############################################################################### 32######################## FACTORIAL and MULTI-FACTORIAL ######################## 33############################################################################### 34 35 36class factorial(CombinatorialFunction): 37 r"""Implementation of factorial function over nonnegative integers. 38 By convention (consistent with the gamma function and the binomial 39 coefficients), factorial of a negative integer is complex infinity. 40 41 The factorial is very important in combinatorics where it gives 42 the number of ways in which `n` objects can be permuted. It also 43 arises in calculus, probability, number theory, etc. 44 45 There is strict relation of factorial with gamma function. In 46 fact `n! = gamma(n+1)` for nonnegative integers. Rewrite of this 47 kind is very useful in case of combinatorial simplification. 48 49 Computation of the factorial is done using two algorithms. For 50 small arguments a precomputed look up table is used. However for bigger 51 input algorithm Prime-Swing is used. It is the fastest algorithm 52 known and computes `n!` via prime factorization of special class 53 of numbers, called here the 'Swing Numbers'. 54 55 Examples 56 ======== 57 58 >>> from sympy import Symbol, factorial, S 59 >>> n = Symbol('n', integer=True) 60 61 >>> factorial(0) 62 1 63 64 >>> factorial(7) 65 5040 66 67 >>> factorial(-2) 68 zoo 69 70 >>> factorial(n) 71 factorial(n) 72 73 >>> factorial(2*n) 74 factorial(2*n) 75 76 >>> factorial(S(1)/2) 77 factorial(1/2) 78 79 See Also 80 ======== 81 82 factorial2, RisingFactorial, FallingFactorial 83 """ 84 85 def fdiff(self, argindex=1): 86 from sympy import gamma, polygamma 87 if argindex == 1: 88 return gamma(self.args[0] + 1)*polygamma(0, self.args[0] + 1) 89 else: 90 raise ArgumentIndexError(self, argindex) 91 92 _small_swing = [ 93 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429, 6435, 6435, 109395, 94 12155, 230945, 46189, 969969, 88179, 2028117, 676039, 16900975, 1300075, 95 35102025, 5014575, 145422675, 9694845, 300540195, 300540195 96 ] 97 98 _small_factorials = [] # type: List[int] 99 100 @classmethod 101 def _swing(cls, n): 102 if n < 33: 103 return cls._small_swing[n] 104 else: 105 N, primes = int(_sqrt(n)), [] 106 107 for prime in sieve.primerange(3, N + 1): 108 p, q = 1, n 109 110 while True: 111 q //= prime 112 113 if q > 0: 114 if q & 1 == 1: 115 p *= prime 116 else: 117 break 118 119 if p > 1: 120 primes.append(p) 121 122 for prime in sieve.primerange(N + 1, n//3 + 1): 123 if (n // prime) & 1 == 1: 124 primes.append(prime) 125 126 L_product = R_product = 1 127 128 for prime in sieve.primerange(n//2 + 1, n + 1): 129 L_product *= prime 130 131 for prime in primes: 132 R_product *= prime 133 134 return L_product*R_product 135 136 @classmethod 137 def _recursive(cls, n): 138 if n < 2: 139 return 1 140 else: 141 return (cls._recursive(n//2)**2)*cls._swing(n) 142 143 @classmethod 144 def eval(cls, n): 145 n = sympify(n) 146 147 if n.is_Number: 148 if n.is_zero: 149 return S.One 150 elif n is S.Infinity: 151 return S.Infinity 152 elif n.is_Integer: 153 if n.is_negative: 154 return S.ComplexInfinity 155 else: 156 n = n.p 157 158 if n < 20: 159 if not cls._small_factorials: 160 result = 1 161 for i in range(1, 20): 162 result *= i 163 cls._small_factorials.append(result) 164 result = cls._small_factorials[n-1] 165 166 # GMPY factorial is faster, use it when available 167 elif HAS_GMPY: 168 from sympy.core.compatibility import gmpy 169 result = gmpy.fac(n) 170 171 else: 172 bits = bin(n).count('1') 173 result = cls._recursive(n)*2**(n - bits) 174 175 return Integer(result) 176 177 def _facmod(self, n, q): 178 res, N = 1, int(_sqrt(n)) 179 180 # Exponent of prime p in n! is e_p(n) = [n/p] + [n/p**2] + ... 181 # for p > sqrt(n), e_p(n) < sqrt(n), the primes with [n/p] = m, 182 # occur consecutively and are grouped together in pw[m] for 183 # simultaneous exponentiation at a later stage 184 pw = [1]*N 185 186 m = 2 # to initialize the if condition below 187 for prime in sieve.primerange(2, n + 1): 188 if m > 1: 189 m, y = 0, n // prime 190 while y: 191 m += y 192 y //= prime 193 if m < N: 194 pw[m] = pw[m]*prime % q 195 else: 196 res = res*pow(prime, m, q) % q 197 198 for ex, bs in enumerate(pw): 199 if ex == 0 or bs == 1: 200 continue 201 if bs == 0: 202 return 0 203 res = res*pow(bs, ex, q) % q 204 205 return res 206 207 def _eval_Mod(self, q): 208 n = self.args[0] 209 if n.is_integer and n.is_nonnegative and q.is_integer: 210 aq = abs(q) 211 d = aq - n 212 if d.is_nonpositive: 213 return S.Zero 214 else: 215 isprime = aq.is_prime 216 if d == 1: 217 # Apply Wilson's theorem (if a natural number n > 1 218 # is a prime number, then (n-1)! = -1 mod n) and 219 # its inverse (if n > 4 is a composite number, then 220 # (n-1)! = 0 mod n) 221 if isprime: 222 return S(-1 % q) 223 elif isprime is False and (aq - 6).is_nonnegative: 224 return S.Zero 225 elif n.is_Integer and q.is_Integer: 226 n, d, aq = map(int, (n, d, aq)) 227 if isprime and (d - 1 < n): 228 fc = self._facmod(d - 1, aq) 229 fc = pow(fc, aq - 2, aq) 230 if d%2: 231 fc = -fc 232 else: 233 fc = self._facmod(n, aq) 234 235 return S(fc % q) 236 237 def _eval_rewrite_as_gamma(self, n, piecewise=True, **kwargs): 238 from sympy import gamma 239 return gamma(n + 1) 240 241 def _eval_rewrite_as_Product(self, n, **kwargs): 242 from sympy import Product 243 if n.is_nonnegative and n.is_integer: 244 i = Dummy('i', integer=True) 245 return Product(i, (i, 1, n)) 246 247 def _eval_is_integer(self): 248 if self.args[0].is_integer and self.args[0].is_nonnegative: 249 return True 250 251 def _eval_is_positive(self): 252 if self.args[0].is_integer and self.args[0].is_nonnegative: 253 return True 254 255 def _eval_is_even(self): 256 x = self.args[0] 257 if x.is_integer and x.is_nonnegative: 258 return (x - 2).is_nonnegative 259 260 def _eval_is_composite(self): 261 x = self.args[0] 262 if x.is_integer and x.is_nonnegative: 263 return (x - 3).is_nonnegative 264 265 def _eval_is_real(self): 266 x = self.args[0] 267 if x.is_nonnegative or x.is_noninteger: 268 return True 269 270 def _eval_as_leading_term(self, x, logx=None, cdir=0): 271 arg = self.args[0].as_leading_term(x) 272 arg0 = arg.subs(x, 0) 273 if arg0.is_zero: 274 return S.One 275 elif not arg0.is_infinite: 276 return self.func(arg) 277 raise PoleError("Cannot expand %s around 0" % (self)) 278 279class MultiFactorial(CombinatorialFunction): 280 pass 281 282 283class subfactorial(CombinatorialFunction): 284 r"""The subfactorial counts the derangements of n items and is 285 defined for non-negative integers as: 286 287 .. math:: !n = \begin{cases} 1 & n = 0 \\ 0 & n = 1 \\ 288 (n-1)(!(n-1) + !(n-2)) & n > 1 \end{cases} 289 290 It can also be written as ``int(round(n!/exp(1)))`` but the 291 recursive definition with caching is implemented for this function. 292 293 An interesting analytic expression is the following [2]_ 294 295 .. math:: !x = \Gamma(x + 1, -1)/e 296 297 which is valid for non-negative integers `x`. The above formula 298 is not very useful incase of non-integers. :math:`\Gamma(x + 1, -1)` is 299 single-valued only for integral arguments `x`, elsewhere on the positive 300 real axis it has an infinite number of branches none of which are real. 301 302 References 303 ========== 304 305 .. [1] https://en.wikipedia.org/wiki/Subfactorial 306 .. [2] http://mathworld.wolfram.com/Subfactorial.html 307 308 Examples 309 ======== 310 311 >>> from sympy import subfactorial 312 >>> from sympy.abc import n 313 >>> subfactorial(n + 1) 314 subfactorial(n + 1) 315 >>> subfactorial(5) 316 44 317 318 See Also 319 ======== 320 321 sympy.functions.combinatorial.factorials.factorial, 322 sympy.utilities.iterables.generate_derangements, 323 sympy.functions.special.gamma_functions.uppergamma 324 """ 325 326 @classmethod 327 @cacheit 328 def _eval(self, n): 329 if not n: 330 return S.One 331 elif n == 1: 332 return S.Zero 333 else: 334 z1, z2 = 1, 0 335 for i in range(2, n + 1): 336 z1, z2 = z2, (i - 1)*(z2 + z1) 337 return z2 338 339 @classmethod 340 def eval(cls, arg): 341 if arg.is_Number: 342 if arg.is_Integer and arg.is_nonnegative: 343 return cls._eval(arg) 344 elif arg is S.NaN: 345 return S.NaN 346 elif arg is S.Infinity: 347 return S.Infinity 348 349 def _eval_is_even(self): 350 if self.args[0].is_odd and self.args[0].is_nonnegative: 351 return True 352 353 def _eval_is_integer(self): 354 if self.args[0].is_integer and self.args[0].is_nonnegative: 355 return True 356 357 def _eval_rewrite_as_factorial(self, arg, **kwargs): 358 from sympy import summation 359 i = Dummy('i') 360 f = S.NegativeOne**i / factorial(i) 361 return factorial(arg) * summation(f, (i, 0, arg)) 362 363 def _eval_rewrite_as_gamma(self, arg, piecewise=True, **kwargs): 364 from sympy import exp, gamma, I, lowergamma 365 return ((-1)**(arg + 1)*exp(-I*pi*arg)*lowergamma(arg + 1, -1) + gamma(arg + 1))*exp(-1) 366 367 def _eval_rewrite_as_uppergamma(self, arg, **kwargs): 368 from sympy import uppergamma 369 return uppergamma(arg + 1, -1)/S.Exp1 370 371 def _eval_is_nonnegative(self): 372 if self.args[0].is_integer and self.args[0].is_nonnegative: 373 return True 374 375 def _eval_is_odd(self): 376 if self.args[0].is_even and self.args[0].is_nonnegative: 377 return True 378 379 380class factorial2(CombinatorialFunction): 381 r"""The double factorial `n!!`, not to be confused with `(n!)!` 382 383 The double factorial is defined for nonnegative integers and for odd 384 negative integers as: 385 386 .. math:: n!! = \begin{cases} 1 & n = 0 \\ 387 n(n-2)(n-4) \cdots 1 & n\ \text{positive odd} \\ 388 n(n-2)(n-4) \cdots 2 & n\ \text{positive even} \\ 389 (n+2)!!/(n+2) & n\ \text{negative odd} \end{cases} 390 391 References 392 ========== 393 394 .. [1] https://en.wikipedia.org/wiki/Double_factorial 395 396 Examples 397 ======== 398 399 >>> from sympy import factorial2, var 400 >>> n = var('n') 401 >>> n 402 n 403 >>> factorial2(n + 1) 404 factorial2(n + 1) 405 >>> factorial2(5) 406 15 407 >>> factorial2(-1) 408 1 409 >>> factorial2(-5) 410 1/3 411 412 See Also 413 ======== 414 415 factorial, RisingFactorial, FallingFactorial 416 """ 417 418 @classmethod 419 def eval(cls, arg): 420 # TODO: extend this to complex numbers? 421 422 if arg.is_Number: 423 if not arg.is_Integer: 424 raise ValueError("argument must be nonnegative integer " 425 "or negative odd integer") 426 427 # This implementation is faster than the recursive one 428 # It also avoids "maximum recursion depth exceeded" runtime error 429 if arg.is_nonnegative: 430 if arg.is_even: 431 k = arg / 2 432 return 2**k * factorial(k) 433 return factorial(arg) / factorial2(arg - 1) 434 435 436 if arg.is_odd: 437 return arg*(S.NegativeOne)**((1 - arg)/2) / factorial2(-arg) 438 raise ValueError("argument must be nonnegative integer " 439 "or negative odd integer") 440 441 442 def _eval_is_even(self): 443 # Double factorial is even for every positive even input 444 n = self.args[0] 445 if n.is_integer: 446 if n.is_odd: 447 return False 448 if n.is_even: 449 if n.is_positive: 450 return True 451 if n.is_zero: 452 return False 453 454 def _eval_is_integer(self): 455 # Double factorial is an integer for every nonnegative input, and for 456 # -1 and -3 457 n = self.args[0] 458 if n.is_integer: 459 if (n + 1).is_nonnegative: 460 return True 461 if n.is_odd: 462 return (n + 3).is_nonnegative 463 464 def _eval_is_odd(self): 465 # Double factorial is odd for every odd input not smaller than -3, and 466 # for 0 467 n = self.args[0] 468 if n.is_odd: 469 return (n + 3).is_nonnegative 470 if n.is_even: 471 if n.is_positive: 472 return False 473 if n.is_zero: 474 return True 475 476 def _eval_is_positive(self): 477 # Double factorial is positive for every nonnegative input, and for 478 # every odd negative input which is of the form -1-4k for an 479 # nonnegative integer k 480 n = self.args[0] 481 if n.is_integer: 482 if (n + 1).is_nonnegative: 483 return True 484 if n.is_odd: 485 return ((n + 1) / 2).is_even 486 487 def _eval_rewrite_as_gamma(self, n, piecewise=True, **kwargs): 488 from sympy import gamma, Piecewise, sqrt 489 return 2**(n/2)*gamma(n/2 + 1) * Piecewise((1, Eq(Mod(n, 2), 0)), 490 (sqrt(2/pi), Eq(Mod(n, 2), 1))) 491 492 493############################################################################### 494######################## RISING and FALLING FACTORIALS ######################## 495############################################################################### 496 497 498class RisingFactorial(CombinatorialFunction): 499 r""" 500 Rising factorial (also called Pochhammer symbol) is a double valued 501 function arising in concrete mathematics, hypergeometric functions 502 and series expansions. It is defined by: 503 504 .. math:: rf(x,k) = x \cdot (x+1) \cdots (x+k-1) 505 506 where `x` can be arbitrary expression and `k` is an integer. For 507 more information check "Concrete mathematics" by Graham, pp. 66 508 or visit http://mathworld.wolfram.com/RisingFactorial.html page. 509 510 When `x` is a Poly instance of degree >= 1 with a single variable, 511 `rf(x,k) = x(y) \cdot x(y+1) \cdots x(y+k-1)`, where `y` is the 512 variable of `x`. This is as described in Peter Paule, "Greatest 513 Factorial Factorization and Symbolic Summation", Journal of 514 Symbolic Computation, vol. 20, pp. 235-268, 1995. 515 516 Examples 517 ======== 518 519 >>> from sympy import rf, Poly 520 >>> from sympy.abc import x 521 >>> rf(x, 0) 522 1 523 >>> rf(1, 5) 524 120 525 >>> rf(x, 5) == x*(1 + x)*(2 + x)*(3 + x)*(4 + x) 526 True 527 >>> rf(Poly(x**3, x), 2) 528 Poly(x**6 + 3*x**5 + 3*x**4 + x**3, x, domain='ZZ') 529 530 Rewriting is complicated unless the relationship between 531 the arguments is known, but rising factorial can 532 be rewritten in terms of gamma, factorial and binomial 533 and falling factorial. 534 535 >>> from sympy import Symbol, factorial, ff, binomial, gamma 536 >>> n = Symbol('n', integer=True, positive=True) 537 >>> R = rf(n, n + 2) 538 >>> for i in (rf, ff, factorial, binomial, gamma): 539 ... R.rewrite(i) 540 ... 541 RisingFactorial(n, n + 2) 542 FallingFactorial(2*n + 1, n + 2) 543 factorial(2*n + 1)/factorial(n - 1) 544 binomial(2*n + 1, n + 2)*factorial(n + 2) 545 gamma(2*n + 2)/gamma(n) 546 547 See Also 548 ======== 549 550 factorial, factorial2, FallingFactorial 551 552 References 553 ========== 554 555 .. [1] https://en.wikipedia.org/wiki/Pochhammer_symbol 556 557 """ 558 559 @classmethod 560 def eval(cls, x, k): 561 x = sympify(x) 562 k = sympify(k) 563 564 if x is S.NaN or k is S.NaN: 565 return S.NaN 566 elif x is S.One: 567 return factorial(k) 568 elif k.is_Integer: 569 if k.is_zero: 570 return S.One 571 else: 572 if k.is_positive: 573 if x is S.Infinity: 574 return S.Infinity 575 elif x is S.NegativeInfinity: 576 if k.is_odd: 577 return S.NegativeInfinity 578 else: 579 return S.Infinity 580 else: 581 if isinstance(x, Poly): 582 gens = x.gens 583 if len(gens)!= 1: 584 raise ValueError("rf only defined for " 585 "polynomials on one generator") 586 else: 587 return reduce(lambda r, i: 588 r*(x.shift(i)), 589 range(0, int(k)), 1) 590 else: 591 return reduce(lambda r, i: r*(x + i), 592 range(0, int(k)), 1) 593 594 else: 595 if x is S.Infinity: 596 return S.Infinity 597 elif x is S.NegativeInfinity: 598 return S.Infinity 599 else: 600 if isinstance(x, Poly): 601 gens = x.gens 602 if len(gens)!= 1: 603 raise ValueError("rf only defined for " 604 "polynomials on one generator") 605 else: 606 return 1/reduce(lambda r, i: 607 r*(x.shift(-i)), 608 range(1, abs(int(k)) + 1), 1) 609 else: 610 return 1/reduce(lambda r, i: 611 r*(x - i), 612 range(1, abs(int(k)) + 1), 1) 613 614 if k.is_integer == False: 615 if x.is_integer and x.is_negative: 616 return S.Zero 617 618 def _eval_rewrite_as_gamma(self, x, k, piecewise=True, **kwargs): 619 from sympy import gamma, Piecewise 620 if not piecewise: 621 if (x <= 0) == True: 622 return (-1)**k*gamma(1 - x) / gamma(-k - x + 1) 623 return gamma(x + k) / gamma(x) 624 return Piecewise( 625 (gamma(x + k) / gamma(x), x > 0), 626 ((-1)**k*gamma(1 - x) / gamma(-k - x + 1), True)) 627 628 def _eval_rewrite_as_FallingFactorial(self, x, k, **kwargs): 629 return FallingFactorial(x + k - 1, k) 630 631 def _eval_rewrite_as_factorial(self, x, k, **kwargs): 632 from sympy import Piecewise 633 if x.is_integer and k.is_integer: 634 return Piecewise( 635 (factorial(k + x - 1)/factorial(x - 1), x > 0), 636 ((-1)**k*factorial(-x)/factorial(-k - x), True)) 637 638 def _eval_rewrite_as_binomial(self, x, k, **kwargs): 639 if k.is_integer: 640 return factorial(k) * binomial(x + k - 1, k) 641 642 def _eval_rewrite_as_tractable(self, x, k, limitvar=None, **kwargs): 643 from sympy import gamma 644 if limitvar: 645 k_lim = k.subs(limitvar, S.Infinity) 646 if k_lim is S.Infinity: 647 return (gamma(x + k).rewrite('tractable', deep=True) / gamma(x)) 648 elif k_lim is S.NegativeInfinity: 649 return ((-1)**k*gamma(1 - x) / gamma(-k - x + 1).rewrite('tractable', deep=True)) 650 return self.rewrite(gamma).rewrite('tractable', deep=True) 651 652 def _eval_is_integer(self): 653 return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer, 654 self.args[1].is_nonnegative)) 655 656 657class FallingFactorial(CombinatorialFunction): 658 r""" 659 Falling factorial (related to rising factorial) is a double valued 660 function arising in concrete mathematics, hypergeometric functions 661 and series expansions. It is defined by 662 663 .. math:: ff(x,k) = x \cdot (x-1) \cdots (x-k+1) 664 665 where `x` can be arbitrary expression and `k` is an integer. For 666 more information check "Concrete mathematics" by Graham, pp. 66 667 or visit http://mathworld.wolfram.com/FallingFactorial.html page. 668 669 When `x` is a Poly instance of degree >= 1 with single variable, 670 `ff(x,k) = x(y) \cdot x(y-1) \cdots x(y-k+1)`, where `y` is the 671 variable of `x`. This is as described in Peter Paule, "Greatest 672 Factorial Factorization and Symbolic Summation", Journal of 673 Symbolic Computation, vol. 20, pp. 235-268, 1995. 674 675 >>> from sympy import ff, Poly, Symbol 676 >>> from sympy.abc import x 677 >>> n = Symbol('n', integer=True) 678 679 >>> ff(x, 0) 680 1 681 >>> ff(5, 5) 682 120 683 >>> ff(x, 5) == x*(x - 1)*(x - 2)*(x - 3)*(x - 4) 684 True 685 >>> ff(Poly(x**2, x), 2) 686 Poly(x**4 - 2*x**3 + x**2, x, domain='ZZ') 687 >>> ff(n, n) 688 factorial(n) 689 690 Rewriting is complicated unless the relationship between 691 the arguments is known, but falling factorial can 692 be rewritten in terms of gamma, factorial and binomial 693 and rising factorial. 694 695 >>> from sympy import factorial, rf, gamma, binomial, Symbol 696 >>> n = Symbol('n', integer=True, positive=True) 697 >>> F = ff(n, n - 2) 698 >>> for i in (rf, ff, factorial, binomial, gamma): 699 ... F.rewrite(i) 700 ... 701 RisingFactorial(3, n - 2) 702 FallingFactorial(n, n - 2) 703 factorial(n)/2 704 binomial(n, n - 2)*factorial(n - 2) 705 gamma(n + 1)/2 706 707 See Also 708 ======== 709 710 factorial, factorial2, RisingFactorial 711 712 References 713 ========== 714 715 .. [1] http://mathworld.wolfram.com/FallingFactorial.html 716 717 """ 718 719 @classmethod 720 def eval(cls, x, k): 721 x = sympify(x) 722 k = sympify(k) 723 724 if x is S.NaN or k is S.NaN: 725 return S.NaN 726 elif k.is_integer and x == k: 727 return factorial(x) 728 elif k.is_Integer: 729 if k.is_zero: 730 return S.One 731 else: 732 if k.is_positive: 733 if x is S.Infinity: 734 return S.Infinity 735 elif x is S.NegativeInfinity: 736 if k.is_odd: 737 return S.NegativeInfinity 738 else: 739 return S.Infinity 740 else: 741 if isinstance(x, Poly): 742 gens = x.gens 743 if len(gens)!= 1: 744 raise ValueError("ff only defined for " 745 "polynomials on one generator") 746 else: 747 return reduce(lambda r, i: 748 r*(x.shift(-i)), 749 range(0, int(k)), 1) 750 else: 751 return reduce(lambda r, i: r*(x - i), 752 range(0, int(k)), 1) 753 else: 754 if x is S.Infinity: 755 return S.Infinity 756 elif x is S.NegativeInfinity: 757 return S.Infinity 758 else: 759 if isinstance(x, Poly): 760 gens = x.gens 761 if len(gens)!= 1: 762 raise ValueError("rf only defined for " 763 "polynomials on one generator") 764 else: 765 return 1/reduce(lambda r, i: 766 r*(x.shift(i)), 767 range(1, abs(int(k)) + 1), 1) 768 else: 769 return 1/reduce(lambda r, i: r*(x + i), 770 range(1, abs(int(k)) + 1), 1) 771 772 def _eval_rewrite_as_gamma(self, x, k, piecewise=True, **kwargs): 773 from sympy import gamma, Piecewise 774 if not piecewise: 775 if (x < 0) == True: 776 return (-1)**k*gamma(k - x) / gamma(-x) 777 return gamma(x + 1) / gamma(x - k + 1) 778 return Piecewise( 779 (gamma(x + 1) / gamma(x - k + 1), x >= 0), 780 ((-1)**k*gamma(k - x) / gamma(-x), True)) 781 782 def _eval_rewrite_as_RisingFactorial(self, x, k, **kwargs): 783 return rf(x - k + 1, k) 784 785 def _eval_rewrite_as_binomial(self, x, k, **kwargs): 786 if k.is_integer: 787 return factorial(k) * binomial(x, k) 788 789 def _eval_rewrite_as_factorial(self, x, k, **kwargs): 790 from sympy import Piecewise 791 if x.is_integer and k.is_integer: 792 return Piecewise( 793 (factorial(x)/factorial(-k + x), x >= 0), 794 ((-1)**k*factorial(k - x - 1)/factorial(-x - 1), True)) 795 796 def _eval_rewrite_as_tractable(self, x, k, limitvar=None, **kwargs): 797 from sympy import gamma 798 if limitvar: 799 k_lim = k.subs(limitvar, S.Infinity) 800 if k_lim is S.Infinity: 801 return ((-1)**k*gamma(k - x).rewrite('tractable', deep=True) / gamma(-x)) 802 elif k_lim is S.NegativeInfinity: 803 return (gamma(x + 1) / gamma(x - k + 1).rewrite('tractable', deep=True)) 804 return self.rewrite(gamma).rewrite('tractable', deep=True) 805 806 def _eval_is_integer(self): 807 return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer, 808 self.args[1].is_nonnegative)) 809 810 811rf = RisingFactorial 812ff = FallingFactorial 813 814############################################################################### 815########################### BINOMIAL COEFFICIENTS ############################# 816############################################################################### 817 818 819class binomial(CombinatorialFunction): 820 r"""Implementation of the binomial coefficient. It can be defined 821 in two ways depending on its desired interpretation: 822 823 .. math:: \binom{n}{k} = \frac{n!}{k!(n-k)!}\ \text{or}\ 824 \binom{n}{k} = \frac{ff(n, k)}{k!} 825 826 First, in a strict combinatorial sense it defines the 827 number of ways we can choose `k` elements from a set of 828 `n` elements. In this case both arguments are nonnegative 829 integers and binomial is computed using an efficient 830 algorithm based on prime factorization. 831 832 The other definition is generalization for arbitrary `n`, 833 however `k` must also be nonnegative. This case is very 834 useful when evaluating summations. 835 836 For the sake of convenience for negative integer `k` this function 837 will return zero no matter what valued is the other argument. 838 839 To expand the binomial when `n` is a symbol, use either 840 ``expand_func()`` or ``expand(func=True)``. The former will keep 841 the polynomial in factored form while the latter will expand the 842 polynomial itself. See examples for details. 843 844 Examples 845 ======== 846 847 >>> from sympy import Symbol, Rational, binomial, expand_func 848 >>> n = Symbol('n', integer=True, positive=True) 849 850 >>> binomial(15, 8) 851 6435 852 853 >>> binomial(n, -1) 854 0 855 856 Rows of Pascal's triangle can be generated with the binomial function: 857 858 >>> for N in range(8): 859 ... print([binomial(N, i) for i in range(N + 1)]) 860 ... 861 [1] 862 [1, 1] 863 [1, 2, 1] 864 [1, 3, 3, 1] 865 [1, 4, 6, 4, 1] 866 [1, 5, 10, 10, 5, 1] 867 [1, 6, 15, 20, 15, 6, 1] 868 [1, 7, 21, 35, 35, 21, 7, 1] 869 870 As can a given diagonal, e.g. the 4th diagonal: 871 872 >>> N = -4 873 >>> [binomial(N, i) for i in range(1 - N)] 874 [1, -4, 10, -20, 35] 875 876 >>> binomial(Rational(5, 4), 3) 877 -5/128 878 >>> binomial(Rational(-5, 4), 3) 879 -195/128 880 881 >>> binomial(n, 3) 882 binomial(n, 3) 883 884 >>> binomial(n, 3).expand(func=True) 885 n**3/6 - n**2/2 + n/3 886 887 >>> expand_func(binomial(n, 3)) 888 n*(n - 2)*(n - 1)/6 889 890 References 891 ========== 892 893 .. [1] https://www.johndcook.com/blog/binomial_coefficients/ 894 895 """ 896 897 def fdiff(self, argindex=1): 898 from sympy import polygamma 899 if argindex == 1: 900 # http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/01/ 901 n, k = self.args 902 return binomial(n, k)*(polygamma(0, n + 1) - \ 903 polygamma(0, n - k + 1)) 904 elif argindex == 2: 905 # http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/02/ 906 n, k = self.args 907 return binomial(n, k)*(polygamma(0, n - k + 1) - \ 908 polygamma(0, k + 1)) 909 else: 910 raise ArgumentIndexError(self, argindex) 911 912 @classmethod 913 def _eval(self, n, k): 914 # n.is_Number and k.is_Integer and k != 1 and n != k 915 916 if k.is_Integer: 917 if n.is_Integer and n >= 0: 918 n, k = int(n), int(k) 919 920 if k > n: 921 return S.Zero 922 elif k > n // 2: 923 k = n - k 924 925 if HAS_GMPY: 926 from sympy.core.compatibility import gmpy 927 return Integer(gmpy.bincoef(n, k)) 928 929 d, result = n - k, 1 930 for i in range(1, k + 1): 931 d += 1 932 result = result * d // i 933 return Integer(result) 934 else: 935 d, result = n - k, 1 936 for i in range(1, k + 1): 937 d += 1 938 result *= d 939 result /= i 940 return result 941 942 @classmethod 943 def eval(cls, n, k): 944 n, k = map(sympify, (n, k)) 945 d = n - k 946 n_nonneg, n_isint = n.is_nonnegative, n.is_integer 947 if k.is_zero or ((n_nonneg or n_isint is False) 948 and d.is_zero): 949 return S.One 950 if (k - 1).is_zero or ((n_nonneg or n_isint is False) 951 and (d - 1).is_zero): 952 return n 953 if k.is_integer: 954 if k.is_negative or (n_nonneg and n_isint and d.is_negative): 955 return S.Zero 956 elif n.is_number: 957 res = cls._eval(n, k) 958 return res.expand(basic=True) if res else res 959 elif n_nonneg is False and n_isint: 960 # a special case when binomial evaluates to complex infinity 961 return S.ComplexInfinity 962 elif k.is_number: 963 from sympy import gamma 964 return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1)) 965 966 def _eval_Mod(self, q): 967 n, k = self.args 968 969 if any(x.is_integer is False for x in (n, k, q)): 970 raise ValueError("Integers expected for binomial Mod") 971 972 if all(x.is_Integer for x in (n, k, q)): 973 n, k = map(int, (n, k)) 974 aq, res = abs(q), 1 975 976 # handle negative integers k or n 977 if k < 0: 978 return S.Zero 979 if n < 0: 980 n = -n + k - 1 981 res = -1 if k%2 else 1 982 983 # non negative integers k and n 984 if k > n: 985 return S.Zero 986 987 isprime = aq.is_prime 988 aq = int(aq) 989 if isprime: 990 if aq < n: 991 # use Lucas Theorem 992 N, K = n, k 993 while N or K: 994 res = res*binomial(N % aq, K % aq) % aq 995 N, K = N // aq, K // aq 996 997 else: 998 # use Factorial Modulo 999 d = n - k 1000 if k > d: 1001 k, d = d, k 1002 kf = 1 1003 for i in range(2, k + 1): 1004 kf = kf*i % aq 1005 df = kf 1006 for i in range(k + 1, d + 1): 1007 df = df*i % aq 1008 res *= df 1009 for i in range(d + 1, n + 1): 1010 res = res*i % aq 1011 1012 res *= pow(kf*df % aq, aq - 2, aq) 1013 res %= aq 1014 1015 else: 1016 # Binomial Factorization is performed by calculating the 1017 # exponents of primes <= n in `n! /(k! (n - k)!)`, 1018 # for non-negative integers n and k. As the exponent of 1019 # prime in n! is e_p(n) = [n/p] + [n/p**2] + ... 1020 # the exponent of prime in binomial(n, k) would be 1021 # e_p(n) - e_p(k) - e_p(n - k) 1022 M = int(_sqrt(n)) 1023 for prime in sieve.primerange(2, n + 1): 1024 if prime > n - k: 1025 res = res*prime % aq 1026 elif prime > n // 2: 1027 continue 1028 elif prime > M: 1029 if n % prime < k % prime: 1030 res = res*prime % aq 1031 else: 1032 N, K = n, k 1033 exp = a = 0 1034 1035 while N > 0: 1036 a = int((N % prime) < (K % prime + a)) 1037 N, K = N // prime, K // prime 1038 exp += a 1039 1040 if exp > 0: 1041 res *= pow(prime, exp, aq) 1042 res %= aq 1043 1044 return S(res % q) 1045 1046 def _eval_expand_func(self, **hints): 1047 """ 1048 Function to expand binomial(n, k) when m is positive integer 1049 Also, 1050 n is self.args[0] and k is self.args[1] while using binomial(n, k) 1051 """ 1052 n = self.args[0] 1053 if n.is_Number: 1054 return binomial(*self.args) 1055 1056 k = self.args[1] 1057 if (n-k).is_Integer: 1058 k = n - k 1059 1060 if k.is_Integer: 1061 if k.is_zero: 1062 return S.One 1063 elif k.is_negative: 1064 return S.Zero 1065 else: 1066 n, result = self.args[0], 1 1067 for i in range(1, k + 1): 1068 result *= n - k + i 1069 result /= i 1070 return result 1071 else: 1072 return binomial(*self.args) 1073 1074 def _eval_rewrite_as_factorial(self, n, k, **kwargs): 1075 return factorial(n)/(factorial(k)*factorial(n - k)) 1076 1077 def _eval_rewrite_as_gamma(self, n, k, piecewise=True, **kwargs): 1078 from sympy import gamma 1079 return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1)) 1080 1081 def _eval_rewrite_as_tractable(self, n, k, limitvar=None, **kwargs): 1082 return self._eval_rewrite_as_gamma(n, k).rewrite('tractable') 1083 1084 def _eval_rewrite_as_FallingFactorial(self, n, k, **kwargs): 1085 if k.is_integer: 1086 return ff(n, k) / factorial(k) 1087 1088 def _eval_is_integer(self): 1089 n, k = self.args 1090 if n.is_integer and k.is_integer: 1091 return True 1092 elif k.is_integer is False: 1093 return False 1094 1095 def _eval_is_nonnegative(self): 1096 n, k = self.args 1097 if n.is_integer and k.is_integer: 1098 if n.is_nonnegative or k.is_negative or k.is_even: 1099 return True 1100 elif k.is_even is False: 1101 return False 1102 1103 def _eval_as_leading_term(self, x, logx=None, cdir=0): 1104 from sympy import gamma 1105 return self.rewrite(gamma)._eval_as_leading_term(x, logx=logx, cdir=cdir) 1106