Polynomial interface

[asbuch]

I have been looking at the polynomial part of AbstractAlgebra and Nemo. The speed is fantastic! But I find that the interface could be more consistent. There are multiple types of polynomial implementations with similar capabilities, but their APIs are not consistent. It would be great if one could replace PolynomialRing(ZZ, "x") with PolynomialRing(ZZ, ["x"]) or LaurentPolynomialRing(ZZ, "x") or LaurentPolynomialRing(ZZ, ["x"]), etc. without breaking any code, say if one suspected that sparse polynomials might do the job faster or if 1/x is needed for new functionality. Maybe some functions don't make sense for all rings, but when they do make sense, it would be great to have them available with consistent names.

I also wonder why shifting of Laurent polynomials and series is split into shift_left and shift_right. How about shift(p, s), where s can have either sign?

Here is a start of some function signatures that ought to make sense for all types of polynomials and series:

PS, s = PolynomialRing(ZZ, "s")
P1, (z,) = PolynomialRing(ZZ, ["z"])
P2, (x,y) = PolynomialRing(ZZ, ["x", "y"])
PU = UniversalPolynomialRing(ZZ)
u = gen(PU, "u")
LS, t = LaurentPolynomialRing(ZZ, "t")
L1, (c,) = LaurentPolynomialRing(ZZ, ["c"])
L2, (a,b) = LaurentPolynomialRing(ZZ, ["a", "b"])

gen(polyring)

gen(PS) == s
gen(P1) == z
gen(P2) # error
gen(PU) # u
gen(LS) == t
gen(L1) == c
gen(L2) # error

gen(polyring, i::Int)

gen(PS, 1) == s
gen(PS, 2) # error
gen(P1, 1) == z
gen(P2, 2) == y
gen(PU, 1) == u
gen(LS, 1) == t
gen(L1, 1) == c
gen(L2, 1) == a

gen(polyring, v::Union{String,Symbol})

gen(PS, "s") == s
gen(PS, :s) == s
gen(P1, "z") == z
gen(P2, :x) == x
gen(P2, "y") == y
gen(PU, "u") == u
gen(LS, :t) == t
gen(L2, "b") == b

gens(polyring, vars::Vector)

gens(PS, ["s"]) == (s,)
gens(P1, ["z"]) == (z,)
gens(P2, ["y"]) == (y,)
gens(P2, ["y", "x"]) == (y, x)
gens(PU, ["u"]) == (u,)
gens(LS, ["t"]) == (t,)
gens(L1, ["c"]) == (c,)
gens(L2, ["a"]) == (a,)
gens(L2, ["a", "b"]) == (a, b)

coeff(poly, e::Int) # get single coefficient

coeff((s+2)^3, 2) == 6
coeff(z, 1) == 1
coeff(x, 1) # error
coeff(u, 1) == 1 # could also be an error, since number of variables can change.
coeff((t - 2t^-1)^3, -1) == 12
coeff((c - 2c^-1)^3, -1) == 12
coeff(a, 1) # error

coeff(poly, exp::Vector{Int}) # get single coefficient

coeff(s, [1]) == 1
coeff(z, [1]) == 1
coeff(x, [1]) # error
coeff(x, [1,0]) == 1
coeff(u, [1]) == 1
coeff(t, [1]) == 1
coeff((a+b)^2, [1,1]) == 2

coeff(poly, vars::Vector, exp::Vector{Int})

coeff(s, [s], [1]) == 1
coeff(u, [u], [1]) == 1
coeff(xy, [x], [1]) == y
coeff(u, [u], [1]) == 1
coeff((c-c^-1)^3, [c], -1) ==
coeff((2a + 3b^-1)^3, [a,b], [1,-2]) == 54
coeff((2a + 3b^-1)^3, [a], [1]) == 54a*b^-2

shift(poly, s::Int)

shift(s, 1) == s^2
shift(s, -1) == 1
shift(s, -2) # error
shift(z, 1) == z^2
shift(x, 1) # error
shift(u, 1) == u^2
shift(t, 1) == t^2
shift(t, -3) == t^-2
shift(c, -2) == c^-1
shift(a, 1) # error

shift(poly, s::Vector{Int})

shift(s, [1]) == s^2
shift(z, [1]) == z^2
shift(x, [1]) # error
shift(x, [-1,2]) == y^2
shift(u, [1]) == u^2
shift(t, [1]) == t^2
shift(c, [1]) == c^2
shift(b, [1,-3]) == a*b^-2

shift(poly, vars::Vector, s::vector{Int})

shift(s, [s], [1]) == s^2
shift(z, [z], [1]) == z^2
shift(x, [x], [1]) == x^2
shift(x, [y,x], [3,-1]) == y^3
shift(u, [u], [1]) == u^2
shift(t, [t], [1]) == t^2
shift(c, [c], [1]) == c^2
shift(a, [b], [-1]) == ab^-1
shift(a, [a,b], [-3,3]) == a^-2b^3

[thofma]

gen(P2) # error

I don't know what it should do except throwing an error?

[fingolfin]

I think this was not a list of how things works, but rather of how @asbuch things they should work; in this case, the current implementation and what they think is appropriate simply happen to match.

Where we get divergence is e.g. for `gen(P1), were currently one gets the exact same error:

julia> gen(P1)
ERROR: MethodError: no method matching gen(::FmpzMPolyRing)

And dually, for the univariate ring:

julia> gens(PS, 1)
ERROR: MethodError: no method matching gens(::FmpzPolyRing, ::Int64)

They are arguing now that multivariate polynomial rings that happen to be in just one variable should "emulate" univariate polynomial rings as much as possible, so that one can quickly switch between a univariate and multivariate implementation, without changing code.

So if I understand right, in essence they are suggesting something like this be added:

function AbstractAlgebra.gen(R::MPolyRing)
  nvars(R) == 1 && return gen(R, 1)
  error("gen is not supported for multivariate rings in more than one variable")
end

or perhaps even this (but IMHO this one is a really bad idea as it is too easy to use it incorrectly by mistake):

AbstractAlgebra.gen(R::MPolyRing) = gen(R, 1)

Dually, one could add:

function AbstractAlgebra.gen(R::PolyRing, i::Int)
  i == 1 && return gen(R)
  error("univariate polynomial rings only have one variable, not $i")
end

And so on

[fieker]

To me the diskussion is in the wrong direction. If there is demand, we can add a sparse-univariate type, following the univariate interface. However, "just swapping in and out" is always tricky: in the dense case, iteration is usually over all coefficients (by degree), while in the sparse case, it is usually only over the non-zero ones. Asymptotically fast methods that apply to the dense case are at best not harmful in the sparse case, ...
To randomly support the univariate interface if the number of variables happens to be one is also tricky. I expect to have ambiguties with e.g. evaluation.
All of this is fine if I create the ring directly, interactively as then I know that I am univariate od not. However programmatically, i might just happen to hit the univariate case...

[asbuch]

AbstractAlgebra does a great job of polymorphism as long as we speak about operations that apply to all rings, such as addition, multiplication, comparison, etc. But for operations that only apply to polynomials, it suddenly seems like each polynomial ring implementation has its own notation, so polymorphism is non-existent in that context. This is a shame!

Lately I have been working on some code that uses multivariate polynomials, but I had the idea that things might be faster if I changed the algorithm to work with Laurent polynomials (because then some tricks are possible that will probably justify the small overhead of Laurent polynomials). Unfortunately it is not so easy to try this out, since the operations I use for polynomials are not available or have different names for Laurant polynomials.

It would be great if somebody would think through all the operations that apply to polynomials and similar objects, and then ensure that they are available under consistent names in all the cases where they make sense. @fieker mentions evaluation, this function is certainly also among those where I wish the notation was more consistent. I don't see any problem with fixing this. If x is the only generator of a polynomial ring, then I suggest that this should work:

evaluate(x^3, 2) == evaluate(x^3, [x], [2]) == 8

I would also like this to work if x is the generator of a Laurent polynomial ring, with an error thrown if the evaluation results in division by zero. Same for elements of the fraction field of a polynomial ring for that matter.

Regarding the comment, "asymptotically fast methods that apply to the dense case are at best not harmful in the sparse case". Is the intention to prevent people from possibly writing inefficient code? I think there will always be lots of opportunities to do this. For example, one could use Z[x,y]/(x*y-1) instead of Z[x,1/x], probably with much worse results. The user should be empowered to make good choices about which implementations to use, not prevented from experimenting. I am sure there are many border cases where it is not so clear if dense or sparse polynomials will be faster, and where the simplest approach is to try both options. This is analogous to trying a different monomial order if a Groebner operation is too slow, I hope nobody would argue that this should be difficult to do?

[thofma]

I agree that we should empower the users to experiment with the available tools. I am also in favour of most of your suggestions, in particular since we have no sparse univariate polynomials and I doubt this will change any time soon. What would be your suggestion to something like iterating over coefficients? At the moment (and this won't change) we have this:

julia> collect(coefficients(f))
3-element Vector{fmpq}:
 0
 0
 1

julia> Qx, (x,) = @PolynomialRing(QQ, x[1:1])
(Multivariate Polynomial Ring in x[1] over Rational Field, fmpq_mpoly[x[1]])

julia> f = x^2
x[1]^2

julia> collect(coefficients(f))
1-element Vector{fmpq}:
 1

[asbuch]

I think this behavior is ok, as long as one can call exponent_vectors(f) when f is a dense polynomial, and its output is consistent with coefficients(f). One could consider giving both iterators an optional "all" argument to control if zero coefficients should be skipped:

coefficients(poly, all::Bool)
exponent_vectors(poly, all::Bool)

Both should result in errors if parent(poly) has more than one variable and all==true. Or maybe they could be allowed if is_univariate(poly)==true.

Let me continue with some suggestion that ignore breaking existing behavior. For multivariate polynomials there are four iterators:

coefficients(p::MPoly)
monomials(p::MPoly)
terms(p::MPoly)
exponent_vectors(a::MPoly)

I would be unlikely to ever use monomials() and terms(), since they have to create a new polynomial object for each iteration. I think it would be better to let terms(f) return the same as zip(coefficients(f), exponent_vectors(f)). Maybe this would be slightly faster than the zip method, depending on how ingenious Julia's optimizer is? Of course, this iterator could also be called something else to avoid breaking anything.

The existing behavior of terms could be recovered if a polynomial ring got a "term" constructor, for example:

R,(x,y,z) = PolynomialRing(ZZ, ["x","y","z"])
term(R, 7, [1,2,3]) == 7*x*y^2*z^3
term(R, 7, [x,z], [1,2]) == 7*x*z^2

Then the following might be as fast as using the current terms() iterator:

for (c, e) in terms(f) 
    println(term(R, c, e))
end

This terms() iterator would also make perfect sense for a dense polynomial ring, where the natural behavior would be to skip zero coefficients.

[tthsqe12]

for what it is worth, there is an incomplete sparse univariate in AA.

[fieker]

I agree it would be nice to have

• sparse and dense univariates
• sparse multivariates
• sparse and dense univariate Laurent
• sparse multivariate Laurent
• sparse universal
  I don't see any good reason for adding dense multivariate objects at this point.

But I think it is a bad idea, unless we absolutely have to, to treat implicitly mulitivariate types where there is only one variable, as univariate.
I fear this is going to mask programming errors - when the wrong interface is used by accident. It is also going to add many signatures to multivariates that intrinsically do not apply to that type (derivative(f) is not (well) defined in general for mpoly types, gen(R) is not defined, ...)
I see less problems in extending the mpoly interface to u-polys: gens(Zx) -> [x], gen(Zx, 1) -> x, gen(Zx, 2) -> error

Discussion is also about access to non-existing coefficients: in all(?) dense versions, non-existing coeffs are always accessible and yield zero. In the multivariate world, non-existing coeffs are errors. If you want a generic interface, this has to be along coeffs + exponents, however, I have a feeling that using this interface for dense polys means one is doing things badly. A user wants this, their choice.

[asbuch]

The issue with many signatures would be improved if some abstract type was inserted between Ring and all the rings with polynomial behavior. Maybe GeneratedRing or VariableRing or PolynomialLikeRing? Hopefully there is a better name! I already suggested this in the discussion of some small pull request. This would also make it easier to write new methods that apply to many types of polynomials and series. And it would make it easier to keep the APIs consistent in the future as they keep evolving, by encouraging anybody introducing a new function to at least consider which types of polynomials it makes sense for.

Regarding the fear of programming errors, the same argument can be applied to polymorphism in general. Polymorphism is about giving the user choices, and users with choices may make the wrong choices. So polymorphism should be avoided or made complicated???

About derivative: If f is a polynomial in a single variable x, then derivative(f), derivative(f, 1), and derivative(f, x) should return the same result. I believe this is consistent with both the dense and sparse APIs.

I don't think non-existent coefficients result in errors anywhere, at least this works (thankfully!):

R, (x,y) = PolynomialRing(ZZ, ["x", "y"])
coeff(x*y, [1,1]) == 1
coeff(x*y, [1,0]) == 0

Finally, for iterating over coefficients, I think users of univariate polynomials (dense or sparse) might be better served if they got an exponents iterator that returns the integer exponents that could also be obtained from exponent_vectors (wrapped in vectors). For a dense univariate polynomial, this could be as simple as this:

exponents(p::PolynomialElem) = 0:degree(p) 

Any application of coefficients(p) that cares about the relationship between coefficients and exponents would probably need to deal with something like 0:degree(p) anyway, so in many cases it might be possible to write code that works for both sparse and dense polynomials without any penalty in the dense case. Of course, exponent_vectors(p) should also be available, just in case it may be desirable to add more variables to the ring another day.