Coverage for /usr/lib/python3/dist-packages/sympy/concrete/products.py: 18%

1from typing import Tuple as tTuple

3from .expr_with_intlimits import ExprWithIntLimits

4from .summations import Sum, summation, _dummy_with_inherited_properties_concrete

5from sympy.core.expr import Expr

6from sympy.core.exprtools import factor_terms

7from sympy.core.function import Derivative

8from sympy.core.mul import Mul

9from sympy.core.singleton import S

10from sympy.core.symbol import Dummy, Symbol

11from sympy.functions.combinatorial.factorials import RisingFactorial

12from sympy.functions.elementary.exponential import exp, log

13from sympy.functions.special.tensor_functions import KroneckerDelta

14from sympy.polys import quo, roots

17class Product(ExprWithIntLimits):

18 r"""

19 Represents unevaluated products.

21 Explanation

22 ===========

24 ``Product`` represents a finite or infinite product, with the first

25 argument being the general form of terms in the series, and the second

26 argument being ``(dummy_variable, start, end)``, with ``dummy_variable``

27 taking all integer values from ``start`` through ``end``. In accordance

28 with long-standing mathematical convention, the end term is included in

29 the product.

31 Finite products

32 ===============

34 For finite products (and products with symbolic limits assumed to be finite)

35 we follow the analogue of the summation convention described by Karr [1],

36 especially definition 3 of section 1.4. The product:

38 .. math::

40 \prod_{m \leq i < n} f(i)

42 has *the obvious meaning* for `m < n`, namely:

44 .. math::

46 \prod_{m \leq i < n} f(i) = f(m) f(m+1) \cdot \ldots \cdot f(n-2) f(n-1)

48 with the upper limit value `f(n)` excluded. The product over an empty set is

49 one if and only if `m = n`:

51 .. math::

53 \prod_{m \leq i < n} f(i) = 1 \quad \mathrm{for} \quad m = n

55 Finally, for all other products over empty sets we assume the following

56 definition:

58 .. math::

60 \prod_{m \leq i < n} f(i) = \frac{1}{\prod_{n \leq i < m} f(i)} \quad \mathrm{for} \quad m > n

62 It is important to note that above we define all products with the upper

63 limit being exclusive. This is in contrast to the usual mathematical notation,

64 but does not affect the product convention. Indeed we have:

66 .. math::

68 \prod_{m \leq i < n} f(i) = \prod_{i = m}^{n - 1} f(i)

70 where the difference in notation is intentional to emphasize the meaning,

71 with limits typeset on the top being inclusive.

73 Examples

74 ========

76 >>> from sympy.abc import a, b, i, k, m, n, x

77 >>> from sympy import Product, oo

78 >>> Product(k, (k, 1, m))

79 Product(k, (k, 1, m))

80 >>> Product(k, (k, 1, m)).doit()

81 factorial(m)

82 >>> Product(k**2,(k, 1, m))

83 Product(k**2, (k, 1, m))

84 >>> Product(k**2,(k, 1, m)).doit()

85 factorial(m)**2

87 Wallis' product for pi:

89 >>> W = Product(2*i/(2*i-1) * 2*i/(2*i+1), (i, 1, oo))

90 >>> W

91 Product(4*i**2/((2*i - 1)*(2*i + 1)), (i, 1, oo))

93 Direct computation currently fails:

95 >>> W.doit()

96 Product(4*i**2/((2*i - 1)*(2*i + 1)), (i, 1, oo))

98 But we can approach the infinite product by a limit of finite products:

100 >>> from sympy import limit

101 >>> W2 = Product(2*i/(2*i-1)*2*i/(2*i+1), (i, 1, n))

102 >>> W2

103 Product(4*i**2/((2*i - 1)*(2*i + 1)), (i, 1, n))

104 >>> W2e = W2.doit()

105 >>> W2e

106 4**n*factorial(n)**2/(2**(2*n)*RisingFactorial(1/2, n)*RisingFactorial(3/2, n))

107 >>> limit(W2e, n, oo)

108 pi/2

109

110 By the same formula we can compute sin(pi/2):

111

112 >>> from sympy import combsimp, pi, gamma, simplify

113 >>> P = pi * x * Product(1 - x**2/k**2, (k, 1, n))

114 >>> P = P.subs(x, pi/2)

115 >>> P

116 pi**2*Product(1 - pi**2/(4*k**2), (k, 1, n))/2

117 >>> Pe = P.doit()

118 >>> Pe

119 pi**2*RisingFactorial(1 - pi/2, n)*RisingFactorial(1 + pi/2, n)/(2*factorial(n)**2)

120 >>> limit(Pe, n, oo).gammasimp()

121 sin(pi**2/2)

122 >>> Pe.rewrite(gamma)

123 (-1)**n*pi**2*gamma(pi/2)*gamma(n + 1 + pi/2)/(2*gamma(1 + pi/2)*gamma(-n + pi/2)*gamma(n + 1)**2)

124

125 Products with the lower limit being larger than the upper one:

126

127 >>> Product(1/i, (i, 6, 1)).doit()

128 120

129 >>> Product(i, (i, 2, 5)).doit()

130 120

131

132 The empty product:

133

134 >>> Product(i, (i, n, n-1)).doit()

135 1

136

137 An example showing that the symbolic result of a product is still

138 valid for seemingly nonsensical values of the limits. Then the Karr

139 convention allows us to give a perfectly valid interpretation to

140 those products by interchanging the limits according to the above rules:

141

142 >>> P = Product(2, (i, 10, n)).doit()

143 >>> P

144 2**(n - 9)

145 >>> P.subs(n, 5)

146 1/16

147 >>> Product(2, (i, 10, 5)).doit()

148 1/16

149 >>> 1/Product(2, (i, 6, 9)).doit()

150 1/16

151

152 An explicit example of the Karr summation convention applied to products:

153

154 >>> P1 = Product(x, (i, a, b)).doit()

155 >>> P1

156 x**(-a + b + 1)

157 >>> P2 = Product(x, (i, b+1, a-1)).doit()

158 >>> P2

159 x**(a - b - 1)

160 >>> simplify(P1 * P2)

161 1

162

163 And another one:

164

165 >>> P1 = Product(i, (i, b, a)).doit()

166 >>> P1

167 RisingFactorial(b, a - b + 1)

168 >>> P2 = Product(i, (i, a+1, b-1)).doit()

169 >>> P2

170 RisingFactorial(a + 1, -a + b - 1)

171 >>> P1 * P2

172 RisingFactorial(b, a - b + 1)*RisingFactorial(a + 1, -a + b - 1)

173 >>> combsimp(P1 * P2)

174 1

175

176 See Also

177 ========

178

179 Sum, summation

180 product

181

182 References

183 ==========

184

185 .. [1] Michael Karr, "Summation in Finite Terms", Journal of the ACM,

186 Volume 28 Issue 2, April 1981, Pages 305-350

187 https://dl.acm.org/doi/10.1145/322248.322255

188 .. [2] https://en.wikipedia.org/wiki/Multiplication#Capital_Pi_notation

189 .. [3] https://en.wikipedia.org/wiki/Empty_product

190 """

191

192 __slots__ = ()

193

194 limits: tTuple[tTuple[Symbol, Expr, Expr]]

195

196 def __new__(cls, function, *symbols, **assumptions):

197 obj = ExprWithIntLimits.__new__(cls, function, *symbols, **assumptions)

198 return obj

199

200 def _eval_rewrite_as_Sum(self, *args, **kwargs):

201 return exp(Sum(log(self.function), *self.limits))

202

203 @property

204 def term(self):

205 return self._args[0]

206 function = term

207

208 def _eval_is_zero(self):

209 if self.has_empty_sequence:

210 return False

211

212 z = self.term.is_zero

213 if z is True:

214 return True

215 if self.has_finite_limits:

216 # A Product is zero only if its term is zero assuming finite limits.

217 return z

218

219 def _eval_is_extended_real(self):

220 if self.has_empty_sequence:

221 return True

222

223 return self.function.is_extended_real

224

225 def _eval_is_positive(self):

226 if self.has_empty_sequence:

227 return True

228 if self.function.is_positive and self.has_finite_limits:

229 return True

230

231 def _eval_is_nonnegative(self):

232 if self.has_empty_sequence:

233 return True

234 if self.function.is_nonnegative and self.has_finite_limits:

235 return True

236

237 def _eval_is_extended_nonnegative(self):

238 if self.has_empty_sequence:

239 return True

240 if self.function.is_extended_nonnegative:

241 return True

242

243 def _eval_is_extended_nonpositive(self):

244 if self.has_empty_sequence:

245 return True

246

247 def _eval_is_finite(self):

248 if self.has_finite_limits and self.function.is_finite:

249 return True

250

251 def doit(self, **hints):

252 # first make sure any definite limits have product

253 # variables with matching assumptions

254 reps = {}

255 for xab in self.limits:

256 d = _dummy_with_inherited_properties_concrete(xab)

257 if d:

258 reps[xab[0]] = d

259 if reps:

260 undo = {v: k for k, v in reps.items()}

261 did = self.xreplace(reps).doit(**hints)

262 if isinstance(did, tuple): # when separate=True

263 did = tuple([i.xreplace(undo) for i in did])

264 else:

265 did = did.xreplace(undo)

266 return did

267

268 from sympy.simplify.powsimp import powsimp

269 f = self.function

270 for index, limit in enumerate(self.limits):

271 i, a, b = limit

272 dif = b - a

273 if dif.is_integer and dif.is_negative:

274 a, b = b + 1, a - 1

275 f = 1 / f

276

277 g = self._eval_product(f, (i, a, b))

278 if g in (None, S.NaN):

279 return self.func(powsimp(f), *self.limits[index:])

280 else:

281 f = g

282

283 if hints.get('deep', True):

284 return f.doit(**hints)

285 else:

286 return powsimp(f)

287

288 def _eval_adjoint(self):

289 if self.is_commutative:

290 return self.func(self.function.adjoint(), *self.limits)

291 return None

292

293 def _eval_conjugate(self):

294 return self.func(self.function.conjugate(), *self.limits)

295

296 def _eval_product(self, term, limits):

297

298 (k, a, n) = limits

299

300 if k not in term.free_symbols:

301 if (term - 1).is_zero:

302 return S.One

303 return term**(n - a + 1)

304

305 if a == n:

306 return term.subs(k, a)

307

308 from .delta import deltaproduct, _has_simple_delta

309 if term.has(KroneckerDelta) and _has_simple_delta(term, limits[0]):

310 return deltaproduct(term, limits)

311

312 dif = n - a

313 definite = dif.is_Integer

314 if definite and (dif < 100):

315 return self._eval_product_direct(term, limits)

316

317 elif term.is_polynomial(k):

318 poly = term.as_poly(k)

319

320 A = B = Q = S.One

321

322 all_roots = roots(poly)

323

324 M = 0

325 for r, m in all_roots.items():

326 M += m

327 A *= RisingFactorial(a - r, n - a + 1)**m

328 Q *= (n - r)**m

329

330 if M < poly.degree():

331 arg = quo(poly, Q.as_poly(k))

332 B = self.func(arg, (k, a, n)).doit()

333

334 return poly.LC()**(n - a + 1) * A * B

335

336 elif term.is_Add:

337 factored = factor_terms(term, fraction=True)

338 if factored.is_Mul:

339 return self._eval_product(factored, (k, a, n))

340

341 elif term.is_Mul:

342 # Factor in part without the summation variable and part with

343 without_k, with_k = term.as_coeff_mul(k)

344

345 if len(with_k) >= 2:

346 # More than one term including k, so still a multiplication

347 exclude, include = [], []

348 for t in with_k:

349 p = self._eval_product(t, (k, a, n))

350

351 if p is not None:

352 exclude.append(p)

353 else:

354 include.append(t)

355

356 if not exclude:

357 return None

358 else:

359 arg = term._new_rawargs(*include)

360 A = Mul(*exclude)

361 B = self.func(arg, (k, a, n)).doit()

362 return without_k**(n - a + 1)*A * B

363 else:

364 # Just a single term

365 p = self._eval_product(with_k[0], (k, a, n))

366 if p is None:

367 p = self.func(with_k[0], (k, a, n)).doit()

368 return without_k**(n - a + 1)*p

369

370

371 elif term.is_Pow:

372 if not term.base.has(k):

373 s = summation(term.exp, (k, a, n))

374

375 return term.base**s

376 elif not term.exp.has(k):

377 p = self._eval_product(term.base, (k, a, n))

378

379 if p is not None:

380 return p**term.exp

381

382 elif isinstance(term, Product):

383 evaluated = term.doit()

384 f = self._eval_product(evaluated, limits)

385 if f is None:

386 return self.func(evaluated, limits)

387 else:

388 return f

389

390 if definite:

391 return self._eval_product_direct(term, limits)

392

393 def _eval_simplify(self, **kwargs):

394 from sympy.simplify.simplify import product_simplify

395 rv = product_simplify(self, **kwargs)

396 return rv.doit() if kwargs['doit'] else rv

397

398 def _eval_transpose(self):

399 if self.is_commutative:

400 return self.func(self.function.transpose(), *self.limits)

401 return None

402

403 def _eval_product_direct(self, term, limits):

404 (k, a, n) = limits

405 return Mul(*[term.subs(k, a + i) for i in range(n - a + 1)])

406

407 def _eval_derivative(self, x):

408 if isinstance(x, Symbol) and x not in self.free_symbols:

409 return S.Zero

410 f, limits = self.function, list(self.limits)

411 limit = limits.pop(-1)

412 if limits:

413 f = self.func(f, *limits)

414 i, a, b = limit

415 if x in a.free_symbols or x in b.free_symbols:

416 return None

417 h = Dummy()

418 rv = Sum( Product(f, (i, a, h - 1)) * Product(f, (i, h + 1, b)) * Derivative(f, x, evaluate=True).subs(i, h), (h, a, b))

419 return rv

420

421 def is_convergent(self):

422 r"""

423 See docs of :obj:`.Sum.is_convergent()` for explanation of convergence

424 in SymPy.

425

426 Explanation

427 ===========

428

429 The infinite product:

430

431 .. math::

432

433 \prod_{1 \leq i < \infty} f(i)

434

435 is defined by the sequence of partial products:

436

437 .. math::

438

439 \prod_{i=1}^{n} f(i) = f(1) f(2) \cdots f(n)

440

441 as n increases without bound. The product converges to a non-zero

442 value if and only if the sum:

443

444 .. math::

445

446 \sum_{1 \leq i < \infty} \log{f(n)}

447

448 converges.

449

450 Examples

451 ========

452

453 >>> from sympy import Product, Symbol, cos, pi, exp, oo

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

455 >>> Product(n/(n + 1), (n, 1, oo)).is_convergent()

456 False

457 >>> Product(1/n**2, (n, 1, oo)).is_convergent()

458 False

459 >>> Product(cos(pi/n), (n, 1, oo)).is_convergent()

460 True

461 >>> Product(exp(-n**2), (n, 1, oo)).is_convergent()

462 False

463

464 References

465 ==========

466

467 .. [1] https://en.wikipedia.org/wiki/Infinite_product

468 """

469 sequence_term = self.function

470 log_sum = log(sequence_term)

471 lim = self.limits

472 try:

473 is_conv = Sum(log_sum, *lim).is_convergent()

474 except NotImplementedError:

475 if Sum(sequence_term - 1, *lim).is_absolutely_convergent() is S.true:

476 return S.true

477 raise NotImplementedError("The algorithm to find the product convergence of %s "

478 "is not yet implemented" % (sequence_term))

479 return is_conv

480

481 def reverse_order(expr, *indices):

482 """

483 Reverse the order of a limit in a Product.

484

485 Explanation

486 ===========

487

488 ``reverse_order(expr, *indices)`` reverses some limits in the expression

489 ``expr`` which can be either a ``Sum`` or a ``Product``. The selectors in

490 the argument ``indices`` specify some indices whose limits get reversed.

491 These selectors are either variable names or numerical indices counted

492 starting from the inner-most limit tuple.

493

494 Examples

495 ========

496

497 >>> from sympy import gamma, Product, simplify, Sum

498 >>> from sympy.abc import x, y, a, b, c, d

499 >>> P = Product(x, (x, a, b))

500 >>> Pr = P.reverse_order(x)

501 >>> Pr

502 Product(1/x, (x, b + 1, a - 1))

503 >>> Pr = Pr.doit()

504 >>> Pr

505 1/RisingFactorial(b + 1, a - b - 1)

506 >>> simplify(Pr.rewrite(gamma))

507 Piecewise((gamma(b + 1)/gamma(a), b > -1), ((-1)**(-a + b + 1)*gamma(1 - a)/gamma(-b), True))

508 >>> P = P.doit()

509 >>> P

510 RisingFactorial(a, -a + b + 1)

511 >>> simplify(P.rewrite(gamma))

512 Piecewise((gamma(b + 1)/gamma(a), a > 0), ((-1)**(-a + b + 1)*gamma(1 - a)/gamma(-b), True))

513

514 While one should prefer variable names when specifying which limits

515 to reverse, the index counting notation comes in handy in case there

516 are several symbols with the same name.

517

518 >>> S = Sum(x*y, (x, a, b), (y, c, d))

519 >>> S

520 Sum(x*y, (x, a, b), (y, c, d))

521 >>> S0 = S.reverse_order(0)

522 >>> S0

523 Sum(-x*y, (x, b + 1, a - 1), (y, c, d))

524 >>> S1 = S0.reverse_order(1)

525 >>> S1

526 Sum(x*y, (x, b + 1, a - 1), (y, d + 1, c - 1))

527

528 Of course we can mix both notations:

529

530 >>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(x, 1)

531 Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))

532 >>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(y, x)

533 Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))

534

535 See Also

536 ========

537

538 sympy.concrete.expr_with_intlimits.ExprWithIntLimits.index,

539 reorder_limit,

540 sympy.concrete.expr_with_intlimits.ExprWithIntLimits.reorder

541

542 References

543 ==========

544

545 .. [1] Michael Karr, "Summation in Finite Terms", Journal of the ACM,

546 Volume 28 Issue 2, April 1981, Pages 305-350

547 https://dl.acm.org/doi/10.1145/322248.322255

548

549 """

550 l_indices = list(indices)

551

552 for i, indx in enumerate(l_indices):

553 if not isinstance(indx, int):

554 l_indices[i] = expr.index(indx)

555

556 e = 1

557 limits = []

558 for i, limit in enumerate(expr.limits):

559 l = limit

560 if i in l_indices:

561 e = -e

562 l = (limit[0], limit[2] + 1, limit[1] - 1)

563 limits.append(l)

564

565 return Product(expr.function ** e, *limits)

566

567

568def product(*args, **kwargs):

569 r"""

570 Compute the product.

571

572 Explanation

573 ===========

574

575 The notation for symbols is similar to the notation used in Sum or

576 Integral. product(f, (i, a, b)) computes the product of f with

577 respect to i from a to b, i.e.,

578

579 ::

580

581 b

582 _____

583 product(f(n), (i, a, b)) = | | f(n)

584 | |

585 i = a

586

587 If it cannot compute the product, it returns an unevaluated Product object.

588 Repeated products can be computed by introducing additional symbols tuples::

589

590 Examples

591 ========

592

593 >>> from sympy import product, symbols

594 >>> i, n, m, k = symbols('i n m k', integer=True)

595

596 >>> product(i, (i, 1, k))

597 factorial(k)

598 >>> product(m, (i, 1, k))

599 m**k

600 >>> product(i, (i, 1, k), (k, 1, n))

601 Product(factorial(k), (k, 1, n))

602

603 """

604

605 prod = Product(*args, **kwargs)

606

607 if isinstance(prod, Product):

608 return prod.doit(deep=False)

609 else:

610 return prod