qfunctions.py 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
  1. from .functions import defun, defun_wrapped
  2. @defun
  3. def qp(ctx, a, q=None, n=None, **kwargs):
  4. r"""
  5. Evaluates the q-Pochhammer symbol (or q-rising factorial)
  6. .. math ::
  7. (a; q)_n = \prod_{k=0}^{n-1} (1-a q^k)
  8. where `n = \infty` is permitted if `|q| < 1`. Called with two arguments,
  9. ``qp(a,q)`` computes `(a;q)_{\infty}`; with a single argument, ``qp(q)``
  10. computes `(q;q)_{\infty}`. The special case
  11. .. math ::
  12. \phi(q) = (q; q)_{\infty} = \prod_{k=1}^{\infty} (1-q^k) =
  13. \sum_{k=-\infty}^{\infty} (-1)^k q^{(3k^2-k)/2}
  14. is also known as the Euler function, or (up to a factor `q^{-1/24}`)
  15. the Dedekind eta function.
  16. **Examples**
  17. If `n` is a positive integer, the function amounts to a finite product::
  18. >>> from mpmath import *
  19. >>> mp.dps = 25; mp.pretty = True
  20. >>> qp(2,3,5)
  21. -725305.0
  22. >>> fprod(1-2*3**k for k in range(5))
  23. -725305.0
  24. >>> qp(2,3,0)
  25. 1.0
  26. Complex arguments are allowed::
  27. >>> qp(2-1j, 0.75j)
  28. (0.4628842231660149089976379 + 4.481821753552703090628793j)
  29. The regular Pochhammer symbol `(a)_n` is obtained in the
  30. following limit as `q \to 1`::
  31. >>> a, n = 4, 7
  32. >>> limit(lambda q: qp(q**a,q,n) / (1-q)**n, 1)
  33. 604800.0
  34. >>> rf(a,n)
  35. 604800.0
  36. The Taylor series of the reciprocal Euler function gives
  37. the partition function `P(n)`, i.e. the number of ways of writing
  38. `n` as a sum of positive integers::
  39. >>> taylor(lambda q: 1/qp(q), 0, 10)
  40. [1.0, 1.0, 2.0, 3.0, 5.0, 7.0, 11.0, 15.0, 22.0, 30.0, 42.0]
  41. Special values include::
  42. >>> qp(0)
  43. 1.0
  44. >>> findroot(diffun(qp), -0.4) # location of maximum
  45. -0.4112484791779547734440257
  46. >>> qp(_)
  47. 1.228348867038575112586878
  48. The q-Pochhammer symbol is related to the Jacobi theta functions.
  49. For example, the following identity holds::
  50. >>> q = mpf(0.5) # arbitrary
  51. >>> qp(q)
  52. 0.2887880950866024212788997
  53. >>> root(3,-2)*root(q,-24)*jtheta(2,pi/6,root(q,6))
  54. 0.2887880950866024212788997
  55. """
  56. a = ctx.convert(a)
  57. if n is None:
  58. n = ctx.inf
  59. else:
  60. n = ctx.convert(n)
  61. if n < 0:
  62. raise ValueError("n cannot be negative")
  63. if q is None:
  64. q = a
  65. else:
  66. q = ctx.convert(q)
  67. if n == 0:
  68. return ctx.one + 0*(a+q)
  69. infinite = (n == ctx.inf)
  70. same = (a == q)
  71. if infinite:
  72. if abs(q) >= 1:
  73. if same and (q == -1 or q == 1):
  74. return ctx.zero * q
  75. raise ValueError("q-function only defined for |q| < 1")
  76. elif q == 0:
  77. return ctx.one - a
  78. maxterms = kwargs.get('maxterms', 50*ctx.prec)
  79. if infinite and same:
  80. # Euler's pentagonal theorem
  81. def terms():
  82. t = 1
  83. yield t
  84. k = 1
  85. x1 = q
  86. x2 = q**2
  87. while 1:
  88. yield (-1)**k * x1
  89. yield (-1)**k * x2
  90. x1 *= q**(3*k+1)
  91. x2 *= q**(3*k+2)
  92. k += 1
  93. if k > maxterms:
  94. raise ctx.NoConvergence
  95. return ctx.sum_accurately(terms)
  96. # return ctx.nprod(lambda k: 1-a*q**k, [0,n-1])
  97. def factors():
  98. k = 0
  99. r = ctx.one
  100. while 1:
  101. yield 1 - a*r
  102. r *= q
  103. k += 1
  104. if k >= n:
  105. return
  106. if k > maxterms:
  107. raise ctx.NoConvergence
  108. return ctx.mul_accurately(factors)
  109. @defun_wrapped
  110. def qgamma(ctx, z, q, **kwargs):
  111. r"""
  112. Evaluates the q-gamma function
  113. .. math ::
  114. \Gamma_q(z) = \frac{(q; q)_{\infty}}{(q^z; q)_{\infty}} (1-q)^{1-z}.
  115. **Examples**
  116. Evaluation for real and complex arguments::
  117. >>> from mpmath import *
  118. >>> mp.dps = 25; mp.pretty = True
  119. >>> qgamma(4,0.75)
  120. 4.046875
  121. >>> qgamma(6,6)
  122. 121226245.0
  123. >>> qgamma(3+4j, 0.5j)
  124. (0.1663082382255199834630088 + 0.01952474576025952984418217j)
  125. The q-gamma function satisfies a functional equation similar
  126. to that of the ordinary gamma function::
  127. >>> q = mpf(0.25)
  128. >>> z = mpf(2.5)
  129. >>> qgamma(z+1,q)
  130. 1.428277424823760954685912
  131. >>> (1-q**z)/(1-q)*qgamma(z,q)
  132. 1.428277424823760954685912
  133. """
  134. if abs(q) > 1:
  135. return ctx.qgamma(z,1/q)*q**((z-2)*(z-1)*0.5)
  136. return ctx.qp(q, q, None, **kwargs) / \
  137. ctx.qp(q**z, q, None, **kwargs) * (1-q)**(1-z)
  138. @defun_wrapped
  139. def qfac(ctx, z, q, **kwargs):
  140. r"""
  141. Evaluates the q-factorial,
  142. .. math ::
  143. [n]_q! = (1+q)(1+q+q^2)\cdots(1+q+\cdots+q^{n-1})
  144. or more generally
  145. .. math ::
  146. [z]_q! = \frac{(q;q)_z}{(1-q)^z}.
  147. **Examples**
  148. >>> from mpmath import *
  149. >>> mp.dps = 25; mp.pretty = True
  150. >>> qfac(0,0)
  151. 1.0
  152. >>> qfac(4,3)
  153. 2080.0
  154. >>> qfac(5,6)
  155. 121226245.0
  156. >>> qfac(1+1j, 2+1j)
  157. (0.4370556551322672478613695 + 0.2609739839216039203708921j)
  158. """
  159. if ctx.isint(z) and ctx._re(z) > 0:
  160. n = int(ctx._re(z))
  161. return ctx.qp(q, q, n, **kwargs) / (1-q)**n
  162. return ctx.qgamma(z+1, q, **kwargs)
  163. @defun
  164. def qhyper(ctx, a_s, b_s, q, z, **kwargs):
  165. r"""
  166. Evaluates the basic hypergeometric series or hypergeometric q-series
  167. .. math ::
  168. \,_r\phi_s \left[\begin{matrix}
  169. a_1 & a_2 & \ldots & a_r \\
  170. b_1 & b_2 & \ldots & b_s
  171. \end{matrix} ; q,z \right] =
  172. \sum_{n=0}^\infty
  173. \frac{(a_1;q)_n, \ldots, (a_r;q)_n}
  174. {(b_1;q)_n, \ldots, (b_s;q)_n}
  175. \left((-1)^n q^{n\choose 2}\right)^{1+s-r}
  176. \frac{z^n}{(q;q)_n}
  177. where `(a;q)_n` denotes the q-Pochhammer symbol (see :func:`~mpmath.qp`).
  178. **Examples**
  179. Evaluation works for real and complex arguments::
  180. >>> from mpmath import *
  181. >>> mp.dps = 25; mp.pretty = True
  182. >>> qhyper([0.5], [2.25], 0.25, 4)
  183. -0.1975849091263356009534385
  184. >>> qhyper([0.5], [2.25], 0.25-0.25j, 4)
  185. (2.806330244925716649839237 + 3.568997623337943121769938j)
  186. >>> qhyper([1+j], [2,3+0.5j], 0.25, 3+4j)
  187. (9.112885171773400017270226 - 1.272756997166375050700388j)
  188. Comparing with a summation of the defining series, using
  189. :func:`~mpmath.nsum`::
  190. >>> b, q, z = 3, 0.25, 0.5
  191. >>> qhyper([], [b], q, z)
  192. 0.6221136748254495583228324
  193. >>> nsum(lambda n: z**n / qp(q,q,n)/qp(b,q,n) * q**(n*(n-1)), [0,inf])
  194. 0.6221136748254495583228324
  195. """
  196. #a_s = [ctx._convert_param(a)[0] for a in a_s]
  197. #b_s = [ctx._convert_param(b)[0] for b in b_s]
  198. #q = ctx._convert_param(q)[0]
  199. a_s = [ctx.convert(a) for a in a_s]
  200. b_s = [ctx.convert(b) for b in b_s]
  201. q = ctx.convert(q)
  202. z = ctx.convert(z)
  203. r = len(a_s)
  204. s = len(b_s)
  205. d = 1+s-r
  206. maxterms = kwargs.get('maxterms', 50*ctx.prec)
  207. def terms():
  208. t = ctx.one
  209. yield t
  210. qk = 1
  211. k = 0
  212. x = 1
  213. while 1:
  214. for a in a_s:
  215. p = 1 - a*qk
  216. t *= p
  217. for b in b_s:
  218. p = 1 - b*qk
  219. if not p:
  220. raise ValueError
  221. t /= p
  222. t *= z
  223. x *= (-1)**d * qk ** d
  224. qk *= q
  225. t /= (1 - qk)
  226. k += 1
  227. yield t * x
  228. if k > maxterms:
  229. raise ctx.NoConvergence
  230. return ctx.sum_accurately(terms)