exact division in polynomial rings

[mantepse]

Steps To Reproduce

sage: Q.<z> = Frac(QQ['z'])
sage: R.<x,y> = Q[]
sage: r = x*y - (2*z-1)/(z^2+z+1) * x + y/z
sage: p = r * (x + z*y - 1/z^2)
sage: p // p

Expected Behavior

Should give 1.

Actual Behavior

0*x^2*y + (0*z)*x*y^2

Additional Information

No response

Environment

• OS: Ubuntu 20.04
• Sage Version: 10.5.beta8

Checklist

• I have searched the existing issues for a bug report that matches the one I want to file, without success.
• I have read the documentation and troubleshoot guide

[vincentmacri]

I tried responding the sage-devel post about this but I haven't posted there before so I think it needs moderator approval before going through.

I can reproduce this in Sage 10.4, and I found some additional weird behaviour:

• Doing p // p a second time gives 0.
• Doing p / p causes a segmentation fault and crashes Sage.

Do you get the same thing?

[mantepse]

Yes :-( also in sage cell, which has Sage 10.2

[mantepse]

Possibly related: #35428 and #35886

[saraedum]

Copying what I just wrote to sage-devel here:

Versions 9.5 - 10.4 produce the output you posted. Versions 9.0 - 9.4 take about 6 minutes and produce "1". Versions 7.0 - 8.8 do not terminate within 12 minutes.

Fwiw, you can test something like this "interactively" in bash or zsh (assuming you have a good amount of disk space):

cat - | eval `echo -n 'tee -p'; for major in $(seq 7 10); do for minor in $(seq 0 9); do echo -n " >(docker run --rm -i sagemath/sagemath:$major.$minor sage -q | sed --unbuffered s/^/\[$major.$minor\]/ )"; done; done`

[mantepse]

I inserted some print statements in MPolynomial_libsingular._floordiv_. Here is the code:

sage/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

Lines 4155 to 4204 in 209ae4c

	cdef MPolynomialRing_libsingular parent = self._parent
	cdef MPolynomial_libsingular _right = <MPolynomial_libsingular>right
	cdef ring *r = self._parent_ring
	cdef poly *quo
	cdef poly *temp
	cdef poly *p
	if _right._poly == NULL:
	raise ZeroDivisionError
	elif p_IsOne(_right._poly, r):
	return self
	if n_GetChar(r.cf) > 1<<29:
	raise NotImplementedError("Division of multivariate polynomials over prime fields with characteristic > 2^29 is not implemented.")
	if r.cf.type != n_unknown:
	if (singular_polynomial_length_bounded(_right._poly, 2) == 1
	and r.cf.cfIsOne(p_GetCoeff(_right._poly, r), r.cf)):
	p = self._poly
	quo = p_ISet(0,r)
	while p:
	if p_DivisibleBy(_right._poly, p, r):
	temp = p_MDivide(p, _right._poly, r)
	p_SetCoeff0(temp, n_Copy(p_GetCoeff(p, r), r.cf), r)
	quo = p_Add_q(quo, temp, r)
	p = pNext(p)
	return new_MP(parent, quo)
	if r.cf.type == n_Znm or r.cf.type == n_Zn or r.cf.type == n_Z2m:
	raise NotImplementedError("Division of multivariate polynomials over non fields by non-monomials not implemented.")
	count = singular_polynomial_length_bounded(self._poly, 15)
	# fast in the most common case where the division is exact; returns zero otherwise
	if count >= 15: # note that _right._poly must be of shorter length than self._poly for us to care about this call
	sig_on()
	quo = p_Divide(p_Copy(self._poly, r), p_Copy(_right._poly, r), r)
	if count >= 15:
	sig_off()
	if quo == NULL:
	if r.cf.type == n_Z:
	P = parent.change_ring(QQ)
	f = (<MPolynomial_libsingular>P(self))._floordiv_(P(right))
	return parent(sum([c.floor() * m for c, m in f]))
	else:
	sig_on()
	quo = singclap_pdivide(self._poly, _right._poly, r)
	sig_off()
	return new_MP(parent, quo)

Apparently, the result of p_Divide in line 4190 returns something "!= NULL" (the branch in 4194 is not executed). Curiously, with the print statements I don't get a crash anymore, but always 0.

[mantepse]

Copying what I just wrote to sage-devel here:

Versions 9.5 - 10.4 produce the output you posted. Versions 9.0 - 9.4 take about 6 minutes and produce "1". Versions 7.0 - 8.8 do not terminate within 12 minutes.

Fwiw, you can test something like this "interactively" in bash or zsh (assuming you have a good amount of disk space):

cat - | eval `echo -n 'tee -p'; for major in $(seq 7 10); do for minor in $(seq 0 9); do echo -n " >(docker run --rm -i sagemath/sagemath:$major.$minor sage -q | sed --unbuffered s/^/\[$major.$minor\]/ )"; done; done`

Cool! I don't have docker (yet), but that information is extremely valuable.

I'm going to see whether floordiv was robust in 9.4 and see what changed since then. In floordiv itself I don't see much change, though. r has been replaced by r.cf in several places.

I suspect that 9.4 did not use singular for these rings.

[mantepse]

Unfortunately, I cannot build 9.4 on my computer (I'm on Ubuntu 22.04.5 LTS; cysignals-1.10.3 and linbox-1.6.3.p1 fail to build) and I'm not sure whether I can afford to get docker right now.

However, I do have a suspect, which is #16567 and was merged in 9.5.beta7.

@saraedum, could you check in 9.4:

sage: Q.<z> = Frac(QQ['z'])
sage: R.<x,y> = Q[]
sage: type(R)

@miguelmarco, can you help? I have never used singular, so I really have no idea how to fix this.

[mantepse]

More info: after p_Divide in line 4190, errorreported equals 1 and sage.libs.singular.singular.error_messages is ['conversion error: denominator!= 1'].

Does that mean that singular cannot do this division?

[maxale]

Looks like duplicate of #38555

[saraedum]

I'm not sure whether I can afford to get docker right now.

I'm on Zulip if you need any assistance with that.

On 9.4:

sage: sage: Q.<z> = Frac(QQ['z'])
....: sage: R.<x,y> = Q[]
....: sage: type(R)
....:
<class 'sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict_domain_with_category'>

[miguelmarco]

This seems related (and might be easier to debug because it crashes and gives a trace):

sage: Q.<z> = Frac(QQ['z'])
sage: R.<x,y> = Q[]
sage: r = x*y - (2*z-1)/(z^2+z+1) * x + y/z
sage: p = r * (x + z*y - 1/z^2)
sage: p.quo_rem(p)