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.

Sign up for GitHub

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

Refactor subs() of multivariate polynomials for readability and efficiency #35210

Merged
merged 3 commits into from
Apr 1, 2023
Merged

Refactor subs() of multivariate polynomials for readability and efficiency #35210

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1c436bd
Refactor subs method
kwankyu Feb 27, 2023
d49f43b
Restore the old comment
kwankyu Feb 28, 2023
8e03691
Merge branch 'develop' into pr/35210/refactor-subs-multivariate
kwankyu Mar 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
213 changes: 85 additions & 128 deletions src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2090,20 +2090,19 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
res_parent = coercion_model.common_parent(parent._base, *x)
cdef poly *res # ownership will be transferred to us in the else block
try:
# Attempt evaluation via singular.
coerced_x = [parent.coerce(e) for e in x]
except TypeError:
# give up, evaluate functional
sage_res = parent.base_ring().zero()
for (m,c) in self.dict().iteritems():
for m, c in self.dict().iteritems():
sage_res += c * mul([x[i] ** m[i] for i in m.nonzero_positions()])
else:
singular_polynomial_call(&res, self._poly, _ring, coerced_x, MPolynomial_libsingular_get_element)

if res == NULL:
return res_parent(0)
if p_LmIsConstant(res, _ring):
sage_res = si2sa( p_GetCoeff(res, _ring), _ring, parent._base )
sage_res = si2sa(p_GetCoeff(res, _ring), _ring, parent._base)
p_Delete(&res, _ring) # sage_res contains copy
else:
sage_res = new_MP(parent, res) # pass on ownership of res to sage_res
Expand Down Expand Up @@ -3401,8 +3400,7 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
- ``fixed`` - (optional) dict with variable:value pairs
- ``**kw`` - names parameters

OUTPUT:
a new multivariate polynomial
OUTPUT: a new multivariate polynomial

EXAMPLES::

Expand Down Expand Up @@ -3462,7 +3460,7 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):

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

We are catching overflows::
Expand Down Expand Up @@ -3513,163 +3511,122 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
c*x^2*y + d*x^2 + (2*c)*x*y + (2*d)*x + c*y + d

"""
cdef int mi, i, need_map, try_symbolic
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)
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 *to_id = idInit(_ring.N, 1)
cdef ideal *from_id
cdef ideal *res_id
need_map = 0
try_symbolic = 0
cdef dict gd

if _p == NULL:
# the polynomial is 0. There is nothing to do except to change the
# ring
try_symbolic = 1

if not try_symbolic and 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")
try:
v = parent.coerce(v)
except TypeError:
try_symbolic = 1
break
_f = (<MPolynomial_libsingular>v)._poly
if p_IsConstant(_f, _ring):
singular_polynomial_subst(&_p, mi-1, _f, _ring)
else:
need_map = 1
degree = <unsigned long>p_GetExp(_p, mi, _ring) * <unsigned long>p_GetMaxExp(_f, _ring)
if degree > _ring.bitmask:
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 OverflowError("exponent overflow (%d)"%(degree))
to_id.m[mi-1] = p_Copy(_f, _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:
# polynomial becomes 0 after some substitution
try_symbolic = 1
change_ring = True
break

cdef dict gd

if not try_symbolic:
gd = parent.gens_dict(copy=False)
for m,v in kw.iteritems():
m = gd[m]
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")
i = g.index(m)
try:
v = parent.coerce(v)
except TypeError:
try_symbolic = 1
break
_f = (<MPolynomial_libsingular>v)._poly
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, mi-1, _f, _ring)
singular_polynomial_subst(&_p, i, _f, _ring)
else:
if to_id.m[mi-1] != NULL:
p_Delete(&to_id.m[mi-1],_ring)
to_id.m[mi-1] = p_Copy(_f, _ring)
degree = <unsigned long>p_GetExp(_p, mi, _ring) * <unsigned long>p_GetMaxExp(_f, _ring)
if degree > _ring.bitmask:
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))
need_map = 1

if _p == NULL:
# the polynomial is 0
try_symbolic = 1
break
raise OverflowError("exponent overflow (%d)" % (degree,))
to_id.m[i] = p_Copy(_f, _ring)

if 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)
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
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]
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
from_id.m[0] = NULL
res_id.m[0] = NULL

id_Delete(&from_id, _ring)
id_Delete(&res_id, _ring)
id_Delete(&from_id, _ring)
id_Delete(&res_id, _ring)

id_Delete(&to_id, _ring)

if not try_symbolic:
if not change_ring:
return new_MP(parent,_p)

# now as everything else failed, try to do it symbolically with call

cdef list g = list(parent.gens())

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:
raise TypeError("key does not match")
else:
raise TypeError("keys do not match self's parent")

g[mi-1] = v
# finally change the parent of the result

gd = parent.gens_dict(copy=False)
for m,v in kw.iteritems():
m = gd[m]
for i from 0 < i <= _ring.N:
if p_GetExp((<MPolynomial_libsingular>m)._poly, i, _ring) != 0:
mi = i
break
if i > _ring.N:
raise TypeError("key does not match")
res_parent = coercion_model.common_parent(parent._base, *[v for _, v in items])

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This appears to be more general than the old code that tries to coerce every v into parent._base.
Is this just incidental, or is there a use case for this greater generality?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps there is no difference in this respect. Both first try to coerce every v into parent._base. If fails, the old code finds the common parent in return self(*g) (the last line) via the constructor __init__ while the new code does the same in the place. However, for some paths (substitutions with constant values), the new code does not make a duplicate call to singular_polynomial_call by avoiding going through the constructor. (if I remember correctly...)

Copy link
Collaborator Author

@kwankyu kwankyu Feb 28, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps it is the substitution with constant values that results in zero that makes the duplicate call.

g[mi-1] = v
if _p == NULL:
return res_parent(0)

return self(*g)
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

def monomials(self):
"""
Expand Down