Skip to content

Commit

Permalink
Optimize map_coefficients, change_base_ring for some inputs
Browse files Browse the repository at this point in the history
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)
  • Loading branch information
fingolfin committed Jan 21, 2025
1 parent fb7f84d commit 550e2e6
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 18 deletions.
15 changes: 0 additions & 15 deletions src/HeckeMoreStuff.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
69 changes: 66 additions & 3 deletions src/flint/fmpq_mpoly.jl
Original file line number Diff line number Diff line change
Expand Up @@ -115,11 +115,16 @@ function content(a::QQMPolyRingElem)
end

function denominator(a::QQMPolyRingElem)
c = ZZRingElem()
@ccall libflint.fmpq_mpoly_get_denominator(c::Ref{ZZRingElem}, a::Ref{QQMPolyRingElem}, parent(a)::Ref{QQMPolyRing})::Nothing
return c
z = ZZRingElem()
return denominator!(z, a)
end

function denominator!(z::ZZRingElemOrPtr, a::QQMPolyRingElem)
@ccall libflint.fmpq_mpoly_get_denominator(z::Ref{ZZRingElem}, a::Ref{QQMPolyRingElem}, parent(a)::Ref{QQMPolyRing})::Nothing
return z
end


function fit!(a::QQMPolyRingElem, n::Int)
# needs to exist for the MPoly interface
return nothing
Expand Down Expand Up @@ -657,6 +662,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
Expand Down Expand Up @@ -760,6 +770,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
Expand Down Expand Up @@ -792,6 +818,43 @@ function combine_like_terms!(a::QQMPolyRingElem)
return a
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"
fp = _zpoly_ptr(f)
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

###############################################################################
#
# Manipulating terms and monomials
Expand Down
17 changes: 17 additions & 0 deletions test/flint/fmpq_mpoly-test.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 550e2e6

Please sign in to comment.