AbstractAlgebra.jl icon indicating copy to clipboard operation
AbstractAlgebra.jl copied to clipboard

Performance regression comparing universal polynomial rings after PR #2028

Open fingolfin opened this issue 7 months ago • 4 comments

@karlwessel made the following comment on PR #2028 after it was merged:

But while this improves performance for UnivPoly{Rational{BigInt}}, performance for UnivPoly{QQFieldElem} has decreased by a lot:

julia> using Nemo

Welcome to Nemo version 0.49.1

Nemo comes with absolutely no warranty whatsoever

julia> R = universal_polynomial_ring(QQ)
Universal Polynomial Ring over Rational field

julia> xs200 = gen.(Ref(R), [Symbol("x$i") for i in 1:200])
200-element Vector{AbstractAlgebra.Generic.UnivPoly{QQFieldElem}}:
 x1
 x2
 x3
 x4
 x5
 x6
 x7
 x8
 x9
 x10
 x11
 x12
 x13
 ⋮
 x188
 x189
 x190
 x191
 x192
 x193
 x194
 x195
 x196
 x197
 x198
 x199
 x200

julia> y = gen(R, :y)
y

julia> using BenchmarkTools

julia> @benchmark xs200[end] == y
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  342.801 μs … 727.543 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     343.088 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   353.182 μs ±  31.159 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █ ▄▂▃▂  ▅▄▁                                                   ▁
  █▅█████████▇▇██▆▆▆▅▄▅▅▆▅▄▁▄▅▄▅▅▃▄▃▃▄▁▄▁▃▃▃▁▁▁▃▁▃▄▁▁▃▄▁▁▃▁▁▄▃▅ █
  343 μs        Histogram: log(frequency) by time        544 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

For reference, the same MWE without this PR:

julia> using Nemo

Welcome to Nemo version 0.49.1

Nemo comes with absolutely no warranty whatsoever

julia> R = universal_polynomial_ring(QQ)
Universal Polynomial Ring over Rational field

julia> xs200 = gen.(Ref(R), [Symbol("x$i") for i in 1:200])
200-element Vector{AbstractAlgebra.Generic.UnivPoly{QQFieldElem}}:
 x1
 x2
 x3
 x4
 x5
 x6
 x7
 x8
 x9
 x10
 x11
 x12
 x13
 ⋮
 x188
 x189
 x190
 x191
 x192
 x193
 x194
 x195
 x196
 x197
 x198
 x199
 x200

julia> y = gen(R, :y)
y

julia> using BenchmarkTools

julia> @benchmark xs200[end] == y
BenchmarkTools.Trial: 10000 samples with 196 evaluations per sample.
 Range (min … max):  469.362 ns …  40.940 μs  ┊ GC (min … max):  0.00% … 96.81%
 Time  (median):     548.643 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   703.006 ns ± 918.173 ns  ┊ GC (mean ± σ):  14.93% ± 16.03%

   ▆█▆▃▂▃▁                                   ▁        ▁▁        ▁
  ▇████████▆▇█▇▇▅▆▃▃▁▄▁▁▃▁▁▁▁▁▃▁▃▁▁▁▃▁▁▁▁▁▄▆███▇▆▇▆▁▃███▇▇██▇▅▆ █
  469 ns        Histogram: log(frequency) by time       2.55 μs <

 Memory estimate: 4.84 KiB, allocs estimate: 5.

julia> 

Sorry for not realising this earlier :(

Originally posted by @karlwessel in https://github.com/Nemocas/AbstractAlgebra.jl/issues/2028#issuecomment-2717146939

fingolfin avatar May 01 '25 10:05 fingolfin

My guess is that we need to add specialized _is_gen_with_index methods to Nemo for Nemo types to avoid this?

While at it, should ie perhaps be made official under a name like is_gen_with_index? (Or alternatively, have a var_index argument resp. variant that that makes it now throw but instead return nothing it the argument is not a variable?)

See also issue #2063 which is about a missing generic is_gen(::MPolyRingElem).

fingolfin avatar May 01 '25 10:05 fingolfin

flint does not have this interface function. At the moment we do an is_gen(x, i) check for all i, but this is too slow when repeated 200 times. Not sure one can do something easy directly in Nemo. (The mpoly_is_gen is a tad more complicated due to the unpacking.)

thofma avatar May 01 '25 11:05 thofma

Why can't we an is_monomial check followed by a look at the exponent vector (in Nemo)? That should be way faster? (I think with some effort we could get that allocation free, too)

I could whip something up (no today though, out with family in this holiday)

fingolfin avatar May 01 '25 12:05 fingolfin

Yes, one can do this. It is the same logic as https://github.com/flintlib/flint/blob/726ee2d3f8e85cee82627a45a8a85d4bff482025/src/mpoly/is_gen.c. To call mpoly_get_monomial_ffmpz, one would have to pass fmpz *, but this is a bit cumbersome to do from within Nemo/julia.

thofma avatar May 01 '25 12:05 thofma