Series <-> Poly

[wbhart]

It's often convenient to be able to quickly convert from a dense univariate polynomial p(x) to the absolute series expression p(x) + O(x^(deg(p) + 1)) and back again, for use internally (e.g. polynomial_to_power_sum and power_sum_to_polynomial uses this).

I propose we add functions series(::AbstractAlgebra.PolyElem{T}) and poly(::AbstractAlgebra.AbsSeriesElem{T}) to convert between the two. The variable symbol would be kept the same between the two, though there could also be an optional argument for setting that.

Similarly it would be useful to have constructors series(A::Array{T}, n::Int) and poly(A::Array{T}) for constructing absolute series and polynomials from a given array of coefficients.

Implementation of these would probably depend on having something like dense_poly_type and dense_series_type.

See https://github.com/Nemocas/AbstractAlgebra.jl/issues/837

[fieker]

On Fri, Apr 09, 2021 at 01:55:13AM -0700, wbhart wrote:

It's often convenient to be able to quickly convert from a dense univariate polynomial p(x) to the absolute series expression p(x) + O(x^(deg(p) + 1)) and back again, for use internally (e.g. polynomial_to_power_sum and power_sum_to_polynomial uses this).

I propose we add functions series(::AbstractAlgebra.PolyElem{T}) and poly(::AbstractAlgebra.AbsSeriesElem{T}) to convert between the two. The variable symbol would be kept the same between the two, though there could also be an optional argument for setting that. Slight counter proposal: lift(PolyRing, Series) or lift(Series) similar to our (ab)use of lift for padics I'd also propose lift(FracField(Poly), Series) for RelSeries rational_reconstruction(Series) map_coeffs

poly already "exists":

julia> methods(polynomial)

5 methods for generic function "polynomial":

[1] polynomial(A::Vector{T}) where T<:RingElem in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1111 [2] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:RingElem in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1118 [3] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:Integer in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1122 [4] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:Rational in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1126 [5] polynomial(f::Hecke.PseudoPoly) in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/PseudoPolynomial.jl:20

Problem are RelSeries for both directions.

Similarly it would be useful to have constructors series(A::Array{T}, n::Int) and poly(A::Array{T}) for constructing polynomials from a given array of coefficients.

Why not S(poly) giving the SeriesRing? This would also bypass the problem with Rel/Abs series

Do make this useful, I think we need to stop caching the series rings based on precision

Implementation of these would probably depend on having something like dense_poly_type and dense_series_type. Possibly...

-- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836

[wbhart]

On Fri, Apr 09, 2021 at 01:55:13AM -0700, wbhart wrote: It's often convenient to be able to quickly convert from a dense univariate polynomial p(x) to the absolute series expression p(x) + O(x^(deg(p) + 1)) and back again, for use internally (e.g. polynomial_to_power_sum and power_sum_to_polynomial uses this). I propose we add functions series(::AbstractAlgebra.PolyElem{T}) and poly(::AbstractAlgebra.AbsSeriesElem{T}) to convert between the two. The variable symbol would be kept the same between the two, though there could also be an optional argument for setting that. Slight counter proposal: lift(PolyRing, Series) or lift(Series) similar to our (ab)use of lift for padics

I'm ok with implementing that as well, but this requires one to have constructed the polynomial ring already, and what I am proposing avoids doing that.

I'd also propose lift(FracField(Poly), Series) for RelSeries rational_reconstruction(Series) map_coeffs poly already "exists": julia> methods(polynomial) # 5 methods for generic function "polynomial": [1] polynomial(A::Vector{T}) where T<:RingElem in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1111 [2] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:RingElem in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1118 [3] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:Integer in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1122 [4] polynomial(R::AbstractAlgebra.Ring, A::Vector{T}) where T<:Rational in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/Poly.jl:1126 [5] polynomial(f::Hecke.PseudoPoly) in Hecke at /home/fieker/.julia/dev/Hecke/src/Misc/PseudoPolynomial.jl:20 Problem are RelSeries for both directions.

Why RelSeries. I was going to suggest we not bother with RelSeries. It can't in any sense be thought of as a residue ring so "lift" doesn't exactly make sense here. I was just going to implement conversions in each direction to AbsSeries.

If there's a genuine need for RelSeries instead, perhaps we ought to have functions relative_series(::PolyElem{T}) and absolute_series(::PolyElem{T})

Similarly it would be useful to have constructors series(A::Array{T}, n::Int) and poly(A::Array{T}) for constructing polynomials from a given array of coefficients. Why not S(poly) giving the SeriesRing?

We already have this, I think, but it requires you to have constructed the series ring in advance, and I just don't think that is as useful in algorithms. It's messy to say the least, and there's always the temptation to leave out the AbstractAlgebra prefix to the parent constructor, not to mention the fact that stuff will slow down over time as people put more and more precomputation in the parent objects.

This would also bypass the problem with Rel/Abs series

Could you say something about this. I don't understand what problem you are referring to. Do you just mean the fact that you have to specify whether you want a relative or absolute series?

Do make this useful, I think we need to stop caching the series rings based on precision

Fair enough. But perhaps that should be a separate ticket?

Implementation of these would probably depend on having something like dense_poly_type and dense_series_type. Possibly...

[wbhart]

Why not S(poly) giving the SeriesRing?

Or do you mean a function called S that returns the corresponding series ring of a polynomial ring. I don't like the name, since S is used throughout our code as a type parameter and it is not exactly a very extensible pattern.

This is obviously divergent from how we are handling matrices. Then again, matrices can be constructed without parents, so maybe you have a point. I'll have to think about this.

[thofma]

I gather from the above comments that the idea is to create the polynomial object without populating the parent?

[wbhart]

Perhaps without fully populating it. But I'm not sure how feasible that is. But at least copying as much as possible.

The main aim is to stop the user having to explicitly create the parent itself, which is an error prone and messy business. In my opinion, the parent constructors were always meant to be called by the user, not by library code.

[wbhart]

E.g. I think all Nemo fmpz_mod_blah types use the same fmpz_mod_ctx_t so that need not be generated de nuovo, but could be copied. I'm sure there would be other savings we could realise. (Though performance is really a secondary consideration here.)

[wbhart]

I'm also happy to call the poly function polynomial. However the other day I thought this would run us into problems eventually when we start doing sparse univariates, mpolys and others. Keeping the name the same as the abstract type seems most consistent. However, I can't for the life of me remember what the problem was any more.

[thofma]

Sorry, I got confused. The polynomial(x::AbsSeries) thing is mainly for users or library code? For users to be useful it would have to cache the parent (otherwise one cannot add to polynomials constructed in this way). But this would make it useless for library code, because library functions should not populate the cache.

[wbhart]

@thofma For library code.

@fieker Rather than opening a ticket for not caching on the series precision let me say here that this wouldn't work. In the capped model, the precision is part of the definition of the ring. If a user creates a series ring with precision cap 20, they would not be best pleased if the cache looked it up and returned for them a series ring with precision cap 100! I recall now the discussion we had, and what I think you were advocating for was an absolute series ring that operated precisely like R[x]/(x^n). This would be a different model from the AbsSeries model we currently have. I'd be happy to implement this, but it would have to raise an exception for all operations which cannot return a meaningful result to exactly precision O(x^n). If you really want this, let me know and I'll open a ticket for it. I still don't see how you get around caching on n though. I know you want this for lifting algorithms, and I am just not sure how to do that. What seems more workable is just an absolute series model like the one we have, but with no precision cap (therefore none to cache), but raising an exception when any operation leads to an infinite number of terms. That would also be workable, so let me know if you want this and I'll open a ticket for that.

[fieker]

On Fri, Apr 09, 2021 at 02:28:25AM -0700, wbhart wrote:

Why not S(poly) giving the SeriesRing?

I mean if S is a series ring, then S(poly) can be made to work Or do you mean a function called S that returns the corresponding series ring of a polynomial ring. I don't like the name, since S is used throughout our code as a type parameter and it is not exactly a very extensible pattern.

This is obviously divergent from how we are handling matrices. Then again, matrices can be constructed without parents, so maybe you have a point. I'll have to think about this.

-- You are receiving this because you commented. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-816552453

[thofma]

So the polynomial(...) function you propose will not create something that is cached. I am happy with that, but maybe it could have another name? unsafe_polynomial? It is not really safe to use within the Nemo/AA world as a user.

[fieker]

On Fri, Apr 09, 2021 at 03:11:27AM -0700, thofma wrote:

So the polynomial(...) function you propose will not create something that is cached. I am happy with that, but maybe it could have another name? unsafe_polynomial? It is not really safe to use within the Nemo/AA world as a user.

What can you do with a parent free poly? -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-816577847

[fieker]

On Fri, Apr 09, 2021 at 03:03:24AM -0700, wbhart wrote:

@thofma For library code.

@fieker Rather than opening a ticket for not caching on the series precision let me say here that this wouldn't work. In the capped model, the precision is part of the definition of the ring. If a user creates a series ring with precision cap 20, they would not be best pleased if the cache looked it up and returned for them a series ring with precision cap 100! I recall now the discussion we had, and what I think you were advocating for was an absolute series ring that operated precisely like R[x]/(x^n). This would be a different model from the AbsSeries model we currently have. I'd be happy to implement this, but it would have to raise and exception for all operations which cannot return a meaningful result to precision O(x^n). If you really want this, let me know and I'll open a ticket for it. I still don't see how you get around caching on n though. I know you want this for lifting algorithms, and I am just not sure how to do that. What seems more workable is just an absolute series model like the one we have, but with no precision cap (therefore none to cache), but raising an exception when any operation leads to an infinite number of terms. That would also be workable, so let me know if you want this and I'll open a ticket for that.

It really depends on your application. Try to implement an efficient Hensel/ Newton lifting. You really do not want to create new rings for all lifting steps. If the ring is a completion, you also do not know the final precision in advance, so basically you need to be able to change precisions - at least on elements. As long as arithmetic between elements that are not exact is ignoring the ring precision, this works. If you can change the ring precision, then caching based on the precision is stupid....

Suggestion, seriously, implement a toy Newton lift and see that you need. Maybe I am not seeing how to use the interface, but lifting without changing the precision is a waste of time

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-816573189

[wbhart]

@fieker Yes, S(poly) can be made to work as well. Not what I'm proposing, but fine.

Both: forget the parentless thing. That wasn't the intention. The intention is just that the developer doesn't have to create the parent in library code. It's just easier to write:

series(p, "x") than

R = AbsSeriesRing(base_ring(p), degree(p) + 1, 0)[1]
s = R(p)

It's even worse at the moment because the second line there doesn't even exist.

You'd currently have to write:

A = [coeff(p, i) for i in 0:degree(p)]
R = AbsSeriesRing(base_ring(p), degree(p) + 1, 0)[1]
s = R(A)

@fieker Yes, I see what you want re lifting. Of course that works just fine at present so long as you create the ring with the maximum precision you want to use as cap in advance. You change the precision of the series as we go, not of the parent! We already do this in our code in lots of places. Naturally it would be nice to make this even easier, but I don't see a way around actually raising the precision of the series as you go. The only thing we could make simpler is to not have a precision cap at all. That's exactly what I proposed as one of my suggestions above.

[thofma]

OK. We just have to make sure that the library code is not caching the parents. Library code should not fill the cache automatically.

[wbhart]

@thfoma I agree with that principle. Another reason why having the developer create parent objects in library code is error prone. They frequently forget to turn off caching.

In short I'd like to make it as easy to make polynomials and series (including from one another) as it is to currently create matrices.

[fieker]

On Fri, Apr 09, 2021 at 04:14:11AM -0700, wbhart wrote:

@fieker Yes, S(poly) can be made to work as well. Not what I'm proposing, but fine.

Both: forget the parentless thing. That wasn't the intention. The intention is just that the developer doesn't have to create the parent in library code. It's just easier to write:

series(p, "x") than

R = AbsSeriesRing(base_ring(p), degree(p) + 1, 0)[1]
s = R(p)

It's even worse at the moment because the second line there doesn't even exist.

@fieker Yes, I see what you want re lifting. Of course that works just fine at present so long as you create the ring with the maximum precision you want to use as cap in advance. You change the precision of the series as we go, not of the parent! We already do this in our code in lots of places. Naturally it would be nice to make this even easier, but I don't see a way around actually raising the precision of the series as you go. The only thing we could make simpler is to not have a precision cap at all. That's exactly what I proposed as one of my suggestions above.

I cannot set the ring precision in advance, the ring is a completion.... (And I illegally increase the ring precisino as needed) -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-816609631

[wbhart]

Right, so perhaps a model with no cap is better for what you want to do. We need to give this some thought, given that it has an important application re the project. Others will likely also use what we come up with if it is good.

[fieker]

On Fri, Apr 09, 2021 at 06:38:14AM -0700, wbhart wrote:

Right, so perhaps a model with no cap is better for what you want to do. We need to give this some thought, given that it has an important application re the project. Others will likely also use what we come up with if it is good.

THe "cap" is used to map exact objects in.... -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-816688224

[wbhart]

Yes, for example S(1), zero(S), one(S) etc all use the cap. We can't keep such a series ring as a ring in our system without zero and one and if we make them infinite precision objects it just pushes the can down the road, e.g. 1/f needs a "cap".

But obviously you want to work with one, zero (sometimes two, I see in code in Hecke) without it making the precision of things less just because they were added and had the precision of the cap which happened to be lower than the precision you wanted to work with.

Something that may work is to have infinite precision objects, like S(n), one, zero and just define the precision of every operation in terms of the operands. So 1/f would have the same precision as f. The only problem would be things like S(1)/S(1) which would have to raise an exception. Such a model might work without a cap, but I'm really not sure without sitting down and working it all through.

The difficulty as usual is coming up with something consistent. E.g. some identities should hold, e.g. f * (1/f) should probably be 1 and not 1 + 37658683495863948x^99 + O(x^100).

[wbhart]

Ok, I've had a longer think about this ticket, and here is what I think we should do to start with:

• Introduce polynomial(R, [a0, a1, ....], var::String="x") for constructing a polynomial from an array, similar to what we have for matrices
• Introduce abs_series(R, prec::Int, [a0, a1, ...], var::String="x") for constructing an absolute series from an array
• Introduce rel_series(R, prec::Int, val::Int, [a0, a1, ...], var::String="x") for constructing a relative series
• Introduce similar(p, R::Ring=base_ring(p)) and zero(p, R::Ring=base_ring(p)) for polynomials p (similar does not initialise)
• Introduce similar(p, R::Ring, precision::Int) and zero(p, R::Ring, precision::Int) for absolute series p
• Introduce similar(p, R::Ring, precision::Int, valuation::Int) and zero(p, R::Ring, precision::Int, valuation::Int) for relative series p
• Introduce dense_poly_type, abs_series_type and rel_series_type for type asserts and other uses

I believe all of these should be uncontroversial (let me know if not) and I can implement Series <-> Poly conversions in the polynomial_to_power_sum and power_sum_to_polynomial functions in terms of these. These should alleviate most cases where devs want to construct polynomials or series in library code without using the PolynomialRing or PowerSeriesRing constructors explicitly (which is messy and error prone).

For now, the functions above can just call PolynomialRing or PowerSeriesRing internally (in a safe way) until such time as we sort out the other issues @fieker has raised above. That will also give us time to evaluate what we are still really missing.

How does this all sound, @fieker @thofma and anyone else following this issue?

[wbhart]

Naturally some polynomial stuff already exists in Hecke and will have to be moved over. I didn't do so already as I felt that functions which automatically set the variable name are not candidates for inclusion in AbstractAlgebra/Nemo. But I now see they are useful in library code, so long as the functions never return those objects to the user.

I guess most if not all of the above functions will also need to take an optional variable name which will by default be copied from the existing object.

[thofma]

Looks fine to me. I just have the usual question of caching. For example, will polynomial(R, [a0, a1, ....], var::String="x") create a cached parent?

[wbhart]

I think we have to add an option to allow caching, but maybe make it off by default? What do we do elsewhere for default?

If someone wants to create two polynomials that are compatible, they should use polynomial for one of them perhaps and similar for the others. The latter could copy the parent over. If they want to create multiple compatible polynomials from arrays then I guess they have to turn caching on, or pull the parent out of the first polynomial they created and use that.

I see the problem, but I really don't want to get rid of check_parent throughout the library. Type checkers are also a pain in the neck but they save a lot of unintended bugs.

[wbhart]

Actually, why don't we just leave caching off and not offer an option to turn it on in these functions. We can simply explain the options in our documentation.

f = polynomial(ZZ, [2, 3, 4])
R = parent(f)
g = R([5, 6, 7])
h = f + g

[fieker]

On Mon, Apr 12, 2021 at 03:20:21AM -0700, wbhart wrote:

I think we have to add an option to allow caching, but maybe make it off by default? What do we do elsewhere for default?

If someone wants to create two polynomials that are compatible, they should use polynomial for one of them perhaps and similar for the others. The latter could copy the parent over. If they want to create multiple compatible polynomials from arrays then I guess they have to turn caching on.

I see the problem, but I really don't want to get rid of check_parent throughout the library. Type checkers are also a pain in the neck but they save a lot of unintended bugs.

in map_coeffs and lots of other functions we use parent as an optional kw parameter. If set, you get an object in THIS parent, otherwise an non-cached bla["x"] is returned.

Having said that, due to this being too frequent, in Hecke.Globals.Zx and Hecke.Globals.Qx are two default rings...

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/836#issuecomment-817687268

[wbhart]

Maybe I'm wrong. Maybe the polynomial function is too useful for the user as well and so should take an optional argument which is on by default to specify whether it is cached.

Thus in library code:

f = polynomial(ZZ, [3, 4, 5], false)
R = parent(f)
g = R([4, 5, 6])

This is one thing the developer has to get right and will therefore commonly get wrong. But I think it more closely matches the pattern we use elsewhere.

[wbhart]

Actually, maybe it needs to be an optional named argument "cached", since I don't think it makes sense to have two optional arguments, one for caching and the other for the variable. So it probably has to be the following in library code:

f = polynomial(ZZ, [3, 4, 5], cached=false)
R = parent(f)
g = R([4, 5, 6])

That's a bit irritating, but I don't see a way around it.

[wbhart]

Unless we make the variable name a named argument "var" say.

Let me know what you prefer.

[thofma]

I like that for the polynomial function. I am happy with that.

We usually have cached = true as the default.

[wbhart]

@thofma Which do you like. There were multiple possibilities. Do you want cached or var to be the named parameter?

[thofma]

I guess for consistency I would prefer cached to be named/keyword.

[wbhart]

Ok, then I will start implementing the above over the next few days or so.

[wbhart]

Currently polynomial(R, [a, b, c]) has to call AbstractAlgebra.PolynomialRing internally to get the right type of polynomial, however, I am thinking we should make it return a generic polynomial by default and then overload it in Nemo for all the other polynomial types.

The PolynomialRing function also creates the generators which we don't need, so it is not very efficient as things currently are. It's actually possible to avoid this overhead if we change the interface as I just described.

The only downside I see is that this would mean making it part of the optional interface for polynomial rings and we'd have to implement it for all our polynomial types (not that this is a big job), unless they already create generic polynomials.

What do you guys think @fieker @thofma @tthsqe12 @rfourquet ?

[thofma]

No objection from my part. But what about not returning generic polynomials by default? Then it would truly be optional?

[wbhart]

We don't do this for matrices, so I decided we will simply do the same thing for polynomials for consistency. I've implemented it for all Nemo polynomials now anyway.

[wbhart]

I believe I've implemented everything that was eventually proposed on this ticket. I'm not sold on direct conversions from polys to series and back just now, so will close out the ticket.