Coverage for /usr/lib/python3/dist-packages/sympy/functions/combinatorial/factorials.py: 17%

1from __future__ import annotations

2from functools import reduce

4from sympy.core import S, sympify, Dummy, Mod

5from sympy.core.cache import cacheit

6from sympy.core.function import Function, ArgumentIndexError, PoleError

7from sympy.core.logic import fuzzy_and

8from sympy.core.numbers import Integer, pi, I

9from sympy.core.relational import Eq

10from sympy.external.gmpy import HAS_GMPY, gmpy

11from sympy.ntheory import sieve

12from sympy.polys.polytools import Poly

14from math import factorial as _factorial, prod, sqrt as _sqrt

16class CombinatorialFunction(Function):

17 """Base class for combinatorial functions. """

19 def _eval_simplify(self, **kwargs):

20 from sympy.simplify.combsimp import combsimp

21 # combinatorial function with non-integer arguments is

22 # automatically passed to gammasimp

23 expr = combsimp(self)

24 measure = kwargs['measure']

25 if measure(expr) <= kwargs['ratio']*measure(self):

26 return expr

27 return self

30###############################################################################

31######################## FACTORIAL and MULTI-FACTORIAL ########################

32###############################################################################

35class factorial(CombinatorialFunction):

36 r"""Implementation of factorial function over nonnegative integers.

37 By convention (consistent with the gamma function and the binomial

38 coefficients), factorial of a negative integer is complex infinity.

40 The factorial is very important in combinatorics where it gives

41 the number of ways in which `n` objects can be permuted. It also

42 arises in calculus, probability, number theory, etc.

44 There is strict relation of factorial with gamma function. In

45 fact `n! = gamma(n+1)` for nonnegative integers. Rewrite of this

46 kind is very useful in case of combinatorial simplification.

48 Computation of the factorial is done using two algorithms. For

49 small arguments a precomputed look up table is used. However for bigger

50 input algorithm Prime-Swing is used. It is the fastest algorithm

51 known and computes `n!` via prime factorization of special class

52 of numbers, called here the 'Swing Numbers'.

54 Examples

55 ========

57 >>> from sympy import Symbol, factorial, S

58 >>> n = Symbol('n', integer=True)

60 >>> factorial(0)

61 1

63 >>> factorial(7)

64 5040

66 >>> factorial(-2)

67 zoo

69 >>> factorial(n)

70 factorial(n)

72 >>> factorial(2*n)

73 factorial(2*n)

75 >>> factorial(S(1)/2)

76 factorial(1/2)

78 See Also

79 ========

81 factorial2, RisingFactorial, FallingFactorial

82 """

84 def fdiff(self, argindex=1):

85 from sympy.functions.special.gamma_functions import (gamma, polygamma)

86 if argindex == 1:

87 return gamma(self.args[0] + 1)*polygamma(0, self.args[0] + 1)

88 else:

89 raise ArgumentIndexError(self, argindex)

91 _small_swing = [

92 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429, 6435, 6435, 109395,

93 12155, 230945, 46189, 969969, 88179, 2028117, 676039, 16900975, 1300075,

94 35102025, 5014575, 145422675, 9694845, 300540195, 300540195

95 ]

97 _small_factorials: list[int] = []

99 @classmethod

100 def _swing(cls, n):

101 if n < 33:

102 return cls._small_swing[n]

103 else:

104 N, primes = int(_sqrt(n)), []

105

106 for prime in sieve.primerange(3, N + 1):

107 p, q = 1, n

108

109 while True:

110 q //= prime

111

112 if q > 0:

113 if q & 1 == 1:

114 p *= prime

115 else:

116 break

117

118 if p > 1:

119 primes.append(p)

120

121 for prime in sieve.primerange(N + 1, n//3 + 1):

122 if (n // prime) & 1 == 1:

123 primes.append(prime)

124

125 L_product = prod(sieve.primerange(n//2 + 1, n + 1))

126 R_product = prod(primes)

127

128 return L_product*R_product

129

130 @classmethod

131 def _recursive(cls, n):

132 if n < 2:

133 return 1

134 else:

135 return (cls._recursive(n//2)**2)*cls._swing(n)

136

137 @classmethod

138 def eval(cls, n):

139 n = sympify(n)

140

141 if n.is_Number:

142 if n.is_zero:

143 return S.One

144 elif n is S.Infinity:

145 return S.Infinity

146 elif n.is_Integer:

147 if n.is_negative:

148 return S.ComplexInfinity

149 else:

150 n = n.p

151

152 if n < 20:

153 if not cls._small_factorials:

154 result = 1

155 for i in range(1, 20):

156 result *= i

157 cls._small_factorials.append(result)

158 result = cls._small_factorials[n-1]

159

160 # GMPY factorial is faster, use it when available

161 elif HAS_GMPY:

162 result = gmpy.fac(n)

163

164 else:

165 bits = bin(n).count('1')

166 result = cls._recursive(n)*2**(n - bits)

167

168 return Integer(result)

169

170 def _facmod(self, n, q):

171 res, N = 1, int(_sqrt(n))

172

173 # Exponent of prime p in n! is e_p(n) = [n/p] + [n/p**2] + ...

174 # for p > sqrt(n), e_p(n) < sqrt(n), the primes with [n/p] = m,

175 # occur consecutively and are grouped together in pw[m] for

176 # simultaneous exponentiation at a later stage

177 pw = [1]*N

178

179 m = 2 # to initialize the if condition below

180 for prime in sieve.primerange(2, n + 1):

181 if m > 1:

182 m, y = 0, n // prime

183 while y:

184 m += y

185 y //= prime

186 if m < N:

187 pw[m] = pw[m]*prime % q

188 else:

189 res = res*pow(prime, m, q) % q

190

191 for ex, bs in enumerate(pw):

192 if ex == 0 or bs == 1:

193 continue

194 if bs == 0:

195 return 0

196 res = res*pow(bs, ex, q) % q

197

198 return res

199

200 def _eval_Mod(self, q):

201 n = self.args[0]

202 if n.is_integer and n.is_nonnegative and q.is_integer:

203 aq = abs(q)

204 d = aq - n

205 if d.is_nonpositive:

206 return S.Zero

207 else:

208 isprime = aq.is_prime

209 if d == 1:

210 # Apply Wilson's theorem (if a natural number n > 1

211 # is a prime number, then (n-1)! = -1 mod n) and

212 # its inverse (if n > 4 is a composite number, then

213 # (n-1)! = 0 mod n)

214 if isprime:

215 return -1 % q

216 elif isprime is False and (aq - 6).is_nonnegative:

217 return S.Zero

218 elif n.is_Integer and q.is_Integer:

219 n, d, aq = map(int, (n, d, aq))

220 if isprime and (d - 1 < n):

221 fc = self._facmod(d - 1, aq)

222 fc = pow(fc, aq - 2, aq)

223 if d%2:

224 fc = -fc

225 else:

226 fc = self._facmod(n, aq)

227

228 return fc % q

229

230 def _eval_rewrite_as_gamma(self, n, piecewise=True, **kwargs):

231 from sympy.functions.special.gamma_functions import gamma

232 return gamma(n + 1)

233

234 def _eval_rewrite_as_Product(self, n, **kwargs):

235 from sympy.concrete.products import Product

236 if n.is_nonnegative and n.is_integer:

237 i = Dummy('i', integer=True)

238 return Product(i, (i, 1, n))

239

240 def _eval_is_integer(self):

241 if self.args[0].is_integer and self.args[0].is_nonnegative:

242 return True

243

244 def _eval_is_positive(self):

245 if self.args[0].is_integer and self.args[0].is_nonnegative:

246 return True

247

248 def _eval_is_even(self):

249 x = self.args[0]

250 if x.is_integer and x.is_nonnegative:

251 return (x - 2).is_nonnegative

252

253 def _eval_is_composite(self):

254 x = self.args[0]

255 if x.is_integer and x.is_nonnegative:

256 return (x - 3).is_nonnegative

257

258 def _eval_is_real(self):

259 x = self.args[0]

260 if x.is_nonnegative or x.is_noninteger:

261 return True

262

263 def _eval_as_leading_term(self, x, logx=None, cdir=0):

264 arg = self.args[0].as_leading_term(x)

265 arg0 = arg.subs(x, 0)

266 if arg0.is_zero:

267 return S.One

268 elif not arg0.is_infinite:

269 return self.func(arg)

270 raise PoleError("Cannot expand %s around 0" % (self))

271

272class MultiFactorial(CombinatorialFunction):

273 pass

274

275

276class subfactorial(CombinatorialFunction):

277 r"""The subfactorial counts the derangements of $n$ items and is

278 defined for non-negative integers as:

279

280 .. math:: !n = \begin{cases} 1 & n = 0 \\ 0 & n = 1 \\

281 (n-1)(!(n-1) + !(n-2)) & n > 1 \end{cases}

282

283 It can also be written as ``int(round(n!/exp(1)))`` but the

284 recursive definition with caching is implemented for this function.

285

286 An interesting analytic expression is the following [2]_

287

288 .. math:: !x = \Gamma(x + 1, -1)/e

289

290 which is valid for non-negative integers `x`. The above formula

291 is not very useful in case of non-integers. `\Gamma(x + 1, -1)` is

292 single-valued only for integral arguments `x`, elsewhere on the positive

293 real axis it has an infinite number of branches none of which are real.

294

295 References

296 ==========

297

298 .. [1] https://en.wikipedia.org/wiki/Subfactorial

299 .. [2] https://mathworld.wolfram.com/Subfactorial.html

300

301 Examples

302 ========

303

304 >>> from sympy import subfactorial

305 >>> from sympy.abc import n

306 >>> subfactorial(n + 1)

307 subfactorial(n + 1)

308 >>> subfactorial(5)

309 44

310

311 See Also

312 ========

313

314 factorial, uppergamma,

315 sympy.utilities.iterables.generate_derangements

316 """

317

318 @classmethod

319 @cacheit

320 def _eval(self, n):

321 if not n:

322 return S.One

323 elif n == 1:

324 return S.Zero

325 else:

326 z1, z2 = 1, 0

327 for i in range(2, n + 1):

328 z1, z2 = z2, (i - 1)*(z2 + z1)

329 return z2

330

331 @classmethod

332 def eval(cls, arg):

333 if arg.is_Number:

334 if arg.is_Integer and arg.is_nonnegative:

335 return cls._eval(arg)

336 elif arg is S.NaN:

337 return S.NaN

338 elif arg is S.Infinity:

339 return S.Infinity

340

341 def _eval_is_even(self):

342 if self.args[0].is_odd and self.args[0].is_nonnegative:

343 return True

344

345 def _eval_is_integer(self):

346 if self.args[0].is_integer and self.args[0].is_nonnegative:

347 return True

348

349 def _eval_rewrite_as_factorial(self, arg, **kwargs):

350 from sympy.concrete.summations import summation

351 i = Dummy('i')

352 f = S.NegativeOne**i / factorial(i)

353 return factorial(arg) * summation(f, (i, 0, arg))

354

355 def _eval_rewrite_as_gamma(self, arg, piecewise=True, **kwargs):

356 from sympy.functions.elementary.exponential import exp

357 from sympy.functions.special.gamma_functions import (gamma, lowergamma)

358 return (S.NegativeOne**(arg + 1)*exp(-I*pi*arg)*lowergamma(arg + 1, -1)

359 + gamma(arg + 1))*exp(-1)

360

361 def _eval_rewrite_as_uppergamma(self, arg, **kwargs):

362 from sympy.functions.special.gamma_functions import uppergamma

363 return uppergamma(arg + 1, -1)/S.Exp1

364

365 def _eval_is_nonnegative(self):

366 if self.args[0].is_integer and self.args[0].is_nonnegative:

367 return True

368

369 def _eval_is_odd(self):

370 if self.args[0].is_even and self.args[0].is_nonnegative:

371 return True

372

373

374class factorial2(CombinatorialFunction):

375 r"""The double factorial `n!!`, not to be confused with `(n!)!`

376

377 The double factorial is defined for nonnegative integers and for odd

378 negative integers as:

379

380 .. math:: n!! = \begin{cases} 1 & n = 0 \\

381 n(n-2)(n-4) \cdots 1 & n\ \text{positive odd} \\

382 n(n-2)(n-4) \cdots 2 & n\ \text{positive even} \\

383 (n+2)!!/(n+2) & n\ \text{negative odd} \end{cases}

384

385 References

386 ==========

387

388 .. [1] https://en.wikipedia.org/wiki/Double_factorial

389

390 Examples

391 ========

392

393 >>> from sympy import factorial2, var

394 >>> n = var('n')

395 >>> n

396 n

397 >>> factorial2(n + 1)

398 factorial2(n + 1)

399 >>> factorial2(5)

400 15

401 >>> factorial2(-1)

402 1

403 >>> factorial2(-5)

404 1/3

405

406 See Also

407 ========

408

409 factorial, RisingFactorial, FallingFactorial

410 """

411

412 @classmethod

413 def eval(cls, arg):

414 # TODO: extend this to complex numbers?

415

416 if arg.is_Number:

417 if not arg.is_Integer:

418 raise ValueError("argument must be nonnegative integer "

419 "or negative odd integer")

420

421 # This implementation is faster than the recursive one

422 # It also avoids "maximum recursion depth exceeded" runtime error

423 if arg.is_nonnegative:

424 if arg.is_even:

425 k = arg / 2

426 return 2**k * factorial(k)

427 return factorial(arg) / factorial2(arg - 1)

428

429

430 if arg.is_odd:

431 return arg*(S.NegativeOne)**((1 - arg)/2) / factorial2(-arg)

432 raise ValueError("argument must be nonnegative integer "

433 "or negative odd integer")

434

435

436 def _eval_is_even(self):

437 # Double factorial is even for every positive even input

438 n = self.args[0]

439 if n.is_integer:

440 if n.is_odd:

441 return False

442 if n.is_even:

443 if n.is_positive:

444 return True

445 if n.is_zero:

446 return False

447

448 def _eval_is_integer(self):

449 # Double factorial is an integer for every nonnegative input, and for

450 # -1 and -3

451 n = self.args[0]

452 if n.is_integer:

453 if (n + 1).is_nonnegative:

454 return True

455 if n.is_odd:

456 return (n + 3).is_nonnegative

457

458 def _eval_is_odd(self):

459 # Double factorial is odd for every odd input not smaller than -3, and

460 # for 0

461 n = self.args[0]

462 if n.is_odd:

463 return (n + 3).is_nonnegative

464 if n.is_even:

465 if n.is_positive:

466 return False

467 if n.is_zero:

468 return True

469

470 def _eval_is_positive(self):

471 # Double factorial is positive for every nonnegative input, and for

472 # every odd negative input which is of the form -1-4k for an

473 # nonnegative integer k

474 n = self.args[0]

475 if n.is_integer:

476 if (n + 1).is_nonnegative:

477 return True

478 if n.is_odd:

479 return ((n + 1) / 2).is_even

480

481 def _eval_rewrite_as_gamma(self, n, piecewise=True, **kwargs):

482 from sympy.functions.elementary.miscellaneous import sqrt

483 from sympy.functions.elementary.piecewise import Piecewise

484 from sympy.functions.special.gamma_functions import gamma

485 return 2**(n/2)*gamma(n/2 + 1) * Piecewise((1, Eq(Mod(n, 2), 0)),

486 (sqrt(2/pi), Eq(Mod(n, 2), 1)))

487

488

489###############################################################################

490######################## RISING and FALLING FACTORIALS ########################

491###############################################################################

492

493

494class RisingFactorial(CombinatorialFunction):

495 r"""

496 Rising factorial (also called Pochhammer symbol [1]_) is a double valued

497 function arising in concrete mathematics, hypergeometric functions

498 and series expansions. It is defined by:

499

500 .. math:: \texttt{rf(y, k)} = (x)^k = x \cdot (x+1) \cdots (x+k-1)

501

502 where `x` can be arbitrary expression and `k` is an integer. For

503 more information check "Concrete mathematics" by Graham, pp. 66

504 or visit https://mathworld.wolfram.com/RisingFactorial.html page.

505

506 When `x` is a `~.Poly` instance of degree $\ge 1$ with a single variable,

507 `(x)^k = x(y) \cdot x(y+1) \cdots x(y+k-1)`, where `y` is the

508 variable of `x`. This is as described in [2]_.

509

510 Examples

511 ========

512

513 >>> from sympy import rf, Poly

514 >>> from sympy.abc import x

515 >>> rf(x, 0)

516 1

517 >>> rf(1, 5)

518 120

519 >>> rf(x, 5) == x*(1 + x)*(2 + x)*(3 + x)*(4 + x)

520 True

521 >>> rf(Poly(x**3, x), 2)

522 Poly(x**6 + 3*x**5 + 3*x**4 + x**3, x, domain='ZZ')

523

524 Rewriting is complicated unless the relationship between

525 the arguments is known, but rising factorial can

526 be rewritten in terms of gamma, factorial, binomial,

527 and falling factorial.

528

529 >>> from sympy import Symbol, factorial, ff, binomial, gamma

530 >>> n = Symbol('n', integer=True, positive=True)

531 >>> R = rf(n, n + 2)

532 >>> for i in (rf, ff, factorial, binomial, gamma):

533 ... R.rewrite(i)

534 ...

535 RisingFactorial(n, n + 2)

536 FallingFactorial(2*n + 1, n + 2)

537 factorial(2*n + 1)/factorial(n - 1)

538 binomial(2*n + 1, n + 2)*factorial(n + 2)

539 gamma(2*n + 2)/gamma(n)

540

541 See Also

542 ========

543

544 factorial, factorial2, FallingFactorial

545

546 References

547 ==========

548

549 .. [1] https://en.wikipedia.org/wiki/Pochhammer_symbol

550 .. [2] Peter Paule, "Greatest Factorial Factorization and Symbolic

551 Summation", Journal of Symbolic Computation, vol. 20, pp. 235-268,

552 1995.

553

554 """

555

556 @classmethod

557 def eval(cls, x, k):

558 x = sympify(x)

559 k = sympify(k)

560

561 if x is S.NaN or k is S.NaN:

562 return S.NaN

563 elif x is S.One:

564 return factorial(k)

565 elif k.is_Integer:

566 if k.is_zero:

567 return S.One

568 else:

569 if k.is_positive:

570 if x is S.Infinity:

571 return S.Infinity

572 elif x is S.NegativeInfinity:

573 if k.is_odd:

574 return S.NegativeInfinity

575 else:

576 return S.Infinity

577 else:

578 if isinstance(x, Poly):

579 gens = x.gens

580 if len(gens)!= 1:

581 raise ValueError("rf only defined for "

582 "polynomials on one generator")

583 else:

584 return reduce(lambda r, i:

585 r*(x.shift(i)),

586 range(int(k)), 1)

587 else:

588 return reduce(lambda r, i: r*(x + i),

589 range(int(k)), 1)

590

591 else:

592 if x is S.Infinity:

593 return S.Infinity

594 elif x is S.NegativeInfinity:

595 return S.Infinity

596 else:

597 if isinstance(x, Poly):

598 gens = x.gens

599 if len(gens)!= 1:

600 raise ValueError("rf only defined for "

601 "polynomials on one generator")

602 else:

603 return 1/reduce(lambda r, i:

604 r*(x.shift(-i)),

605 range(1, abs(int(k)) + 1), 1)

606 else:

607 return 1/reduce(lambda r, i:

608 r*(x - i),

609 range(1, abs(int(k)) + 1), 1)

610

611 if k.is_integer == False:

612 if x.is_integer and x.is_negative:

613 return S.Zero

614

615 def _eval_rewrite_as_gamma(self, x, k, piecewise=True, **kwargs):

616 from sympy.functions.elementary.piecewise import Piecewise

617 from sympy.functions.special.gamma_functions import gamma

618 if not piecewise:

619 if (x <= 0) == True:

620 return S.NegativeOne**k*gamma(1 - x) / gamma(-k - x + 1)

621 return gamma(x + k) / gamma(x)

622 return Piecewise(

623 (gamma(x + k) / gamma(x), x > 0),

624 (S.NegativeOne**k*gamma(1 - x) / gamma(-k - x + 1), True))

625

626 def _eval_rewrite_as_FallingFactorial(self, x, k, **kwargs):

627 return FallingFactorial(x + k - 1, k)

628

629 def _eval_rewrite_as_factorial(self, x, k, **kwargs):

630 from sympy.functions.elementary.piecewise import Piecewise

631 if x.is_integer and k.is_integer:

632 return Piecewise(

633 (factorial(k + x - 1)/factorial(x - 1), x > 0),

634 (S.NegativeOne**k*factorial(-x)/factorial(-k - x), True))

635

636 def _eval_rewrite_as_binomial(self, x, k, **kwargs):

637 if k.is_integer:

638 return factorial(k) * binomial(x + k - 1, k)

639

640 def _eval_rewrite_as_tractable(self, x, k, limitvar=None, **kwargs):

641 from sympy.functions.special.gamma_functions import gamma

642 if limitvar:

643 k_lim = k.subs(limitvar, S.Infinity)

644 if k_lim is S.Infinity:

645 return (gamma(x + k).rewrite('tractable', deep=True) / gamma(x))

646 elif k_lim is S.NegativeInfinity:

647 return (S.NegativeOne**k*gamma(1 - x) / gamma(-k - x + 1).rewrite('tractable', deep=True))

648 return self.rewrite(gamma).rewrite('tractable', deep=True)

649

650 def _eval_is_integer(self):

651 return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer,

652 self.args[1].is_nonnegative))

653

654

655class FallingFactorial(CombinatorialFunction):

656 r"""

657 Falling factorial (related to rising factorial) is a double valued

658 function arising in concrete mathematics, hypergeometric functions

659 and series expansions. It is defined by

660

661 .. math:: \texttt{ff(x, k)} = (x)_k = x \cdot (x-1) \cdots (x-k+1)

662

663 where `x` can be arbitrary expression and `k` is an integer. For

664 more information check "Concrete mathematics" by Graham, pp. 66

665 or [1]_.

666

667 When `x` is a `~.Poly` instance of degree $\ge 1$ with single variable,

668 `(x)_k = x(y) \cdot x(y-1) \cdots x(y-k+1)`, where `y` is the

669 variable of `x`. This is as described in

670

671 >>> from sympy import ff, Poly, Symbol

672 >>> from sympy.abc import x

673 >>> n = Symbol('n', integer=True)

674

675 >>> ff(x, 0)

676 1

677 >>> ff(5, 5)

678 120

679 >>> ff(x, 5) == x*(x - 1)*(x - 2)*(x - 3)*(x - 4)

680 True

681 >>> ff(Poly(x**2, x), 2)

682 Poly(x**4 - 2*x**3 + x**2, x, domain='ZZ')

683 >>> ff(n, n)

684 factorial(n)

685

686 Rewriting is complicated unless the relationship between

687 the arguments is known, but falling factorial can

688 be rewritten in terms of gamma, factorial and binomial

689 and rising factorial.

690

691 >>> from sympy import factorial, rf, gamma, binomial, Symbol

692 >>> n = Symbol('n', integer=True, positive=True)

693 >>> F = ff(n, n - 2)

694 >>> for i in (rf, ff, factorial, binomial, gamma):

695 ... F.rewrite(i)

696 ...

697 RisingFactorial(3, n - 2)

698 FallingFactorial(n, n - 2)

699 factorial(n)/2

700 binomial(n, n - 2)*factorial(n - 2)

701 gamma(n + 1)/2

702

703 See Also

704 ========

705

706 factorial, factorial2, RisingFactorial

707

708 References

709 ==========

710

711 .. [1] https://mathworld.wolfram.com/FallingFactorial.html

712 .. [2] Peter Paule, "Greatest Factorial Factorization and Symbolic

713 Summation", Journal of Symbolic Computation, vol. 20, pp. 235-268,

714 1995.

715

716 """

717

718 @classmethod

719 def eval(cls, x, k):

720 x = sympify(x)

721 k = sympify(k)

722

723 if x is S.NaN or k is S.NaN:

724 return S.NaN

725 elif k.is_integer and x == k:

726 return factorial(x)

727 elif k.is_Integer:

728 if k.is_zero:

729 return S.One

730 else:

731 if k.is_positive:

732 if x is S.Infinity:

733 return S.Infinity

734 elif x is S.NegativeInfinity:

735 if k.is_odd:

736 return S.NegativeInfinity

737 else:

738 return S.Infinity

739 else:

740 if isinstance(x, Poly):

741 gens = x.gens

742 if len(gens)!= 1:

743 raise ValueError("ff only defined for "

744 "polynomials on one generator")

745 else:

746 return reduce(lambda r, i:

747 r*(x.shift(-i)),

748 range(int(k)), 1)

749 else:

750 return reduce(lambda r, i: r*(x - i),

751 range(int(k)), 1)

752 else:

753 if x is S.Infinity:

754 return S.Infinity

755 elif x is S.NegativeInfinity:

756 return S.Infinity

757 else:

758 if isinstance(x, Poly):

759 gens = x.gens

760 if len(gens)!= 1:

761 raise ValueError("rf only defined for "

762 "polynomials on one generator")

763 else:

764 return 1/reduce(lambda r, i:

765 r*(x.shift(i)),

766 range(1, abs(int(k)) + 1), 1)

767 else:

768 return 1/reduce(lambda r, i: r*(x + i),

769 range(1, abs(int(k)) + 1), 1)

770

771 def _eval_rewrite_as_gamma(self, x, k, piecewise=True, **kwargs):

772 from sympy.functions.elementary.piecewise import Piecewise

773 from sympy.functions.special.gamma_functions import gamma

774 if not piecewise:

775 if (x < 0) == True:

776 return S.NegativeOne**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 (S.NegativeOne**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.functions.elementary.piecewise import Piecewise

791 if x.is_integer and k.is_integer:

792 return Piecewise(

793 (factorial(x)/factorial(-k + x), x >= 0),

794 (S.NegativeOne**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.functions.special.gamma_functions import gamma

798 if limitvar:

799 k_lim = k.subs(limitvar, S.Infinity)

800 if k_lim is S.Infinity:

801 return (S.NegativeOne**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{(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 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.functions.special.gamma_functions import polygamma

899 if argindex == 1:

900 # https://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 # https://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 return Integer(gmpy.bincoef(n, k))

927

928 d, result = n - k, 1

929 for i in range(1, k + 1):

930 d += 1

931 result = result * d // i

932 return Integer(result)

933 else:

934 d, result = n - k, 1

935 for i in range(1, k + 1):

936 d += 1

937 result *= d

938 return result / _factorial(k)

939

940 @classmethod

941 def eval(cls, n, k):

942 n, k = map(sympify, (n, k))

943 d = n - k

944 n_nonneg, n_isint = n.is_nonnegative, n.is_integer

945 if k.is_zero or ((n_nonneg or n_isint is False)

946 and d.is_zero):

947 return S.One

948 if (k - 1).is_zero or ((n_nonneg or n_isint is False)

949 and (d - 1).is_zero):

950 return n

951 if k.is_integer:

952 if k.is_negative or (n_nonneg and n_isint and d.is_negative):

953 return S.Zero

954 elif n.is_number:

955 res = cls._eval(n, k)

956 return res.expand(basic=True) if res else res

957 elif n_nonneg is False and n_isint:

958 # a special case when binomial evaluates to complex infinity

959 return S.ComplexInfinity

960 elif k.is_number:

961 from sympy.functions.special.gamma_functions import gamma

962 return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1))

963

964 def _eval_Mod(self, q):

965 n, k = self.args

966

967 if any(x.is_integer is False for x in (n, k, q)):

968 raise ValueError("Integers expected for binomial Mod")

969

970 if all(x.is_Integer for x in (n, k, q)):

971 n, k = map(int, (n, k))

972 aq, res = abs(q), 1

973

974 # handle negative integers k or n

975 if k < 0:

976 return S.Zero

977 if n < 0:

978 n = -n + k - 1

979 res = -1 if k%2 else 1

980

981 # non negative integers k and n

982 if k > n:

983 return S.Zero

984

985 isprime = aq.is_prime

986 aq = int(aq)

987 if isprime:

988 if aq < n:

989 # use Lucas Theorem

990 N, K = n, k

991 while N or K:

992 res = res*binomial(N % aq, K % aq) % aq

993 N, K = N // aq, K // aq

994

995 else:

996 # use Factorial Modulo

997 d = n - k

998 if k > d:

999 k, d = d, k

1000 kf = 1

1001 for i in range(2, k + 1):

1002 kf = kf*i % aq

1003 df = kf

1004 for i in range(k + 1, d + 1):

1005 df = df*i % aq

1006 res *= df

1007 for i in range(d + 1, n + 1):

1008 res = res*i % aq

1009

1010 res *= pow(kf*df % aq, aq - 2, aq)

1011 res %= aq

1012

1013 else:

1014 # Binomial Factorization is performed by calculating the

1015 # exponents of primes <= n in `n! /(k! (n - k)!)`,

1016 # for non-negative integers n and k. As the exponent of

1017 # prime in n! is e_p(n) = [n/p] + [n/p**2] + ...

1018 # the exponent of prime in binomial(n, k) would be

1019 # e_p(n) - e_p(k) - e_p(n - k)

1020 M = int(_sqrt(n))

1021 for prime in sieve.primerange(2, n + 1):

1022 if prime > n - k:

1023 res = res*prime % aq

1024 elif prime > n // 2:

1025 continue

1026 elif prime > M:

1027 if n % prime < k % prime:

1028 res = res*prime % aq

1029 else:

1030 N, K = n, k

1031 exp = a = 0

1032

1033 while N > 0:

1034 a = int((N % prime) < (K % prime + a))

1035 N, K = N // prime, K // prime

1036 exp += a

1037

1038 if exp > 0:

1039 res *= pow(prime, exp, aq)

1040 res %= aq

1041

1042 return S(res % q)

1043

1044 def _eval_expand_func(self, **hints):

1045 """

1046 Function to expand binomial(n, k) when m is positive integer

1047 Also,

1048 n is self.args[0] and k is self.args[1] while using binomial(n, k)

1049 """

1050 n = self.args[0]

1051 if n.is_Number:

1052 return binomial(*self.args)

1053

1054 k = self.args[1]

1055 if (n-k).is_Integer:

1056 k = n - k

1057

1058 if k.is_Integer:

1059 if k.is_zero:

1060 return S.One

1061 elif k.is_negative:

1062 return S.Zero

1063 else:

1064 n, result = self.args[0], 1

1065 for i in range(1, k + 1):

1066 result *= n - k + i

1067 return result / _factorial(k)

1068 else:

1069 return binomial(*self.args)

1070

1071 def _eval_rewrite_as_factorial(self, n, k, **kwargs):

1072 return factorial(n)/(factorial(k)*factorial(n - k))

1073

1074 def _eval_rewrite_as_gamma(self, n, k, piecewise=True, **kwargs):

1075 from sympy.functions.special.gamma_functions import gamma

1076 return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1))

1077

1078 def _eval_rewrite_as_tractable(self, n, k, limitvar=None, **kwargs):

1079 return self._eval_rewrite_as_gamma(n, k).rewrite('tractable')

1080

1081 def _eval_rewrite_as_FallingFactorial(self, n, k, **kwargs):

1082 if k.is_integer:

1083 return ff(n, k) / factorial(k)

1084

1085 def _eval_is_integer(self):

1086 n, k = self.args

1087 if n.is_integer and k.is_integer:

1088 return True

1089 elif k.is_integer is False:

1090 return False

1091

1092 def _eval_is_nonnegative(self):

1093 n, k = self.args

1094 if n.is_integer and k.is_integer:

1095 if n.is_nonnegative or k.is_negative or k.is_even:

1096 return True

1097 elif k.is_even is False:

1098 return False

1099

1100 def _eval_as_leading_term(self, x, logx=None, cdir=0):

1101 from sympy.functions.special.gamma_functions import gamma

1102 return self.rewrite(gamma)._eval_as_leading_term(x, logx=logx, cdir=cdir)