"""Advanced tools for dense recursive polynomials in ``K[x]`` or ``K[X]``. """
from __future__ import print_function, division
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.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.polyerrors import (
MultivariatePolynomialError,
DomainError
)
from sympy.utilities import variations
from math import ceil as _ceil, log as _log
from sympy.core.compatibility import range
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
[docs]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)
[docs]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 %s" % (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)
[docs]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)
[docs]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 dup_TC(f, K)
result = K.zero
for c in f:
result *= a
result += c
return result
[docs]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)
[docs]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)
[docs]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)
[docs]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)
[docs]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)
[docs]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)
[docs]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)
[docs]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)
[docs]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
[docs]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
[docs]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)
[docs]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)
[docs]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
[docs]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
[docs]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
[docs]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_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
[docs]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
[docs]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
[docs]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 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)
[docs]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
[docs]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
[docs]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)
