dispersion.py 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212
  1. from sympy.core import S
  2. from sympy.polys import Poly
  3. def dispersionset(p, q=None, *gens, **args):
  4. r"""Compute the *dispersion set* of two polynomials.
  5. For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
  6. and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as:
  7. .. math::
  8. \operatorname{J}(f, g)
  9. & := \{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\} \\
  10. & = \{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\}
  11. For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`.
  12. Examples
  13. ========
  14. >>> from sympy import poly
  15. >>> from sympy.polys.dispersion import dispersion, dispersionset
  16. >>> from sympy.abc import x
  17. Dispersion set and dispersion of a simple polynomial:
  18. >>> fp = poly((x - 3)*(x + 3), x)
  19. >>> sorted(dispersionset(fp))
  20. [0, 6]
  21. >>> dispersion(fp)
  22. 6
  23. Note that the definition of the dispersion is not symmetric:
  24. >>> fp = poly(x**4 - 3*x**2 + 1, x)
  25. >>> gp = fp.shift(-3)
  26. >>> sorted(dispersionset(fp, gp))
  27. [2, 3, 4]
  28. >>> dispersion(fp, gp)
  29. 4
  30. >>> sorted(dispersionset(gp, fp))
  31. []
  32. >>> dispersion(gp, fp)
  33. -oo
  34. Computing the dispersion also works over field extensions:
  35. >>> from sympy import sqrt
  36. >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
  37. >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
  38. >>> sorted(dispersionset(fp, gp))
  39. [2]
  40. >>> sorted(dispersionset(gp, fp))
  41. [1, 4]
  42. We can even perform the computations for polynomials
  43. having symbolic coefficients:
  44. >>> from sympy.abc import a
  45. >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
  46. >>> sorted(dispersionset(fp))
  47. [0, 1]
  48. See Also
  49. ========
  50. dispersion
  51. References
  52. ==========
  53. .. [1] [ManWright94]_
  54. .. [2] [Koepf98]_
  55. .. [3] [Abramov71]_
  56. .. [4] [Man93]_
  57. """
  58. # Check for valid input
  59. same = False if q is not None else True
  60. if same:
  61. q = p
  62. p = Poly(p, *gens, **args)
  63. q = Poly(q, *gens, **args)
  64. if not p.is_univariate or not q.is_univariate:
  65. raise ValueError("Polynomials need to be univariate")
  66. # The generator
  67. if not p.gen == q.gen:
  68. raise ValueError("Polynomials must have the same generator")
  69. gen = p.gen
  70. # We define the dispersion of constant polynomials to be zero
  71. if p.degree() < 1 or q.degree() < 1:
  72. return {0}
  73. # Factor p and q over the rationals
  74. fp = p.factor_list()
  75. fq = q.factor_list() if not same else fp
  76. # Iterate over all pairs of factors
  77. J = set()
  78. for s, unused in fp[1]:
  79. for t, unused in fq[1]:
  80. m = s.degree()
  81. n = t.degree()
  82. if n != m:
  83. continue
  84. an = s.LC()
  85. bn = t.LC()
  86. if not (an - bn).is_zero:
  87. continue
  88. # Note that the roles of `s` and `t` below are switched
  89. # w.r.t. the original paper. This is for consistency
  90. # with the description in the book of W. Koepf.
  91. anm1 = s.coeff_monomial(gen**(m-1))
  92. bnm1 = t.coeff_monomial(gen**(n-1))
  93. alpha = (anm1 - bnm1) / S(n*bn)
  94. if not alpha.is_integer:
  95. continue
  96. if alpha < 0 or alpha in J:
  97. continue
  98. if n > 1 and not (s - t.shift(alpha)).is_zero:
  99. continue
  100. J.add(alpha)
  101. return J
  102. def dispersion(p, q=None, *gens, **args):
  103. r"""Compute the *dispersion* of polynomials.
  104. For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
  105. and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as:
  106. .. math::
  107. \operatorname{dis}(f, g)
  108. & := \max\{ J(f,g) \cup \{0\} \} \\
  109. & = \max\{ \{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\} \cup \{0\} \}
  110. and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`.
  111. Note that we make the definition `\max\{\} := -\infty`.
  112. Examples
  113. ========
  114. >>> from sympy import poly
  115. >>> from sympy.polys.dispersion import dispersion, dispersionset
  116. >>> from sympy.abc import x
  117. Dispersion set and dispersion of a simple polynomial:
  118. >>> fp = poly((x - 3)*(x + 3), x)
  119. >>> sorted(dispersionset(fp))
  120. [0, 6]
  121. >>> dispersion(fp)
  122. 6
  123. Note that the definition of the dispersion is not symmetric:
  124. >>> fp = poly(x**4 - 3*x**2 + 1, x)
  125. >>> gp = fp.shift(-3)
  126. >>> sorted(dispersionset(fp, gp))
  127. [2, 3, 4]
  128. >>> dispersion(fp, gp)
  129. 4
  130. >>> sorted(dispersionset(gp, fp))
  131. []
  132. >>> dispersion(gp, fp)
  133. -oo
  134. The maximum of an empty set is defined to be `-\infty`
  135. as seen in this example.
  136. Computing the dispersion also works over field extensions:
  137. >>> from sympy import sqrt
  138. >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
  139. >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
  140. >>> sorted(dispersionset(fp, gp))
  141. [2]
  142. >>> sorted(dispersionset(gp, fp))
  143. [1, 4]
  144. We can even perform the computations for polynomials
  145. having symbolic coefficients:
  146. >>> from sympy.abc import a
  147. >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
  148. >>> sorted(dispersionset(fp))
  149. [0, 1]
  150. See Also
  151. ========
  152. dispersionset
  153. References
  154. ==========
  155. .. [1] [ManWright94]_
  156. .. [2] [Koepf98]_
  157. .. [3] [Abramov71]_
  158. .. [4] [Man93]_
  159. """
  160. J = dispersionset(p, q, *gens, **args)
  161. if not J:
  162. # Definition for maximum of empty set
  163. j = S.NegativeInfinity
  164. else:
  165. j = max(J)
  166. return j