polyfuncs.py 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321
  1. """High-level polynomials manipulation functions. """
  2. from sympy.core import S, Basic, symbols, Dummy
  3. from sympy.polys.polyerrors import (
  4. PolificationFailed, ComputationFailed,
  5. MultivariatePolynomialError, OptionError)
  6. from sympy.polys.polyoptions import allowed_flags, build_options
  7. from sympy.polys.polytools import poly_from_expr, Poly
  8. from sympy.polys.specialpolys import (
  9. symmetric_poly, interpolating_poly)
  10. from sympy.polys.rings import sring
  11. from sympy.utilities import numbered_symbols, take, public
  12. @public
  13. def symmetrize(F, *gens, **args):
  14. r"""
  15. Rewrite a polynomial in terms of elementary symmetric polynomials.
  16. A symmetric polynomial is a multivariate polynomial that remains invariant
  17. under any variable permutation, i.e., if `f = f(x_1, x_2, \dots, x_n)`,
  18. then `f = f(x_{i_1}, x_{i_2}, \dots, x_{i_n})`, where
  19. `(i_1, i_2, \dots, i_n)` is a permutation of `(1, 2, \dots, n)` (an
  20. element of the group `S_n`).
  21. Returns a tuple of symmetric polynomials ``(f1, f2, ..., fn)`` such that
  22. ``f = f1 + f2 + ... + fn``.
  23. Examples
  24. ========
  25. >>> from sympy.polys.polyfuncs import symmetrize
  26. >>> from sympy.abc import x, y
  27. >>> symmetrize(x**2 + y**2)
  28. (-2*x*y + (x + y)**2, 0)
  29. >>> symmetrize(x**2 + y**2, formal=True)
  30. (s1**2 - 2*s2, 0, [(s1, x + y), (s2, x*y)])
  31. >>> symmetrize(x**2 - y**2)
  32. (-2*x*y + (x + y)**2, -2*y**2)
  33. >>> symmetrize(x**2 - y**2, formal=True)
  34. (s1**2 - 2*s2, -2*y**2, [(s1, x + y), (s2, x*y)])
  35. """
  36. allowed_flags(args, ['formal', 'symbols'])
  37. iterable = True
  38. if not hasattr(F, '__iter__'):
  39. iterable = False
  40. F = [F]
  41. R, F = sring(F, *gens, **args)
  42. gens = R.symbols
  43. opt = build_options(gens, args)
  44. symbols = opt.symbols
  45. symbols = [next(symbols) for i in range(len(gens))]
  46. result = []
  47. for f in F:
  48. p, r, m = f.symmetrize()
  49. result.append((p.as_expr(*symbols), r.as_expr(*gens)))
  50. polys = [(s, g.as_expr()) for s, (_, g) in zip(symbols, m)]
  51. if not opt.formal:
  52. for i, (sym, non_sym) in enumerate(result):
  53. result[i] = (sym.subs(polys), non_sym)
  54. if not iterable:
  55. result, = result
  56. if not opt.formal:
  57. return result
  58. else:
  59. if iterable:
  60. return result, polys
  61. else:
  62. return result + (polys,)
  63. @public
  64. def horner(f, *gens, **args):
  65. """
  66. Rewrite a polynomial in Horner form.
  67. Among other applications, evaluation of a polynomial at a point is optimal
  68. when it is applied using the Horner scheme ([1]).
  69. Examples
  70. ========
  71. >>> from sympy.polys.polyfuncs import horner
  72. >>> from sympy.abc import x, y, a, b, c, d, e
  73. >>> horner(9*x**4 + 8*x**3 + 7*x**2 + 6*x + 5)
  74. x*(x*(x*(9*x + 8) + 7) + 6) + 5
  75. >>> horner(a*x**4 + b*x**3 + c*x**2 + d*x + e)
  76. e + x*(d + x*(c + x*(a*x + b)))
  77. >>> f = 4*x**2*y**2 + 2*x**2*y + 2*x*y**2 + x*y
  78. >>> horner(f, wrt=x)
  79. x*(x*y*(4*y + 2) + y*(2*y + 1))
  80. >>> horner(f, wrt=y)
  81. y*(x*y*(4*x + 2) + x*(2*x + 1))
  82. References
  83. ==========
  84. [1] - https://en.wikipedia.org/wiki/Horner_scheme
  85. """
  86. allowed_flags(args, [])
  87. try:
  88. F, opt = poly_from_expr(f, *gens, **args)
  89. except PolificationFailed as exc:
  90. return exc.expr
  91. form, gen = S.Zero, F.gen
  92. if F.is_univariate:
  93. for coeff in F.all_coeffs():
  94. form = form*gen + coeff
  95. else:
  96. F, gens = Poly(F, gen), gens[1:]
  97. for coeff in F.all_coeffs():
  98. form = form*gen + horner(coeff, *gens, **args)
  99. return form
  100. @public
  101. def interpolate(data, x):
  102. """
  103. Construct an interpolating polynomial for the data points
  104. evaluated at point x (which can be symbolic or numeric).
  105. Examples
  106. ========
  107. >>> from sympy.polys.polyfuncs import interpolate
  108. >>> from sympy.abc import a, b, x
  109. A list is interpreted as though it were paired with a range starting
  110. from 1:
  111. >>> interpolate([1, 4, 9, 16], x)
  112. x**2
  113. This can be made explicit by giving a list of coordinates:
  114. >>> interpolate([(1, 1), (2, 4), (3, 9)], x)
  115. x**2
  116. The (x, y) coordinates can also be given as keys and values of a
  117. dictionary (and the points need not be equispaced):
  118. >>> interpolate([(-1, 2), (1, 2), (2, 5)], x)
  119. x**2 + 1
  120. >>> interpolate({-1: 2, 1: 2, 2: 5}, x)
  121. x**2 + 1
  122. If the interpolation is going to be used only once then the
  123. value of interest can be passed instead of passing a symbol:
  124. >>> interpolate([1, 4, 9], 5)
  125. 25
  126. Symbolic coordinates are also supported:
  127. >>> [(i,interpolate((a, b), i)) for i in range(1, 4)]
  128. [(1, a), (2, b), (3, -a + 2*b)]
  129. """
  130. n = len(data)
  131. if isinstance(data, dict):
  132. if x in data:
  133. return S(data[x])
  134. X, Y = list(zip(*data.items()))
  135. else:
  136. if isinstance(data[0], tuple):
  137. X, Y = list(zip(*data))
  138. if x in X:
  139. return S(Y[X.index(x)])
  140. else:
  141. if x in range(1, n + 1):
  142. return S(data[x - 1])
  143. Y = list(data)
  144. X = list(range(1, n + 1))
  145. try:
  146. return interpolating_poly(n, x, X, Y).expand()
  147. except ValueError:
  148. d = Dummy()
  149. return interpolating_poly(n, d, X, Y).expand().subs(d, x)
  150. @public
  151. def rational_interpolate(data, degnum, X=symbols('x')):
  152. """
  153. Returns a rational interpolation, where the data points are element of
  154. any integral domain.
  155. The first argument contains the data (as a list of coordinates). The
  156. ``degnum`` argument is the degree in the numerator of the rational
  157. function. Setting it too high will decrease the maximal degree in the
  158. denominator for the same amount of data.
  159. Examples
  160. ========
  161. >>> from sympy.polys.polyfuncs import rational_interpolate
  162. >>> data = [(1, -210), (2, -35), (3, 105), (4, 231), (5, 350), (6, 465)]
  163. >>> rational_interpolate(data, 2)
  164. (105*x**2 - 525)/(x + 1)
  165. Values do not need to be integers:
  166. >>> from sympy import sympify
  167. >>> x = [1, 2, 3, 4, 5, 6]
  168. >>> y = sympify("[-1, 0, 2, 22/5, 7, 68/7]")
  169. >>> rational_interpolate(zip(x, y), 2)
  170. (3*x**2 - 7*x + 2)/(x + 1)
  171. The symbol for the variable can be changed if needed:
  172. >>> from sympy import symbols
  173. >>> z = symbols('z')
  174. >>> rational_interpolate(data, 2, X=z)
  175. (105*z**2 - 525)/(z + 1)
  176. References
  177. ==========
  178. .. [1] Algorithm is adapted from:
  179. http://axiom-wiki.newsynthesis.org/RationalInterpolation
  180. """
  181. from sympy.matrices.dense import ones
  182. xdata, ydata = list(zip(*data))
  183. k = len(xdata) - degnum - 1
  184. if k < 0:
  185. raise OptionError("Too few values for the required degree.")
  186. c = ones(degnum + k + 1, degnum + k + 2)
  187. for j in range(max(degnum, k)):
  188. for i in range(degnum + k + 1):
  189. c[i, j + 1] = c[i, j]*xdata[i]
  190. for j in range(k + 1):
  191. for i in range(degnum + k + 1):
  192. c[i, degnum + k + 1 - j] = -c[i, k - j]*ydata[i]
  193. r = c.nullspace()[0]
  194. return (sum(r[i] * X**i for i in range(degnum + 1))
  195. / sum(r[i + degnum + 1] * X**i for i in range(k + 1)))
  196. @public
  197. def viete(f, roots=None, *gens, **args):
  198. """
  199. Generate Viete's formulas for ``f``.
  200. Examples
  201. ========
  202. >>> from sympy.polys.polyfuncs import viete
  203. >>> from sympy import symbols
  204. >>> x, a, b, c, r1, r2 = symbols('x,a:c,r1:3')
  205. >>> viete(a*x**2 + b*x + c, [r1, r2], x)
  206. [(r1 + r2, -b/a), (r1*r2, c/a)]
  207. """
  208. allowed_flags(args, [])
  209. if isinstance(roots, Basic):
  210. gens, roots = (roots,) + gens, None
  211. try:
  212. f, opt = poly_from_expr(f, *gens, **args)
  213. except PolificationFailed as exc:
  214. raise ComputationFailed('viete', 1, exc)
  215. if f.is_multivariate:
  216. raise MultivariatePolynomialError(
  217. "multivariate polynomials are not allowed")
  218. n = f.degree()
  219. if n < 1:
  220. raise ValueError(
  221. "Cannot derive Viete's formulas for a constant polynomial")
  222. if roots is None:
  223. roots = numbered_symbols('r', start=1)
  224. roots = take(roots, n)
  225. if n != len(roots):
  226. raise ValueError("required %s roots, got %s" % (n, len(roots)))
  227. lc, coeffs = f.LC(), f.all_coeffs()
  228. result, sign = [], -1
  229. for i, coeff in enumerate(coeffs[1:]):
  230. poly = symmetric_poly(i + 1, roots)
  231. coeff = sign*(coeff/lc)
  232. result.append((poly, coeff))
  233. sign = -sign
  234. return result