Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

inconsistencies of .subs() for multivariate polynomials #19130

Open
videlec opened this issue Sep 4, 2015 · 10 comments
Open

inconsistencies of .subs() for multivariate polynomials #19130

videlec opened this issue Sep 4, 2015 · 10 comments

Comments

@videlec
Copy link
Contributor

videlec commented Sep 4, 2015

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
Copy link
Contributor

mantepse commented Apr 5, 2024

@mantepse
Copy link
Contributor

mantepse commented Apr 5, 2024

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
Copy link
Contributor

mantepse commented Apr 5, 2024

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
Copy link
Contributor

mantepse commented Apr 5, 2024

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
Copy link
Collaborator

tscrim commented Apr 5, 2024

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
Copy link
Contributor

nbruin commented Apr 5, 2024

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
Copy link
Contributor

mantepse commented Apr 5, 2024

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
Copy link
Member

novoselt commented Apr 5, 2024

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
Copy link
Contributor

nbruin commented Apr 5, 2024

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
Copy link
Contributor

mantepse commented Apr 5, 2024

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

cdef int mi, i
cdef bint change_ring = False # indicates the need to change the parent ring
cdef bint need_map = False
cdef unsigned long degree = 0
cdef MPolynomialRing_libsingular parent = self._parent
cdef ring *_ring = parent._ring
if (_ring != currRing):
rChangeCurrRing(_ring)
cdef poly *_p = p_Copy(self._poly, _ring)
cdef poly *_f
cdef ideal *to_id = idInit(_ring.N, 1)
cdef ideal *from_id
cdef ideal *res_id
cdef dict gd
cdef list g = list(parent.gens())
cdef list items = []
if not change_ring:
if fixed is not None:
for m, v in fixed.items():
if isinstance(m, (int, Integer)):
mi = m + 1
elif isinstance(m, MPolynomial_libsingular) and m.parent() is parent:
for i from 0 < i <= _ring.N:
if p_GetExp((<MPolynomial_libsingular> m)._poly, i, _ring) != 0:
mi = i
break
if i > _ring.N:
id_Delete(&to_id, _ring)
p_Delete(&_p, _ring)
raise TypeError("key does not match")
else:
id_Delete(&to_id, _ring)
p_Delete(&_p, _ring)
raise TypeError("keys do not match self's parent")
items.append((g[mi - 1], v))
if kw:
gd = parent.gens_dict(copy=False)
for m, v in kw.items():
items.append((gd[m], v))
for m, v in items:
if _p == NULL:
change_ring = True
break
i = g.index(m)
try:
v = parent.coerce(v)
except TypeError: # give up, evaluate symbolically
id_Delete(&to_id, _ring)
p_Delete(&_p, _ring)
gg = list(parent.gens())
for m, v in items:
gg[g.index(m)] = v
y = parent.base_ring().zero()
for (m,c) in self.dict().items():
y += c*mul([ gg[i]**m[i] for i in m.nonzero_positions()])
return y
_f = (<MPolynomial_libsingular> v)._poly
if p_IsConstant(_f, _ring):
singular_polynomial_subst(&_p, i, _f, _ring)
else:
need_map = True
degree = (<unsigned long> p_GetExp(_p, i + 1, _ring)) * (<unsigned long> p_GetMaxExp(_f, _ring))
if degree > _ring.bitmask:
id_Delete(&to_id, _ring)
p_Delete(&_p, _ring)
raise OverflowError("exponent overflow (%d)" % (degree,))
to_id.m[i] = p_Copy(_f, _ring)
if not change_ring and need_map:
for mi from 0 <= mi < _ring.N:
if to_id.m[mi] == NULL:
to_id.m[mi] = p_ISet(1,_ring)
p_SetExp(to_id.m[mi], mi + 1, 1, _ring)
p_Setm(to_id.m[mi], _ring)
from_id = idInit(1, 1)
from_id.m[0] = _p
rChangeCurrRing(_ring)
res_id = fast_map_common_subexp(from_id, _ring, to_id, _ring)
_p = res_id.m[0]
from_id.m[0] = NULL
res_id.m[0] = NULL
id_Delete(&from_id, _ring)
id_Delete(&res_id, _ring)
id_Delete(&to_id, _ring)
if not change_ring:
return new_MP(parent,_p)
# finally change the parent of the result
res_parent = coercion_model.common_parent(parent._base, *[v for _, v in items])
if _p == NULL:
return res_parent(0)
if p_LmIsConstant(_p, _ring):
res = si2sa(p_GetCoeff(_p, _ring), _ring, parent._base)
p_Delete(&_p, _ring) # safe to delete; res contains copy
else:
res = new_MP(parent, _p)
if parent(res) is not res_parent:
res = res_parent(res)
return res

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

6 participants