multinomial.py 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188
  1. from sympy.utilities.misc import as_int
  2. def binomial_coefficients(n):
  3. """Return a dictionary containing pairs :math:`{(k1,k2) : C_kn}` where
  4. :math:`C_kn` are binomial coefficients and :math:`n=k1+k2`.
  5. Examples
  6. ========
  7. >>> from sympy.ntheory import binomial_coefficients
  8. >>> binomial_coefficients(9)
  9. {(0, 9): 1, (1, 8): 9, (2, 7): 36, (3, 6): 84,
  10. (4, 5): 126, (5, 4): 126, (6, 3): 84, (7, 2): 36, (8, 1): 9, (9, 0): 1}
  11. See Also
  12. ========
  13. binomial_coefficients_list, multinomial_coefficients
  14. """
  15. n = as_int(n)
  16. d = {(0, n): 1, (n, 0): 1}
  17. a = 1
  18. for k in range(1, n//2 + 1):
  19. a = (a * (n - k + 1))//k
  20. d[k, n - k] = d[n - k, k] = a
  21. return d
  22. def binomial_coefficients_list(n):
  23. """ Return a list of binomial coefficients as rows of the Pascal's
  24. triangle.
  25. Examples
  26. ========
  27. >>> from sympy.ntheory import binomial_coefficients_list
  28. >>> binomial_coefficients_list(9)
  29. [1, 9, 36, 84, 126, 126, 84, 36, 9, 1]
  30. See Also
  31. ========
  32. binomial_coefficients, multinomial_coefficients
  33. """
  34. n = as_int(n)
  35. d = [1] * (n + 1)
  36. a = 1
  37. for k in range(1, n//2 + 1):
  38. a = (a * (n - k + 1))//k
  39. d[k] = d[n - k] = a
  40. return d
  41. def multinomial_coefficients(m, n):
  42. r"""Return a dictionary containing pairs ``{(k1,k2,..,km) : C_kn}``
  43. where ``C_kn`` are multinomial coefficients such that
  44. ``n=k1+k2+..+km``.
  45. Examples
  46. ========
  47. >>> from sympy.ntheory import multinomial_coefficients
  48. >>> multinomial_coefficients(2, 5) # indirect doctest
  49. {(0, 5): 1, (1, 4): 5, (2, 3): 10, (3, 2): 10, (4, 1): 5, (5, 0): 1}
  50. Notes
  51. =====
  52. The algorithm is based on the following result:
  53. .. math::
  54. \binom{n}{k_1, \ldots, k_m} =
  55. \frac{k_1 + 1}{n - k_1} \sum_{i=2}^m \binom{n}{k_1 + 1, \ldots, k_i - 1, \ldots}
  56. Code contributed to Sage by Yann Laigle-Chapuy, copied with permission
  57. of the author.
  58. See Also
  59. ========
  60. binomial_coefficients_list, binomial_coefficients
  61. """
  62. m = as_int(m)
  63. n = as_int(n)
  64. if not m:
  65. if n:
  66. return {}
  67. return {(): 1}
  68. if m == 2:
  69. return binomial_coefficients(n)
  70. if m >= 2*n and n > 1:
  71. return dict(multinomial_coefficients_iterator(m, n))
  72. t = [n] + [0] * (m - 1)
  73. r = {tuple(t): 1}
  74. if n:
  75. j = 0 # j will be the leftmost nonzero position
  76. else:
  77. j = m
  78. # enumerate tuples in co-lex order
  79. while j < m - 1:
  80. # compute next tuple
  81. tj = t[j]
  82. if j:
  83. t[j] = 0
  84. t[0] = tj
  85. if tj > 1:
  86. t[j + 1] += 1
  87. j = 0
  88. start = 1
  89. v = 0
  90. else:
  91. j += 1
  92. start = j + 1
  93. v = r[tuple(t)]
  94. t[j] += 1
  95. # compute the value
  96. # NB: the initialization of v was done above
  97. for k in range(start, m):
  98. if t[k]:
  99. t[k] -= 1
  100. v += r[tuple(t)]
  101. t[k] += 1
  102. t[0] -= 1
  103. r[tuple(t)] = (v * tj) // (n - t[0])
  104. return r
  105. def multinomial_coefficients_iterator(m, n, _tuple=tuple):
  106. """multinomial coefficient iterator
  107. This routine has been optimized for `m` large with respect to `n` by taking
  108. advantage of the fact that when the monomial tuples `t` are stripped of
  109. zeros, their coefficient is the same as that of the monomial tuples from
  110. ``multinomial_coefficients(n, n)``. Therefore, the latter coefficients are
  111. precomputed to save memory and time.
  112. >>> from sympy.ntheory.multinomial import multinomial_coefficients
  113. >>> m53, m33 = multinomial_coefficients(5,3), multinomial_coefficients(3,3)
  114. >>> m53[(0,0,0,1,2)] == m53[(0,0,1,0,2)] == m53[(1,0,2,0,0)] == m33[(0,1,2)]
  115. True
  116. Examples
  117. ========
  118. >>> from sympy.ntheory.multinomial import multinomial_coefficients_iterator
  119. >>> it = multinomial_coefficients_iterator(20,3)
  120. >>> next(it)
  121. ((3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 1)
  122. """
  123. m = as_int(m)
  124. n = as_int(n)
  125. if m < 2*n or n == 1:
  126. mc = multinomial_coefficients(m, n)
  127. yield from mc.items()
  128. else:
  129. mc = multinomial_coefficients(n, n)
  130. mc1 = {}
  131. for k, v in mc.items():
  132. mc1[_tuple(filter(None, k))] = v
  133. mc = mc1
  134. t = [n] + [0] * (m - 1)
  135. t1 = _tuple(t)
  136. b = _tuple(filter(None, t1))
  137. yield (t1, mc[b])
  138. if n:
  139. j = 0 # j will be the leftmost nonzero position
  140. else:
  141. j = m
  142. # enumerate tuples in co-lex order
  143. while j < m - 1:
  144. # compute next tuple
  145. tj = t[j]
  146. if j:
  147. t[j] = 0
  148. t[0] = tj
  149. if tj > 1:
  150. t[j + 1] += 1
  151. j = 0
  152. else:
  153. j += 1
  154. t[j] += 1
  155. t[0] -= 1
  156. t1 = _tuple(t)
  157. b = _tuple(filter(None, t1))
  158. yield (t1, mc[b])