From cbb00e1fb58c5906a690ea311769f133c36a60b7 Mon Sep 17 00:00:00 2001 From: wbhart Date: Fri, 19 Mar 2021 09:58:44 +0100 Subject: [PATCH] Move divhigh, mulhigh_n and divexact_low from Misc. Add tests. (#801) --- src/generic/Misc/Poly.jl | 62 ---------------------------- src/generic/Poly.jl | 85 ++++++++++++++++++++++++++++++++++++++- src/julia/Integer.jl | 3 ++ test/generic/Poly-test.jl | 51 +++++++++++++++++++++++ 4 files changed, 138 insertions(+), 63 deletions(-) diff --git a/src/generic/Misc/Poly.jl b/src/generic/Misc/Poly.jl index 10dd234ef1..a07a6ff797 100644 --- a/src/generic/Misc/Poly.jl +++ b/src/generic/Misc/Poly.jl @@ -31,68 +31,6 @@ function valence(f::PolyElem) return c end -############################################################################## -# -# Binary arithmetic -# -############################################################################## - -# computes the top n coeffs of the product only -function mulhigh_n(a::PolyElem{T}, b::PolyElem{T}, n::Int) where {T} - # sum a_i t^i and sum b_j t^j - # want (i, j) s.th. i + j >= deg a + deg b - n - r = parent(a)() - for i=max(degree(a)-n, 0):degree(a) - for j = max(degree(a) + degree(b) - n - i, 0):degree(b) - setcoeff!(r, i+j, coeff(r, i+j) + coeff(a, i)*coeff(b, j)) - end - end - return r -end - -# assuming b divides a, compute the last n coeffs of the quotient -# will produce garbage otherwise -# div(a, b) mod x^n -function divexact_low(a::PolyElem{T}, b::PolyElem{T}, n::Int) where {T} - r = parent(a)() - a = truncate(a, n) - b = truncate(b, n) - for i = 0:n - 1 - q = divexact(trail(a), trail(b)) - setcoeff!(r, i, q) - a = shift_right(a - q*b, 1) - b = truncate(b, n - i - 1) - # truncate both a and b to n - i - 1 (for generic polys one could just change the length) - end - return r -end - -#computes the top coeffs starting with x^n -function divhigh(a::PolyElem{T}, b::PolyElem{T}, n0::Int) where {T} - r = parent(a)() - n = degree(a) - degree(b) - n0 - fit!(r, degree(a) - degree(b)) - a = deepcopy(a) - k = degree(a) - n0 - da = degree(a) - for i=0:n - if degree(a) < degree(b) - break - end - q = divexact(coeff(a, da), lead(b)) - setcoeff!(r, da - degree(b), q) - for j=da:-1:max(k, da - degree(b)) - setcoeff!(a, j, coeff(a, j)-q*coeff(b, j-da+degree(b))) - end - da -= 1 - end - if iszero(r) - return r - end - set_length!(r, normalise(r, length(r) - 1)) - return r -end - ################################################################################ # # Deflate and inflate diff --git a/src/generic/Poly.jl b/src/generic/Poly.jl index 65a2a11cfa..379e9ffa97 100644 --- a/src/generic/Poly.jl +++ b/src/generic/Poly.jl @@ -151,7 +151,7 @@ lead(a::PolynomialElem) = length(a) == 0 ? base_ring(a)(0) : coeff(a, length(a) Return the trailing coefficient of the given polynomial. This will be the nonzero coefficient of the term with lowest degree unless the polynomial -in the zero polynomial, in which case a zero coefficient is returned. +is the zero polynomial, in which case a zero coefficient is returned. """ function trail(a::PolynomialElem) if iszero(a) @@ -958,6 +958,89 @@ function mullow(a::AbstractAlgebra.PolyElem{T}, b::AbstractAlgebra.PolyElem{T}, return z end +# computes the terms of the product from degree deg(a) + deg(b) down to +# deg(a) + deg(b) - n inclusive, setting the remaining terms to zero +function mulhigh_n(a::PolyElem{T}, b::PolyElem{T}, n::Int) where T <: RingElement + # if a = sum a_i t^i and b = sum b_j t^j + # want (i, j) such that i + j >= deg a + deg b - n + r = parent(a)() + for i = max(degree(a) - n, 0):degree(a) + for j = max(degree(a) + degree(b) - n - i, 0):degree(b) + r = setcoeff!(r, i + j, coeff(r, i + j) + coeff(a, i)*coeff(b, j)) + end + end + return r +end + +# assuming b divides a (behaviour is undefined otherwise), computes the last n +# terms of the quotient, i.e. computes divexact(a, b) mod x^n +function divexact_low(a::PolyElem{T}, b::PolyElem{T}, n::Int) where T <: RingElement + r = parent(a)() + if iszero(b) + return r + end + shift = 0 + for i = 0:degree(b) + if !iszero(coeff(b, i)) + break + end + shift += 1 + end + if shift != 0 + a = shift_right(a, shift) + b = shift_right(b, shift) + end + a = truncate(a, n) + b = truncate(b, n) + for i = 0:n - 1 + flag, q = divides(coeff(a, 0), coeff(b, 0)) + !flag && error("Not an exact division") + r = setcoeff!(r, i, q) + a = shift_right(a - q*b, 1) + b = truncate(b, n - i - 1) + # truncate both a and b to n - i - 1 + end + return r +end + +# computes the top terms of the quotient of a by b starting with the term of +# degree n0 +# if the division is not exact until at least the term of degree n0, an +# exception may be raised +function divhigh(a::PolyElem{T}, b::PolyElem{T}, n0::Int) where T <: RingElement + r = parent(a)() + n = degree(a) - degree(b) - n0 + fit!(r, degree(a) - degree(b) + 1) + a = deepcopy(a) + da = degree(a) + R = base_ring(a) + t = R() + for i = 0:n + if da < degree(b) + break + end + # negate quotient so we can use addeq! below + q = -divexact(coeff(a, da), lead(b)) + r = setcoeff!(r, da - degree(b), q) + da -= 1 + if i != n + c = coeff(a, da) + for j = 0:min(degree(b) - 1, i) + t = mul!(t, coeff(r, da - degree(b) + j + 1), + coeff(b, degree(b) - j - 1)) + c = addeq!(c, t) + end + a = setcoeff!(a, da, c) + end + end + if iszero(r) + return r + end + r = set_length!(r, normalise(r, length(r))) + # negate r to compensate for negation above + return -r +end + ############################################################################### # # Reversal diff --git a/src/julia/Integer.jl b/src/julia/Integer.jl index 573f1067e7..51dd69e4e5 100644 --- a/src/julia/Integer.jl +++ b/src/julia/Integer.jl @@ -104,6 +104,9 @@ end ############################################################################### function divides(a::T, b::T) where T <: Integer + if b == 0 + return a == 0, b + end q, r = divrem(a, b) return r == 0, q end diff --git a/test/generic/Poly-test.jl b/test/generic/Poly-test.jl index 4dd46eba23..38f04f9cdd 100644 --- a/test/generic/Poly-test.jl +++ b/test/generic/Poly-test.jl @@ -609,6 +609,57 @@ end @test truncate(f*g, n) == mullow(f, g, n) end + for iter = 1:100 + lena = rand(0:10) + lenb = rand(0:10) + a = rand(R, lena:lena, -10:10) + b = rand(R, lenb:lenb, -10:10) + c = a*b + for i = 0:length(c) - 1 + d = mulhigh_n(a, b, i) + f = c - d + @test length(f) < length(c) - i + for j = length(c):length(c) - i + @test coeff(d, j) == coeff(c, j) + end + end + end + + for iter = 1:100 + lena = rand(0:10) + lenb = rand(1:10) + a = rand(R, lena:lena, -10:10) + b = R() + while iszero(b) + b = rand(R, lenb:lenb, -10:10) + end + c = a*b + for i = 0:length(b) + d = divhigh(c, b, i) + f = a - d + @test degree(f) < i + for j = 0:i - 1 + @test coeff(f, j) == coeff(a, j) + end + end + end + + for iter = 1:100 + lena = rand(0:10) + lenb = rand(1:10) + a = rand(R, lena:lena, -10:10) + b = R() + while iszero(b) + b = rand(R, lenb:lenb, -10:10) + end + c = a*b + for i = 0:length(a) + d = divexact_low(c, b, i) + f = truncate(a, i) + @test f == d + end + end + # Fake finite field of char 7, degree 2 S, y = PolynomialRing(GF(7), "y") F = ResidueField(S, y^2 + 6y + 3)