inconsistencies of .subs() for multivariate polynomials

[videlec]

sage: R.<x,y,z> = ZZ[]
sage: x.subs({x:1}).parent()
Multivariate Polynomial Ring in x, y, z over Integer Ring
sage: x.subs({x:1, y:1, z:1}).parent()
Multivariate Polynomial Ring in x, y, z over Integer Ring

compared to the univariate case

sage: R.<t> = ZZ[]
sage: t.subs({t:1}).parent()
Integer Ring

Note also that we have

sage: x.subs({x:RDF(1)}).parent()
Real Double Field

CC: @sagetrac-tmonteil @bgrenet

Component: algebra

Issue created by migration from https://trac.sagemath.org/ticket/19130

[mantepse]

See https://groups.google.com/g/sage-devel/c/vGzNJKAQWbs/m/YBRN4XOjAwAJ for some discussion.

[mantepse]

There seems to be consensus that, ideally, the parent of the output should only depend on the parents of the input.

A further question is, whether one should make a distinction between a partial substitution, i.e., where only some of the generators are substituted, and a full substitution.

One idea is to regard any partial substitution as a full substitution, where unspecified generators are substituted by themselves.

Then, the parent of the output could be the common parent of base ring and parents of the values. Doing so,

sage: R.<x,y,z> = ZZ[]
sage: x.subs({x:1}).parent()
Multivariate Polynomial Ring in x, y, z over Integer Ring

would be correct. By contrast, we would have

sage: R.<x,y,z> = ZZ[]
sage: x.subs({x:1, y:0, z:0}).parent()
Integer Ring

By analogy, substitutions in the InfinitePolynomialRing would always be partial.

@tscrim, @nbruin, @videlec, @novoselt, what do you think?

[mantepse]

However, the following doctests violates the proposed rule:

            sage: P = QQ['x,y']
            sage: x = var('x')                                                          # needs sage.symbolic
            sage: parent(P.zero().subs(x=x))                                            # needs sage.symbolic
            Symbolic Ring

Here, the substitution is partial, so given the above the result should be in Multivariate Polynomial Ring in x, y over Symbolic Ring. (Maybe replacing SR with QQbar would make a better example.)

One might argue that constants are special, but it seems that only zero is treated specially:

sage: R.<x,y> = QQ[]
sage: R(0).subs({x:QQbar(sqrt(2))}).parent()
Algebraic Field
sage: R(1).subs({x:QQbar(sqrt(2))}).parent()
Rational Field

[mantepse]

Possibly this question should be discussed together with what the parent of the result of division should be. Currently, we have

sage: R.<x,y> = QQbar[]
sage: (x/R(2)).parent()
Multivariate Polynomial Ring in x, y over Algebraic Field

sage: R.<x> = QQbar[]
sage: (x/R(2)).parent()
Fraction Field of Univariate Polynomial Ring in x over Algebraic Field

[tscrim]

In some ways this is more clear cut, but we do have some exceptions for very important reasons (mainly Laurent polynomials when dividing by a monomial). Whichever way we decide to do it, we should 100% be consistent between the uni- and multi-variate cases.

I don’t like that 0 is treated differently than constants. At the very least I expect those to behave the same.

[nbruin]

I like the idea of filling up a partial substitution to a full one. Hopefully that is not incompatible with other cases. However, isn't that a bit impractical for infinite polynomial rings? How would you "fully evaluate" an infinite polynomial ring member? Are these just not evaluable? (perhaps they are not. Perhaps one should first coerce them into a polynomial ring with all the variables that actually occur and evaluate them from there ... or can you only get an outright full evaluation if you pass it an infinite sequence?)

[mantepse]

Yes, there would be no full evaluation of InfinitePolynomial. And I should also have said that this would only be the definition, not the implementation.

[novoselt]

My main point on such issues is that it is extremely frustrating when there are discrepancies between univariate and multivariate behaviour, as well as different "special" cases like 0 and 1 demonstrated here. When the code is written to treat one case and then suddenly it breaks on another, which is logically the same, it is unpleasant. Since you can substitute elements that will expand the ring, it is natural to expect some changes. But when such a change is happening for a particular substitution, there is no reason to depend on the polynomial. In particular, maybe there is usually a substitution into "a full polynomial", but in some case it happens to be 0. That zero should work in the same way after the substitution, without any surprises with smaller rings.

[nbruin]

I don’t like that 0 is treated differently than constants. At the very least I expect those to behave the same.

I expect it's an implementation issue: with other polynomials you can just let the coercion system figure out common parents as you go along, but for the zero polynomial there is no arithmetic with the arguments involved (because there are no monomials!), so you have to synthesize the parent for 0 yourself.

One would probably get more uniform behaviour if polynomial evaluation/substitution code would analyze its evaluation arguments and construct a target parent to begin with and then coerce everything into there before doing arithmetic. I note that actually doing that may decrease performance: imagine polynomials over some very complicated ring that gets evaluated at integers. The values of the monomials could potentially be computed very efficiently. On the other hand, over a finite field that could blow up quite badly.

[mantepse]

I started to decipher the code of MPolynomial_libsingular.subs (note that, in the other cases subs actually seems to behave as discussed).

https://github.com/sagemath/sage/blob/1cd49900cda6c443df024ca79216a2716fe67710/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx#L3563-L3678

Let me give a brief summary of what I believe to understand so far: 1) line 3585 (if not change_ring:) has no effect 2) in line 3626 we return if we were unable to coerce all values into self.parent(). 3) once we are in line 3661 (if not change_ring:) we are essentially done with computing the result, and at this point change_ring is false if and only if the result is non-zero. 4) if the result is non-zero, we return it in line 3662 (return new_MP(parent, _p)) 5) Thus, the only purpose of the code after line 3664 is to determine the correct parent of zero.

What I'd like to do (for experimentation) is to return in 3662 only if len(items) < len(g), that is, if we were not doing a full substitution. Unfortunately, this makes sage crash - I do not understand enough of that code to create a polynomial with the data in _p and res_parent as parent.

I think that the loop leading to 2) is also quite inefficient - we are possibly doing unnecessary work if only the last v cannot be coerced. But that's for later.