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 4 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
17 changes: 8 additions & 9 deletions src/MPoly.jl
View file
Open in desktop
Copy link
Collaborator

@JohnAAbbott JohnAAbbott Jan 6, 2025

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

@fingolfin fingolfin Jan 6, 2025

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

@fingolfin fingolfin Jan 6, 2025

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

@JohnAAbbott JohnAAbbott Jan 7, 2025

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 @@ -426,19 +426,18 @@ end
iszero(x::MPolyRingElem{T}) where T <: RingElement = length(x) == 0

function is_unit(f::T) where {T <: MPolyRingElem}
# for constant polynomials we delegate to the coefficient ring:
is_constant(f) && return is_unit(constant_coefficient(f))
# Here deg(f) > 0; over an integral domain, non-constant polynomials are never units:
is_domain_type(T) && return false
# 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.
!is_unit(constant_coefficient(f)) && return false
# for constant polynomials we are done now
is_constant(f) && return true
# over a domain, non-constant polynomials are never units
is_domain_type(T) && return false
# remains to check the non-constant non-zero terms have nilpotent coefficients
for (c, expv) in zip(coefficients(f), exponent_vectors(f))
if is_zero(expv)
is_unit(c) || return false
else
is_nilpotent(c) || return false
end
is_zero(expv) && continue # skip constant term
is_nilpotent(c) || return false
end
return true
end
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
Loading