Slow arithmetic in polynomial rings?

[wdecker]

ping @fieker

The following computation came up in discussions with @pbelmans:

m = 13;
k = m-1;
n = m+1;
var = "c";
R, c = graded_polynomial_ring(QQ, var => 1:k; weights = 1:k);
@time inv_c = sum((-sum(c))^i for i in 1:n);
@time I = ideal(inv_c[n-k+1:n]);

The timings suggest that this can be done much faster in other systems.

Any thoughts?

[lgoettgens]

For anyone looking at this: The bottleneck is already sum(c)^14.

[fieker]

Can I have the MC2 code? Singular is same as Oscar, Magma, on nenepapa takes 12sec.
I can make the computation faster - by being slight more careful, e.g. if I want [c^i for i=1:10], so all powers, then this is one of the slowest ways to compute this. Also sum(c) is evaluated k times...

[ederc]

Isn't it also strange that when looking at the last line, accessing the vector entries takes most of the time?

julia> @time inv_c[n-k+1:n];
  4.005624 seconds (28.98 M allocations: 2.042 GiB, 10.34% gc time)

julia> @time I = ideal(inv_c[n-k+1:n]);
  4.234645 seconds (31.18 M allocations: 2.147 GiB, 9.85% gc time, 5.19% compilation time)

Accessing single elements from inv_c takes even longer for me:

julia> @time inv_c[1];
 40.816626 seconds (550.49 M allocations: 17.269 GiB, 18.42% gc time)

julia> @time inv_c[13];
 40.830624 seconds (550.49 M allocations: 17.269 GiB, 18.24% gc time)

[JohnAAbbott]

Isn't it also strange that when looking at the last line, accessing the vector entries takes most of the time?

julia> @time inv_c[n-k+1:n];
  4.005624 seconds (28.98 M allocations: 2.042 GiB, 10.34% gc time)

julia> @time I = ideal(inv_c[n-k+1:n]);
  4.234645 seconds (31.18 M allocations: 2.147 GiB, 9.85% gc time, 5.19% compilation time)

Accessing single elements from inv_c takes even longer for me:

julia> @time inv_c[1];
 40.816626 seconds (550.49 M allocations: 17.269 GiB, 18.42% gc time)

julia> @time inv_c[13];
 40.830624 seconds (550.49 M allocations: 17.269 GiB, 18.24% gc time)

@fieker explained to me that the apparent vector accesses are actually extracting homogeneous components; still this does not explain why a single "index" takes so much longer than a range...

[joschmitt]

Apparently, homogeneous_component (single index) just adds the terms of the polynomial whereas homogeneous_components (plural!, index range) uses the MPolyBuildCtx.

[fieker]

On Wed, Oct 22, 2025 at 05:18:11AM -0700, ederc wrote: ederc left a comment (oscar-system/Oscar.jl#5483) Isn't it also strange that when looking at the last line, accessing the vector entries takes most of the time?

Not strange at all... the code does - zip coeffs, expo vectors - expo vectors creates a Vector{Int} of the correct length - checks that the exponents fit into Int (they could be HUGE) - then calls an ! function to copy into the Vector{Int} - the Vector{Int} is then copied into a ZZMatrix (row) - the matrix is multiplied (inplace) - ...

…

``` julia> @time inv_c[n-k+1:n]; 4.005624 seconds (28.98 M allocations: 2.042 GiB, 10.34% gc time) julia> @time I = ideal(inv_c[n-k+1:n]); 4.234645 seconds (31.18 M allocations: 2.147 GiB, 9.85% gc time, 5.19% compilation time) ``` Accessing single elements from `inv_c` takes even longer for me: ``` julia> @time inv_c[1]; 40.816626 seconds (550.49 M allocations: 17.269 GiB, 18.42% gc time) julia> @time inv_c[13]; 40.830624 seconds (550.49 M allocations: 17.269 GiB, 18.24% gc time) ``` -- Reply to this email directly or view it on GitHub: #5483 (comment) You are receiving this because you were mentioned. Message ID: ***@***.***>

[JohnAAbbott]

It seems to be wasteful to compute the powers independently; a potentially more efficient approach is below.

m = 13;
k = m-1;
n = m+1;
var = "c";
R, c = graded_polynomial_ring(QQ, var => 1:k; weights = 1:k);
L = [zero(R)  for i in 1:n];
s = -sum(c);
L[1] = s;
@time for i in 2:n
  L[i] = s*L[i-1];
end
@time inv_c = sum(L);
@time I = ideal(inv_c[n-k+1:n]);

On my machine I get the following timings: 0.93s, 0.28s, 4.62s

[JohnAAbbott]

Apparently, homogeneous_component (single index) just adds the terms of the polynomial whereas homogeneous_components (plural!, index range) uses the MPolyBuildCtx.

If we have the functions, then why also the weird "vector syntax"? Of course, I don't have to use it if I don't want to, but if I have to read someone else's code which does use it then it's easy to be perplexed -- the function names are indeed very long, but they also make the code easier to comprehend. Simple hack/fix would be to make the single index version just call the range version on a singleton range...

[fieker]

The vector syntax above is in IntersectionTheory, the homogenous_comonents is in Rings, it returns all components. One could conceivable do the filtering immediately...

[joschmitt]

Regarding the line

@time I = ideal(inv_c[n-k+1:n]);

I get the following timings:

On master:

julia> @btime I = ideal(inv_c[n-k+1:n]);
  4.441 s (28980075 allocations: 2.04 GiB)

With PR #5487:

julia> @btime I = ideal(inv_c[n-k+1:n]);
  2.764 s (28974497 allocations: 1.73 GiB)

With PR #5487 and "faking" an in-place iterator over the exponents by using Nemo.exponent_vector!:

julia> @btime I = ideal(inv_c[n-k+1:n]);
  955.683 ms (4408 allocations: 379.76 KiB)

[fieker]

OK: inv_c has > 1e6 terms, just iterating over the coefficients create thus >1e6 rationals (QQFieldElem). This alone takes seconds. To improve this we need

• inplace iterators for both the exponents and the coefficients
• possibly also a combined one (to not use zip)
  To make life harder: the exponent iterator should (needs to) also deal with
• into Vector{Int}
• Vector{ZZRingElem}
• into a ZZMatrix (row, possibly via view)
  Exponents can be > 64 bits