Bug report: Kernel dies after 1 hour while dividing polynomials

edit

The division method is generally implemented, it may work or not work in special circumstances. Here are some simple examples:

sage: R.<U> = PolynomialRing(ZZ)
sage: (U^2-1)/(U+1)
U - 1
sage: _.parent()
Fraction Field of Univariate Polynomial Ring in U over Integer Ring
sage: (U^2-1)._div_(U+1)
U - 1
sage: _.parent()
Fraction Field of Univariate Polynomial Ring in U over Integer Ring

and

sage: parent(_)
Rational Field
sage: 4._div_(2)
2
sage: parent(_)
Rational Field

This is unrelated, but there is an aspect that can be made transparent in these situations. The divisors, the "denominators" are not invertible in the given ring, so sage passes automatically to the quotient field. Then it uses the simplifications in this new context. (We do not get the fraction 4/2, instead the better representation 2. But not as an element of the ring ZZ.) In our case, there is no obvious simplification algorithm, and we would not be happy with the fraction representation. We really want $poly_1/poly_2$ in the quotient ring, not in the fraction field of it! Above, this is not the case, so we already expect "more" (in a special situation).

Let us consider now an example with a quotient ring.

sage: R.<U> = PolynomialRing(QQ)
sage: S.<u> = R.quotient( U^20-1 )
sage: one = S(1)
sage: one / u^5
u^15
sage: one._div_??
Signature: one._div_(right)
Source:   
    def _div_(self, right):
        """
        Return the quotient of two polynomial ring quotient elements.

        EXAMPLES::

            sage: R.<x> = PolynomialRing(QQ)
            sage: S.<a> = R.quotient(x^3-2)
            sage: (a^2 - 4) / (a+2)
            a - 2
        """
        return self * ~right
File:      /usr/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_quotient_ring_element.py
Type:      instancemethod
sage:

So the used method _div_ needs the inverse element. Let us then construct an example without such an inverse in S. The simplest:

sage: S(u-1)/S(u-1)

runs into a

ZeroDivisionError: element u - 1 of quotient polynomial ring not invertible

I do not know the reason for the long computation in our new, more involved setting, but the implementation seems to be done so that in some cases it gives an answer, the answer is then correct, but it is not written to cover all cases. (It would be hard to do this in all possible cases / classes, and to propagate / inherit the _div_ method in a coherent way.)

So let us try to give a valid alternative simple solution, that works in this case, by slightly changing mathematically the point of view.

First of all, let us copy+paste the given situation, and slightly change it. The change will be that we do not pass to the quotient ring.

Q.<E,F,X,Y> = PolynomialRing( QQ )

i1 = F^4 + (E^3+E^2-E+2)*F^3 + (E^3-3*E+1)*F^2 - (E^4+2*E)*F + (E^3+E^2);

i2 = Y^2 \
     + (E^3+E^2*(3*F+2)-E*(F^2-2*F-1)-F*(F^2+3*F+1))*X*Y \
     + (F*(E+1)*(E-F)*(E+F+1)^2*(E^2+E-F)*(E^2+E*F+E-F^2-F))*Y \
     - X^3 - (F*(E+1)*(E-F)*(E+F+1)*(E^2+E-F))*X^2

# I = Q.ideal(i1, i2);

e,f,x,y = E,F,X,Y

poly1 = ( x^4
          +(-e^2*f^4+e^4*f+2*e^3*f^2-e*f^4-2*e*f^3+2*f^4+2*e^2*f
            -6*e*f^2+f^3+5*e^2-5*e*f+11*f^2-17*e+18*f+16)*x^3
          +(-2*e^3-2*e^2*f+e*f^2+f^3-3*e^2+2*f^2-e+f)*x^2*y
          +(-6*e^2*f^5+11*e*f^6-3*f^7+3*e^6-6*e^5*f+17*e^4*f^2
            -69*e^2*f^4+5*e*f^5+17*f^6-5*e^5+26*e^4*f+83*e^3*f^2
            -106*e^2*f^3-104*e*f^4+71*f^5-28*e^4+66*e^3*f+57*e^2*f^2
            -149*e*f^3+54*f^4-55*e^3+78*e^2*f+9*e*f^2-32*f^3-35*e^2
            +70*e*f-35*f^2)*x^2
          +(-e^2*f^7-e^7*f-3*e^6*f^2-2*e^2*f^6-e*f^7-5*e^6*f
            -e^5*f^2-6*e^2*f^5-6*e*f^6+2*f^7+3*e^5*f-2*e^4*f^2
            -2*e^2*f^4-15*e*f^5+2*f^6-2*e^5-7*e^4*f-4*e^3*f^2-9*e^2*f^3
            -4*e*f^4+3*e^4+20*e^3*f-43*e^2*f^2-21*e*f^3+28*f^4-e^3
            -32*e^2*f-61*e*f^2+61*f^3-22*e^2-53*e*f+59*f^2-16*e+16*f)*x*y
          +(e^6+2*e^5*f+3*e^5-4*e^3*f^2-2*e^2*f^3+2*e*f^4+f^5+3*e^4-5*e^2*f^2
            +2*f^4+3*e^3-3*e^2*f-e*f^2+4*f^3-3*e^2+e*f-9*f^2+17*e-18*f-16)*y^2
          +(11*e*f^9-3*f^10-3*e^9-5*e^7*f^2-101*e^2*f^7+30*e*f^8+11*f^9
            -e^8-32*e^7*f-74*e^6*f^2-422*e^2*f^6-117*e*f^7+121*f^8+35*e^7
            -39*e^6*f-134*e^5*f^2-486*e^2*f^5-660*e*f^6+276*f^7+141*e^6
            -106*e^5*f-226*e^4*f^2+384*e^2*f^4-1019*e*f^5+36*f^6+189*e^5
            -367*e^4*f-681*e^3*f^2+1382*e^2*f^3+130*e*f^4-884*f^5+61*e^4
            -1019*e^3*f+565*e^2*f^2+1914*e*f^3-1521*f^4-251*e^3
            -499*e^2*f+1751*e*f^2-1001*f^3-231*e^2+462*e*f-231*f^2)*y
          )

poly2 = ( x+e^5*f+e^4*f^2-e^2*f^4+2*e^4*f+e^3*f^2
          -3*e^2*f^3+f^5+3*e^3*f-3*e^2*f^2-3*e*f^3
          +3*f^4+e^3+e^2*f-5*e*f^2+3*f^3+e^2-2*e*f+f^2 )

J = Q.ideal( [ i1, i2, poly2 ] )
poly1 in J

So instead of working in $$(Q\text{ mod } i_1,i_2)\text{ mod } poly_2\ ,$$ we work directly in $$Q/J\ =\ Q\text{ mod } i_1,i_2, poly_2\ .$$ Then indeed, $poly_1$ lies in $J$,

sage: J = Q.ideal( [ i1, i2, poly2 ] )
....: poly1 in J
....: 
True

and we can explicitly ask for its representation in terms of the generators of $J$:

poly1.lift(J)

The coefficient corresponding to poly2 corresponds to the quotient modulo $i_1,i_2$. It is:

sage: poly1.lift(J)[2].factor()
(-1/6) * (3*E^8*F^2*Y + 3*E^7*F^3*Y - 6*E^6*F^4*Y - 3*E^5*F^5*Y - 3*E^4*F^6*Y + 6*E^3*F^7*Y + 6*E^9*Y - 21*E^7*F^2*Y + 3*E^6*F^3*Y + 72*E^5*F^4*Y - 60*E^4*F^5*Y - 15*E^3*F^6*Y + 9*E^2*F^7*Y - 3*E*F^8*Y - 3*E^5*F^2*X^2 + 30*E^8*Y - 54*E^7*F*Y - 165*E^6*F^2*Y + 272*E^5*F^3*Y + 127*E^4*F^4*Y - 265*E^3*F^5*Y - 18*E^2*F^6*Y + 3*E*F^7*Y - 6*E^6*X^2 + 15*E^5*F*X^2 - 39*E^3*F^3*X^2 + 3*E^2*F^4*X^2 - 3*E*F^5*X^2 + 57*E^7*Y - 233*E^6*F*Y - 50*E^5*F^2*Y + 636*E^4*F^3*Y - 408*E^3*F^4*Y - 178*E^2*F^5*Y - 18*E*F^6*Y - 18*E^5*X^2 + 90*E^4*F*X^2 - 54*E^3*F^2*X^2 - 57*E^2*F^3*X^2 + 18*E*F^4*X^2 - 6*F^5*X^2 + 60*E^6*Y - 287*E^5*F*Y + 383*E^4*F^2*Y - 103*E^3*F^3*Y - 161*E^2*F^4*Y + 26*E*F^5*Y - 103*F^6*Y - 18*E^4*X^2 + 78*E^3*F*X^2 - 75*E^2*F^2*X^2 + 30*E*F^3*X^2 - 15*F^4*X^2 + 89*E^5*Y - 55*E^4*F*Y - 281*E^3*F^2*Y + 70*E^2*F^3*Y + 559*E*F^4*Y - 434*F^5*Y + 3*E^3*F*X*Y + 9*E^2*F^2*X*Y - 3*E*F^3*X*Y - 3*F^4*X*Y - 33*E^3*X^2 + 30*E^2*F*X^2 + 39*E*F^2*X^2 - 36*F^3*X^2 + 157*E^4*Y - 372*E^3*F*Y - 286*E^2*F^2*Y + 1112*E*F^3*Y - 611*F^4*Y - 9*E^3*X*Y - 45*E^2*F*X*Y + 21*E*F^2*X*Y + 6*F^3*X*Y - 27*E^2*X^2 + 54*E*F*X^2 - 27*F^2*X^2 - 3*F*X^3 + 43*E^3*Y - 418*E^2*F*Y + 707*E*F^2*Y - 332*F^3*Y - 24*E^2*X*Y - 39*E*F*X*Y + 48*F^2*X*Y + 15*X^3 - 52*E^2*Y + 104*E*F*Y - 52*F^2*Y - 15*E*X*Y + 15*F*X*Y + 3*F*Y^2 - 21*Y^2)

We can check this:

sage: q = poly1.lift(J)[2]
sage: S.<e,f,x,y> = Q.quotient( [i1, i2] )
sage: S( poly1 - poly2*q )
0

In my opinion, the situation is not a bug, but a (rather special) case where the algorithms can not go in details. (In the same manner, working with polynomial rings over polynomial rings produce problems. But of course, this is not our case, just a parallel.)