Skip to content

Add conformance tests for MPolyRing, and fix bugs in ^ and is_unit for MPoly #1950

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Jan 10, 2025
Merged

Add conformance tests for MPolyRing, and fix bugs in ^ and is_unit for MPoly #1950

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d51ff09
Add conformance tests for MPolyRing
fingolfin Dec 25, 2024
dc96201
Fix bug in ^(f::Generic.MPoly, n::Int)
fingolfin Dec 26, 2024
b90fb4e
Fix is_unit for MPolyRingElem
fingolfin Dec 26, 2024
516a952
fixup! Add conformance tests for MPolyRing
fingolfin Dec 26, 2024
0ed1bf9
Ensure residue ring is quotient of euclidean ring
fingolfin Jan 6, 2025
04625b5
fixup! Fix is_unit for MPolyRingElem
fingolfin Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion src/MPoly.jl
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I avoided calling constant_coefficient because it might iterate over the terms in f, but then we iterate again over the terms in f in the main loop. Here is a partial demonstration:

julia> P,(x,y) = polynomial_ring(QQ, ["x","y"]);
julia> f = (2*x+3*y-4)^3 + (3*x-2*y-6)^3;
julia> f999 = f^999;
julia> @time constant_coefficient(f999);
  3.127883 seconds (17.98 M allocations: 5.133 GiB, 15.79% gc time)

I have just tried a more realistic test, which actually calls the code I wrote:

julia> P,(x,y) = polynomial_ring(ZZmod720, ["x","y"]);
julia> f = (2*x+3*y-4)^3 + (3*x-2*y-6)^3;
julia> f9999 = f^9999; # takes a long time
julia> @time constant_coefficient(f9999);
  0.052131 seconds (1.00 M allocations: 76.614 MiB, 18.73% gc time)
julia> @time is_unit(f9999)
  0.000009 seconds (1 allocation: 80 bytes)
false

So in this case, calling constant_coefficient will be a bit slower. That said, @fingolfin code is neater.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

More importantly, my new code is less wrong ;-).

julia> R,_ = residue_ring(ZZ,ZZ(720));

julia> S,(x,y)=R[:x,:y];

julia> is_unit(30*x)
true

I am happy if somebody adjusts the code to be faster again, but to me correctness trumps speed.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK I thought about it and I think I found a better way, we'll see.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The version with constant_term_is_unit looks fine to me: i.e. should be both fast and correct

Original file line number Diff line number Diff line change
Expand Up @@ -433,14 +433,16 @@ function is_unit(f::T) where {T <: MPolyRingElem}
# A polynomial over a commutative ring is a unit iff its
# constant term is a unit and all other coefficients are nilpotent:
# see e.g. <https://kconrad.math.uconn.edu/blurbs/ringtheory/polynomial-properties.pdf> for a proof.
constant_term_is_unit = false
for (c, expv) in zip(coefficients(f), exponent_vectors(f))
if is_zero(expv)
is_unit(c) || return false
constant_term_is_unit = true
else
is_nilpotent(c) || return false
end
end
return true
return constant_term_is_unit
end

function content(a::MPolyRingElem{T}) where T <: RingElement
Expand Down
4 changes: 3 additions & 1 deletion src/generic/MPoly.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2316,6 +2316,8 @@ function ^(a::MPoly{T}, b::Int) where {T <: RingElement}
return zero(a)
end
elseif length(a) == 1
c = coeff(a, 1)^b
is_zero(c) && return zero(a)
N = size(a.exps, 1)
exps = zeros(UInt, N, 1)
monomial_mul!(exps, 1, a.exps, 1, b, N)
Expand All @@ -2324,7 +2326,7 @@ function ^(a::MPoly{T}, b::Int) where {T <: RingElement}
error("Exponent overflow in powering")
end
end
return parent(a)([coeff(a, 1)^b], exps)
return parent(a)([c], exps)
elseif b == 0
return one(a)
elseif b == 1
Expand Down
142 changes: 139 additions & 3 deletions test/Rings-conformance-tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -719,15 +719,15 @@ function test_Poly_interface(Rx::AbstractAlgebra.PolyRing; reps = 30)
@test symbols(Rx) isa Vector{Symbol}
@test length(symbols(Rx)) == 1
@test is_gen(gen(Rx))
@test is_gen(x)
@test is_monic(x)
@test is_trivial(Rx) || !is_gen(x^2)
for i in 1:reps
a = test_elem(Rx)
@test iszero(a) || degree(a) >= 0
@test equality(a, leading_coefficient(a)*x^max(0, degree(a)) + tail(a))
@test constant_coefficient(a) isa elem_type(R)
@test trailing_coefficient(a) isa elem_type(R)
@test is_gen(x)
@test iszero(one(Rx)) || !is_gen(x^2)
@test is_monic(a) == isone(leading_coefficient(a))
end
end
Expand All @@ -737,6 +737,140 @@ function test_Poly_interface(Rx::AbstractAlgebra.PolyRing; reps = 30)
end


function test_MPoly_interface(Rxy::AbstractAlgebra.MPolyRing; reps = 30)

# for simplicity, these tests for now assume exactly two generators
@assert ngens(Rxy) == 2

T = elem_type(Rxy)

@testset "MPoly interface for $(Rxy) of type $(typeof(Rxy))" begin

test_Ring_interface(Rxy; reps = reps)

@testset "Basic functionality" begin
@test symbols(Rxy) isa Vector{Symbol}
@test length(symbols(Rxy)) == ngens(Rxy)
@test length(gens(Rxy)) == ngens(Rxy)
@test gens(Rxy) == [gen(Rxy, i) for i in 1:ngens(Rxy)]
@test all(is_gen, gens(Rxy)) || is_trivial(Rxy)
end

@testset "Polynomial Constructors" begin
for i in 1:reps
a = test_elem(Rxy)::T
for b in coefficients(a)
@assert Rxy(b) isa T
end

# test MPolyBuildCtx
B = MPolyBuildCtx(Rxy)
for (c, e) in zip(AbstractAlgebra.coefficients(a), AbstractAlgebra.exponent_vectors(a))
push_term!(B, c, e)
end
@test finish(B) == a
end
x, y = gens(Rxy)
f = 13*x^3*y^4 + 2*x - 7
#@test Rxy([2,-7,13], [[1,0],[0,0],[3,4]]) == f # FIXME: interface spec does not say this is required?

R = base_ring(Rxy)
@test Rxy(R.([2,-7,13]), [[1,0],[0,0],[3,4]]) == f
end

# skip trivial rings after this, it is not worth the bother
is_trivial(Rxy) && return

@testset "Element properties" begin
R = base_ring(Rxy)
x, y = gens(Rxy)

a = zero(Rxy)
@test !is_monomial(a)
@test !is_term(a)
@test is_constant(a)
@test !is_gen(a)
@test !is_unit(a)
@test is_nilpotent(a)
@test length(a) == 0
@test total_degree(a) < 0
@test all(is_negative, degrees(a))

a = one(Rxy)
@test is_monomial(a)
@test is_term(a)
@test is_constant(a)
@test !is_gen(a)
@test is_unit(a)
@test !is_nilpotent(a)
@test length(a) == 1
@test total_degree(a) == 0
@test degrees(a) == [0, 0]

a = x
@test is_monomial(a)
@test is_term(a)
@test !is_constant(a)
@test is_gen(a)
@test !is_unit(a)
@test !is_nilpotent(a)
@test length(a) == 1
@test total_degree(a) == 1
@test degrees(a) == [1, 0]

a = x^2
@test is_monomial(a)
@test is_term(a)
@test !is_constant(a)
@test !is_gen(a)
@test !is_unit(a)
@test !is_nilpotent(a)
@test length(a) == 1
@test total_degree(a) == 2
@test degrees(a) == [2, 0]

if !is_zero(R(2))
a = 2*x
@test !is_monomial(a)
@test is_term(a)
@test !is_constant(a)
@test !is_gen(a)
@test !is_unit(a)
@test is_nilpotent(a) == is_nilpotent(R(2))
@test length(a) == 1
@test total_degree(a) == 1
@test degrees(a) == [1, 0]
end

a = x^3 + y^4
@test !is_monomial(a)
@test !is_term(a)
@test !is_constant(a)
@test !is_gen(a)
@test !is_unit(a)
@test !is_nilpotent(a)
@test length(a) == 2
@test total_degree(a) == 4
@test degrees(a) == [3, 4]

for i in 1:reps
a = test_elem(Rxy)
iszero(a) && continue
@test length(a) >= 0
@test sum(degrees(a)) >= total_degree(a)
end

end

# TODO: add more tests, covering everything described in the manual, see
# https://nemocas.github.io/AbstractAlgebra.jl/dev/mpoly_interface/
# https://nemocas.github.io/AbstractAlgebra.jl/dev/mpolynomial/
end

return nothing
end


function test_MatSpace_interface(S::MatSpace; reps = 20)

ST = elem_type(S)
Expand Down Expand Up @@ -891,8 +1025,10 @@ end

function test_Ring_interface_recursive(R::AbstractAlgebra.Ring; reps = 50)
test_Ring_interface(R; reps = reps)
Rx, _ = polynomial_ring(R, "x")
Rx, _ = polynomial_ring(R, :x)
test_Poly_interface(Rx, reps = 2 + fld(reps, 2))
Rxy, _ = polynomial_ring(R, [:x, :y])
test_MPoly_interface(Rxy, reps = 2 + fld(reps, 2))
S = matrix_ring(R, rand(0:3))
test_MatAlgebra_interface(S, reps = 2 + fld(reps, 2))
S = matrix_space(R, rand(0:3), rand(0:3))
Expand Down
2 changes: 1 addition & 1 deletion test/generic/Residue-test.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
test_Ring_interface_recursive(T)

#
S, x = polynomial_ring(ZZ, "x")
S, x = polynomial_ring(QQ, "x")
T, = residue_ring(S, x^2 + 1)
test_Ring_interface_recursive(T)

Expand Down
Loading