densetools.py 25 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309
  1. """Advanced tools for dense recursive polynomials in ``K[x]`` or ``K[X]``. """
  2. from sympy.polys.densearith import (
  3. dup_add_term, dmp_add_term,
  4. dup_lshift,
  5. dup_add, dmp_add,
  6. dup_sub, dmp_sub,
  7. dup_mul, dmp_mul,
  8. dup_sqr,
  9. dup_div,
  10. dup_rem, dmp_rem,
  11. dmp_expand,
  12. dup_mul_ground, dmp_mul_ground,
  13. dup_quo_ground, dmp_quo_ground,
  14. dup_exquo_ground, dmp_exquo_ground,
  15. )
  16. from sympy.polys.densebasic import (
  17. dup_strip, dmp_strip,
  18. dup_convert, dmp_convert,
  19. dup_degree, dmp_degree,
  20. dmp_to_dict,
  21. dmp_from_dict,
  22. dup_LC, dmp_LC, dmp_ground_LC,
  23. dup_TC, dmp_TC,
  24. dmp_zero, dmp_ground,
  25. dmp_zero_p,
  26. dup_to_raw_dict, dup_from_raw_dict,
  27. dmp_zeros
  28. )
  29. from sympy.polys.polyerrors import (
  30. MultivariatePolynomialError,
  31. DomainError
  32. )
  33. from sympy.utilities import variations
  34. from math import ceil as _ceil, log as _log
  35. def dup_integrate(f, m, K):
  36. """
  37. Computes the indefinite integral of ``f`` in ``K[x]``.
  38. Examples
  39. ========
  40. >>> from sympy.polys import ring, QQ
  41. >>> R, x = ring("x", QQ)
  42. >>> R.dup_integrate(x**2 + 2*x, 1)
  43. 1/3*x**3 + x**2
  44. >>> R.dup_integrate(x**2 + 2*x, 2)
  45. 1/12*x**4 + 1/3*x**3
  46. """
  47. if m <= 0 or not f:
  48. return f
  49. g = [K.zero]*m
  50. for i, c in enumerate(reversed(f)):
  51. n = i + 1
  52. for j in range(1, m):
  53. n *= i + j + 1
  54. g.insert(0, K.exquo(c, K(n)))
  55. return g
  56. def dmp_integrate(f, m, u, K):
  57. """
  58. Computes the indefinite integral of ``f`` in ``x_0`` in ``K[X]``.
  59. Examples
  60. ========
  61. >>> from sympy.polys import ring, QQ
  62. >>> R, x,y = ring("x,y", QQ)
  63. >>> R.dmp_integrate(x + 2*y, 1)
  64. 1/2*x**2 + 2*x*y
  65. >>> R.dmp_integrate(x + 2*y, 2)
  66. 1/6*x**3 + x**2*y
  67. """
  68. if not u:
  69. return dup_integrate(f, m, K)
  70. if m <= 0 or dmp_zero_p(f, u):
  71. return f
  72. g, v = dmp_zeros(m, u - 1, K), u - 1
  73. for i, c in enumerate(reversed(f)):
  74. n = i + 1
  75. for j in range(1, m):
  76. n *= i + j + 1
  77. g.insert(0, dmp_quo_ground(c, K(n), v, K))
  78. return g
  79. def _rec_integrate_in(g, m, v, i, j, K):
  80. """Recursive helper for :func:`dmp_integrate_in`."""
  81. if i == j:
  82. return dmp_integrate(g, m, v, K)
  83. w, i = v - 1, i + 1
  84. return dmp_strip([ _rec_integrate_in(c, m, w, i, j, K) for c in g ], v)
  85. def dmp_integrate_in(f, m, j, u, K):
  86. """
  87. Computes the indefinite integral of ``f`` in ``x_j`` in ``K[X]``.
  88. Examples
  89. ========
  90. >>> from sympy.polys import ring, QQ
  91. >>> R, x,y = ring("x,y", QQ)
  92. >>> R.dmp_integrate_in(x + 2*y, 1, 0)
  93. 1/2*x**2 + 2*x*y
  94. >>> R.dmp_integrate_in(x + 2*y, 1, 1)
  95. x*y + y**2
  96. """
  97. if j < 0 or j > u:
  98. raise IndexError("0 <= j <= u expected, got u = %d, j = %d" % (u, j))
  99. return _rec_integrate_in(f, m, u, 0, j, K)
  100. def dup_diff(f, m, K):
  101. """
  102. ``m``-th order derivative of a polynomial in ``K[x]``.
  103. Examples
  104. ========
  105. >>> from sympy.polys import ring, ZZ
  106. >>> R, x = ring("x", ZZ)
  107. >>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 1)
  108. 3*x**2 + 4*x + 3
  109. >>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 2)
  110. 6*x + 4
  111. """
  112. if m <= 0:
  113. return f
  114. n = dup_degree(f)
  115. if n < m:
  116. return []
  117. deriv = []
  118. if m == 1:
  119. for coeff in f[:-m]:
  120. deriv.append(K(n)*coeff)
  121. n -= 1
  122. else:
  123. for coeff in f[:-m]:
  124. k = n
  125. for i in range(n - 1, n - m, -1):
  126. k *= i
  127. deriv.append(K(k)*coeff)
  128. n -= 1
  129. return dup_strip(deriv)
  130. def dmp_diff(f, m, u, K):
  131. """
  132. ``m``-th order derivative in ``x_0`` of a polynomial in ``K[X]``.
  133. Examples
  134. ========
  135. >>> from sympy.polys import ring, ZZ
  136. >>> R, x,y = ring("x,y", ZZ)
  137. >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
  138. >>> R.dmp_diff(f, 1)
  139. y**2 + 2*y + 3
  140. >>> R.dmp_diff(f, 2)
  141. 0
  142. """
  143. if not u:
  144. return dup_diff(f, m, K)
  145. if m <= 0:
  146. return f
  147. n = dmp_degree(f, u)
  148. if n < m:
  149. return dmp_zero(u)
  150. deriv, v = [], u - 1
  151. if m == 1:
  152. for coeff in f[:-m]:
  153. deriv.append(dmp_mul_ground(coeff, K(n), v, K))
  154. n -= 1
  155. else:
  156. for coeff in f[:-m]:
  157. k = n
  158. for i in range(n - 1, n - m, -1):
  159. k *= i
  160. deriv.append(dmp_mul_ground(coeff, K(k), v, K))
  161. n -= 1
  162. return dmp_strip(deriv, u)
  163. def _rec_diff_in(g, m, v, i, j, K):
  164. """Recursive helper for :func:`dmp_diff_in`."""
  165. if i == j:
  166. return dmp_diff(g, m, v, K)
  167. w, i = v - 1, i + 1
  168. return dmp_strip([ _rec_diff_in(c, m, w, i, j, K) for c in g ], v)
  169. def dmp_diff_in(f, m, j, u, K):
  170. """
  171. ``m``-th order derivative in ``x_j`` of a polynomial in ``K[X]``.
  172. Examples
  173. ========
  174. >>> from sympy.polys import ring, ZZ
  175. >>> R, x,y = ring("x,y", ZZ)
  176. >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
  177. >>> R.dmp_diff_in(f, 1, 0)
  178. y**2 + 2*y + 3
  179. >>> R.dmp_diff_in(f, 1, 1)
  180. 2*x*y + 2*x + 4*y + 3
  181. """
  182. if j < 0 or j > u:
  183. raise IndexError("0 <= j <= %s expected, got %s" % (u, j))
  184. return _rec_diff_in(f, m, u, 0, j, K)
  185. def dup_eval(f, a, K):
  186. """
  187. Evaluate a polynomial at ``x = a`` in ``K[x]`` using Horner scheme.
  188. Examples
  189. ========
  190. >>> from sympy.polys import ring, ZZ
  191. >>> R, x = ring("x", ZZ)
  192. >>> R.dup_eval(x**2 + 2*x + 3, 2)
  193. 11
  194. """
  195. if not a:
  196. return dup_TC(f, K)
  197. result = K.zero
  198. for c in f:
  199. result *= a
  200. result += c
  201. return result
  202. def dmp_eval(f, a, u, K):
  203. """
  204. Evaluate a polynomial at ``x_0 = a`` in ``K[X]`` using the Horner scheme.
  205. Examples
  206. ========
  207. >>> from sympy.polys import ring, ZZ
  208. >>> R, x,y = ring("x,y", ZZ)
  209. >>> R.dmp_eval(2*x*y + 3*x + y + 2, 2)
  210. 5*y + 8
  211. """
  212. if not u:
  213. return dup_eval(f, a, K)
  214. if not a:
  215. return dmp_TC(f, K)
  216. result, v = dmp_LC(f, K), u - 1
  217. for coeff in f[1:]:
  218. result = dmp_mul_ground(result, a, v, K)
  219. result = dmp_add(result, coeff, v, K)
  220. return result
  221. def _rec_eval_in(g, a, v, i, j, K):
  222. """Recursive helper for :func:`dmp_eval_in`."""
  223. if i == j:
  224. return dmp_eval(g, a, v, K)
  225. v, i = v - 1, i + 1
  226. return dmp_strip([ _rec_eval_in(c, a, v, i, j, K) for c in g ], v)
  227. def dmp_eval_in(f, a, j, u, K):
  228. """
  229. Evaluate a polynomial at ``x_j = a`` in ``K[X]`` using the Horner scheme.
  230. Examples
  231. ========
  232. >>> from sympy.polys import ring, ZZ
  233. >>> R, x,y = ring("x,y", ZZ)
  234. >>> f = 2*x*y + 3*x + y + 2
  235. >>> R.dmp_eval_in(f, 2, 0)
  236. 5*y + 8
  237. >>> R.dmp_eval_in(f, 2, 1)
  238. 7*x + 4
  239. """
  240. if j < 0 or j > u:
  241. raise IndexError("0 <= j <= %s expected, got %s" % (u, j))
  242. return _rec_eval_in(f, a, u, 0, j, K)
  243. def _rec_eval_tail(g, i, A, u, K):
  244. """Recursive helper for :func:`dmp_eval_tail`."""
  245. if i == u:
  246. return dup_eval(g, A[-1], K)
  247. else:
  248. h = [ _rec_eval_tail(c, i + 1, A, u, K) for c in g ]
  249. if i < u - len(A) + 1:
  250. return h
  251. else:
  252. return dup_eval(h, A[-u + i - 1], K)
  253. def dmp_eval_tail(f, A, u, K):
  254. """
  255. Evaluate a polynomial at ``x_j = a_j, ...`` in ``K[X]``.
  256. Examples
  257. ========
  258. >>> from sympy.polys import ring, ZZ
  259. >>> R, x,y = ring("x,y", ZZ)
  260. >>> f = 2*x*y + 3*x + y + 2
  261. >>> R.dmp_eval_tail(f, [2])
  262. 7*x + 4
  263. >>> R.dmp_eval_tail(f, [2, 2])
  264. 18
  265. """
  266. if not A:
  267. return f
  268. if dmp_zero_p(f, u):
  269. return dmp_zero(u - len(A))
  270. e = _rec_eval_tail(f, 0, A, u, K)
  271. if u == len(A) - 1:
  272. return e
  273. else:
  274. return dmp_strip(e, u - len(A))
  275. def _rec_diff_eval(g, m, a, v, i, j, K):
  276. """Recursive helper for :func:`dmp_diff_eval`."""
  277. if i == j:
  278. return dmp_eval(dmp_diff(g, m, v, K), a, v, K)
  279. v, i = v - 1, i + 1
  280. return dmp_strip([ _rec_diff_eval(c, m, a, v, i, j, K) for c in g ], v)
  281. def dmp_diff_eval_in(f, m, a, j, u, K):
  282. """
  283. Differentiate and evaluate a polynomial in ``x_j`` at ``a`` in ``K[X]``.
  284. Examples
  285. ========
  286. >>> from sympy.polys import ring, ZZ
  287. >>> R, x,y = ring("x,y", ZZ)
  288. >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
  289. >>> R.dmp_diff_eval_in(f, 1, 2, 0)
  290. y**2 + 2*y + 3
  291. >>> R.dmp_diff_eval_in(f, 1, 2, 1)
  292. 6*x + 11
  293. """
  294. if j > u:
  295. raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
  296. if not j:
  297. return dmp_eval(dmp_diff(f, m, u, K), a, u, K)
  298. return _rec_diff_eval(f, m, a, u, 0, j, K)
  299. def dup_trunc(f, p, K):
  300. """
  301. Reduce a ``K[x]`` polynomial modulo a constant ``p`` in ``K``.
  302. Examples
  303. ========
  304. >>> from sympy.polys import ring, ZZ
  305. >>> R, x = ring("x", ZZ)
  306. >>> R.dup_trunc(2*x**3 + 3*x**2 + 5*x + 7, ZZ(3))
  307. -x**3 - x + 1
  308. """
  309. if K.is_ZZ:
  310. g = []
  311. for c in f:
  312. c = c % p
  313. if c > p // 2:
  314. g.append(c - p)
  315. else:
  316. g.append(c)
  317. else:
  318. g = [ c % p for c in f ]
  319. return dup_strip(g)
  320. def dmp_trunc(f, p, u, K):
  321. """
  322. Reduce a ``K[X]`` polynomial modulo a polynomial ``p`` in ``K[Y]``.
  323. Examples
  324. ========
  325. >>> from sympy.polys import ring, ZZ
  326. >>> R, x,y = ring("x,y", ZZ)
  327. >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
  328. >>> g = (y - 1).drop(x)
  329. >>> R.dmp_trunc(f, g)
  330. 11*x**2 + 11*x + 5
  331. """
  332. return dmp_strip([ dmp_rem(c, p, u - 1, K) for c in f ], u)
  333. def dmp_ground_trunc(f, p, u, K):
  334. """
  335. Reduce a ``K[X]`` polynomial modulo a constant ``p`` in ``K``.
  336. Examples
  337. ========
  338. >>> from sympy.polys import ring, ZZ
  339. >>> R, x,y = ring("x,y", ZZ)
  340. >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
  341. >>> R.dmp_ground_trunc(f, ZZ(3))
  342. -x**2 - x*y - y
  343. """
  344. if not u:
  345. return dup_trunc(f, p, K)
  346. v = u - 1
  347. return dmp_strip([ dmp_ground_trunc(c, p, v, K) for c in f ], u)
  348. def dup_monic(f, K):
  349. """
  350. Divide all coefficients by ``LC(f)`` in ``K[x]``.
  351. Examples
  352. ========
  353. >>> from sympy.polys import ring, ZZ, QQ
  354. >>> R, x = ring("x", ZZ)
  355. >>> R.dup_monic(3*x**2 + 6*x + 9)
  356. x**2 + 2*x + 3
  357. >>> R, x = ring("x", QQ)
  358. >>> R.dup_monic(3*x**2 + 4*x + 2)
  359. x**2 + 4/3*x + 2/3
  360. """
  361. if not f:
  362. return f
  363. lc = dup_LC(f, K)
  364. if K.is_one(lc):
  365. return f
  366. else:
  367. return dup_exquo_ground(f, lc, K)
  368. def dmp_ground_monic(f, u, K):
  369. """
  370. Divide all coefficients by ``LC(f)`` in ``K[X]``.
  371. Examples
  372. ========
  373. >>> from sympy.polys import ring, ZZ, QQ
  374. >>> R, x,y = ring("x,y", ZZ)
  375. >>> f = 3*x**2*y + 6*x**2 + 3*x*y + 9*y + 3
  376. >>> R.dmp_ground_monic(f)
  377. x**2*y + 2*x**2 + x*y + 3*y + 1
  378. >>> R, x,y = ring("x,y", QQ)
  379. >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
  380. >>> R.dmp_ground_monic(f)
  381. x**2*y + 8/3*x**2 + 5/3*x*y + 2*x + 2/3*y + 1
  382. """
  383. if not u:
  384. return dup_monic(f, K)
  385. if dmp_zero_p(f, u):
  386. return f
  387. lc = dmp_ground_LC(f, u, K)
  388. if K.is_one(lc):
  389. return f
  390. else:
  391. return dmp_exquo_ground(f, lc, u, K)
  392. def dup_content(f, K):
  393. """
  394. Compute the GCD of coefficients of ``f`` in ``K[x]``.
  395. Examples
  396. ========
  397. >>> from sympy.polys import ring, ZZ, QQ
  398. >>> R, x = ring("x", ZZ)
  399. >>> f = 6*x**2 + 8*x + 12
  400. >>> R.dup_content(f)
  401. 2
  402. >>> R, x = ring("x", QQ)
  403. >>> f = 6*x**2 + 8*x + 12
  404. >>> R.dup_content(f)
  405. 2
  406. """
  407. from sympy.polys.domains import QQ
  408. if not f:
  409. return K.zero
  410. cont = K.zero
  411. if K == QQ:
  412. for c in f:
  413. cont = K.gcd(cont, c)
  414. else:
  415. for c in f:
  416. cont = K.gcd(cont, c)
  417. if K.is_one(cont):
  418. break
  419. return cont
  420. def dmp_ground_content(f, u, K):
  421. """
  422. Compute the GCD of coefficients of ``f`` in ``K[X]``.
  423. Examples
  424. ========
  425. >>> from sympy.polys import ring, ZZ, QQ
  426. >>> R, x,y = ring("x,y", ZZ)
  427. >>> f = 2*x*y + 6*x + 4*y + 12
  428. >>> R.dmp_ground_content(f)
  429. 2
  430. >>> R, x,y = ring("x,y", QQ)
  431. >>> f = 2*x*y + 6*x + 4*y + 12
  432. >>> R.dmp_ground_content(f)
  433. 2
  434. """
  435. from sympy.polys.domains import QQ
  436. if not u:
  437. return dup_content(f, K)
  438. if dmp_zero_p(f, u):
  439. return K.zero
  440. cont, v = K.zero, u - 1
  441. if K == QQ:
  442. for c in f:
  443. cont = K.gcd(cont, dmp_ground_content(c, v, K))
  444. else:
  445. for c in f:
  446. cont = K.gcd(cont, dmp_ground_content(c, v, K))
  447. if K.is_one(cont):
  448. break
  449. return cont
  450. def dup_primitive(f, K):
  451. """
  452. Compute content and the primitive form of ``f`` in ``K[x]``.
  453. Examples
  454. ========
  455. >>> from sympy.polys import ring, ZZ, QQ
  456. >>> R, x = ring("x", ZZ)
  457. >>> f = 6*x**2 + 8*x + 12
  458. >>> R.dup_primitive(f)
  459. (2, 3*x**2 + 4*x + 6)
  460. >>> R, x = ring("x", QQ)
  461. >>> f = 6*x**2 + 8*x + 12
  462. >>> R.dup_primitive(f)
  463. (2, 3*x**2 + 4*x + 6)
  464. """
  465. if not f:
  466. return K.zero, f
  467. cont = dup_content(f, K)
  468. if K.is_one(cont):
  469. return cont, f
  470. else:
  471. return cont, dup_quo_ground(f, cont, K)
  472. def dmp_ground_primitive(f, u, K):
  473. """
  474. Compute content and the primitive form of ``f`` in ``K[X]``.
  475. Examples
  476. ========
  477. >>> from sympy.polys import ring, ZZ, QQ
  478. >>> R, x,y = ring("x,y", ZZ)
  479. >>> f = 2*x*y + 6*x + 4*y + 12
  480. >>> R.dmp_ground_primitive(f)
  481. (2, x*y + 3*x + 2*y + 6)
  482. >>> R, x,y = ring("x,y", QQ)
  483. >>> f = 2*x*y + 6*x + 4*y + 12
  484. >>> R.dmp_ground_primitive(f)
  485. (2, x*y + 3*x + 2*y + 6)
  486. """
  487. if not u:
  488. return dup_primitive(f, K)
  489. if dmp_zero_p(f, u):
  490. return K.zero, f
  491. cont = dmp_ground_content(f, u, K)
  492. if K.is_one(cont):
  493. return cont, f
  494. else:
  495. return cont, dmp_quo_ground(f, cont, u, K)
  496. def dup_extract(f, g, K):
  497. """
  498. Extract common content from a pair of polynomials in ``K[x]``.
  499. Examples
  500. ========
  501. >>> from sympy.polys import ring, ZZ
  502. >>> R, x = ring("x", ZZ)
  503. >>> R.dup_extract(6*x**2 + 12*x + 18, 4*x**2 + 8*x + 12)
  504. (2, 3*x**2 + 6*x + 9, 2*x**2 + 4*x + 6)
  505. """
  506. fc = dup_content(f, K)
  507. gc = dup_content(g, K)
  508. gcd = K.gcd(fc, gc)
  509. if not K.is_one(gcd):
  510. f = dup_quo_ground(f, gcd, K)
  511. g = dup_quo_ground(g, gcd, K)
  512. return gcd, f, g
  513. def dmp_ground_extract(f, g, u, K):
  514. """
  515. Extract common content from a pair of polynomials in ``K[X]``.
  516. Examples
  517. ========
  518. >>> from sympy.polys import ring, ZZ
  519. >>> R, x,y = ring("x,y", ZZ)
  520. >>> R.dmp_ground_extract(6*x*y + 12*x + 18, 4*x*y + 8*x + 12)
  521. (2, 3*x*y + 6*x + 9, 2*x*y + 4*x + 6)
  522. """
  523. fc = dmp_ground_content(f, u, K)
  524. gc = dmp_ground_content(g, u, K)
  525. gcd = K.gcd(fc, gc)
  526. if not K.is_one(gcd):
  527. f = dmp_quo_ground(f, gcd, u, K)
  528. g = dmp_quo_ground(g, gcd, u, K)
  529. return gcd, f, g
  530. def dup_real_imag(f, K):
  531. """
  532. Return bivariate polynomials ``f1`` and ``f2``, such that ``f = f1 + f2*I``.
  533. Examples
  534. ========
  535. >>> from sympy.polys import ring, ZZ
  536. >>> R, x,y = ring("x,y", ZZ)
  537. >>> R.dup_real_imag(x**3 + x**2 + x + 1)
  538. (x**3 + x**2 - 3*x*y**2 + x - y**2 + 1, 3*x**2*y + 2*x*y - y**3 + y)
  539. """
  540. if not K.is_ZZ and not K.is_QQ:
  541. raise DomainError("computing real and imaginary parts is not supported over %s" % K)
  542. f1 = dmp_zero(1)
  543. f2 = dmp_zero(1)
  544. if not f:
  545. return f1, f2
  546. g = [[[K.one, K.zero]], [[K.one], []]]
  547. h = dmp_ground(f[0], 2)
  548. for c in f[1:]:
  549. h = dmp_mul(h, g, 2, K)
  550. h = dmp_add_term(h, dmp_ground(c, 1), 0, 2, K)
  551. H = dup_to_raw_dict(h)
  552. for k, h in H.items():
  553. m = k % 4
  554. if not m:
  555. f1 = dmp_add(f1, h, 1, K)
  556. elif m == 1:
  557. f2 = dmp_add(f2, h, 1, K)
  558. elif m == 2:
  559. f1 = dmp_sub(f1, h, 1, K)
  560. else:
  561. f2 = dmp_sub(f2, h, 1, K)
  562. return f1, f2
  563. def dup_mirror(f, K):
  564. """
  565. Evaluate efficiently the composition ``f(-x)`` in ``K[x]``.
  566. Examples
  567. ========
  568. >>> from sympy.polys import ring, ZZ
  569. >>> R, x = ring("x", ZZ)
  570. >>> R.dup_mirror(x**3 + 2*x**2 - 4*x + 2)
  571. -x**3 + 2*x**2 + 4*x + 2
  572. """
  573. f = list(f)
  574. for i in range(len(f) - 2, -1, -2):
  575. f[i] = -f[i]
  576. return f
  577. def dup_scale(f, a, K):
  578. """
  579. Evaluate efficiently composition ``f(a*x)`` in ``K[x]``.
  580. Examples
  581. ========
  582. >>> from sympy.polys import ring, ZZ
  583. >>> R, x = ring("x", ZZ)
  584. >>> R.dup_scale(x**2 - 2*x + 1, ZZ(2))
  585. 4*x**2 - 4*x + 1
  586. """
  587. f, n, b = list(f), len(f) - 1, a
  588. for i in range(n - 1, -1, -1):
  589. f[i], b = b*f[i], b*a
  590. return f
  591. def dup_shift(f, a, K):
  592. """
  593. Evaluate efficiently Taylor shift ``f(x + a)`` in ``K[x]``.
  594. Examples
  595. ========
  596. >>> from sympy.polys import ring, ZZ
  597. >>> R, x = ring("x", ZZ)
  598. >>> R.dup_shift(x**2 - 2*x + 1, ZZ(2))
  599. x**2 + 2*x + 1
  600. """
  601. f, n = list(f), len(f) - 1
  602. for i in range(n, 0, -1):
  603. for j in range(0, i):
  604. f[j + 1] += a*f[j]
  605. return f
  606. def dup_transform(f, p, q, K):
  607. """
  608. Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.
  609. Examples
  610. ========
  611. >>> from sympy.polys import ring, ZZ
  612. >>> R, x = ring("x", ZZ)
  613. >>> R.dup_transform(x**2 - 2*x + 1, x**2 + 1, x - 1)
  614. x**4 - 2*x**3 + 5*x**2 - 4*x + 4
  615. """
  616. if not f:
  617. return []
  618. n = len(f) - 1
  619. h, Q = [f[0]], [[K.one]]
  620. for i in range(0, n):
  621. Q.append(dup_mul(Q[-1], q, K))
  622. for c, q in zip(f[1:], Q[1:]):
  623. h = dup_mul(h, p, K)
  624. q = dup_mul_ground(q, c, K)
  625. h = dup_add(h, q, K)
  626. return h
  627. def dup_compose(f, g, K):
  628. """
  629. Evaluate functional composition ``f(g)`` in ``K[x]``.
  630. Examples
  631. ========
  632. >>> from sympy.polys import ring, ZZ
  633. >>> R, x = ring("x", ZZ)
  634. >>> R.dup_compose(x**2 + x, x - 1)
  635. x**2 - x
  636. """
  637. if len(g) <= 1:
  638. return dup_strip([dup_eval(f, dup_LC(g, K), K)])
  639. if not f:
  640. return []
  641. h = [f[0]]
  642. for c in f[1:]:
  643. h = dup_mul(h, g, K)
  644. h = dup_add_term(h, c, 0, K)
  645. return h
  646. def dmp_compose(f, g, u, K):
  647. """
  648. Evaluate functional composition ``f(g)`` in ``K[X]``.
  649. Examples
  650. ========
  651. >>> from sympy.polys import ring, ZZ
  652. >>> R, x,y = ring("x,y", ZZ)
  653. >>> R.dmp_compose(x*y + 2*x + y, y)
  654. y**2 + 3*y
  655. """
  656. if not u:
  657. return dup_compose(f, g, K)
  658. if dmp_zero_p(f, u):
  659. return f
  660. h = [f[0]]
  661. for c in f[1:]:
  662. h = dmp_mul(h, g, u, K)
  663. h = dmp_add_term(h, c, 0, u, K)
  664. return h
  665. def _dup_right_decompose(f, s, K):
  666. """Helper function for :func:`_dup_decompose`."""
  667. n = len(f) - 1
  668. lc = dup_LC(f, K)
  669. f = dup_to_raw_dict(f)
  670. g = { s: K.one }
  671. r = n // s
  672. for i in range(1, s):
  673. coeff = K.zero
  674. for j in range(0, i):
  675. if not n + j - i in f:
  676. continue
  677. if not s - j in g:
  678. continue
  679. fc, gc = f[n + j - i], g[s - j]
  680. coeff += (i - r*j)*fc*gc
  681. g[s - i] = K.quo(coeff, i*r*lc)
  682. return dup_from_raw_dict(g, K)
  683. def _dup_left_decompose(f, h, K):
  684. """Helper function for :func:`_dup_decompose`."""
  685. g, i = {}, 0
  686. while f:
  687. q, r = dup_div(f, h, K)
  688. if dup_degree(r) > 0:
  689. return None
  690. else:
  691. g[i] = dup_LC(r, K)
  692. f, i = q, i + 1
  693. return dup_from_raw_dict(g, K)
  694. def _dup_decompose(f, K):
  695. """Helper function for :func:`dup_decompose`."""
  696. df = len(f) - 1
  697. for s in range(2, df):
  698. if df % s != 0:
  699. continue
  700. h = _dup_right_decompose(f, s, K)
  701. if h is not None:
  702. g = _dup_left_decompose(f, h, K)
  703. if g is not None:
  704. return g, h
  705. return None
  706. def dup_decompose(f, K):
  707. """
  708. Computes functional decomposition of ``f`` in ``K[x]``.
  709. Given a univariate polynomial ``f`` with coefficients in a field of
  710. characteristic zero, returns list ``[f_1, f_2, ..., f_n]``, where::
  711. f = f_1 o f_2 o ... f_n = f_1(f_2(... f_n))
  712. and ``f_2, ..., f_n`` are monic and homogeneous polynomials of at
  713. least second degree.
  714. Unlike factorization, complete functional decompositions of
  715. polynomials are not unique, consider examples:
  716. 1. ``f o g = f(x + b) o (g - b)``
  717. 2. ``x**n o x**m = x**m o x**n``
  718. 3. ``T_n o T_m = T_m o T_n``
  719. where ``T_n`` and ``T_m`` are Chebyshev polynomials.
  720. Examples
  721. ========
  722. >>> from sympy.polys import ring, ZZ
  723. >>> R, x = ring("x", ZZ)
  724. >>> R.dup_decompose(x**4 - 2*x**3 + x**2)
  725. [x**2, x**2 - x]
  726. References
  727. ==========
  728. .. [1] [Kozen89]_
  729. """
  730. F = []
  731. while True:
  732. result = _dup_decompose(f, K)
  733. if result is not None:
  734. f, h = result
  735. F = [h] + F
  736. else:
  737. break
  738. return [f] + F
  739. def dmp_lift(f, u, K):
  740. """
  741. Convert algebraic coefficients to integers in ``K[X]``.
  742. Examples
  743. ========
  744. >>> from sympy.polys import ring, QQ
  745. >>> from sympy import I
  746. >>> K = QQ.algebraic_field(I)
  747. >>> R, x = ring("x", K)
  748. >>> f = x**2 + K([QQ(1), QQ(0)])*x + K([QQ(2), QQ(0)])
  749. >>> R.dmp_lift(f)
  750. x**8 + 2*x**6 + 9*x**4 - 8*x**2 + 16
  751. """
  752. if K.is_GaussianField:
  753. K1 = K.as_AlgebraicField()
  754. f = dmp_convert(f, u, K, K1)
  755. K = K1
  756. if not K.is_Algebraic:
  757. raise DomainError(
  758. 'computation can be done only in an algebraic domain')
  759. F, monoms, polys = dmp_to_dict(f, u), [], []
  760. for monom, coeff in F.items():
  761. if not coeff.is_ground:
  762. monoms.append(monom)
  763. perms = variations([-1, 1], len(monoms), repetition=True)
  764. for perm in perms:
  765. G = dict(F)
  766. for sign, monom in zip(perm, monoms):
  767. if sign == -1:
  768. G[monom] = -G[monom]
  769. polys.append(dmp_from_dict(G, u, K))
  770. return dmp_convert(dmp_expand(polys, u, K), u, K, K.dom)
  771. def dup_sign_variations(f, K):
  772. """
  773. Compute the number of sign variations of ``f`` in ``K[x]``.
  774. Examples
  775. ========
  776. >>> from sympy.polys import ring, ZZ
  777. >>> R, x = ring("x", ZZ)
  778. >>> R.dup_sign_variations(x**4 - x**2 - x + 1)
  779. 2
  780. """
  781. prev, k = K.zero, 0
  782. for coeff in f:
  783. if K.is_negative(coeff*prev):
  784. k += 1
  785. if coeff:
  786. prev = coeff
  787. return k
  788. def dup_clear_denoms(f, K0, K1=None, convert=False):
  789. """
  790. Clear denominators, i.e. transform ``K_0`` to ``K_1``.
  791. Examples
  792. ========
  793. >>> from sympy.polys import ring, QQ
  794. >>> R, x = ring("x", QQ)
  795. >>> f = QQ(1,2)*x + QQ(1,3)
  796. >>> R.dup_clear_denoms(f, convert=False)
  797. (6, 3*x + 2)
  798. >>> R.dup_clear_denoms(f, convert=True)
  799. (6, 3*x + 2)
  800. """
  801. if K1 is None:
  802. if K0.has_assoc_Ring:
  803. K1 = K0.get_ring()
  804. else:
  805. K1 = K0
  806. common = K1.one
  807. for c in f:
  808. common = K1.lcm(common, K0.denom(c))
  809. if not K1.is_one(common):
  810. f = dup_mul_ground(f, common, K0)
  811. if not convert:
  812. return common, f
  813. else:
  814. return common, dup_convert(f, K0, K1)
  815. def _rec_clear_denoms(g, v, K0, K1):
  816. """Recursive helper for :func:`dmp_clear_denoms`."""
  817. common = K1.one
  818. if not v:
  819. for c in g:
  820. common = K1.lcm(common, K0.denom(c))
  821. else:
  822. w = v - 1
  823. for c in g:
  824. common = K1.lcm(common, _rec_clear_denoms(c, w, K0, K1))
  825. return common
  826. def dmp_clear_denoms(f, u, K0, K1=None, convert=False):
  827. """
  828. Clear denominators, i.e. transform ``K_0`` to ``K_1``.
  829. Examples
  830. ========
  831. >>> from sympy.polys import ring, QQ
  832. >>> R, x,y = ring("x,y", QQ)
  833. >>> f = QQ(1,2)*x + QQ(1,3)*y + 1
  834. >>> R.dmp_clear_denoms(f, convert=False)
  835. (6, 3*x + 2*y + 6)
  836. >>> R.dmp_clear_denoms(f, convert=True)
  837. (6, 3*x + 2*y + 6)
  838. """
  839. if not u:
  840. return dup_clear_denoms(f, K0, K1, convert=convert)
  841. if K1 is None:
  842. if K0.has_assoc_Ring:
  843. K1 = K0.get_ring()
  844. else:
  845. K1 = K0
  846. common = _rec_clear_denoms(f, u, K0, K1)
  847. if not K1.is_one(common):
  848. f = dmp_mul_ground(f, common, u, K0)
  849. if not convert:
  850. return common, f
  851. else:
  852. return common, dmp_convert(f, u, K0, K1)
  853. def dup_revert(f, n, K):
  854. """
  855. Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.
  856. This function computes first ``2**n`` terms of a polynomial that
  857. is a result of inversion of a polynomial modulo ``x**n``. This is
  858. useful to efficiently compute series expansion of ``1/f``.
  859. Examples
  860. ========
  861. >>> from sympy.polys import ring, QQ
  862. >>> R, x = ring("x", QQ)
  863. >>> f = -QQ(1,720)*x**6 + QQ(1,24)*x**4 - QQ(1,2)*x**2 + 1
  864. >>> R.dup_revert(f, 8)
  865. 61/720*x**6 + 5/24*x**4 + 1/2*x**2 + 1
  866. """
  867. g = [K.revert(dup_TC(f, K))]
  868. h = [K.one, K.zero, K.zero]
  869. N = int(_ceil(_log(n, 2)))
  870. for i in range(1, N + 1):
  871. a = dup_mul_ground(g, K(2), K)
  872. b = dup_mul(f, dup_sqr(g, K), K)
  873. g = dup_rem(dup_sub(a, b, K), h, K)
  874. h = dup_lshift(h, dup_degree(h), K)
  875. return g
  876. def dmp_revert(f, g, u, K):
  877. """
  878. Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.
  879. Examples
  880. ========
  881. >>> from sympy.polys import ring, QQ
  882. >>> R, x,y = ring("x,y", QQ)
  883. """
  884. if not u:
  885. return dup_revert(f, g, K)
  886. else:
  887. raise MultivariatePolynomialError(f, g)