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

Sage crashes when inverting/dividing large number field elements #20693

Closed
sagetrac-ehlen mannequin opened this issue May 27, 2016 · 61 comments
Closed

Sage crashes when inverting/dividing large number field elements #20693

sagetrac-ehlen mannequin opened this issue May 27, 2016 · 61 comments

Comments

@sagetrac-ehlen
Copy link
Mannequin

sagetrac-ehlen mannequin commented May 27, 2016

This ticket used to be about a crash that occurred when computing newforms for a certain character of modulus 23 in sage 7.2.
Here's how to reproduce it (you have to wait 10 minutes or so until the crash happens):

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: N=Newforms(c,6, names='a')

It turned out that this was due to NTL running out of FFT primes when inverting number field elements with humongous denominators. Moreover, it turned out that we only ran into this problem in the example (and other examples in the comments) because the function _invert_c_() of a number field element was doing unnecessary work.

Component: number fields

Author: Stephan Ehlen

Branch/Commit: eb3da68

Reviewer: Peter Bruin, Fredrik Stromberg

Issue created by migration from https://trac.sagemath.org/ticket/20693

@sagetrac-ehlen sagetrac-ehlen mannequin added this to the sage-7.3 milestone May 27, 2016
@sagetrac-ehlen

This comment has been minimized.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented May 27, 2016

comment:2

Let me add something that might be helpful in finding the bug:
I get the same / a similar crash when doing the following, which is how I discovered it in the first place.
It confirms what the crash looks like above, that the real bug lies somewhere in the linear algebra and/or number field code.

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: M=ModularSymbols(c,6,sign=1)
sage: S=M.cuspidal_subspace().new_subspace()
sage: A=D[0]
sage: v  = A.dual_eigenvector(names='a', lift=False)

And crash...

------------------------------------------------------------------------
0   signals.so                          0x00000001067855c5 print_backtrace + 37
------------------------------------------------------------------------
Unhandled SIGABRT: An abort() occurred.
This probably occurred because a *compiled* module has a bug
in it and is not properly wrapped with sig_on(), sig_off().
Python will now terminate.
------------------------------------------------------------------------
[1]    7238 abort      ~/Documents/math/devel/sage-git/sage

And let me also add that I tried a binary distribution of sage 7.1 on yet another machine as well with the Newforms(...) function and got the same crash.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented May 29, 2016

comment:4

I would like to add a remark: something really extremely bad happened with the implementation of relative Number Fields.
It was always slow but now, and this might be related to this crash here, spaces that I computed using older versions of sage (6.1 to be concrete and give a pointer) within minutes or up to an hour or a bit more, do not finish to be computed within a day now. A concrete example is Gamma0(19), character [zeta18^2], weight 14. I had the modular symbols space computed with sage 6.1 and stored the cputime used to compute it: 3780s. And now it didn't finish within 7 hours. I believe this happens in the very same function computing the dual_eigenvector. To be precise, it happens when coercing the elements of the base ring into the extension created within that code. I'm not sure if it is related to the crash or if it is a different problem but both of the problems are in fact problems with relative number fields, it seems

@nbruin
Copy link
Contributor

nbruin commented May 30, 2016

comment:5

I confirm the crash and I'm getting the same traceback.

For regressions in relative number field performance: it would be nice to have a smaller example where both 6.1 and 7.2 run in reasonable time. We can then just profile the code. There's a good chance that something will be sticking out there, leading to the place of the regression.

@jdemeyer
Copy link

comment:6

@sagetrac-ehlen: are you writing code or intend to write code to fix this? For me, seeing an Author filled in is a good reason not to investigate the bug.

@jdemeyer
Copy link

comment:7

Replying to @sagetrac-ehlen:

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: M=ModularSymbols(c,6,sign=1)
sage: S=M.cuspidal_subspace().new_subspace()
sage: A=D[0]
sage: v  = A.dual_eigenvector(names='a', lift=False)

For me, this gives

AttributeError: 'DirichletGroup_class_with_category.element_class' object has no attribute 'dual_eigenvector'

and no crash...

@jdemeyer
Copy link

Replying to @sagetrac-ehlen:

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: N=Newforms(c,6, names='a')

I let this run for a few minutes and didn't get a crash. Do you have a simpler crashing example?

@pjbruin
Copy link
Contributor

pjbruin commented May 30, 2016

comment:9

To investigate the crash, it may help to add sig_on()...sig_off() around the NTL calls in sage.rings.number_field.number_field_element.NumberFieldElement._invert_c_().

For the slowdown, I tried the following with increasing weights:

sage: G = DirichletGroup(17)
sage: %prun Newforms(G[2], 8, names='a')

Indeed, the method dual_eigenvector() takes up more and more of the time; more precisely, most of it is spent in a number of calls to the PARI function eltreltoabs():

      403   21.366    0.053   21.366    0.053 {method '_eltreltoabs' of 'sage.libs.pari.gen.gen' objects}
       17    1.843    0.108    2.290    0.135 {method 'echelon_form' of 'sage.matrix.matrix_cyclo_dense.Matrix_cyclo_dense' objects}
        1    1.706    1.706    1.963    1.963 {method 'nonpivots' of 'sage.matrix.matrix0.Matrix' objects}
        1    1.549    1.549    4.381    4.381 {method 'height' of 'sage.matrix.matrix_cyclo_dense.Matrix_cyclo_dense' objects}
        1    1.135    1.135   24.522   24.522 module.py:1076(dual_eigenvector)
     4056    0.603    0.000   14.743    0.004 matrix_space.py:1270(matrix)
72532/72531    0.477    0.000    0.928    0.000 number_field.py:9303(_element_constructor_)
65119/65111    0.457    0.000    1.202    0.000 polynomial_ring_constructor.py:50(PolynomialRing)
[...]

Perhaps #18727, #18740 and/or #252 could be relevant for this problem?

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented May 30, 2016

comment:10

@jdemeyer

I'm sorry, I guess I misinterpreted the "Author" field. I probably won't write code for this as I think the bugs and performance problems come from relative extensions of number fields and I don't know much about the code (and most of it is pari in some way, I guess).

My example in #20693 comment:2 was missing lines, sorry again.
I guess what I meant was

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: M=ModularSymbols(c,6,sign=1)
sage: S=M.cuspidal_subspace().new_subspace()
sage: Dec = S.decomposition()
sage: A=Dec[0]
sage: v =A.dual_eigenvector(names='a', lift=False)

To get a crash you have to let it run quite some time (I don't remember how long exactly it was, maybe 15 minutes, I can restart it and let you know).
I'm not sure if there much simpler/faster examples that crash but I can check.

I can come up with more examples for sure. Indeed you have to wait quite a bit until you see a crash sometimes, I had a computation running for Gamma0(19) and weight 14 for 2 days or so until I saw the crash (but I didn't run out of memory or anything else trivial and it also used to work on the same machine with the same space (and sage should really be able to do it!).

@nbruin I'll try to find some time to check out sage 6.1 again and compute some spaces that also run in sage 7.2 and then get back to you with some profiling results.

I did similar profiling tests as you did, @pjbruin and discovered that the main reason for this being slow is that elements in the base ring (a cyclotomic field) get coerced into a relative extension (obtained by adjoining the roots of a hecke polynomial). This is damn slow for some reason (note that it is stated in the documentation of relative number fields that doing arithmetic in number field towers is very slow but it seems that it got worse than it used to be). I discovered that such operations can be accelerated by very stupid means. For instance, instead of coercing an element e of K into an extension L directly, write it in terms of a power basis (e.list()), coerce the generator of K into L and create the element corresponding to e in L directly. I can document some examples later.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented May 30, 2016

Changed author from Stephan Ehlen to none

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented May 30, 2016

comment:12

@jdemeyer, OK, so the corrected example I gave in #20693 comment:10 crashes after at most 43 minutes for me. I can probably give examples that crash faster, but not sure. The runtime is most certainly related to the degree of the relative extension (the dimension of the modular symbols space) of the cyclotomic field and it seems that we need to have this large enough to really cause a crash.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 1, 2016

comment:13

Replying to @sagetrac-ehlen:

I did similar profiling tests as you did, @pjbruin and discovered that the main reason for this being slow is that elements in the base ring (a cyclotomic field) get coerced into a relative extension (obtained by adjoining the roots of a hecke polynomial). This is damn slow for some reason (note that it is stated in the documentation of relative number fields that doing arithmetic in number field towers is very slow but it seems that it got worse than it used to be). I discovered that such operations can be accelerated by very stupid means. For instance, instead of coercing an element e of K into an extension L directly, write it in terms of a power basis (e.list()), coerce the generator of K into L and create the element corresponding to e in L directly. I can document some examples later.

It turns out that NumberField_relative.__base_inclusion() uses a PARI function that is too general (and hence too slow) for this purpose. I opened #20749 to address this.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 1, 2016

comment:14

Replying to @pjbruin:

For the slowdown, I tried the following with increasing weights:

sage: G = DirichletGroup(17)
sage: %prun Newforms(G[2], 8, names='a')

Indeed, the method dual_eigenvector() takes up more and more of the time; more precisely, most of it is spent in a number of calls to the PARI function eltreltoabs():

After #20749, dual_eigenvector() only takes up 28% of the time; before, it was 69%.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 1, 2016

comment:15

My old examples still crash with #20749 applied.
Unfortunately, it does not seem like these were very easy to compute with sage 6.1; everything I tested so far runs in fact much slower in sage 6.1 (which is good, I guess). [Let me mention, although a bit off-topic, that magma computes the newforms in weight 4 in the example above in 0.5s and in weight 6 in 7.5s.]

Now, as for the crash, first of all there is no crash with #20749 when running

N=Newforms(DirichletGroup(23).gen()^2,6, names='a')

However, I now get the more informative

PariError: the PARI stack overflows (current size: 2147483648; maximum size: 2147483648)
You can use pari.allocatemem() to change the stack size and try again.

After doing so:

sage: pari.allocatemem(2147483648*4)

I get a crash again. For some reason this crash is not informative at all on OSX, I will have to run it on Linux to see the backtrace, I guess.

Also, your level 17 example is quite good because by increasing the weight only by one to 9, I the same behaviour (PARI stack overflow, crash when increasing it). Try:

sage: c=DirichletGroup(17).gen()
sage: %time Newforms(c,9, names='a')

@pjbruin
Copy link
Contributor

pjbruin commented Jun 2, 2016

comment:16

From the traceback it appears that NTL runs out of primes for its FFT-based ZZX_XGCD function. We should catch the NTL error in NumberFieldElement._invert_c_() and use an alternative extended GCD implementation if that happens.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 3, 2016

comment:17

Replying to @pjbruin:

From the traceback it appears that NTL runs out of primes for its FFT-based ZZX_XGCD function. We should catch the NTL error in NumberFieldElement._invert_c_() and use an alternative extended GCD implementation if that happens.

Indeed, after finally managing to get sage run in gdb on OSX, I get

#4  0x00000001076cc353 in NTL::UseFFTPrime(long) () from /Users/stephan/sage/local/lib/libntl.25.dylib
#5  0x00000001077bfaf8 in NTL::zz_pContext::zz_pContext(NTL::INIT_FFT_STRUCT const&, long) ()
   from /Users/stephan/sage/local/lib/libntl.25.dylib
#6  0x00000001077bfd07 in NTL::zz_p::FFTInit(long) () from /Users/stephan/sage/local/lib/libntl.25.dylib
#7  0x0000000107743dba in NTL::resultant(NTL::ZZ&, NTL::ZZX const&, NTL::ZZX const&, long) ()
   from /Users/stephan/sage/local/lib/libntl.25.dylib
#8  0x0000000107749594 in NTL::XGCD(NTL::ZZ&, NTL::ZZX&, NTL::ZZX&, NTL::ZZX const&, NTL::ZZX const&, long) ()
   from /Users/stephan/sage/local/lib/libntl.25.dylib
#9  0x000000018b253179 in __pyx_f_4sage_5rings_12number_field_20number_field_element_18NumberFieldElement__invert_c_(__pyx_obj_4sage_5rings_12number_field_20number_field_element_NumberFieldElement*, NTL::ZZX*, NTL::ZZ*) ()

Is this a bug in NTL? Is it documented anywhere that this might fail? Can we know in advance that this will happen so that we don't have to run into _invert_c_c()?

@pjbruin
Copy link
Contributor

pjbruin commented Jun 9, 2016

comment:18

After applying #20749, #20759 and #20791, the examples in the ticket description and in comment:10 crash already after about 1-2 minutes; on this ticket we can now focus on the crash instead of the slowness.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 9, 2016

comment:19

Replying to @sagetrac-ehlen:

I would like to add a remark: something really extremely bad happened with the implementation of relative Number Fields.
It was always slow but now, and this might be related to this crash here, spaces that I computed using older versions of sage (6.1 to be concrete and give a pointer) within minutes or up to an hour or a bit more, do not finish to be computed within a day now. A concrete example is Gamma0(19), character [zeta18^2], weight 14. I had the modular symbols space computed with sage 6.1 and stored the cputime used to compute it: 3780s. And now it didn't finish within 7 hours.

With the same branches applied as in comment:18, I tried

%prun Newforms(DirichletGroup(19)[2], 14, names='a')

and got a similar NTL crash as in comment:17 after about 1.5 hours.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 10, 2016

comment:20

I tried to debug this and catch the NTL exception but I failed.
Can someone help me? How do I properly catch the exception that NTL is probably raising?
I tried adding except + to the declarations of ZZX_XGCD in ZZX.pxd and then subsequently to _invert_c and wrapped the call to _invert_c in try:...:except: but that didn't work... I'm probably doing something stupid.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 15, 2016

comment:21

Replying to @sagetrac-ehlen:

I tried to debug this and catch the NTL exception but I failed.
Can someone help me? How do I properly catch the exception that NTL is probably raising?
I tried adding except + to the declarations of ZZX_XGCD in ZZX.pxd and then subsequently to _invert_c and wrapped the call to _invert_c in try:...:except: but that didn't work... I'm probably doing something stupid.

The method _invert_c() should be declared except * instead of except +, and you have to wrap the NTL calls in sig_on()...sig_off(). I made the following changes (not sure if the except + on ZZX_XGCD is actually necessary):

--- a/src/sage/libs/ntl/ZZX.pxd
+++ b/src/sage/libs/ntl/ZZX.pxd
@@ -35,7 +35,7 @@ cdef extern from "sage/libs/ntl/ntlwrap.cpp":
     void ZZX_div_ZZ "div"( ZZX_c x, ZZX_c a, ZZ_c b)
     long ZZX_deg "deg"( ZZX_c x )
     void ZZX_rem "rem"(ZZX_c r, ZZX_c a, ZZX_c b)
-    void ZZX_XGCD "XGCD"(ZZ_c r, ZZX_c s, ZZX_c t, ZZX_c a, ZZX_c b, long deterministic)
+    void ZZX_XGCD "XGCD"(ZZ_c r, ZZX_c s, ZZX_c t, ZZX_c a, ZZX_c b, long deterministic) except +
     void ZZX_content "content"(ZZ_c d, ZZX_c f)
     void ZZX_factor "factor"(ZZ_c c, vec_pair_ZZX_long_c factors, ZZX_c f, long verbose, long bnd)
 
--- a/src/sage/rings/number_field/number_field_element.pxd
+++ b/src/sage/rings/number_field/number_field_element.pxd
@@ -26,7 +26,7 @@ cdef class NumberFieldElement(FieldElement):
     cdef void _ntl_coeff_as_mpz(self, mpz_t z, long i)
     cdef void _ntl_denom_as_mpz(self, mpz_t z)
 
-    cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den)
+    cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) except *
     cdef void _reduce_c_(self)
     cpdef ModuleElement _add_(self, ModuleElement right)
     cpdef ModuleElement _sub_(self, ModuleElement right)
--- a/src/sage/rings/number_field/number_field_element.pyx
+++ b/src/sage/rings/number_field/number_field_element.pyx
@@ -2258,7 +2258,7 @@ cdef class NumberFieldElement(FieldElement):
         """
         return long(self.polynomial())
 
-    cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den):
+    cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) except *:
         """
         Computes the numerator and denominator of the multiplicative
         inverse of this element.
@@ -2276,11 +2276,13 @@ cdef class NumberFieldElement(FieldElement):
         """
         cdef ZZX_c t # unneeded except to be there
         cdef ZZX_c a, b
+        sig_on()
         ZZX_mul_ZZ( a, self.__numerator, self.__fld_denominator.x )
         ZZX_mul_ZZ( b, self.__fld_numerator.x, self.__denominator )
         ZZX_XGCD( den[0], num[0],  t, a, b, 1 )
         ZZX_mul_ZZ( num[0], num[0], self.__fld_denominator.x )
         ZZX_mul_ZZ( num[0], num[0], self.__denominator )
+        sig_off()
 
     def __invert__(self):
         """

Then I get

sage: D=DirichletGroup(23)
sage: c=D.gen()^2
sage: N=Newforms(c,6, names='a')
---------------------------------------------------------------------------
NTLError                                  Traceback (most recent call last)
<ipython-input-3-e9e314cc3f45> in <module>()
----> 1 N=Newforms(c,Integer(6), names='a')

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/modform/constructor.pyc in Newforms(group, weight, base_ring, names)
    452 
    453     """
--> 454     return CuspForms(group, weight, base_ring).newforms(names)
    455 
    456 

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/modform/space.pyc in newforms(self, names)
   1680             names = 'a'
   1681         return [ element.Newform(self, factors[i], names=(names+str(i)) )
-> 1682                  for i in range(len(factors)) ]
   1683 
   1684     def eisenstein_submodule(self):

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/modform/element.pyc in __init__(self, parent, component, names, check)
   1070             if not component.is_simple():
   1071                 raise ValueError("component must be simple")
-> 1072         extension_field = component.eigenvalue(1,name=names).parent()
   1073         if extension_field != parent.base_ring(): # .degree() != 1 and rings.is_NumberField(extension_field):
   1074             assert extension_field.base_field() == parent.base_ring()

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/hecke/module.pyc in eigenvalue(self, n, name)
   1304         if (arith.is_prime(n) or n==1):
   1305             Tn_e = self._eigen_nonzero_element(n)
-> 1306             an = self._element_eigenvalue(Tn_e, name=name)
   1307             _dict_set(ev, n, name, an)
   1308             return an

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/hecke/module.pyc in _element_eigenvalue(self, x, name)
    694         if not x in self.ambient_hecke_module():
    695             raise ArithmeticError("x must be in the ambient Hecke module.")
--> 696         v = self.dual_eigenvector(names=name)
    697         return v.dot_product(x.element())
    698 

/home/bruinpj/src/sage/local/lib/python2.7/site-packages/sage/modular/hecke/module.pyc in dual_eigenvector(self, names, lift, nz)
   1197             x = self._eigen_nonzero_element()
   1198         alpha = w_lift.dot_product(x.element())
-> 1199         beta = ~alpha
   1200         w_lift = w_lift * beta
   1201         w = w * beta

/home/bruinpj/src/sage/src/sage/rings/number_field/number_field_element.pyx in sage.rings.number_field.number_field_element.NumberFieldElement.__invert__ (/home/bruinpj/src/sage/src/build/cythonized/sage/rings/number_field/number_field_element.cpp:22816)()
   2301         cdef NumberFieldElement x
   2302         x = self._new()
-> 2303         self._invert_c_(&x.__numerator, &x.__denominator)
   2304         x._reduce_c_()
   2305         return x

/home/bruinpj/src/sage/src/sage/rings/number_field/number_field_element.pyx in sage.rings.number_field.number_field_element.NumberFieldElement._invert_c_ (/home/bruinpj/src/sage/src/build/cythonized/sage/rings/number_field/number_field_element.cpp:22657)()
   2277         cdef ZZX_c t # unneeded except to be there
   2278         cdef ZZX_c a, b
-> 2279         sig_on()
   2280         ZZX_mul_ZZ( a, self.__numerator, self.__fld_denominator.x )
   2281         ZZX_mul_ZZ( b, self.__fld_numerator.x, self.__denominator )

/home/bruinpj/src/sage/src/sage/libs/ntl/error.pyx in sage.libs.ntl.error.NTL_error_callback (/home/bruinpj/src/sage/src/build/cythonized/sage/libs/ntl/error.cpp:794)()
     39 
     40 cdef void NTL_error_callback(const char* s) except *:
---> 41     raise NTLError(s)
     42 
     43 

NTLError: FFT prime index too large

@jdemeyer
Copy link

comment:22

Instead of cdef void ... except *, it's more efficient to use cdef int ... except -1. It behaves the same in practice but it doesn't need a call to PyErr_Occurred().

@pjbruin
Copy link
Contributor

pjbruin commented Jun 15, 2016

comment:23

Another solution, using PARI as a fallback:

--- a/src/sage/rings/number_field/number_field_element.pyx
+++ b/src/sage/rings/number_field/number_field_element.pyx
@@ -40,6 +40,8 @@ from sage.libs.gmp.mpz cimport *
 from sage.libs.gmp.mpq cimport *
 from sage.libs.mpfi cimport mpfi_t, mpfi_init, mpfi_set, mpfi_clear, mpfi_div_z, mpfi_init2, mpfi_get_prec, mpfi_set_prec
 from sage.libs.mpfr cimport mpfr_less_p, mpfr_greater_p, mpfr_greaterequal_p
+from sage.libs.ntl.error import NTLError
+
 from cpython.object cimport Py_EQ, Py_NE, Py_LT, Py_GT, Py_LE, Py_GE
 from sage.structure.sage_object cimport rich_to_bool
 
@@ -2297,9 +2299,14 @@ cdef class NumberFieldElement(FieldElement):
         if IsZero_ZZX(self.__numerator):
             raise ZeroDivisionError
         cdef NumberFieldElement x
-        x = self._new()
-        self._invert_c_(&x.__numerator, &x.__denominator)
-        x._reduce_c_()
+        try:
+            x = self._new()
+            sig_on()
+            self._invert_c_(&x.__numerator, &x.__denominator)
+            x._reduce_c_()
+            sig_off()
+        except NTLError:
+            x = self._parent(~self._pari_())
         return x
 
     def _integer_(self, Z=None):

It no longer crashes and returns an answer after about 20 minutes...

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 16, 2016

comment:24

@pjbruin: Great! Thanks a lot for your work!

However, the pari alternative seems to be kind of slow (which is much better than crashing, of course).
With the example I gave in #20749, the pari inversion takes about 7 minutes but the code below, which mimics the _invert_c function, runs in about 1.5 minutes for me:

d = alpha.polynomial().denominator()
D = alpha.parent().absolute_polynomial().denominator()
r,s,tt = xgcd(alpha.polynomial().numerator()*D,alpha.parent().absolute_polynomial().numerator()*d)
b = s*d/D
F = alpha.parent()
beta = F(F.polynomial_quotient_ring()(b))

I don't like the last line of the code and I still have to make sure I understood everything correctly and it works in all cases but then I could submit a patch for the number field element to use a variant of this code as an alternative. What do you think?

@pjbruin
Copy link
Contributor

pjbruin commented Jun 16, 2016

comment:25

Replying to @sagetrac-ehlen:

However, the pari alternative seems to be kind of slow (which is much better than crashing, of course).

I would have hoped for PARI to be quite fast, but it seems this is not particularly optimised in PARI. In any case I didn't try hard to make it fast; it was just the easiest solution that didn't crash.

With the example I gave in #20749, the pari inversion takes about 7 minutes but the code below, which mimics the _invert_c function, runs in about 1.5 minutes for me:

d = alpha.polynomial().denominator()
D = alpha.parent().absolute_polynomial().denominator()
r,s,tt = xgcd(alpha.polynomial().numerator()*D,alpha.parent().absolute_polynomial().numerator()*d)
b = s*d/D
F = alpha.parent()
beta = F(F.polynomial_quotient_ring()(b))

I don't like the last line of the code and I still have to make sure I understood everything correctly and it works in all cases but then I could submit a patch for the number field element to use a variant of this code as an alternative. What do you think?

This is definitely a good start. The xgcd is implemented using FLINT, which suggests a slightly more direct approach: use the functions from sage.libs.flint.ntl_interface to convert alpha.__numerator and alpha.__denominator to FLINT fmpz_poly and fmpz objects, respectively (and similarly for the defining polynomial of the number field), then call the FLINT function fmpz_poly_xgcd, and then convert back to NTL. It could even be reasonable to simply always use FLINT instead of NTL in __invert__.

@pjbruin pjbruin changed the title Sage crashes when computing newforms Sage crashes when inverting/dividing large number field elements Jun 20, 2016
@jdemeyer
Copy link

comment:39

Replying to @pjbruin:

Changing the summary and component of this ticket to reflect the fact that this is purely a bug in the number field code.

Please also update the ticket description.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

comment:40

Replying to @jdemeyer:

Replying to @pjbruin:

Changing the summary and component of this ticket to reflect the fact that this is purely a bug in the number field code.

Please also update the ticket description.

This appears to be impossible at the moment due to Trac server problems...

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

Changed author from Stephan Ehlen to Stephan Ehlen, Peter Bruin

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

Changed commit from a26dd5f to 1620359

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

comment:41

Pushing my commit and changing the branch does seem to work.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 20, 2016

comment:42

I did some cleaning up of _div_() and __invert__(). There seems to be no reason anymore to have a separate _invert_c_() method; I removed this. Note that _div_() is just __invert__() followed by _mul_(), with an intermediate reduction removed.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 20, 2016

comment:43

@pjbruin Very good, I didn't look at _div_(), which was stupid ;-) Anyway, do we really want to have the code for inverting and multiplying duplicated in _div_()? I guess this is for performance reasons?
I.e. would there be anything bad about just doing

return self._mul_(right.__invert__())

in _div_()?

@sagetrac-ehlen

This comment has been minimized.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 20, 2016

comment:45

@pjbruin I can remove the doctest.
What I really miss in sage is some more unit tests that would check many many cases and maybe even some random cases to work that should be separate from the doctests but this should be discussed elsewhere.

If you agree that simplifying _div_() makes sense, I could push both changes at any time.
(I made some randomized tests and I don't see a performance advantage of the currently duplicated code versus the suggested simplification.)

Also, did you check very carefully that removing those additional multiplications that used to be done in _invert_c() do not cause any trouble anywhere? I wonder how someone could even come up with this.... could it really be that there was really no reason at all for this??? (This is exactly why I think there should be (even) more (complex) testing in sage.)
Maybe we should ask the author Joel B. Mohler? However, I don't see activity by him on the github repo and the code has been sitting there since 2009... It seems really trivial to remove it and I can't find any problems but it just seems so hard for me to believe that someone did such a complicated thing without any reason ;-)

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 20, 2016

comment:46

PS: my branch u/ehlen/sage_crashes_when_computing_newforms includes these changes now (I didn't change the branch of the ticket since I wanted to wait for your answer because you might have had a good reason for not simplifying div as much as I did).

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Jun 21, 2016

Branch pushed to git repo; I updated commit sha1. New commits:

d089903Removed doctest that takes long and simplified division of number field elements a lot.
eb3da68Trac 20693: reviewer patch

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Jun 21, 2016

Changed commit from 1620359 to eb3da68

@pjbruin
Copy link
Contributor

pjbruin commented Jun 21, 2016

Changed author from Stephan Ehlen, Peter Bruin to Stephan Ehlen

@pjbruin
Copy link
Contributor

pjbruin commented Jun 21, 2016

comment:48

Replying to @sagetrac-ehlen:

@pjbruin I can remove the doctest.
What I really miss in sage is some more unit tests that would check many many cases and maybe even some random cases to work that should be separate from the doctests but this should be discussed elsewhere.

Yes, it might make sense to have a place for tests that are too long to write down or take too much time to run. I agree that this is not the place to discuss this further.

If you agree that simplifying _div_() makes sense, I could push both changes at any time.
(I made some randomized tests and I don't see a performance advantage of the currently duplicated code versus the suggested simplification.)

I included your commit in my branch and added a reviewer commit (typographical fixes and using x.__invert__ by the shorter and faster ~x, which does not require a method name lookup). I did not try any timings myself, but can imagine that this simplification makes no practical difference (the running time is probably dominated by the extended GCD computation).

Also, did you check very carefully that removing those additional multiplications that used to be done in _invert_c() do not cause any trouble anywhere?

Here is a proof that the simplified code is correct. Suppose our number field K, is Q[X]/(f) where we may assume without loss of generality that f is in Z[X]. We want to invert x = (g mod f)/d in K, where g is in Z[X] and is coprime to f, and where d is in Z and non-zero. We call XGCD(e, h, t, g, f); then e is the resultant of g and f (a non-zero integer) and h, t are in Z[X] satisfying e = hg + tf. Modulo f, this becomes (h mod f)(g mod f) = e in K, and hence 1/x = d/(g mod f) = ((h mod f)d)/e. The right-hand is precisely what is computed by the simplified version of the code. I assume that the previous version resulted in an equivalent, but less often reduced, fraction.

I wonder how someone could even come up with this.... could it really be that there was really no reason at all for this??? (This is exactly why I think there should be (even) more (complex) testing in sage.)

Yes, this could very well be; it is surprising how often one runs into code that is much more complicated than necessary...

It now seems to makes more sense to me to be the reviewer instead of an author. I am now running doctests. Maybe it is still good if someone else takes a look at this; these number field operations are after all essential and heavily used...

@pjbruin
Copy link
Contributor

pjbruin commented Jun 21, 2016

Reviewer: Peter Bruin

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 21, 2016

comment:49

@pjbruin: I completely agree with the proof and this is essentially what I tried to write in comment:30.
I was only worried I might overlook a detail in the implementation and that my assumptions (for instance the representation of x to begin with) are not entirely correct in every case.

And yes, it would be good if someone else could have a look as well although the code is pretty solid now, I think.
But the impact of changing number field code is indeed so great that it can't hurt to double-check.

@pjbruin
Copy link
Contributor

pjbruin commented Jun 21, 2016

comment:50

Replying to @sagetrac-ehlen:

@pjbruin: I completely agree with the proof and this is essentially what I tried to write in comment:30.
I was only worried I might overlook a detail in the implementation and that my assumptions (for instance the representation of x to begin with) are not entirely correct in every case.

OK, good to see that we agree on this. I don't think there are any hidden assumptions that could be violated...

And yes, it would be good if someone else could have a look as well although the code is pretty solid now, I think.
But the impact of changing number field code is indeed so great that it can't hurt to double-check.

I agree. By the way, doctests still pass.

@sagetrac-ehlen
Copy link
Mannequin Author

sagetrac-ehlen mannequin commented Jun 21, 2016

comment:51

One thing I don't yet understand (and maybe we don't have to) is why it was possible to compute some of the spaces of newforms that caused NTL to run out of FFT primes in older versions of sage. One particular example that I have is level 17, weight 9, character=generator of the dirichlet group. I have a file lying around that says that this was computed with sage 5.12beta4, release date: 2013-08-30. Why did inversion work back then if the cython version number_field_element.pyx exists since 2007 and the inversion has not been changed since then? Somehow I can only imagine that NTL used to be able to handle these large polynomials that the code used to produce. Any other ideas?

@fredstro
Copy link

comment:52

I have gone through the proposed patch and everything seems to work as stated.

  • I ran sage -t --all and the few doctests which failed seems to not be related to the proposed patch (and in particular all tests in sage -t src/sage/rings/number_field/ passes).

  • The example given above for level 23 crash on Sage 7.2 without the patch and works fine with the patch.

  • The new code introduced in src/sage/rings/number_field/number_field_element.pyx is well-formatted and does what is intended.

  • The code which was removed from src/sage/rings/number_field/number_field_element.pyx (i.e. the 'invert_c' routine) does indeed seem to have been unnecessarily complicated.

  • The new code introduce a test which crashes in Sage 7.2 without the patch but passes with the patch.

@fredstro
Copy link

Changed reviewer from Peter Bruin to Peter Bruin, Fredrik Stromberg

@vbraun
Copy link
Member

vbraun commented Jun 30, 2016

Changed branch from u/pbruin/20693-sage_crashes_when_computing_newforms to eb3da68

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

5 participants