Overhaul preprocessing for radical computations.

[HechtiDerLachs]

As discussed last Wednesday with @wdecker, @jankoboehm, @hannes14, and @ederc . Currently this is WIP.

[HechtiDerLachs]

An example taken from #3711 for testing:

First without manual flattening, directly over towers of numberfields.

julia> P1, t1 = QQ[:t1];

julia> kk1, a = extension_field(t1^8 + 1);

julia> P2, t2 = kk1[:t2];

julia> kk2, b = extension_field(-a^6 + a^4 - a^2 + t2^2);

julia> R, (x, y, t) = kk2[:x, :y, :t];

julia> g = [x^3 - x*t^8 - x - y^2, x^2 - 1//3*t^8 - 1//3, y, x*t];

julia> J = ideal(R, g);

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=false, eliminate_variables=false);
  5.282182 seconds (2.37 M allocations: 63.570 MiB)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=true, eliminate_variables=false);
  4.977875 seconds (1.24 M allocations: 52.309 MiB)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=false, eliminate_variables=true);
  5.754451 seconds (2.49 M allocations: 65.264 MiB)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=true, eliminate_variables=true);
  6.275381 seconds (2.14 M allocations: 65.900 MiB, 11.58% gc time)

Now with manual flattening:

julia> R, (a, b, x, y, t) = QQ[:a, :b, :x, :y, :t];

julia> g = QQMPolyRingElem[x^3 - x*t^8 - x - y^2, x^2 - 1//3*t^8 - 1//3, y, x*t, a^8 + 1, -a^6 + a^4 - a^2 + b^2];

julia> J = ideal(R, g);

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=false, eliminate_variables=false);
 11.450386 seconds (1.52 M allocations: 98.770 MiB, 2.42% gc time)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=true, eliminate_variables=false);
  4.992918 seconds (653.11 k allocations: 27.977 MiB, 5.76% gc time)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=false, eliminate_variables=true);
 12.019289 seconds (4.06 M allocations: 206.570 MiB, 6.61% gc time)

julia> I = J^5;

julia> @time radical(I; use_squarefree_parts_of_generators=true, eliminate_variables=true);
  4.479714 seconds (689.18 k allocations: 29.212 MiB)

Test code for comparison with native Singular, started with Singular --ticks-per-sec 1000:

LIB "primdec.lib";
ring R = 0, (a, b, x, y, t), dp;
ideal I = x^3 - x*t^8 - x - y^2, x^2 - 1/3*t^8 - 1/3, y, x*t, a^8 + 1, -a^6 + a^4 - a^2 + b^2;
int t0 = timer;
radical(I^5);
timer - t0;

gives between 3600-3800 ms for me.

[HechtiDerLachs]

Test code for absolute primary decomposition:

julia> P1, t1 = QQ[:t1];

julia> kk1, a = extension_field(t1^8 + 1);

julia> P2, t2 = kk1[:t2];

julia> kk2, b = extension_field(-a^6 + a^4 - a^2 + t2^2);

julia> R, (x,) = polynomial_ring(kk2, [:x]);

julia> I = ideal(R, x^2 + 1)
Ideal generated by
  x^2 + 1

julia> @time dec = absolute_primary_decomposition(I)
 10.923452 seconds (21.09 M allocations: 8.204 GiB, 27.92% gc time)
2-element Vector{Any}:
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)

julia> I = ideal(R, x^2 + 1)
Ideal generated by
  x^2 + 1

julia> @time dec = absolute_primary_decomposition(I)
  5.495546 seconds (20.96 M allocations: 4.406 GiB, 31.38% gc time)
2-element Vector{Any}:
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)

julia> I = ideal(R, x^2 + 1)
Ideal generated by
  x^2 + 1

julia> @time dec = absolute_primary_decomposition(I)
  7.633494 seconds (20.93 M allocations: 5.316 GiB, 33.11% gc time)
2-element Vector{Any}:
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)
 (Ideal with 3 generators, Ideal with 3 generators, Ideal with 3 generators, 1)

In comparison with direct computations over QQ:

julia> R, (a, b, x) = QQ[:a, :b, :x];

julia> I = ideal(R, [x^2 + 1, a^8 + 1, -a^6 + a^4 - a^2 + b^2]);

julia> @time dec = absolute_primary_decomposition(I)
  9.097868 seconds (20.19 M allocations: 6.628 GiB, 35.20% gc time)
2-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}, Int64}}:
 (Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal with 3 generators, 16)
 (Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal with 3 generators, 16)

julia> I = ideal(R, [x^2 + 1, a^8 + 1, -a^6 + a^4 - a^2 + b^2]);

julia> @time dec = absolute_primary_decomposition(I)
  6.472232 seconds (20.28 M allocations: 5.211 GiB, 31.73% gc time)
2-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}, Int64}}:
 (Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal with 3 generators, 16)
 (Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal with 3 generators, 16)

julia> I = ideal(R, [x^2 + 1, a^8 + 1, -a^6 + a^4 - a^2 + b^2]);

julia> @time dec = absolute_primary_decomposition(I)
  7.240318 seconds (20.19 M allocations: 5.603 GiB, 30.65% gc time)
2-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}, Int64}}:
 (Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal (x^2 + 1, b^4 + 2*b^2*x + 1, 2*a^2 + b^2*x - b^2 - x - 1), Ideal with 3 generators, 16)
 (Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal (x^2 + 1, b^4 - 2*b^2*x + 1, 2*a^2 - b^2*x - b^2 + x - 1), Ideal with 3 generators, 16)

[HechtiDerLachs]

@hannes14 : When testing the above code, I found that actually a lot of processes were spawned on my system again. Even worse: They keep running, even after I quit julia and my laptop sounds like a vacuum cleaner. So the absolute primary decomposition might also be affected by this spawning problem.

[HechtiDerLachs]

What should we do about the failing booktests?

We should decide about the renaming of the kwarg factor_generators to use_squarefree....

Other than that, this should be good to go.

[lgoettgens]

What should we do about the failing booktests?

as far as I remember, you can adapt this one if it is just a reordering of the triples. but pinging @joschmitt to double-check

[joschmitt]

What should we do about the failing booktests?

as far as I remember, you can adapt this one if it is just a reordering of the triples. but pinging @joschmitt to double-check

Yes, there is no fixed order. Just update it.

[HechtiDerLachs]

@wdecker : Please let me know whether you're happy with the changes and the examples provided. The timings still suggest that there is room for improvement, but it seems to me that this can only be solved within Singular and taking care of the spawning issue.

[HechtiDerLachs]

@wdecker : Could you please make a suggestion for what the kwarg for the preprocessing via taking squarefree parts of the generators should be called? Thank you.

[wdecker]

@wdecker : Could you please make a suggestion for what the kwarg for the preprocessing via taking squarefree parts of the generators should be called? Thank you.

To keep things short: How about squarefree::Bool=true

[jankoboehm]

@wdecker : Could you please make a suggestion for what the kwarg for the preprocessing via taking squarefree parts of the generators should be called? Thank you.

To keep things short: How about squarefree::Bool=true

I think that is a bit too short, since one could understand that this actually has an impact on the result (which is even more confusing since it is supposed to be the radical). Can we introduce an algorithm keyword which has as standard behaviour the orginal behaviour? (Keeping factor_generators for compatibility)

[wdecker]

I think that is a bit too short, since one could understand that this actually has an impact on the result (which is even more confusing since it is supposed to be the radical). Can we introduce an algorithm keyword which has as standard behaviour the orginal behaviour? (Keeping factor_generators for compatibility)

How about preprocessing::Bool=true and explaining in the doctest what preprocessing means.
Also, keep factor_generators::Bool=true for backwards compatibility, but do not showing this in the docstring.

[fingolfin]

In triage just now @jankoboehm stated that he agrees with @wdecker on preprocessing so @HechtiDerLachs will now implement that.

[codecov]

Codecov Report

Attention: Patch coverage is 93.18182% with 3 lines in your changes missing coverage. Please review.

Project coverage is 84.53%. Comparing base (0d59f06) to head (b4233f8).
Report is 10 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Rings/primary_decomposition_helpers.jl	80.00%	2 Missing ⚠️
experimental/GaloisGrp/src/Solve.jl	50.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #4488      +/-   ##
==========================================
- Coverage   84.55%   84.53%   -0.02%     
==========================================
  Files         672      671       -1     
  Lines       88833    88909      +76     
==========================================
+ Hits        75109    75161      +52     
- Misses      13724    13748      +24     

Files with missing lines	Coverage Δ
...ometry/Schemes/AffineSchemes/Objects/Attributes.jl	91.66% <ø> (ø)
...try/Schemes/CoveredSchemes/Morphisms/Attributes.jl	85.71% <ø> (ø)
...metry/Schemes/CoveredSchemes/Objects/Attributes.jl	98.18% <ø> (ø)
src/AlgebraicGeometry/Schemes/Sheaves/Types.jl	93.58% <ø> (ø)
src/Rings/MPolyQuo.jl	93.47% <100.00%> (ø)
src/Rings/mpoly-graded.jl	90.26% <100.00%> (+0.01%)	⬆️
src/Rings/mpoly-ideals.jl	94.19% <100.00%> (-0.19%)	⬇️
experimental/GaloisGrp/src/Solve.jl	73.71% <50.00%> (ø)
src/Rings/primary_decomposition_helpers.jl	85.60% <80.00%> (-0.91%)	⬇️

... and 16 files with indirect coverage changes

[HechtiDerLachs]

Tests are passing. I think this is fine from my side.