# Source code for diofant.polys.densearith

"""Arithmetics for dense recursive polynomials in K[x] or K[X]. """

from .densebasic import (dmp_degree_in, dmp_LC, dmp_one, dmp_one_p, dmp_slice,
dmp_strip, dmp_zero, dmp_zero_p, dmp_zeros)
from .polyerrors import ExactQuotientFailed, PolynomialDivisionFailed

"""
Add c*x**i to f in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_add_term(x**2 - 1, ZZ(2), 4)
2*x**4 + x**2 - 1
"""
if not c:
return f

n = len(f)
m = n - i - 1

if i == n - 1:
return dmp_strip([f[0] + c] + f[1:], 0)
else:
if i >= n:
return [c] + [K.zero]*(i - n) + f
else:
return f[:m] + [f[m] + c] + f[m + 1:]

[docs]def dmp_add_term(f, c, i, u, K):
"""
Add c(x_2..x_u)*x_0**i to f in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_add_term(x*y + 1, 2, 2)
2*x**2 + x*y + 1
"""
if not u:

v = u - 1

if dmp_zero_p(c, v):
return f

n = len(f)
m = n - i - 1

if i == n - 1:
return dmp_strip([dmp_add(f[0], c, v, K)] + f[1:], u)
else:
if i >= n:
return [c] + dmp_zeros(i - n, v, K) + f
else:
return f[:m] + [dmp_add(f[m], c, v, K)] + f[m + 1:]

[docs]def dmp_sub_term(f, c, i, u, K):
"""
Subtract c(x_2..x_u)*x_0**i from f in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_sub_term(2*x**2 + x*y + 1, 2, 2)
x*y + 1
"""
if not u:

v = u - 1

if dmp_zero_p(c, v):
return f

n = len(f)
m = n - i - 1

if i == n - 1:
return dmp_strip([dmp_sub(f[0], c, v, K)] + f[1:], u)
else:
if i >= n:
return [dmp_neg(c, v, K)] + dmp_zeros(i - n, v, K) + f
else:
return f[:m] + [dmp_sub(f[m], c, v, K)] + f[m + 1:]

[docs]def dup_mul_term(f, c, i, K):
"""
Multiply f by c*x**i in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_mul_term(x**2 - 1, ZZ(3), 2)
3*x**4 - 3*x**2
"""
if not c or not f:
return []
else:
return [ cf * c for cf in f ] + [K.zero]*i

[docs]def dmp_mul_term(f, c, i, u, K):
"""
Multiply f by c(x_2..x_u)*x_0**i in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_mul_term(x**2*y + x, 3*y, 2)
3*x**4*y**2 + 3*x**3*y
"""
if not u:
return dup_mul_term(f, c, i, K)

v = u - 1

if dmp_zero_p(f, u):
return f
if dmp_zero_p(c, v):
return dmp_zero(u)
else:
return [ dmp_mul(cf, c, v, K) for cf in f ] + dmp_zeros(i, v, K)

[docs]def dmp_mul_ground(f, c, u, K):
"""
Multiply f by a constant value in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_mul_ground(2*x + 2*y, ZZ(3))
6*x + 6*y
"""
if not u:
return dmp_strip([coeff * c for coeff in f], u)
else:
v = u - 1
return [dmp_mul_ground(coeff, c, v, K) for coeff in f]

[docs]def dmp_quo_ground(f, c, u, K):
"""
Quotient by a constant in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)
>>> R.dmp_quo_ground(2*x**2*y + 3*x, ZZ(2))
x**2*y + x

>>> R, x, y = ring("x y", QQ)
>>> R.dmp_quo_ground(2*x**2*y + 3*x, QQ(2))
x**2*y + 3/2*x
"""
if not u:
if not c:
raise ZeroDivisionError('polynomial division')
if not f:
return f

if K.is_Field:
return [K.quo(coeff, c) for coeff in f]
else:
return [coeff // c for coeff in f]

v = u - 1
return [dmp_quo_ground(coeff, c, v, K) for coeff in f]

[docs]def dmp_exquo_ground(f, c, u, K):
"""
Exact quotient by a constant in K[X].

Examples
========

>>> R, x, y = ring("x y", QQ)

>>> R.dmp_exquo_ground(x**2*y + 2*x, QQ(2))
1/2*x**2*y + x
"""
if not u:
if not c:
raise ZeroDivisionError('polynomial division')
if not f:
return f

return [K.exquo(coeff, c) for coeff in f]

v = u - 1
return [dmp_exquo_ground(coeff, c, v, K) for coeff in f]

[docs]def dup_lshift(f, n, K):
"""
Efficiently multiply f by x**n in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dup_lshift(x**2 + 1, 2)
x**4 + x**2
"""
if not f:
return f
else:
return f + [K.zero]*n

[docs]def dup_rshift(f, n, K):
"""
Efficiently divide f by x**n in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dup_rshift(x**4 + x**2, 2)
x**2 + 1
>>> R.dup_rshift(x**4 + x**2 + 2, 2)
x**2 + 1
"""
return f[:-n]

[docs]def dmp_abs(f, u, K):
"""
Make all coefficients positive in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_abs(x**2*y - x)
x**2*y + x
"""
if not u:
return [abs(coeff) for coeff in f]
else:
v = u - 1
return [dmp_abs(coeff, v, K) for coeff in f]

[docs]def dmp_neg(f, u, K):
"""
Negate a polynomial in K[X].

Examples
=======

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_neg(x**2*y - x)
-x**2*y + x
"""
if not u:
return [-coeff for coeff in f]
else:
v = u - 1
return [dmp_neg(coeff, v, K) for coeff in f]

"""
Add dense polynomials in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_add(x**2 - 1, x - 2)
x**2 + x - 3
"""
if not f:
return g
if not g:
return f

df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if df == dg:
return dmp_strip([a + b for a, b in zip(f, g)], 0)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = g[:k], g[k:]

return h + [ a + b for a, b in zip(f, g) ]

"""
Add dense polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_add(x**2 + y, x**2*y + x)
x**2*y + x**2 + x + y
"""
if not u:

df = dmp_degree_in(f, 0, u)

if df < 0:
return g

dg = dmp_degree_in(g, 0, u)

if dg < 0:
return f

v = u - 1

if df == dg:
return dmp_strip([ dmp_add(a, b, v, K) for a, b in zip(f, g) ], u)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = g[:k], g[k:]

return h + [ dmp_add(a, b, v, K) for a, b in zip(f, g) ]

[docs]def dup_sub(f, g, K):
"""
Subtract dense polynomials in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_sub(x**2 - 1, x - 2)
x**2 - x + 1
"""
if not f:
return dmp_neg(g, 0, K)
if not g:
return f

df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if df == dg:
return dmp_strip([a - b for a, b in zip(f, g)], 0)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = dmp_neg(g[:k], 0, K), g[k:]

return h + [ a - b for a, b in zip(f, g) ]

[docs]def dmp_sub(f, g, u, K):
"""
Subtract dense polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_sub(x**2 + y, x**2*y + x)
-x**2*y + x**2 - x + y
"""
if not u:
return dup_sub(f, g, K)

df = dmp_degree_in(f, 0, u)

if df < 0:
return dmp_neg(g, u, K)

dg = dmp_degree_in(g, 0, u)

if dg < 0:
return f

v = u - 1

if df == dg:
return dmp_strip([ dmp_sub(a, b, v, K) for a, b in zip(f, g) ], u)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = dmp_neg(g[:k], u, K), g[k:]

return h + [ dmp_sub(a, b, v, K) for a, b in zip(f, g) ]

[docs]def dmp_add_mul(f, g, h, u, K):
"""
Returns f + g*h where f, g, h are in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_add_mul(x**2 + y, x, x + 2)
2*x**2 + 2*x + y

"""
return dmp_add(f, dmp_mul(g, h, u, K), u, K)

[docs]def dmp_sub_mul(f, g, h, u, K):
"""
Returns f - g*h where f, g, h are in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_sub_mul(x**2 + y, x, x + 2)
-2*x + y
"""
return dmp_sub(f, dmp_mul(g, h, u, K), u, K)

[docs]def dup_mul(f, g, K):
"""
Multiply dense polynomials in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_mul(x - 2, x + 2)
x**2 - 4
"""
if f == g:
return dup_sqr(f, K)

if not (f and g):
return []

df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

n = max(df, dg) + 1

if n < 100:
h = []

for i in range(df + dg + 1):
coeff = K.zero

for j in range(max(0, i - dg), min(df, i) + 1):
coeff += f[j]*g[i - j]

h.append(coeff)

return dmp_strip(h, 0)
else:
# Use Karatsuba's algorithm (divide and conquer), see e.g.:
# Joris van der Hoeven, Relax But Don't Be Too Lazy,
# J. Symbolic Computation, 11 (2002), section 3.1.1.
n2 = n//2

fl, gl = dmp_slice(f, 0, n2, 0, K), dmp_slice(g, 0, n2, 0, K)

fh = dup_rshift(dmp_slice(f, n2, n, 0, K), n2, K)
gh = dup_rshift(dmp_slice(g, n2, n, 0, K), n2, K)

lo, hi = dup_mul(fl, gl, K), dup_mul(fh, gh, K)

mid = dup_sub(mid, dup_add(lo, hi, K), K)

dup_lshift(hi, 2*n2, K), K)

[docs]def dmp_mul(f, g, u, K):
"""
Multiply dense polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_mul(x*y + 1, x)
x**2*y + x
"""
if not u:
return dup_mul(f, g, K)

if f == g:
return dmp_sqr(f, u, K)

df = dmp_degree_in(f, 0, u)

if df < 0:
return f

dg = dmp_degree_in(g, 0, u)

if dg < 0:
return g

h, v = [], u - 1

for i in range(df + dg + 1):
coeff = dmp_zero(v)

for j in range(max(0, i - dg), min(df, i) + 1):
coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)

h.append(coeff)

return dmp_strip(h, u)

[docs]def dup_sqr(f, K):
"""
Square dense polynomials in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_sqr(x**2 + 1)
x**4 + 2*x**2 + 1
"""
df, h = len(f) - 1, []

for i in range(2*df + 1):
c = K.zero

jmin = max(0, i - df)
jmax = min(i, df)

n = jmax - jmin + 1

jmax = jmin + n // 2 - 1

for j in range(jmin, jmax + 1):
c += f[j]*f[i - j]

c += c

if n & 1:
elem = f[jmax + 1]
c += elem**2

h.append(c)

return dmp_strip(h, 0)

[docs]def dmp_sqr(f, u, K):
"""
Square dense polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_sqr(x**2 + x*y + y**2)
x**4 + 2*x**3*y + 3*x**2*y**2 + 2*x*y**3 + y**4
"""
if not u:
return dup_sqr(f, K)

df = dmp_degree_in(f, 0, u)

if df < 0:
return f

h, v = [], u - 1

for i in range(2*df + 1):
c = dmp_zero(v)

jmin = max(0, i - df)
jmax = min(i, df)

n = jmax - jmin + 1

jmax = jmin + n // 2 - 1

for j in range(jmin, jmax + 1):
c = dmp_add(c, dmp_mul(f[j], f[i - j], v, K), v, K)

c = dmp_mul_ground(c, K(2), v, K)

if n & 1:
elem = dmp_sqr(f[jmax + 1], v, K)
c = dmp_add(c, elem, v, K)

h.append(c)

return dmp_strip(h, u)

[docs]def dmp_pow(f, n, u, K):
"""
Raise f to the n-th power in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_pow(x*y + 1, 3)
x**3*y**3 + 3*x**2*y**2 + 3*x*y + 1
"""
if not n:
return dmp_one(u, K)
if n < 0:
raise ValueError("can't raise polynomial to a negative power")
if n == 1 or dmp_zero_p(f, u) or dmp_one_p(f, u, K):
return f

g = dmp_one(u, K)

while True:
n, m = n//2, n

if m & 1:
g = dmp_mul(g, f, u, K)

if not n:
break

f = dmp_sqr(f, u, K)

return g

[docs]def dup_pexquo(f, g, K):
"""
Polynomial pseudo-quotient in K[x].

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dup_pexquo(x**2 - 1, 2*x - 2)
2*x + 2

>>> R.dup_pexquo(x**2 + 1, 2*x - 4)
Traceback (most recent call last):
...
ExactQuotientFailed: [2, -4] does not divide [1, 0, 1]
"""
q, r = dmp_pdiv(f, g, 0, K)

if not r:
return q
else:
raise ExactQuotientFailed(f, g)

[docs]def dmp_pdiv(f, g, u, K):
"""
Polynomial pseudo-division in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_pdiv(x**2 + x*y, 2*x + 2)
(2*x + 2*y - 2, -4*y + 4)
"""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if dg < 0:
raise ZeroDivisionError("polynomial division")

q, r, dr = dmp_zero(u), f, df

if df < dg:
return q, r

N = df - dg + 1
lc_g = dmp_LC(g, K)

while True:
lc_r = dmp_LC(r, K)
j, N = dr - dg, N - 1

Q = dmp_mul_term(q, lc_g, 0, u, K)
q = dmp_add_term(Q, lc_r, j, u, K)

R = dmp_mul_term(r, lc_g, 0, u, K)
G = dmp_mul_term(g, lc_r, j, u, K)
r = dmp_sub(R, G, u, K)

_dr, dr = dr, dmp_degree_in(r, 0, u)

if dr < dg:
break
assert dr < _dr

if u:
c = dmp_pow(lc_g, N, u - 1, K)
else:
c = lc_g**N

q = dmp_mul_term(q, c, 0, u, K)
r = dmp_mul_term(r, c, 0, u, K)

return q, r

[docs]def dmp_prem(f, g, u, K):
"""
Polynomial pseudo-remainder in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_prem(x**2 + x*y, 2*x + 2)
-4*y + 4
"""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if dg < 0:
raise ZeroDivisionError("polynomial division")

r, dr = f, df

if df < dg:
return r

N = df - dg + 1
lc_g = dmp_LC(g, K)

while True:
lc_r = dmp_LC(r, K)
j, N = dr - dg, N - 1

R = dmp_mul_term(r, lc_g, 0, u, K)
G = dmp_mul_term(g, lc_r, j, u, K)
r = dmp_sub(R, G, u, K)

_dr, dr = dr, dmp_degree_in(r, 0, u)

if dr < dg:
break
assert dr < _dr

if u:
c = dmp_pow(lc_g, N, u - 1, K)
else:
c = lc_g**N

return dmp_mul_term(r, c, 0, u, K)

[docs]def dmp_pquo(f, g, u, K):
"""
Polynomial exact pseudo-quotient in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + x*y
>>> g = 2*x + 2*y
>>> h = 2*x + 2

>>> R.dmp_pquo(f, g)
2*x

>>> R.dmp_pquo(f, h)
2*x + 2*y - 2
"""
return dmp_pdiv(f, g, u, K)[0]

[docs]def dmp_pexquo(f, g, u, K):
"""
Polynomial pseudo-quotient in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + x*y
>>> g = 2*x + 2*y
>>> h = 2*x + 2

>>> R.dmp_pexquo(f, g)
2*x

>>> R.dmp_pexquo(f, h)
Traceback (most recent call last):
...
ExactQuotientFailed: [[2], [2]] does not divide [[1], [1, 0], []]
"""
q, r = dmp_pdiv(f, g, u, K)

if dmp_zero_p(r, u):
return q
else:
raise ExactQuotientFailed(f, g)

[docs]def dmp_rr_div(f, g, u, K):
"""
Multivariate division with remainder over a ring.

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_rr_div(x**2 + x*y, 2*x + 2)
(0, x**2 + x*y)
"""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if dg < 0:
raise ZeroDivisionError("polynomial division")

q, r, dr = dmp_zero(u), f, df

if df < dg:
return q, r

lc_g, v = dmp_LC(g, K), u - 1

while True:
lc_r = dmp_LC(r, K)

if v >= 0:
c, R = dmp_rr_div(lc_r, lc_g, v, K)
if not dmp_zero_p(R, v):
break
else:
if lc_r % lc_g:
break
c = K.exquo(lc_r, lc_g)

j = dr - dg

q = dmp_add_term(q, c, j, u, K)
h = dmp_mul_term(g, c, j, u, K)
r = dmp_sub(r, h, u, K)

_dr, dr = dr, dmp_degree_in(r, 0, u)

if dr < dg:
break
assert dr < _dr

return q, r

[docs]def dmp_ff_div(f, g, u, K):
"""
Polynomial division with remainder over a field.

Examples
========

>>> R, x, y = ring("x y", QQ)

>>> R.dmp_ff_div(x**2 + x*y, 2*x + 2)
(1/2*x + 1/2*y - 1/2, -y + 1)
"""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if dg < 0:
raise ZeroDivisionError("polynomial division")

q, r, dr = dmp_zero(u), f, df

if df < dg:
return q, r

lc_g, v = dmp_LC(g, K), u - 1

while True:
lc_r = dmp_LC(r, K)

if v >= 0:
c, R = dmp_ff_div(lc_r, lc_g, v, K)
if not dmp_zero_p(R, v):
break
else:
c = K.exquo(lc_r, lc_g)

j = dr - dg

q = dmp_add_term(q, c, j, u, K)
h = dmp_mul_term(g, c, j, u, K)
r = dmp_sub(r, h, u, K)

_dr, dr = dr, dmp_degree_in(r, 0, u)

if dr < dg:
break
elif not (dr < _dr):
raise PolynomialDivisionFailed(f, g, K)

return q, r

[docs]def dmp_div(f, g, u, K):
"""
Polynomial division with remainder in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)
>>> R.dmp_div(x**2 + x*y, 2*x + 2)
(0, x**2 + x*y)

>>> R, x, y = ring("x y", QQ)
>>> R.dmp_div(x**2 + x*y, 2*x + 2)
(1/2*x + 1/2*y - 1/2, -y + 1)
"""
if K.is_Field:
return dmp_ff_div(f, g, u, K)
else:
return dmp_rr_div(f, g, u, K)

[docs]def dmp_rem(f, g, u, K):
"""
Returns polynomial remainder in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)
>>> R.dmp_rem(x**2 + x*y, 2*x + 2)
x**2 + x*y

>>> R, x, y = ring("x y", QQ)
>>> R.dmp_rem(x**2 + x*y, 2*x + 2)
-y + 1
"""
return dmp_div(f, g, u, K)[1]

[docs]def dmp_quo(f, g, u, K):
"""
Returns exact polynomial quotient in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)
>>> R.dmp_quo(x**2 + x*y, 2*x + 2)
0

>>> R, x, y = ring("x y", QQ)
>>> R.dmp_quo(x**2 + x*y, 2*x + 2)
1/2*x + 1/2*y - 1/2
"""
return dmp_div(f, g, u, K)[0]

[docs]def dmp_exquo(f, g, u, K):
"""
Returns polynomial quotient in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + x*y
>>> g = x + y
>>> h = 2*x + 2

>>> R.dmp_exquo(f, g)
x

>>> R.dmp_exquo(f, h)
Traceback (most recent call last):
...
ExactQuotientFailed: [[2], [2]] does not divide [[1], [1, 0], []]
"""
q, r = dmp_div(f, g, u, K)

if dmp_zero_p(r, u):
return q
else:
raise ExactQuotientFailed(f, g)

[docs]def dmp_max_norm(f, u, K):
"""
Returns maximum norm of a polynomial in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_max_norm(2*x*y - x - 3)
3
"""
if not u:
return max(dmp_abs(f, 0, K), default=K.zero)

v = u - 1
return max(dmp_max_norm(c, v, K) for c in f)

[docs]def dmp_l1_norm(f, u, K):
"""
Returns l1 norm of a polynomial in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_l1_norm(2*x*y - x - 3)
6
"""
if not u:
return sum(dmp_abs(f, u, K), K.zero)

v = u - 1
return sum(dmp_l1_norm(c, v, K) for c in f)

[docs]def dmp_expand(polys, u, K):
"""
Multiply together several polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_expand([x**2 + y**2, x + 1])
x**3 + x**2 + x*y**2 + y**2
"""
if not polys:
return dmp_one(u, K)

f = polys[0]

for g in polys[1:]:
f = dmp_mul(f, g, u, K)

return f