Coverage for /usr/lib/python3/dist-packages/sympy/simplify/gammasimp.py: 4%

1from sympy.core import Function, S, Mul, Pow, Add

2from sympy.core.sorting import ordered, default_sort_key

3from sympy.core.function import expand_func

4from sympy.core.symbol import Dummy

5from sympy.functions import gamma, sqrt, sin

6from sympy.polys import factor, cancel

7from sympy.utilities.iterables import sift, uniq

10def gammasimp(expr):

11 r"""

12 Simplify expressions with gamma functions.

14 Explanation

15 ===========

17 This function takes as input an expression containing gamma

18 functions or functions that can be rewritten in terms of gamma

19 functions and tries to minimize the number of those functions and

20 reduce the size of their arguments.

22 The algorithm works by rewriting all gamma functions as expressions

23 involving rising factorials (Pochhammer symbols) and applies

24 recurrence relations and other transformations applicable to rising

25 factorials, to reduce their arguments, possibly letting the resulting

26 rising factorial to cancel. Rising factorials with the second argument

27 being an integer are expanded into polynomial forms and finally all

28 other rising factorial are rewritten in terms of gamma functions.

30 Then the following two steps are performed.

32 1. Reduce the number of gammas by applying the reflection theorem

33 gamma(x)*gamma(1-x) == pi/sin(pi*x).

34 2. Reduce the number of gammas by applying the multiplication theorem

35 gamma(x)*gamma(x+1/n)*...*gamma(x+(n-1)/n) == C*gamma(n*x).

37 It then reduces the number of prefactors by absorbing them into gammas

38 where possible and expands gammas with rational argument.

40 All transformation rules can be found (or were derived from) here:

42 .. [1] https://functions.wolfram.com/GammaBetaErf/Pochhammer/17/01/02/

43 .. [2] https://functions.wolfram.com/GammaBetaErf/Pochhammer/27/01/0005/

45 Examples

46 ========

48 >>> from sympy.simplify import gammasimp

49 >>> from sympy import gamma, Symbol

50 >>> from sympy.abc import x

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

53 >>> gammasimp(gamma(x)/gamma(x - 3))

54 (x - 3)*(x - 2)*(x - 1)

55 >>> gammasimp(gamma(n + 3))

56 gamma(n + 3)

58 """

60 expr = expr.rewrite(gamma)

62 # compute_ST will be looking for Functions and we don't want

63 # it looking for non-gamma functions: issue 22606

64 # so we mask free, non-gamma functions

65 f = expr.atoms(Function)

66 # take out gammas

67 gammas = {i for i in f if isinstance(i, gamma)}

68 if not gammas:

69 return expr # avoid side effects like factoring

70 f -= gammas

71 # keep only those without bound symbols

72 f = f & expr.as_dummy().atoms(Function)

73 if f:

74 dum, fun, simp = zip(*[

75 (Dummy(), fi, fi.func(*[

76 _gammasimp(a, as_comb=False) for a in fi.args]))

77 for fi in ordered(f)])

78 d = expr.xreplace(dict(zip(fun, dum)))

79 return _gammasimp(d, as_comb=False).xreplace(dict(zip(dum, simp)))

81 return _gammasimp(expr, as_comb=False)

84def _gammasimp(expr, as_comb):

85 """

86 Helper function for gammasimp and combsimp.

88 Explanation

89 ===========

91 Simplifies expressions written in terms of gamma function. If

92 as_comb is True, it tries to preserve integer arguments. See

93 docstring of gammasimp for more information. This was part of

94 combsimp() in combsimp.py.

95 """

96 expr = expr.replace(gamma,

97 lambda n: _rf(1, (n - 1).expand()))

99 if as_comb:

100 expr = expr.replace(_rf,

101 lambda a, b: gamma(b + 1))

102 else:

103 expr = expr.replace(_rf,

104 lambda a, b: gamma(a + b)/gamma(a))

105

106 def rule_gamma(expr, level=0):

107 """ Simplify products of gamma functions further. """

108

109 if expr.is_Atom:

110 return expr

111

112 def gamma_rat(x):

113 # helper to simplify ratios of gammas

114 was = x.count(gamma)

115 xx = x.replace(gamma, lambda n: _rf(1, (n - 1).expand()

116 ).replace(_rf, lambda a, b: gamma(a + b)/gamma(a)))

117 if xx.count(gamma) < was:

118 x = xx

119 return x

120

121 def gamma_factor(x):

122 # return True if there is a gamma factor in shallow args

123 if isinstance(x, gamma):

124 return True

125 if x.is_Add or x.is_Mul:

126 return any(gamma_factor(xi) for xi in x.args)

127 if x.is_Pow and (x.exp.is_integer or x.base.is_positive):

128 return gamma_factor(x.base)

129 return False

130

131 # recursion step

132 if level == 0:

133 expr = expr.func(*[rule_gamma(x, level + 1) for x in expr.args])

134 level += 1

135

136 if not expr.is_Mul:

137 return expr

138

139 # non-commutative step

140 if level == 1:

141 args, nc = expr.args_cnc()

142 if not args:

143 return expr

144 if nc:

145 return rule_gamma(Mul._from_args(args), level + 1)*Mul._from_args(nc)

146 level += 1

147

148 # pure gamma handling, not factor absorption

149 if level == 2:

150 T, F = sift(expr.args, gamma_factor, binary=True)

151 gamma_ind = Mul(*F)

152 d = Mul(*T)

153

154 nd, dd = d.as_numer_denom()

155 for ipass in range(2):

156 args = list(ordered(Mul.make_args(nd)))

157 for i, ni in enumerate(args):

158 if ni.is_Add:

159 ni, dd = Add(*[

160 rule_gamma(gamma_rat(a/dd), level + 1) for a in ni.args]

161 ).as_numer_denom()

162 args[i] = ni

163 if not dd.has(gamma):

164 break

165 nd = Mul(*args)

166 if ipass == 0 and not gamma_factor(nd):

167 break

168 nd, dd = dd, nd # now process in reversed order

169 expr = gamma_ind*nd/dd

170 if not (expr.is_Mul and (gamma_factor(dd) or gamma_factor(nd))):

171 return expr

172 level += 1

173

174 # iteration until constant

175 if level == 3:

176 while True:

177 was = expr

178 expr = rule_gamma(expr, 4)

179 if expr == was:

180 return expr

181

182 numer_gammas = []

183 denom_gammas = []

184 numer_others = []

185 denom_others = []

186 def explicate(p):

187 if p is S.One:

188 return None, []

189 b, e = p.as_base_exp()

190 if e.is_Integer:

191 if isinstance(b, gamma):

192 return True, [b.args[0]]*e

193 else:

194 return False, [b]*e

195 else:

196 return False, [p]

197

198 newargs = list(ordered(expr.args))

199 while newargs:

200 n, d = newargs.pop().as_numer_denom()

201 isg, l = explicate(n)

202 if isg:

203 numer_gammas.extend(l)

204 elif isg is False:

205 numer_others.extend(l)

206 isg, l = explicate(d)

207 if isg:

208 denom_gammas.extend(l)

209 elif isg is False:

210 denom_others.extend(l)

211

212 # =========== level 2 work: pure gamma manipulation =========

213

214 if not as_comb:

215 # Try to reduce the number of gamma factors by applying the

216 # reflection formula gamma(x)*gamma(1-x) = pi/sin(pi*x)

217 for gammas, numer, denom in [(

218 numer_gammas, numer_others, denom_others),

219 (denom_gammas, denom_others, numer_others)]:

220 new = []

221 while gammas:

222 g1 = gammas.pop()

223 if g1.is_integer:

224 new.append(g1)

225 continue

226 for i, g2 in enumerate(gammas):

227 n = g1 + g2 - 1

228 if not n.is_Integer:

229 continue

230 numer.append(S.Pi)

231 denom.append(sin(S.Pi*g1))

232 gammas.pop(i)

233 if n > 0:

234 for k in range(n):

235 numer.append(1 - g1 + k)

236 elif n < 0:

237 for k in range(-n):

238 denom.append(-g1 - k)

239 break

240 else:

241 new.append(g1)

242 # /!\ updating IN PLACE

243 gammas[:] = new

244

245 # Try to reduce the number of gammas by using the duplication

246 # theorem to cancel an upper and lower: gamma(2*s)/gamma(s) =

247 # 2**(2*s + 1)/(4*sqrt(pi))*gamma(s + 1/2). Although this could

248 # be done with higher argument ratios like gamma(3*x)/gamma(x),

249 # this would not reduce the number of gammas as in this case.

250 for ng, dg, no, do in [(numer_gammas, denom_gammas, numer_others,

251 denom_others),

252 (denom_gammas, numer_gammas, denom_others,

253 numer_others)]:

254

255 while True:

256 for x in ng:

257 for y in dg:

258 n = x - 2*y

259 if n.is_Integer:

260 break

261 else:

262 continue

263 break

264 else:

265 break

266 ng.remove(x)

267 dg.remove(y)

268 if n > 0:

269 for k in range(n):

270 no.append(2*y + k)

271 elif n < 0:

272 for k in range(-n):

273 do.append(2*y - 1 - k)

274 ng.append(y + S.Half)

275 no.append(2**(2*y - 1))

276 do.append(sqrt(S.Pi))

277

278 # Try to reduce the number of gamma factors by applying the

279 # multiplication theorem (used when n gammas with args differing

280 # by 1/n mod 1 are encountered).

281 #

282 # run of 2 with args differing by 1/2

283 #

284 # >>> gammasimp(gamma(x)*gamma(x+S.Half))

285 # 2*sqrt(2)*2**(-2*x - 1/2)*sqrt(pi)*gamma(2*x)

286 #

287 # run of 3 args differing by 1/3 (mod 1)

288 #

289 # >>> gammasimp(gamma(x)*gamma(x+S(1)/3)*gamma(x+S(2)/3))

290 # 6*3**(-3*x - 1/2)*pi*gamma(3*x)

291 # >>> gammasimp(gamma(x)*gamma(x+S(1)/3)*gamma(x+S(5)/3))

292 # 2*3**(-3*x - 1/2)*pi*(3*x + 2)*gamma(3*x)

293 #

294 def _run(coeffs):

295 # find runs in coeffs such that the difference in terms (mod 1)

296 # of t1, t2, ..., tn is 1/n

297 u = list(uniq(coeffs))

298 for i in range(len(u)):

299 dj = ([((u[j] - u[i]) % 1, j) for j in range(i + 1, len(u))])

300 for one, j in dj:

301 if one.p == 1 and one.q != 1:

302 n = one.q

303 got = [i]

304 get = list(range(1, n))

305 for d, j in dj:

306 m = n*d

307 if m.is_Integer and m in get:

308 get.remove(m)

309 got.append(j)

310 if not get:

311 break

312 else:

313 continue

314 for i, j in enumerate(got):

315 c = u[j]

316 coeffs.remove(c)

317 got[i] = c

318 return one.q, got[0], got[1:]

319

320 def _mult_thm(gammas, numer, denom):

321 # pull off and analyze the leading coefficient from each gamma arg

322 # looking for runs in those Rationals

323

324 # expr -> coeff + resid -> rats[resid] = coeff

325 rats = {}

326 for g in gammas:

327 c, resid = g.as_coeff_Add()

328 rats.setdefault(resid, []).append(c)

329

330 # look for runs in Rationals for each resid

331 keys = sorted(rats, key=default_sort_key)

332 for resid in keys:

333 coeffs = sorted(rats[resid])

334 new = []

335 while True:

336 run = _run(coeffs)

337 if run is None:

338 break

339

340 # process the sequence that was found:

341 # 1) convert all the gamma functions to have the right

342 # argument (could be off by an integer)

343 # 2) append the factors corresponding to the theorem

344 # 3) append the new gamma function

345

346 n, ui, other = run

347

348 # (1)

349 for u in other:

350 con = resid + u - 1

351 for k in range(int(u - ui)):

352 numer.append(con - k)

353

354 con = n*(resid + ui) # for (2) and (3)

355

356 # (2)

357 numer.append((2*S.Pi)**(S(n - 1)/2)*

358 n**(S.Half - con))

359 # (3)

360 new.append(con)

361

362 # restore resid to coeffs

363 rats[resid] = [resid + c for c in coeffs] + new

364

365 # rebuild the gamma arguments

366 g = []

367 for resid in keys:

368 g += rats[resid]

369 # /!\ updating IN PLACE

370 gammas[:] = g

371

372 for l, numer, denom in [(numer_gammas, numer_others, denom_others),

373 (denom_gammas, denom_others, numer_others)]:

374 _mult_thm(l, numer, denom)

375

376 # =========== level >= 2 work: factor absorption =========

377

378 if level >= 2:

379 # Try to absorb factors into the gammas: x*gamma(x) -> gamma(x + 1)

380 # and gamma(x)/(x - 1) -> gamma(x - 1)

381 # This code (in particular repeated calls to find_fuzzy) can be very

382 # slow.

383 def find_fuzzy(l, x):

384 if not l:

385 return

386 S1, T1 = compute_ST(x)

387 for y in l:

388 S2, T2 = inv[y]

389 if T1 != T2 or (not S1.intersection(S2) and

390 (S1 != set() or S2 != set())):

391 continue

392 # XXX we want some simplification (e.g. cancel or

393 # simplify) but no matter what it's slow.

394 a = len(cancel(x/y).free_symbols)

395 b = len(x.free_symbols)

396 c = len(y.free_symbols)

397 # TODO is there a better heuristic?

398 if a == 0 and (b > 0 or c > 0):

399 return y

400

401 # We thus try to avoid expensive calls by building the following

402 # "invariants": For every factor or gamma function argument

403 # - the set of free symbols S

404 # - the set of functional components T

405 # We will only try to absorb if T1==T2 and (S1 intersect S2 != emptyset

406 # or S1 == S2 == emptyset)

407 inv = {}

408

409 def compute_ST(expr):

410 if expr in inv:

411 return inv[expr]

412 return (expr.free_symbols, expr.atoms(Function).union(

413 {e.exp for e in expr.atoms(Pow)}))

414

415 def update_ST(expr):

416 inv[expr] = compute_ST(expr)

417 for expr in numer_gammas + denom_gammas + numer_others + denom_others:

418 update_ST(expr)

419

420 for gammas, numer, denom in [(

421 numer_gammas, numer_others, denom_others),

422 (denom_gammas, denom_others, numer_others)]:

423 new = []

424 while gammas:

425 g = gammas.pop()

426 cont = True

427 while cont:

428 cont = False

429 y = find_fuzzy(numer, g)

430 if y is not None:

431 numer.remove(y)

432 if y != g:

433 numer.append(y/g)

434 update_ST(y/g)

435 g += 1

436 cont = True

437 y = find_fuzzy(denom, g - 1)

438 if y is not None:

439 denom.remove(y)

440 if y != g - 1:

441 numer.append((g - 1)/y)

442 update_ST((g - 1)/y)

443 g -= 1

444 cont = True

445 new.append(g)

446 # /!\ updating IN PLACE

447 gammas[:] = new

448

449 # =========== rebuild expr ==================================

450

451 return Mul(*[gamma(g) for g in numer_gammas]) \

452 / Mul(*[gamma(g) for g in denom_gammas]) \

453 * Mul(*numer_others) / Mul(*denom_others)

454

455 was = factor(expr)

456 # (for some reason we cannot use Basic.replace in this case)

457 expr = rule_gamma(was)

458 if expr != was:

459 expr = factor(expr)

460

461 expr = expr.replace(gamma,

462 lambda n: expand_func(gamma(n)) if n.is_Rational else gamma(n))

463

464 return expr

465

466

467class _rf(Function):

468 @classmethod

469 def eval(cls, a, b):

470 if b.is_Integer:

471 if not b:

472 return S.One

473

474 n = int(b)

475

476 if n > 0:

477 return Mul(*[a + i for i in range(n)])

478 elif n < 0:

479 return 1/Mul(*[a - i for i in range(1, -n + 1)])

480 else:

481 if b.is_Add:

482 c, _b = b.as_coeff_Add()

483

484 if c.is_Integer:

485 if c > 0:

486 return _rf(a, _b)*_rf(a + _b, c)

487 elif c < 0:

488 return _rf(a, _b)/_rf(a + _b + c, -c)

489

490 if a.is_Add:

491 c, _a = a.as_coeff_Add()

492

493 if c.is_Integer:

494 if c > 0:

495 return _rf(_a, b)*_rf(_a + b, c)/_rf(_a, c)

496 elif c < 0:

497 return _rf(_a, b)*_rf(_a + c, -c)/_rf(_a + b + c, -c)