From 1f40354d74a901cdee335f00b37163f4cc5b78ec Mon Sep 17 00:00:00 2001 From: Max Horn Date: Fri, 17 Jan 2025 16:55:10 +0100 Subject: [PATCH 1/3] Optimize map_coefficients, change_base_ring for some inputs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Namely if the input polynomial is a `QQMPolyRingElem` and the transformation "function" is `ZZ` resp. a `fpField`. Before this patch: julia> Qxy,(x,y) = QQ[:x,:y]; Zxy,_ = ZZ[:u,:v]; julia> f = 11*(x^2/2 + ZZ(3)^50*x*y^5/7 + y^6/12); julia> @b $f*84 160.000 ns (4 allocs: 200 bytes) julia> @b change_base_ring($ZZ, $f*84) 1.559 μs (45 allocs: 1.633 KiB) julia> @b change_base_ring($ZZ, $f*84; parent = $Zxy) 1.002 μs (24 allocs: 808 bytes) julia> @b map_coefficients($ZZ, $f*84; parent = $Zxy) 996.680 ns (24 allocs: 808 bytes) After this patch: julia> Qxy,(x,y) = QQ[:x,:y]; Zxy,_ = ZZ[:u,:v]; julia> f = 11*(x^2/2 + ZZ(3)^50*x*y^5/7 + y^6/12); julia> @b $f*84 111.693 ns (4 allocs: 200 bytes) julia> @b change_base_ring($ZZ, $f*84) 733.806 ns (30 allocs: 1.227 KiB) julia> @b change_base_ring($ZZ, $f*84; parent = $Zxy) 242.117 ns (9 allocs: 392 bytes) julia> @b map_coefficients($ZZ, $f*84; parent = $Zxy) 249.607 ns (9 allocs: 392 bytes) As a bonus we now also have a special `mul!` method which is even faster but is not suitable for the faint of heart: julia> @b mul!($Zxy(), $f, 84) 129.041 ns (5 allocs: 192 bytes) --- src/HeckeMoreStuff.jl | 15 --------- src/flint/fmpq_mpoly.jl | 63 +++++++++++++++++++++++++++++++++++ test/flint/fmpq_mpoly-test.jl | 17 ++++++++++ 3 files changed, 80 insertions(+), 15 deletions(-) diff --git a/src/HeckeMoreStuff.jl b/src/HeckeMoreStuff.jl index 262ab8148c..c073aac1c2 100644 --- a/src/HeckeMoreStuff.jl +++ b/src/HeckeMoreStuff.jl @@ -182,21 +182,6 @@ function Base.hash(f::zzModMPolyRingElem, h::UInt) return UInt(1) # TODO: enhance or throw error end -function AbstractAlgebra.map_coefficients(F::fpField, f::QQMPolyRingElem; parent=polynomial_ring(F, nvars(parent(f)), cached=false)[1]) - dF = denominator(f) - d = F(dF) - if iszero(d) - error("Denominator divisible by p!") - end - m = inv(d) - ctx = MPolyBuildCtx(parent) - for x in zip(coefficients(f), exponent_vectors(f)) - el = numerator(x[1] * dF) - push_term!(ctx, F(el) * m, x[2]) - end - return finish(ctx) -end - function tdivpow2!(B::ZZMatrix, t::Int) @ccall libflint.fmpz_mat_scalar_tdiv_q_2exp(B::Ref{ZZMatrix}, B::Ref{ZZMatrix}, t::Cint)::Nothing end diff --git a/src/flint/fmpq_mpoly.jl b/src/flint/fmpq_mpoly.jl index d3bd7d51bc..8db9a253e9 100644 --- a/src/flint/fmpq_mpoly.jl +++ b/src/flint/fmpq_mpoly.jl @@ -657,6 +657,11 @@ end # ############################################################################### +_content_ptr(c::QQMPolyRingElem) = Ptr{QQFieldElem}(pointer_from_objref(c)) +_content_ptr(c::Ptr{QQMPolyRingElem}) = Ptr{QQFieldElem}(c) +_content_ptr(c::Ref{QQMPolyRingElem}) = _content_ptr(c[]) +_zpoly_ptr(c::QQMPolyRingElemOrPtr) = Ptr{ZZMPolyRingElem}(_content_ptr(c) + sizeof(QQFieldElem)) + function zero!(a::QQMPolyRingElem) @ccall libflint.fmpq_mpoly_zero(a::Ref{QQMPolyRingElem}, a.parent::Ref{QQMPolyRing})::Nothing return a @@ -760,6 +765,22 @@ sub!(a::QQMPolyRingElem, b::RationalUnion, c::QQMPolyRingElem) = neg!(sub!(a, c, mul!(a::QQMPolyRingElem, b::QQMPolyRingElem, c::RationalUnion) = mul!(a, b, flintify(c)) mul!(a::QQMPolyRingElem, b::RationalUnion, c::QQMPolyRingElem) = mul!(a, c, b) +# special: multiply a QQMPolyRingElem by an integer and store the result as a ZZMPolyRingElem. +# obviously this only works if the integers is a multiple of the denominator of the polynomial +function mul!(a::ZZMPolyRingElem, b::QQMPolyRingElem, c::IntegerUnion) + @assert ngens(parent(a)) == ngens(parent(b)) + x = QQFieldElem() + GC.@preserve a b x begin + x = mul!(x, _content_ptr(b), c) + @assert isinteger(x) + bp = _zpoly_ptr(b) + xp = _num_ptr(x) + @ccall libflint.fmpz_mpoly_scalar_mul_fmpz(a::Ref{ZZMPolyRingElem}, bp::Ref{ZZMPolyRingElem}, xp::Ref{ZZRingElem}, parent(a)::Ref{ZZMPolyRing})::Nothing + end + return a +end + + divexact!(a::QQMPolyRingElem, b::QQMPolyRingElem, c::RationalUnion) = divexact!(a, b, flintify(c)) # Set the n-th coefficient of a to c. If zero coefficients are inserted, they @@ -1059,3 +1080,45 @@ function (R::QQMPolyRing)(a::Vector{Any}, b::Vector{Vector{T}}) where T return R(newaa, newbb) end + +############################################################################### +# +# Changing base ring, mapping coefficients +# +############################################################################### + +function map_coefficients(R::ZZRing, f::QQMPolyRingElem; + cached::Bool = true, + parent::ZZMPolyRing = AbstractAlgebra._change_mpoly_ring(R, parent(f), cached)) + @req isinteger(_content_ptr(f)) "input polynomial does not have integral coefficients" + return mul!(zero(parent), f, 1) +end + +function map_coefficients(F::fpField, f::QQMPolyRingElem; + cached::Bool = true, + parent::MPolyRing = AbstractAlgebra._change_mpoly_ring(F, parent(f), cached)) + dF = denominator(f) + d = F(dF) + if iszero(d) + error("Denominator divisible by p!") + end + m = inv(d) + ctx = MPolyBuildCtx(parent) + for x in zip(coefficients(f), exponent_vectors(f)) + el = numerator(x[1] * dF) + push_term!(ctx, F(el) * m, x[2]) + end + return finish(ctx) +end + +function change_base_ring(R::ZZRing, f::QQMPolyRingElem; + cached::Bool = true, + parent::ZZMPolyRing = AbstractAlgebra._change_mpoly_ring(R, parent(f), cached)) + return map_coefficients(R, f; cached, parent) +end + +function change_base_ring(F::fpField, f::QQMPolyRingElem; + cached::Bool = true, + parent::fpMPolyRing = AbstractAlgebra._change_mpoly_ring(F, parent(f), cached)) + return map_coefficients(R, f; cached, parent) +end diff --git a/test/flint/fmpq_mpoly-test.jl b/test/flint/fmpq_mpoly-test.jl index 37662000dc..470cea320f 100644 --- a/test/flint/fmpq_mpoly-test.jl +++ b/test/flint/fmpq_mpoly-test.jl @@ -813,3 +813,20 @@ end @test abc - a - b == c @test abc - ab == c end + +@testset "QQMPolyRingElem.convert_to_ZZMPolyRingElem" begin + Qxy, (x,y) = QQ[:x,:y] + f = 11*(x^2/2 + ZZ(3)^50*x*y^5/7 + y^6/12) + + Zxy, (x,y) = ZZ[:u,:v] + g = 462*x^2 + 94762534375324541717672868*x*y^5 + 77*y^6 + + @test g == @inferred mul!(zero(Zxy), f, 84) + @test g == change_base_ring(ZZ, f*84; parent = Zxy) + @test g == map_coefficients(ZZ, f*84; parent = Zxy) + + # test error handling + @test_throws ArgumentError change_base_ring(ZZ, f; parent = Zxy) + @test_throws ArgumentError map_coefficients(ZZ, f; parent = Zxy) + +end From 5a88ed0ca737f704028aa45259831c5fc5771283 Mon Sep 17 00:00:00 2001 From: Max Horn Date: Tue, 21 Jan 2025 11:17:14 +0100 Subject: [PATCH 2/3] input validation --- src/flint/fmpq_mpoly.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/flint/fmpq_mpoly.jl b/src/flint/fmpq_mpoly.jl index 8db9a253e9..16be895ebe 100644 --- a/src/flint/fmpq_mpoly.jl +++ b/src/flint/fmpq_mpoly.jl @@ -1090,7 +1090,9 @@ end function map_coefficients(R::ZZRing, f::QQMPolyRingElem; cached::Bool = true, parent::ZZMPolyRing = AbstractAlgebra._change_mpoly_ring(R, parent(f), cached)) - @req isinteger(_content_ptr(f)) "input polynomial does not have integral coefficients" + @req isinteger(_content_ptr(f)) "input polynomial must have integral coefficients" + @req ngens(parent) == ngens(Nemo.parent(f)) "parents must have matching numbers of generators" + @req internal_ordering(parent) == internal_ordering(Nemo.parent(f)) "parents must have matching internal ordering" return mul!(zero(parent), f, 1) end From 47be76c588945414ba82e22dd80e62771426906e Mon Sep 17 00:00:00 2001 From: Max Horn Date: Mon, 27 Jan 2025 11:45:41 +0100 Subject: [PATCH 3/3] remove @assert --- src/flint/fmpq_mpoly.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/flint/fmpq_mpoly.jl b/src/flint/fmpq_mpoly.jl index 16be895ebe..8d92ee6d63 100644 --- a/src/flint/fmpq_mpoly.jl +++ b/src/flint/fmpq_mpoly.jl @@ -768,11 +768,9 @@ mul!(a::QQMPolyRingElem, b::RationalUnion, c::QQMPolyRingElem) = mul!(a, c, b) # special: multiply a QQMPolyRingElem by an integer and store the result as a ZZMPolyRingElem. # obviously this only works if the integers is a multiple of the denominator of the polynomial function mul!(a::ZZMPolyRingElem, b::QQMPolyRingElem, c::IntegerUnion) - @assert ngens(parent(a)) == ngens(parent(b)) x = QQFieldElem() GC.@preserve a b x begin x = mul!(x, _content_ptr(b), c) - @assert isinteger(x) bp = _zpoly_ptr(b) xp = _num_ptr(x) @ccall libflint.fmpz_mpoly_scalar_mul_fmpz(a::Ref{ZZMPolyRingElem}, bp::Ref{ZZMPolyRingElem}, xp::Ref{ZZRingElem}, parent(a)::Ref{ZZMPolyRing})::Nothing