Ask Your Question
2

Bug report: Kernel dies after 1 hour while dividing polynomials

asked 2018-02-02 11:14:30 -0500

updated 2018-02-07 04:37:55 -0500

I want to submit an error report for SAGE 7.0 and 8.0:

The division of two polynomials in an ideal causes an kernel death after one hour of computation. Singular performs the calculation in about 10 seconds.

Sage source code:

Q.<E,F,X,Y> = QQ['E', 'F', 'X', 'Y'];

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;

J = Q.ideal(i1, i2);

R.<e,f,x,y> = QuotientRing(Q, J);

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;

print (poly1 / poly2);

Singular source code:

In Singular you cannot use the operator "/" to divide polynomials, since "non divisible terms" will be discarded and set to zero, see manual https://www.singular.uni-kl.de/Manual/4-0-3/sing_150.htm#SEC189.

So in Singular we have to use the "lift" command to calculate the quotient "poly 1 / poly2":

ring R = 0,(x,y,e,f),dp;

poly 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;

poly 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;

ideal I = i1,i2; 

ideal J = std(I);

poly 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 + 7 
0*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 - 10 
19*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;

poly 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;

poly t1 = lift(poly2,poly1)[1,1];

Singular output:

t1;

-ye9+5ye7f2-5ye8+7ye7f+15ye6f2+4ye2f6-2yef7+x2e6-10ye7-2x2e5f+14ye6f-x2e4f2+20ye5f2+26ye2f5+yef6-5yf7+3x2e5-11ye6-7x2e4f+22ye5f+3ye4f2+4x2e2f3-2x2ef4+49ye2f4+43yef5-22yf6+3x2e4+4ye5-12x2e3f-7ye4f+10x2e2f2-27ye3f2+4x2ef3+10ye2f3-5x2f4+91yef4-18yf5-x2e3+xye3+29ye4-5x2e2f+7xye2f+11ye3f+13x2ef2-2xyef2-99ye2f2-7x2f3-2xyf3-4yef3+63yf4-2x2e2+3xye2+73ye3+4x2ef+6xyef-29ye2f-2x2f2-7xyf2-161yef2+117yf3-2x3+2xye+53ye2-2xyf-106yef+53yf2+3y2
edit retag flag offensive close merge delete

Comments

Please edit the code above, using the corresponding markdown, for this, copy+paste the code, and as long it is shown as marked text

  • either press Control+K "on it",
  • or press the edit button with that "code-pre-signs" 101 and 010.

This will show the stars for instance, instead of changing something like 1*2222222222*3 into the following: 122222222223, i.e. making the text between the stars italic...

dan_fulea gravatar imagedan_fulea ( 2018-02-03 15:24:08 -0500 )edit

Note that for Sage input you don't need semi-colon ; at the end of each line.

slelievre gravatar imageslelievre ( 2018-02-05 14:57:02 -0500 )edit

Please write "Sage" rather than "SAGE", see the "From SAGE to Sage to SageMath" wiki page.

slelievre gravatar imageslelievre ( 2018-02-05 14:59:15 -0500 )edit

1 answer

Sort by ยป oldest newest most voted
1

answered 2018-02-07 09:32:21 -0500

dan_fulea gravatar image

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.)

edit flag offensive delete link more

Comments

"lift" is exactly the right functionality. It's also a horrible undiscoverable name.

nbruin gravatar imagenbruin ( 2018-02-07 14:55:06 -0500 )edit

Thank you very much for the given approach, that helps me a lot, because I still have some calculations of this kind ahead of me.

To round off this example, I would like to briefly complete the background to this question:

i2 is the general form of an elliptic curve with 17-torsion point P = (0,0) and base point O.

i1 describes the relationship between the parameters e and f and is usually referred as X_1(17).

poly1 and poly2 are meremorphic functions on the elliptic curve with the divisors

div (poly1) = 4[P] + 2[4P] +[-4P] +[-8P] +[-8P] - 8[O]

and

div (poly2) =[4P] +[-4P] - 2[O]

Thus, the quotient has the divisor

div (poly1/poly2) = 4[P] +[4P] +[-8P] - 6[O]

and can therefore be represented as a denomination-free polynomial, since the base point O is the only pole.

Patrick Reichert gravatar imagePatrick Reichert ( 2018-02-07 16:32:16 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2018-02-02 11:14:30 -0500

Seen: 97 times

Last updated: Feb 07