Add puiseux_expansion

[HechtiDerLachs]

This is a first take at wrapping Puiseux expansion from PuiseuxExpansion.lib as requested by @simonbrandhorst .

@jankoboehm : Could you please have a look? I tried to make sense of the output from Singular, but I could only make an educated guess at what we hold in our hands there. There is lots of data being discarded and I'm wondering whether 1. my interpretation of the output is correct so far and 2. whether we want to keep more of the data and make it available in Oscar.

One particular question: Singular's puiseux command takes an integer value max_ord. I am not sure how that relates to the prec argument in the constructor puiseux_series_ring. So far, I just set them to be equal, i.e. I call Singular's puiseux with max_ord=k and wrap the output in a puiseux_series_ring with prec=k for the same k. But I'm afraid that this might not be the correct approach. @jankoboehm : If you remember the implementation in Singular, do you see how max_ord and prec are related?

Pinging also @YueRen for critical input on this topic.

[HechtiDerLachs]

I have a problem here. If I run my code on this polynomial

julia> g = x^8 + x^2 + y^4
x^8 + x^2 + y^4

julia> hg = Oscar.puiseux_expansion(g, 20)
ERROR: BoundsError: attempt to access 2-element Vector{Any} at index [3]
Stacktrace:
 [1] throw_boundserror(A::Vector{Any}, I::Tuple{Int64})
   @ Base ./essentials.jl:14
 [2] getindex
   @ ./essentials.jl:916 [inlined]
 [3] indexed_iterate(a::Vector{Any}, i::Int64, state::Int64)
   @ Base ./tuple.jl:160
 [4] _puiseux(f::QQMPolyRingElem, max_deg::Int64, at_origin::Bool)
   @ Oscar ~/Oscar.jl/src/Rings/puiseux_wrapper.jl:18
 [5] puiseux_expansion(f::QQMPolyRingElem, precision::Int64; parent::AbstractAlgebra.Generic.PuiseuxSeriesField{QQFieldElem})
   @ Oscar ~/Oscar.jl/src/Rings/puiseux_wrapper.jl:62
 [6] puiseux_expansion(f::QQMPolyRingElem, precision::Int64)
   @ Oscar ~/Oscar.jl/src/Rings/puiseux_wrapper.jl:52
 [7] top-level scope
   @ REPL[184]:1

I get an error. The reason is that all of a sudden Singular's puiseux command decides to return a completely different format of output. That seems to depend on the polynomial which is put in. Any ideas why this happens? Maybe @hannes14 ?

Edit: To be more precise about what the issue is, I include the following code snippet.

julia> R, (x, y) = QQ[:x, :y]
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> g = y^4 + x^2 + x^8 # exponent 4 for y
x^8 + x^2 + y^4

julia> SR = Oscar.singular_poly_ring(R)
Singular polynomial ring (QQ),(x,y),(dp(2),C)

julia> gg = SR(g)
x^8 + y^4 + x^2

julia> Singular.LibPuiseuxexpansions.puiseux(gg, 20, 1)
1-element Vector{Vector{Any}}:
 [Singular polynomial ring (0,@a),(x,y),(dp(2),C), Dict{Symbol, Vector}(:erg => Singular.n_algExt{Singular.Rationals}[@a], :PE => Vector{Any}[[(-77*@a)//2048*x^49 + (7*@a)//128*x^37 - (3*@a)//32*x^25 + @a//4*x^13 + @a*x, 2, nothing, nothing, nothing, [1, 1, 1, 1, 1], [1]]], :minPolys => Any[1, y^4 + 1])]

julia> g = y^3 + x^2 + x^8 # exponent 3 for y
x^8 + x^2 + y^3

julia> gg = SR(g)
x^8 + y^3 + x^2

julia> Singular.LibPuiseuxexpansions.puiseux(gg, 20, 1)
1-element Vector{Vector{Any}}:
 [10//243*x^74 - 5//81*x^56 + 1//9*x^38 - 1//3*x^20 - x^2, 3, nothing, nothing, nothing, [1, 1, 1, 1, 1], [2]]

It is unclear for me how Singular decides to return something completely different in the two cases for y^4 and y^3. It looks like it tries to produce some extension field for the result to live in/over. That might be a fair thing to do, but if this is not foreseeable, then how do I make sure to collect the info in the output correctly?

[jankoboehm]

I will have a look!

[jankoboehm]

Here are a few notes about how the Singular library works:

The output of the command puiseux is mainly designed for direct use in the integral basis algorithm, which is why only basic data types are used. One can easily wrap this into a newstruct type to make it type stable, but for Oscar this would currently be more of a disadvantage since we cannot easily transfer newstructs. So we are lucky.

puiseux returns a list in which each entry corresponds to either a single Puiseux expansion or a conjugate group of expansions. In the first case (i.e. a rational expansion) the data is given directly; in the second case a ring is returned, which then contains the Puiseux expansion in PE, which depends on the variable of the minimal polynomial minpoly of the ring. Via type we can distinguish the cases.

The expansions themselves always consist of a list. Only the first two entries should be relevant:
the first is a polynomial, the second is the lcm of the denominators of the exponents. (The remaining entries are for the integral basis algorithm, where one needs information on how to construct polynomials over the base field from Puiseux factors.)

One more remark: the option atOrigin controls whether all expansions at $x=0$ are computed, or only those that pass through $(0,0)$.

@hannes14 is looking into how to properly pack and unpack the information for Oscar.

[HechtiDerLachs]

Thank you, @jankoboehm ! I will meet with @hannes14 to discuss further steps.

[HechtiDerLachs]

I made another attempt at wrapping Singular's output correctly. The price is that the user can not specify parent objects for the output anymore, as they are dynamically generated during the process. Unfortunately one test still fails which is probably due to the inconsistency of precision vs. max_order described above.

[HechtiDerLachs]

I found that I accidentally hardcoded an exponent. That is fixed now and the tests should pass. However, it would be good if e.g. @jankoboehm could comment on the validity of the choice of max_order vs. precision.

[HechtiDerLachs]

I have now made max_ord and precision different arguments of the Puiseux expansion command. It is now up to the user to choose these values with the default of them being equal.

For the time being this should bring the PR to a mergable state as the code is not producing wrong results, unless the user themself can be held accountable for it due to faulty input.

In the long term it would be desirable to have the default set a bit more conciously. But that will eventually come up through people actually using it, knowing what they're doing, and getting back to us with "Hey, you know: In my application, it is useful to set it to ..."

[HechtiDerLachs]

Failing tests do not seem to be my fault. Good to go from my side.

[lgoettgens]

is it a feature that parents are not compatible? If not, I would suggest to use a @test_broken instead to mark that this is the current state of things, but eventually it would be good to fix it

[HechtiDerLachs]

It is not a feature, but it also can not be avoided. Thus I would say it is unfortunate, but expected behavior.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 84.07%. Comparing base (3e1cd0e) to head (39c47df).
⚠️ Report is 16 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #5588      +/-   ##
==========================================
+ Coverage   84.05%   84.07%   +0.02%     
==========================================
  Files         744      745       +1     
  Lines      101320   101403      +83     
==========================================
+ Hits        85168    85259      +91     
+ Misses      16152    16144       -8     

Files with missing lines	Coverage Δ
src/Rings/puiseux_wrapper.jl	100.00% <100.00%> (ø)

... and 13 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jankoboehm]

For maxDeg > 0, the series are computed to all exponents < maxDeg, and if the next exponent would be exactly maxDeg, that term is included and then the process stops. If maxDeg <=0 then the singular part is computed. So positive maxDeg works as the precision, you might get a bit more if this is anyway known, but that does not harm. One could make use of the extra though.

[HechtiDerLachs]

So positive maxDeg works as the precision

Alright, this is the confirmation we needed. Thanks!

[simonbrandhorst]

Many thanks!
Could you please add it to the documentation and export it?
I think that we should also hook it up to the (plane) curves. But if you prefer, then we can do this another time.

[HechtiDerLachs]

I tried to find the correct spot to put it in the documentation; it seems to be here. But those files seem to be imported from AA and Nemo in some ways which I can not really follow.

Does anyone have a good suggestion where to put the docstring for puiseux_expansion?

[simonbrandhorst]

How about AlgebraicGeometry/Curves?