You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1310 lines
25 KiB
1310 lines
25 KiB
5 months ago
|
"""Advanced tools for dense recursive polynomials in ``K[x]`` or ``K[X]``. """
|
||
|
|
||
|
|
||
|
from sympy.polys.densearith import (
|
||
|
dup_add_term, dmp_add_term,
|
||
|
dup_lshift,
|
||
|
dup_add, dmp_add,
|
||
|
dup_sub, dmp_sub,
|
||
|
dup_mul, dmp_mul,
|
||
|
dup_sqr,
|
||
|
dup_div,
|
||
|
dup_rem, dmp_rem,
|
||
|
dmp_expand,
|
||
|
dup_mul_ground, dmp_mul_ground,
|
||
|
dup_quo_ground, dmp_quo_ground,
|
||
|
dup_exquo_ground, dmp_exquo_ground,
|
||
|
)
|
||
|
from sympy.polys.densebasic import (
|
||
|
dup_strip, dmp_strip,
|
||
|
dup_convert, dmp_convert,
|
||
|
dup_degree, dmp_degree,
|
||
|
dmp_to_dict,
|
||
|
dmp_from_dict,
|
||
|
dup_LC, dmp_LC, dmp_ground_LC,
|
||
|
dup_TC, dmp_TC,
|
||
|
dmp_zero, dmp_ground,
|
||
|
dmp_zero_p,
|
||
|
dup_to_raw_dict, dup_from_raw_dict,
|
||
|
dmp_zeros
|
||
|
)
|
||
|
from sympy.polys.polyerrors import (
|
||
|
MultivariatePolynomialError,
|
||
|
DomainError
|
||
|
)
|
||
|
from sympy.utilities import variations
|
||
|
|
||
|
from math import ceil as _ceil, log as _log
|
||
|
|
||
|
def dup_integrate(f, m, K):
|
||
|
"""
|
||
|
Computes the indefinite integral of ``f`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
|
||
|
>>> R.dup_integrate(x**2 + 2*x, 1)
|
||
|
1/3*x**3 + x**2
|
||
|
>>> R.dup_integrate(x**2 + 2*x, 2)
|
||
|
1/12*x**4 + 1/3*x**3
|
||
|
|
||
|
"""
|
||
|
if m <= 0 or not f:
|
||
|
return f
|
||
|
|
||
|
g = [K.zero]*m
|
||
|
|
||
|
for i, c in enumerate(reversed(f)):
|
||
|
n = i + 1
|
||
|
|
||
|
for j in range(1, m):
|
||
|
n *= i + j + 1
|
||
|
|
||
|
g.insert(0, K.exquo(c, K(n)))
|
||
|
|
||
|
return g
|
||
|
|
||
|
|
||
|
def dmp_integrate(f, m, u, K):
|
||
|
"""
|
||
|
Computes the indefinite integral of ``f`` in ``x_0`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
|
||
|
>>> R.dmp_integrate(x + 2*y, 1)
|
||
|
1/2*x**2 + 2*x*y
|
||
|
>>> R.dmp_integrate(x + 2*y, 2)
|
||
|
1/6*x**3 + x**2*y
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_integrate(f, m, K)
|
||
|
|
||
|
if m <= 0 or dmp_zero_p(f, u):
|
||
|
return f
|
||
|
|
||
|
g, v = dmp_zeros(m, u - 1, K), u - 1
|
||
|
|
||
|
for i, c in enumerate(reversed(f)):
|
||
|
n = i + 1
|
||
|
|
||
|
for j in range(1, m):
|
||
|
n *= i + j + 1
|
||
|
|
||
|
g.insert(0, dmp_quo_ground(c, K(n), v, K))
|
||
|
|
||
|
return g
|
||
|
|
||
|
|
||
|
def _rec_integrate_in(g, m, v, i, j, K):
|
||
|
"""Recursive helper for :func:`dmp_integrate_in`."""
|
||
|
if i == j:
|
||
|
return dmp_integrate(g, m, v, K)
|
||
|
|
||
|
w, i = v - 1, i + 1
|
||
|
|
||
|
return dmp_strip([ _rec_integrate_in(c, m, w, i, j, K) for c in g ], v)
|
||
|
|
||
|
|
||
|
def dmp_integrate_in(f, m, j, u, K):
|
||
|
"""
|
||
|
Computes the indefinite integral of ``f`` in ``x_j`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
|
||
|
>>> R.dmp_integrate_in(x + 2*y, 1, 0)
|
||
|
1/2*x**2 + 2*x*y
|
||
|
>>> R.dmp_integrate_in(x + 2*y, 1, 1)
|
||
|
x*y + y**2
|
||
|
|
||
|
"""
|
||
|
if j < 0 or j > u:
|
||
|
raise IndexError("0 <= j <= u expected, got u = %d, j = %d" % (u, j))
|
||
|
|
||
|
return _rec_integrate_in(f, m, u, 0, j, K)
|
||
|
|
||
|
|
||
|
def dup_diff(f, m, K):
|
||
|
"""
|
||
|
``m``-th order derivative of a polynomial in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 1)
|
||
|
3*x**2 + 4*x + 3
|
||
|
>>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 2)
|
||
|
6*x + 4
|
||
|
|
||
|
"""
|
||
|
if m <= 0:
|
||
|
return f
|
||
|
|
||
|
n = dup_degree(f)
|
||
|
|
||
|
if n < m:
|
||
|
return []
|
||
|
|
||
|
deriv = []
|
||
|
|
||
|
if m == 1:
|
||
|
for coeff in f[:-m]:
|
||
|
deriv.append(K(n)*coeff)
|
||
|
n -= 1
|
||
|
else:
|
||
|
for coeff in f[:-m]:
|
||
|
k = n
|
||
|
|
||
|
for i in range(n - 1, n - m, -1):
|
||
|
k *= i
|
||
|
|
||
|
deriv.append(K(k)*coeff)
|
||
|
n -= 1
|
||
|
|
||
|
return dup_strip(deriv)
|
||
|
|
||
|
|
||
|
def dmp_diff(f, m, u, K):
|
||
|
"""
|
||
|
``m``-th order derivative in ``x_0`` of a polynomial in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
|
||
|
|
||
|
>>> R.dmp_diff(f, 1)
|
||
|
y**2 + 2*y + 3
|
||
|
>>> R.dmp_diff(f, 2)
|
||
|
0
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_diff(f, m, K)
|
||
|
if m <= 0:
|
||
|
return f
|
||
|
|
||
|
n = dmp_degree(f, u)
|
||
|
|
||
|
if n < m:
|
||
|
return dmp_zero(u)
|
||
|
|
||
|
deriv, v = [], u - 1
|
||
|
|
||
|
if m == 1:
|
||
|
for coeff in f[:-m]:
|
||
|
deriv.append(dmp_mul_ground(coeff, K(n), v, K))
|
||
|
n -= 1
|
||
|
else:
|
||
|
for coeff in f[:-m]:
|
||
|
k = n
|
||
|
|
||
|
for i in range(n - 1, n - m, -1):
|
||
|
k *= i
|
||
|
|
||
|
deriv.append(dmp_mul_ground(coeff, K(k), v, K))
|
||
|
n -= 1
|
||
|
|
||
|
return dmp_strip(deriv, u)
|
||
|
|
||
|
|
||
|
def _rec_diff_in(g, m, v, i, j, K):
|
||
|
"""Recursive helper for :func:`dmp_diff_in`."""
|
||
|
if i == j:
|
||
|
return dmp_diff(g, m, v, K)
|
||
|
|
||
|
w, i = v - 1, i + 1
|
||
|
|
||
|
return dmp_strip([ _rec_diff_in(c, m, w, i, j, K) for c in g ], v)
|
||
|
|
||
|
|
||
|
def dmp_diff_in(f, m, j, u, K):
|
||
|
"""
|
||
|
``m``-th order derivative in ``x_j`` of a polynomial in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
|
||
|
|
||
|
>>> R.dmp_diff_in(f, 1, 0)
|
||
|
y**2 + 2*y + 3
|
||
|
>>> R.dmp_diff_in(f, 1, 1)
|
||
|
2*x*y + 2*x + 4*y + 3
|
||
|
|
||
|
"""
|
||
|
if j < 0 or j > u:
|
||
|
raise IndexError("0 <= j <= %s expected, got %s" % (u, j))
|
||
|
|
||
|
return _rec_diff_in(f, m, u, 0, j, K)
|
||
|
|
||
|
|
||
|
def dup_eval(f, a, K):
|
||
|
"""
|
||
|
Evaluate a polynomial at ``x = a`` in ``K[x]`` using Horner scheme.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_eval(x**2 + 2*x + 3, 2)
|
||
|
11
|
||
|
|
||
|
"""
|
||
|
if not a:
|
||
|
return K.convert(dup_TC(f, K))
|
||
|
|
||
|
result = K.zero
|
||
|
|
||
|
for c in f:
|
||
|
result *= a
|
||
|
result += c
|
||
|
|
||
|
return result
|
||
|
|
||
|
|
||
|
def dmp_eval(f, a, u, K):
|
||
|
"""
|
||
|
Evaluate a polynomial at ``x_0 = a`` in ``K[X]`` using the Horner scheme.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> R.dmp_eval(2*x*y + 3*x + y + 2, 2)
|
||
|
5*y + 8
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_eval(f, a, K)
|
||
|
|
||
|
if not a:
|
||
|
return dmp_TC(f, K)
|
||
|
|
||
|
result, v = dmp_LC(f, K), u - 1
|
||
|
|
||
|
for coeff in f[1:]:
|
||
|
result = dmp_mul_ground(result, a, v, K)
|
||
|
result = dmp_add(result, coeff, v, K)
|
||
|
|
||
|
return result
|
||
|
|
||
|
|
||
|
def _rec_eval_in(g, a, v, i, j, K):
|
||
|
"""Recursive helper for :func:`dmp_eval_in`."""
|
||
|
if i == j:
|
||
|
return dmp_eval(g, a, v, K)
|
||
|
|
||
|
v, i = v - 1, i + 1
|
||
|
|
||
|
return dmp_strip([ _rec_eval_in(c, a, v, i, j, K) for c in g ], v)
|
||
|
|
||
|
|
||
|
def dmp_eval_in(f, a, j, u, K):
|
||
|
"""
|
||
|
Evaluate a polynomial at ``x_j = a`` in ``K[X]`` using the Horner scheme.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = 2*x*y + 3*x + y + 2
|
||
|
|
||
|
>>> R.dmp_eval_in(f, 2, 0)
|
||
|
5*y + 8
|
||
|
>>> R.dmp_eval_in(f, 2, 1)
|
||
|
7*x + 4
|
||
|
|
||
|
"""
|
||
|
if j < 0 or j > u:
|
||
|
raise IndexError("0 <= j <= %s expected, got %s" % (u, j))
|
||
|
|
||
|
return _rec_eval_in(f, a, u, 0, j, K)
|
||
|
|
||
|
|
||
|
def _rec_eval_tail(g, i, A, u, K):
|
||
|
"""Recursive helper for :func:`dmp_eval_tail`."""
|
||
|
if i == u:
|
||
|
return dup_eval(g, A[-1], K)
|
||
|
else:
|
||
|
h = [ _rec_eval_tail(c, i + 1, A, u, K) for c in g ]
|
||
|
|
||
|
if i < u - len(A) + 1:
|
||
|
return h
|
||
|
else:
|
||
|
return dup_eval(h, A[-u + i - 1], K)
|
||
|
|
||
|
|
||
|
def dmp_eval_tail(f, A, u, K):
|
||
|
"""
|
||
|
Evaluate a polynomial at ``x_j = a_j, ...`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = 2*x*y + 3*x + y + 2
|
||
|
|
||
|
>>> R.dmp_eval_tail(f, [2])
|
||
|
7*x + 4
|
||
|
>>> R.dmp_eval_tail(f, [2, 2])
|
||
|
18
|
||
|
|
||
|
"""
|
||
|
if not A:
|
||
|
return f
|
||
|
|
||
|
if dmp_zero_p(f, u):
|
||
|
return dmp_zero(u - len(A))
|
||
|
|
||
|
e = _rec_eval_tail(f, 0, A, u, K)
|
||
|
|
||
|
if u == len(A) - 1:
|
||
|
return e
|
||
|
else:
|
||
|
return dmp_strip(e, u - len(A))
|
||
|
|
||
|
|
||
|
def _rec_diff_eval(g, m, a, v, i, j, K):
|
||
|
"""Recursive helper for :func:`dmp_diff_eval`."""
|
||
|
if i == j:
|
||
|
return dmp_eval(dmp_diff(g, m, v, K), a, v, K)
|
||
|
|
||
|
v, i = v - 1, i + 1
|
||
|
|
||
|
return dmp_strip([ _rec_diff_eval(c, m, a, v, i, j, K) for c in g ], v)
|
||
|
|
||
|
|
||
|
def dmp_diff_eval_in(f, m, a, j, u, K):
|
||
|
"""
|
||
|
Differentiate and evaluate a polynomial in ``x_j`` at ``a`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1
|
||
|
|
||
|
>>> R.dmp_diff_eval_in(f, 1, 2, 0)
|
||
|
y**2 + 2*y + 3
|
||
|
>>> R.dmp_diff_eval_in(f, 1, 2, 1)
|
||
|
6*x + 11
|
||
|
|
||
|
"""
|
||
|
if j > u:
|
||
|
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
|
||
|
if not j:
|
||
|
return dmp_eval(dmp_diff(f, m, u, K), a, u, K)
|
||
|
|
||
|
return _rec_diff_eval(f, m, a, u, 0, j, K)
|
||
|
|
||
|
|
||
|
def dup_trunc(f, p, K):
|
||
|
"""
|
||
|
Reduce a ``K[x]`` polynomial modulo a constant ``p`` in ``K``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_trunc(2*x**3 + 3*x**2 + 5*x + 7, ZZ(3))
|
||
|
-x**3 - x + 1
|
||
|
|
||
|
"""
|
||
|
if K.is_ZZ:
|
||
|
g = []
|
||
|
|
||
|
for c in f:
|
||
|
c = c % p
|
||
|
|
||
|
if c > p // 2:
|
||
|
g.append(c - p)
|
||
|
else:
|
||
|
g.append(c)
|
||
|
else:
|
||
|
g = [ c % p for c in f ]
|
||
|
|
||
|
return dup_strip(g)
|
||
|
|
||
|
|
||
|
def dmp_trunc(f, p, u, K):
|
||
|
"""
|
||
|
Reduce a ``K[X]`` polynomial modulo a polynomial ``p`` in ``K[Y]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
|
||
|
>>> g = (y - 1).drop(x)
|
||
|
|
||
|
>>> R.dmp_trunc(f, g)
|
||
|
11*x**2 + 11*x + 5
|
||
|
|
||
|
"""
|
||
|
return dmp_strip([ dmp_rem(c, p, u - 1, K) for c in f ], u)
|
||
|
|
||
|
|
||
|
def dmp_ground_trunc(f, p, u, K):
|
||
|
"""
|
||
|
Reduce a ``K[X]`` polynomial modulo a constant ``p`` in ``K``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
|
||
|
|
||
|
>>> R.dmp_ground_trunc(f, ZZ(3))
|
||
|
-x**2 - x*y - y
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_trunc(f, p, K)
|
||
|
|
||
|
v = u - 1
|
||
|
|
||
|
return dmp_strip([ dmp_ground_trunc(c, p, v, K) for c in f ], u)
|
||
|
|
||
|
|
||
|
def dup_monic(f, K):
|
||
|
"""
|
||
|
Divide all coefficients by ``LC(f)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
>>> R.dup_monic(3*x**2 + 6*x + 9)
|
||
|
x**2 + 2*x + 3
|
||
|
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
>>> R.dup_monic(3*x**2 + 4*x + 2)
|
||
|
x**2 + 4/3*x + 2/3
|
||
|
|
||
|
"""
|
||
|
if not f:
|
||
|
return f
|
||
|
|
||
|
lc = dup_LC(f, K)
|
||
|
|
||
|
if K.is_one(lc):
|
||
|
return f
|
||
|
else:
|
||
|
return dup_exquo_ground(f, lc, K)
|
||
|
|
||
|
|
||
|
def dmp_ground_monic(f, u, K):
|
||
|
"""
|
||
|
Divide all coefficients by ``LC(f)`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
>>> f = 3*x**2*y + 6*x**2 + 3*x*y + 9*y + 3
|
||
|
|
||
|
>>> R.dmp_ground_monic(f)
|
||
|
x**2*y + 2*x**2 + x*y + 3*y + 1
|
||
|
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
>>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3
|
||
|
|
||
|
>>> R.dmp_ground_monic(f)
|
||
|
x**2*y + 8/3*x**2 + 5/3*x*y + 2*x + 2/3*y + 1
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_monic(f, K)
|
||
|
|
||
|
if dmp_zero_p(f, u):
|
||
|
return f
|
||
|
|
||
|
lc = dmp_ground_LC(f, u, K)
|
||
|
|
||
|
if K.is_one(lc):
|
||
|
return f
|
||
|
else:
|
||
|
return dmp_exquo_ground(f, lc, u, K)
|
||
|
|
||
|
|
||
|
def dup_content(f, K):
|
||
|
"""
|
||
|
Compute the GCD of coefficients of ``f`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
>>> f = 6*x**2 + 8*x + 12
|
||
|
|
||
|
>>> R.dup_content(f)
|
||
|
2
|
||
|
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
>>> f = 6*x**2 + 8*x + 12
|
||
|
|
||
|
>>> R.dup_content(f)
|
||
|
2
|
||
|
|
||
|
"""
|
||
|
from sympy.polys.domains import QQ
|
||
|
|
||
|
if not f:
|
||
|
return K.zero
|
||
|
|
||
|
cont = K.zero
|
||
|
|
||
|
if K == QQ:
|
||
|
for c in f:
|
||
|
cont = K.gcd(cont, c)
|
||
|
else:
|
||
|
for c in f:
|
||
|
cont = K.gcd(cont, c)
|
||
|
|
||
|
if K.is_one(cont):
|
||
|
break
|
||
|
|
||
|
return cont
|
||
|
|
||
|
|
||
|
def dmp_ground_content(f, u, K):
|
||
|
"""
|
||
|
Compute the GCD of coefficients of ``f`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
>>> f = 2*x*y + 6*x + 4*y + 12
|
||
|
|
||
|
>>> R.dmp_ground_content(f)
|
||
|
2
|
||
|
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
>>> f = 2*x*y + 6*x + 4*y + 12
|
||
|
|
||
|
>>> R.dmp_ground_content(f)
|
||
|
2
|
||
|
|
||
|
"""
|
||
|
from sympy.polys.domains import QQ
|
||
|
|
||
|
if not u:
|
||
|
return dup_content(f, K)
|
||
|
|
||
|
if dmp_zero_p(f, u):
|
||
|
return K.zero
|
||
|
|
||
|
cont, v = K.zero, u - 1
|
||
|
|
||
|
if K == QQ:
|
||
|
for c in f:
|
||
|
cont = K.gcd(cont, dmp_ground_content(c, v, K))
|
||
|
else:
|
||
|
for c in f:
|
||
|
cont = K.gcd(cont, dmp_ground_content(c, v, K))
|
||
|
|
||
|
if K.is_one(cont):
|
||
|
break
|
||
|
|
||
|
return cont
|
||
|
|
||
|
|
||
|
def dup_primitive(f, K):
|
||
|
"""
|
||
|
Compute content and the primitive form of ``f`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
>>> f = 6*x**2 + 8*x + 12
|
||
|
|
||
|
>>> R.dup_primitive(f)
|
||
|
(2, 3*x**2 + 4*x + 6)
|
||
|
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
>>> f = 6*x**2 + 8*x + 12
|
||
|
|
||
|
>>> R.dup_primitive(f)
|
||
|
(2, 3*x**2 + 4*x + 6)
|
||
|
|
||
|
"""
|
||
|
if not f:
|
||
|
return K.zero, f
|
||
|
|
||
|
cont = dup_content(f, K)
|
||
|
|
||
|
if K.is_one(cont):
|
||
|
return cont, f
|
||
|
else:
|
||
|
return cont, dup_quo_ground(f, cont, K)
|
||
|
|
||
|
|
||
|
def dmp_ground_primitive(f, u, K):
|
||
|
"""
|
||
|
Compute content and the primitive form of ``f`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ, QQ
|
||
|
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
>>> f = 2*x*y + 6*x + 4*y + 12
|
||
|
|
||
|
>>> R.dmp_ground_primitive(f)
|
||
|
(2, x*y + 3*x + 2*y + 6)
|
||
|
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
>>> f = 2*x*y + 6*x + 4*y + 12
|
||
|
|
||
|
>>> R.dmp_ground_primitive(f)
|
||
|
(2, x*y + 3*x + 2*y + 6)
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_primitive(f, K)
|
||
|
|
||
|
if dmp_zero_p(f, u):
|
||
|
return K.zero, f
|
||
|
|
||
|
cont = dmp_ground_content(f, u, K)
|
||
|
|
||
|
if K.is_one(cont):
|
||
|
return cont, f
|
||
|
else:
|
||
|
return cont, dmp_quo_ground(f, cont, u, K)
|
||
|
|
||
|
|
||
|
def dup_extract(f, g, K):
|
||
|
"""
|
||
|
Extract common content from a pair of polynomials in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_extract(6*x**2 + 12*x + 18, 4*x**2 + 8*x + 12)
|
||
|
(2, 3*x**2 + 6*x + 9, 2*x**2 + 4*x + 6)
|
||
|
|
||
|
"""
|
||
|
fc = dup_content(f, K)
|
||
|
gc = dup_content(g, K)
|
||
|
|
||
|
gcd = K.gcd(fc, gc)
|
||
|
|
||
|
if not K.is_one(gcd):
|
||
|
f = dup_quo_ground(f, gcd, K)
|
||
|
g = dup_quo_ground(g, gcd, K)
|
||
|
|
||
|
return gcd, f, g
|
||
|
|
||
|
|
||
|
def dmp_ground_extract(f, g, u, K):
|
||
|
"""
|
||
|
Extract common content from a pair of polynomials in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> R.dmp_ground_extract(6*x*y + 12*x + 18, 4*x*y + 8*x + 12)
|
||
|
(2, 3*x*y + 6*x + 9, 2*x*y + 4*x + 6)
|
||
|
|
||
|
"""
|
||
|
fc = dmp_ground_content(f, u, K)
|
||
|
gc = dmp_ground_content(g, u, K)
|
||
|
|
||
|
gcd = K.gcd(fc, gc)
|
||
|
|
||
|
if not K.is_one(gcd):
|
||
|
f = dmp_quo_ground(f, gcd, u, K)
|
||
|
g = dmp_quo_ground(g, gcd, u, K)
|
||
|
|
||
|
return gcd, f, g
|
||
|
|
||
|
|
||
|
def dup_real_imag(f, K):
|
||
|
"""
|
||
|
Return bivariate polynomials ``f1`` and ``f2``, such that ``f = f1 + f2*I``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> R.dup_real_imag(x**3 + x**2 + x + 1)
|
||
|
(x**3 + x**2 - 3*x*y**2 + x - y**2 + 1, 3*x**2*y + 2*x*y - y**3 + y)
|
||
|
|
||
|
"""
|
||
|
if not K.is_ZZ and not K.is_QQ:
|
||
|
raise DomainError("computing real and imaginary parts is not supported over %s" % K)
|
||
|
|
||
|
f1 = dmp_zero(1)
|
||
|
f2 = dmp_zero(1)
|
||
|
|
||
|
if not f:
|
||
|
return f1, f2
|
||
|
|
||
|
g = [[[K.one, K.zero]], [[K.one], []]]
|
||
|
h = dmp_ground(f[0], 2)
|
||
|
|
||
|
for c in f[1:]:
|
||
|
h = dmp_mul(h, g, 2, K)
|
||
|
h = dmp_add_term(h, dmp_ground(c, 1), 0, 2, K)
|
||
|
|
||
|
H = dup_to_raw_dict(h)
|
||
|
|
||
|
for k, h in H.items():
|
||
|
m = k % 4
|
||
|
|
||
|
if not m:
|
||
|
f1 = dmp_add(f1, h, 1, K)
|
||
|
elif m == 1:
|
||
|
f2 = dmp_add(f2, h, 1, K)
|
||
|
elif m == 2:
|
||
|
f1 = dmp_sub(f1, h, 1, K)
|
||
|
else:
|
||
|
f2 = dmp_sub(f2, h, 1, K)
|
||
|
|
||
|
return f1, f2
|
||
|
|
||
|
|
||
|
def dup_mirror(f, K):
|
||
|
"""
|
||
|
Evaluate efficiently the composition ``f(-x)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_mirror(x**3 + 2*x**2 - 4*x + 2)
|
||
|
-x**3 + 2*x**2 + 4*x + 2
|
||
|
|
||
|
"""
|
||
|
f = list(f)
|
||
|
|
||
|
for i in range(len(f) - 2, -1, -2):
|
||
|
f[i] = -f[i]
|
||
|
|
||
|
return f
|
||
|
|
||
|
|
||
|
def dup_scale(f, a, K):
|
||
|
"""
|
||
|
Evaluate efficiently composition ``f(a*x)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_scale(x**2 - 2*x + 1, ZZ(2))
|
||
|
4*x**2 - 4*x + 1
|
||
|
|
||
|
"""
|
||
|
f, n, b = list(f), len(f) - 1, a
|
||
|
|
||
|
for i in range(n - 1, -1, -1):
|
||
|
f[i], b = b*f[i], b*a
|
||
|
|
||
|
return f
|
||
|
|
||
|
|
||
|
def dup_shift(f, a, K):
|
||
|
"""
|
||
|
Evaluate efficiently Taylor shift ``f(x + a)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_shift(x**2 - 2*x + 1, ZZ(2))
|
||
|
x**2 + 2*x + 1
|
||
|
|
||
|
"""
|
||
|
f, n = list(f), len(f) - 1
|
||
|
|
||
|
for i in range(n, 0, -1):
|
||
|
for j in range(0, i):
|
||
|
f[j + 1] += a*f[j]
|
||
|
|
||
|
return f
|
||
|
|
||
|
|
||
|
def dup_transform(f, p, q, K):
|
||
|
"""
|
||
|
Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_transform(x**2 - 2*x + 1, x**2 + 1, x - 1)
|
||
|
x**4 - 2*x**3 + 5*x**2 - 4*x + 4
|
||
|
|
||
|
"""
|
||
|
if not f:
|
||
|
return []
|
||
|
|
||
|
n = len(f) - 1
|
||
|
h, Q = [f[0]], [[K.one]]
|
||
|
|
||
|
for i in range(0, n):
|
||
|
Q.append(dup_mul(Q[-1], q, K))
|
||
|
|
||
|
for c, q in zip(f[1:], Q[1:]):
|
||
|
h = dup_mul(h, p, K)
|
||
|
q = dup_mul_ground(q, c, K)
|
||
|
h = dup_add(h, q, K)
|
||
|
|
||
|
return h
|
||
|
|
||
|
|
||
|
def dup_compose(f, g, K):
|
||
|
"""
|
||
|
Evaluate functional composition ``f(g)`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_compose(x**2 + x, x - 1)
|
||
|
x**2 - x
|
||
|
|
||
|
"""
|
||
|
if len(g) <= 1:
|
||
|
return dup_strip([dup_eval(f, dup_LC(g, K), K)])
|
||
|
|
||
|
if not f:
|
||
|
return []
|
||
|
|
||
|
h = [f[0]]
|
||
|
|
||
|
for c in f[1:]:
|
||
|
h = dup_mul(h, g, K)
|
||
|
h = dup_add_term(h, c, 0, K)
|
||
|
|
||
|
return h
|
||
|
|
||
|
|
||
|
def dmp_compose(f, g, u, K):
|
||
|
"""
|
||
|
Evaluate functional composition ``f(g)`` in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x,y = ring("x,y", ZZ)
|
||
|
|
||
|
>>> R.dmp_compose(x*y + 2*x + y, y)
|
||
|
y**2 + 3*y
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_compose(f, g, K)
|
||
|
|
||
|
if dmp_zero_p(f, u):
|
||
|
return f
|
||
|
|
||
|
h = [f[0]]
|
||
|
|
||
|
for c in f[1:]:
|
||
|
h = dmp_mul(h, g, u, K)
|
||
|
h = dmp_add_term(h, c, 0, u, K)
|
||
|
|
||
|
return h
|
||
|
|
||
|
|
||
|
def _dup_right_decompose(f, s, K):
|
||
|
"""Helper function for :func:`_dup_decompose`."""
|
||
|
n = len(f) - 1
|
||
|
lc = dup_LC(f, K)
|
||
|
|
||
|
f = dup_to_raw_dict(f)
|
||
|
g = { s: K.one }
|
||
|
|
||
|
r = n // s
|
||
|
|
||
|
for i in range(1, s):
|
||
|
coeff = K.zero
|
||
|
|
||
|
for j in range(0, i):
|
||
|
if not n + j - i in f:
|
||
|
continue
|
||
|
|
||
|
if not s - j in g:
|
||
|
continue
|
||
|
|
||
|
fc, gc = f[n + j - i], g[s - j]
|
||
|
coeff += (i - r*j)*fc*gc
|
||
|
|
||
|
g[s - i] = K.quo(coeff, i*r*lc)
|
||
|
|
||
|
return dup_from_raw_dict(g, K)
|
||
|
|
||
|
|
||
|
def _dup_left_decompose(f, h, K):
|
||
|
"""Helper function for :func:`_dup_decompose`."""
|
||
|
g, i = {}, 0
|
||
|
|
||
|
while f:
|
||
|
q, r = dup_div(f, h, K)
|
||
|
|
||
|
if dup_degree(r) > 0:
|
||
|
return None
|
||
|
else:
|
||
|
g[i] = dup_LC(r, K)
|
||
|
f, i = q, i + 1
|
||
|
|
||
|
return dup_from_raw_dict(g, K)
|
||
|
|
||
|
|
||
|
def _dup_decompose(f, K):
|
||
|
"""Helper function for :func:`dup_decompose`."""
|
||
|
df = len(f) - 1
|
||
|
|
||
|
for s in range(2, df):
|
||
|
if df % s != 0:
|
||
|
continue
|
||
|
|
||
|
h = _dup_right_decompose(f, s, K)
|
||
|
|
||
|
if h is not None:
|
||
|
g = _dup_left_decompose(f, h, K)
|
||
|
|
||
|
if g is not None:
|
||
|
return g, h
|
||
|
|
||
|
return None
|
||
|
|
||
|
|
||
|
def dup_decompose(f, K):
|
||
|
"""
|
||
|
Computes functional decomposition of ``f`` in ``K[x]``.
|
||
|
|
||
|
Given a univariate polynomial ``f`` with coefficients in a field of
|
||
|
characteristic zero, returns list ``[f_1, f_2, ..., f_n]``, where::
|
||
|
|
||
|
f = f_1 o f_2 o ... f_n = f_1(f_2(... f_n))
|
||
|
|
||
|
and ``f_2, ..., f_n`` are monic and homogeneous polynomials of at
|
||
|
least second degree.
|
||
|
|
||
|
Unlike factorization, complete functional decompositions of
|
||
|
polynomials are not unique, consider examples:
|
||
|
|
||
|
1. ``f o g = f(x + b) o (g - b)``
|
||
|
2. ``x**n o x**m = x**m o x**n``
|
||
|
3. ``T_n o T_m = T_m o T_n``
|
||
|
|
||
|
where ``T_n`` and ``T_m`` are Chebyshev polynomials.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_decompose(x**4 - 2*x**3 + x**2)
|
||
|
[x**2, x**2 - x]
|
||
|
|
||
|
References
|
||
|
==========
|
||
|
|
||
|
.. [1] [Kozen89]_
|
||
|
|
||
|
"""
|
||
|
F = []
|
||
|
|
||
|
while True:
|
||
|
result = _dup_decompose(f, K)
|
||
|
|
||
|
if result is not None:
|
||
|
f, h = result
|
||
|
F = [h] + F
|
||
|
else:
|
||
|
break
|
||
|
|
||
|
return [f] + F
|
||
|
|
||
|
|
||
|
def dmp_lift(f, u, K):
|
||
|
"""
|
||
|
Convert algebraic coefficients to integers in ``K[X]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> from sympy import I
|
||
|
|
||
|
>>> K = QQ.algebraic_field(I)
|
||
|
>>> R, x = ring("x", K)
|
||
|
|
||
|
>>> f = x**2 + K([QQ(1), QQ(0)])*x + K([QQ(2), QQ(0)])
|
||
|
|
||
|
>>> R.dmp_lift(f)
|
||
|
x**8 + 2*x**6 + 9*x**4 - 8*x**2 + 16
|
||
|
|
||
|
"""
|
||
|
if K.is_GaussianField:
|
||
|
K1 = K.as_AlgebraicField()
|
||
|
f = dmp_convert(f, u, K, K1)
|
||
|
K = K1
|
||
|
|
||
|
if not K.is_Algebraic:
|
||
|
raise DomainError(
|
||
|
'computation can be done only in an algebraic domain')
|
||
|
|
||
|
F, monoms, polys = dmp_to_dict(f, u), [], []
|
||
|
|
||
|
for monom, coeff in F.items():
|
||
|
if not coeff.is_ground:
|
||
|
monoms.append(monom)
|
||
|
|
||
|
perms = variations([-1, 1], len(monoms), repetition=True)
|
||
|
|
||
|
for perm in perms:
|
||
|
G = dict(F)
|
||
|
|
||
|
for sign, monom in zip(perm, monoms):
|
||
|
if sign == -1:
|
||
|
G[monom] = -G[monom]
|
||
|
|
||
|
polys.append(dmp_from_dict(G, u, K))
|
||
|
|
||
|
return dmp_convert(dmp_expand(polys, u, K), u, K, K.dom)
|
||
|
|
||
|
|
||
|
def dup_sign_variations(f, K):
|
||
|
"""
|
||
|
Compute the number of sign variations of ``f`` in ``K[x]``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, ZZ
|
||
|
>>> R, x = ring("x", ZZ)
|
||
|
|
||
|
>>> R.dup_sign_variations(x**4 - x**2 - x + 1)
|
||
|
2
|
||
|
|
||
|
"""
|
||
|
prev, k = K.zero, 0
|
||
|
|
||
|
for coeff in f:
|
||
|
if K.is_negative(coeff*prev):
|
||
|
k += 1
|
||
|
|
||
|
if coeff:
|
||
|
prev = coeff
|
||
|
|
||
|
return k
|
||
|
|
||
|
|
||
|
def dup_clear_denoms(f, K0, K1=None, convert=False):
|
||
|
"""
|
||
|
Clear denominators, i.e. transform ``K_0`` to ``K_1``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
|
||
|
>>> f = QQ(1,2)*x + QQ(1,3)
|
||
|
|
||
|
>>> R.dup_clear_denoms(f, convert=False)
|
||
|
(6, 3*x + 2)
|
||
|
>>> R.dup_clear_denoms(f, convert=True)
|
||
|
(6, 3*x + 2)
|
||
|
|
||
|
"""
|
||
|
if K1 is None:
|
||
|
if K0.has_assoc_Ring:
|
||
|
K1 = K0.get_ring()
|
||
|
else:
|
||
|
K1 = K0
|
||
|
|
||
|
common = K1.one
|
||
|
|
||
|
for c in f:
|
||
|
common = K1.lcm(common, K0.denom(c))
|
||
|
|
||
|
if not K1.is_one(common):
|
||
|
f = dup_mul_ground(f, common, K0)
|
||
|
|
||
|
if not convert:
|
||
|
return common, f
|
||
|
else:
|
||
|
return common, dup_convert(f, K0, K1)
|
||
|
|
||
|
|
||
|
def _rec_clear_denoms(g, v, K0, K1):
|
||
|
"""Recursive helper for :func:`dmp_clear_denoms`."""
|
||
|
common = K1.one
|
||
|
|
||
|
if not v:
|
||
|
for c in g:
|
||
|
common = K1.lcm(common, K0.denom(c))
|
||
|
else:
|
||
|
w = v - 1
|
||
|
|
||
|
for c in g:
|
||
|
common = K1.lcm(common, _rec_clear_denoms(c, w, K0, K1))
|
||
|
|
||
|
return common
|
||
|
|
||
|
|
||
|
def dmp_clear_denoms(f, u, K0, K1=None, convert=False):
|
||
|
"""
|
||
|
Clear denominators, i.e. transform ``K_0`` to ``K_1``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
|
||
|
>>> f = QQ(1,2)*x + QQ(1,3)*y + 1
|
||
|
|
||
|
>>> R.dmp_clear_denoms(f, convert=False)
|
||
|
(6, 3*x + 2*y + 6)
|
||
|
>>> R.dmp_clear_denoms(f, convert=True)
|
||
|
(6, 3*x + 2*y + 6)
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_clear_denoms(f, K0, K1, convert=convert)
|
||
|
|
||
|
if K1 is None:
|
||
|
if K0.has_assoc_Ring:
|
||
|
K1 = K0.get_ring()
|
||
|
else:
|
||
|
K1 = K0
|
||
|
|
||
|
common = _rec_clear_denoms(f, u, K0, K1)
|
||
|
|
||
|
if not K1.is_one(common):
|
||
|
f = dmp_mul_ground(f, common, u, K0)
|
||
|
|
||
|
if not convert:
|
||
|
return common, f
|
||
|
else:
|
||
|
return common, dmp_convert(f, u, K0, K1)
|
||
|
|
||
|
|
||
|
def dup_revert(f, n, K):
|
||
|
"""
|
||
|
Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.
|
||
|
|
||
|
This function computes first ``2**n`` terms of a polynomial that
|
||
|
is a result of inversion of a polynomial modulo ``x**n``. This is
|
||
|
useful to efficiently compute series expansion of ``1/f``.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x = ring("x", QQ)
|
||
|
|
||
|
>>> f = -QQ(1,720)*x**6 + QQ(1,24)*x**4 - QQ(1,2)*x**2 + 1
|
||
|
|
||
|
>>> R.dup_revert(f, 8)
|
||
|
61/720*x**6 + 5/24*x**4 + 1/2*x**2 + 1
|
||
|
|
||
|
"""
|
||
|
g = [K.revert(dup_TC(f, K))]
|
||
|
h = [K.one, K.zero, K.zero]
|
||
|
|
||
|
N = int(_ceil(_log(n, 2)))
|
||
|
|
||
|
for i in range(1, N + 1):
|
||
|
a = dup_mul_ground(g, K(2), K)
|
||
|
b = dup_mul(f, dup_sqr(g, K), K)
|
||
|
g = dup_rem(dup_sub(a, b, K), h, K)
|
||
|
h = dup_lshift(h, dup_degree(h), K)
|
||
|
|
||
|
return g
|
||
|
|
||
|
|
||
|
def dmp_revert(f, g, u, K):
|
||
|
"""
|
||
|
Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.polys import ring, QQ
|
||
|
>>> R, x,y = ring("x,y", QQ)
|
||
|
|
||
|
"""
|
||
|
if not u:
|
||
|
return dup_revert(f, g, K)
|
||
|
else:
|
||
|
raise MultivariatePolynomialError(f, g)
|