Opened 7 years ago
Closed 6 years ago
#20693 closed defect (fixed)
Sage crashes when inverting/dividing large number field elements
Reported by:  Stephan Ehlen  Owned by:  

Priority:  critical  Milestone:  sage7.3 
Component:  number fields  Keywords:  
Cc:  Merged in:  
Authors:  Stephan Ehlen  Reviewers:  Peter Bruin, Fredrik Stromberg 
Report Upstream:  N/A  Work issues:  
Branch:  eb3da68 (Commits, GitHub, GitLab)  Commit:  eb3da684a613c81294dec6eff22972152950790f 
Dependencies:  Stopgaps: 
Description (last modified by )
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.
Change History (53)
comment:1 Changed 7 years ago by
Description:  modified (diff) 

comment:2 followup: 7 Changed 7 years ago by
comment:3 Changed 7 years ago by
Priority:  major → critical 

comment:4 followup: 19 Changed 7 years ago by
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
comment:5 Changed 7 years ago by
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.
comment:6 Changed 7 years ago by
@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.
comment:7 Changed 7 years ago by
Replying to 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...
comment:8 Changed 7 years ago by
Replying to 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?
comment:9 followup: 14 Changed 7 years ago by
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?
comment:10 followup: 13 Changed 7 years ago by
@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 http://trac.sagemath.org/ticket/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, @pbruin 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.
comment:11 Changed 7 years ago by
Authors:  Stephan Ehlen 

comment:12 Changed 7 years ago by
@jdemeyer, OK, so the corrected example I gave in http://trac.sagemath.org/ticket/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.
comment:13 Changed 7 years ago by
Replying to ehlen:
I did similar profiling tests as you did, @pbruin 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.
comment:14 Changed 7 years ago by
Replying to pbruin:
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 functioneltreltoabs()
:
After #20749, dual_eigenvector()
only takes up 28% of the time; before, it was 69%.
comment:15 Changed 7 years ago by
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 offtopic, 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')
comment:16 followup: 17 Changed 7 years ago by
From the traceback it appears that NTL runs out of primes for its FFTbased ZZX_XGCD
function. We should catch the NTL error in NumberFieldElement._invert_c_()
and use an alternative extended GCD implementation if that happens.
comment:17 Changed 7 years ago by
Replying to pbruin:
From the traceback it appears that NTL runs out of primes for its FFTbased
ZZX_XGCD
function. We should catch the NTL error inNumberFieldElement._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()?
comment:18 Changed 7 years ago by
After applying #20749, #20759 and #20791, the examples in the ticket description and in comment:10 crash already after about 12 minutes; on this ticket we can now focus on the crash instead of the slowness.
comment:19 Changed 7 years ago by
Replying to 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.
comment:20 followup: 21 Changed 7 years ago by
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.
comment:21 Changed 6 years ago by
Replying to 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):

src/sage/libs/ntl/ZZX.pxd
a b cdef extern from "sage/libs/ntl/ntlwrap.cpp": 35 35 void ZZX_div_ZZ "div"( ZZX_c x, ZZX_c a, ZZ_c b) 36 36 long ZZX_deg "deg"( ZZX_c x ) 37 37 void ZZX_rem "rem"(ZZX_c r, ZZX_c a, ZZX_c b) 38 void ZZX_XGCD "XGCD"(ZZ_c r, ZZX_c s, ZZX_c t, ZZX_c a, ZZX_c b, long deterministic) 38 void ZZX_XGCD "XGCD"(ZZ_c r, ZZX_c s, ZZX_c t, ZZX_c a, ZZX_c b, long deterministic) except + 39 39 void ZZX_content "content"(ZZ_c d, ZZX_c f) 40 40 void ZZX_factor "factor"(ZZ_c c, vec_pair_ZZX_long_c factors, ZZX_c f, long verbose, long bnd) 41 41 
src/sage/rings/number_field/number_field_element.pxd
a b cdef class NumberFieldElement(FieldElement): 26 26 cdef void _ntl_coeff_as_mpz(self, mpz_t z, long i) 27 27 cdef void _ntl_denom_as_mpz(self, mpz_t z) 28 28 29 cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) 29 cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) except * 30 30 cdef void _reduce_c_(self) 31 31 cpdef ModuleElement _add_(self, ModuleElement right) 32 32 cpdef ModuleElement _sub_(self, ModuleElement right) 
src/sage/rings/number_field/number_field_element.pyx
a b cdef class NumberFieldElement(FieldElement): 2258 2258 """ 2259 2259 return long(self.polynomial()) 2260 2260 2261 cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) :2261 cdef void _invert_c_(self, ZZX_c *num, ZZ_c *den) except *: 2262 2262 """ 2263 2263 Computes the numerator and denominator of the multiplicative 2264 2264 inverse of this element. … … cdef class NumberFieldElement(FieldElement): 2276 2276 """ 2277 2277 cdef ZZX_c t # unneeded except to be there 2278 2278 cdef ZZX_c a, b 2279 sig_on() 2279 2280 ZZX_mul_ZZ( a, self.__numerator, self.__fld_denominator.x ) 2280 2281 ZZX_mul_ZZ( b, self.__fld_numerator.x, self.__denominator ) 2281 2282 ZZX_XGCD( den[0], num[0], t, a, b, 1 ) 2282 2283 ZZX_mul_ZZ( num[0], num[0], self.__fld_denominator.x ) 2283 2284 ZZX_mul_ZZ( num[0], num[0], self.__denominator ) 2285 sig_off() 2284 2286 2285 2287 def __invert__(self): 2286 2288 """
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) <ipythoninput3e9e314cc3f45> in <module>() > 1 N=Newforms(c,Integer(6), names='a') /home/bruinpj/src/sage/local/lib/python2.7/sitepackages/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/sitepackages/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/sitepackages/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/sitepackages/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/sitepackages/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/sitepackages/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
comment:22 Changed 6 years ago by
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()
.
comment:23 Changed 6 years ago by
Another solution, using PARI as a fallback:

src/sage/rings/number_field/number_field_element.pyx
a b from sage.libs.gmp.mpz cimport * 40 40 from sage.libs.gmp.mpq cimport * 41 41 from sage.libs.mpfi cimport mpfi_t, mpfi_init, mpfi_set, mpfi_clear, mpfi_div_z, mpfi_init2, mpfi_get_prec, mpfi_set_prec 42 42 from sage.libs.mpfr cimport mpfr_less_p, mpfr_greater_p, mpfr_greaterequal_p 43 from sage.libs.ntl.error import NTLError 44 43 45 from cpython.object cimport Py_EQ, Py_NE, Py_LT, Py_GT, Py_LE, Py_GE 44 46 from sage.structure.sage_object cimport rich_to_bool 45 47 … … cdef class NumberFieldElement(FieldElement): 2297 2299 if IsZero_ZZX(self.__numerator): 2298 2300 raise ZeroDivisionError 2299 2301 cdef NumberFieldElement x 2300 x = self._new() 2301 self._invert_c_(&x.__numerator, &x.__denominator) 2302 x._reduce_c_() 2302 try: 2303 x = self._new() 2304 sig_on() 2305 self._invert_c_(&x.__numerator, &x.__denominator) 2306 x._reduce_c_() 2307 sig_off() 2308 except NTLError: 2309 x = self._parent(~self._pari_()) 2303 2310 return x 2304 2311 2305 2312 def _integer_(self, Z=None):
It no longer crashes and returns an answer after about 20 minutes...
comment:24 followup: 25 Changed 6 years ago by
@pbruin: 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?
comment:25 Changed 6 years ago by
Replying to 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__
.
comment:27 Changed 6 years ago by
Sorry for the delay.
It's a bit weird but I really must have made a mistake when I timed the implementation in comment:24, although I really don't know what went wrong. Anyway, I cannot reproduce the results and instead all my tests indicate that pari inversion is about as fast as xgcd and also as the following more direct approach using flint which can be used in number_field_element.pyx
cdef void _invert_flint_(self, ZZX_c *num, ZZ_c *den): cdef ZZX_c a, b ZZX_mul_ZZ( a, self.__numerator, self.__fld_denominator.x ) ZZX_mul_ZZ( b, self.__fld_numerator.x, self.__denominator ) cdef fmpz_poly_t a_f, b_f, t_f cdef fmpz_poly_t num_f #flint numerator cdef fmpz_t den_f #flint denominator fmpz_poly_init(a_f); fmpz_poly_init(b_f); fmpz_poly_init(t_f); fmpz_poly_init(num_f); fmpz_init(den_f) fmpz_poly_set_ZZX(a_f, a) #convert to flint fmpz_poly_set_ZZX(b_f, b) #convert to flint fmpz_poly_xgcd_modular(den_f, num_f, t_f, a_f, b_f) fmpz_poly_get_ZZX(num[0], num_f) #convert back to NTL fmpz_get_ZZ(den[0], den_f) #convert back to NTL #finishing up ZZX_mul_ZZ(num[0], num[0], self.__fld_denominator.x ) ZZX_mul_ZZ(num[0], num[0], self.__denominator )
So, I'm not sure if it's worth it to include the flint alternative. But I found something else that even avoids the NTLError to happen and makes the newforms computations way faster in many cases! I'm cleaning up my branch now after a lot of testing and will then push the changes, including your pari alternative.
comment:28 Changed 6 years ago by
Branch:  → u/ehlen/sage_crashes_when_computing_newforms 

comment:29 Changed 6 years ago by
Commit:  → 9218fb3214dbc776e1f914bd68f8fb4df1dd437e 

Branch pushed to git repo; I updated commit sha1. New commits:
9218fb3  removed unnecessary import

comment:30 Changed 6 years ago by
OK, so I really tried a lot of things and I came to the conclusion that what is currently done in _invert_c_
of number_field_element.pyx
is too much work, "causes" the NTLError in the cases at hand and inversion can be made much much faster with a minimal change.
The main point is this: Currently, given a number field element of the form x/d, where (if I understood it correctly), x is a polynomial with integer coeffs and d is an integer, modulo a polynomial M/D, where again M is a poly with integer coeffs and M an integer, the code sets
a = x*D b = M*d
and computes r, s, t
, such that
r = s*a + t*b.
Thus, s/r
is the inverse of a=x*D
modulo b=M*d
(with rational coefficients, in Q[X]/(dM)).
Now, I don't really now why this is done and maybe I misunderstand something but I think taking dM
here is unnecessary and slows everything down. Anyway, the result is that
1 = s/r*D*d*x/d + t*D/(r*d)*d*M/D
And so the code returns s/r*D*d
as the inverse.
What I did was to instead set
a = x*D b = M
and then again compute r,s,t
as above so that now s/r
is the inverse of a
modulo M
with rational coefficients, i.e.
in Q[X]/(M). This means that we get
1 = s/r*D*d*x/d + t*D/r*M/D
which shows that modulo M/D, we get that s/r*D*d
is the inverse.
Now, after I wrote this down, it also seems that multiplying x
by D
is not needed.
I will try this out and push the changes immediately.
Btw, all doctests pass.
I also added some doctests to make sure we always test that this issue here is resolved. However, since the example I had at hand is so large, it adds a lot of junk to the source code which I don't like very much. Is there any way to add an external test, that's not necessary in the dotstring in sage (no nosetests at all???)?
comment:31 Changed 6 years ago by
These changes, together with what was done before have an enormous impact on computing newforms! Example, in the current branch:
sage: %time N=Newforms(DirichletGroup(17).gen(), 7, names='a') CPU times: user 34.1 s, sys: 372 ms, total: 34.5 s Wall time: 34.4 s
use to run for over 6 minutes in sage 7.1.
Moreover, as a remark: In fact, I tested a few cases against magma and it seems we are now able to beat magma even though sage provides somewhat more information since we have the relative extensions over the cyclotomic fields and magma provides absolute number fields. To compute the same space as above in magma takes only less than 2s but to get to to Fourier coefficients takes additional computing time. In sage I get 100 coefficients now in an additional 9s. In magma it takes 44s (although on a different machine which might very well be a bit slower). Anyway, subsequent calls to get more Fourier coefficients scale now very well in sage. To get 200 coefficients with 100 precomputed takes just 9s more but in magma it takes me 90s! So the investment made to compute the space in the beginning with sage pays certainly off if you need more Fourier coefficients.
comment:32 Changed 6 years ago by
Authors:  → Stephan Ehlen 

comment:33 Changed 6 years ago by
I made the change mentioned above but currently I'm unable to push:
STDERR: ssh: connect to host trac.sagemath.org port 22: No route to host
I will try again later. All tests pass (and the computation mentioned above is one second faster now but this might well be within the usual variation).
comment:34 Changed 6 years ago by
Commit:  9218fb3214dbc776e1f914bd68f8fb4df1dd437e → e4ba18d98841d267f0fc9911e319f5f6eebaa91f 

Branch pushed to git repo; I updated commit sha1. New commits:
e4ba18d  remove further unneded operations when inverting

comment:35 Changed 6 years ago by
OK, so pushing worked now again. All tests pass and I did some randomized test for a couple of number fields which all worked (trying to make sure I didn't miss anything with the inversion, but maybe I did and don't see it in the test I do, so it would be very good if someone could check this thoroughly). Anyway, in these tests, I generated in a for loop 1000 times 2 elements r,s in a number field, and test if indeed r/s == ~(s/r) and ~r is the inverse of r and of course ~s the inverse of s... In sage 7.1 the routine usually (using timeit) runs 1.5s and in this branch about a second.
comment:36 Changed 6 years ago by
Status:  new → needs_review 

comment:37 Changed 6 years ago by
Commit:  e4ba18d98841d267f0fc9911e319f5f6eebaa91f → a26dd5fd1a0e12c05e3b02f0715a11705f21e04c 

Branch pushed to git repo; I updated commit sha1. New commits:
a26dd5f  Simpler doctest that crashes in sage 7.1 but works now.

comment:38 followup: 39 Changed 6 years ago by
Component:  modular forms → number fields 

Summary:  Sage crashes when computing newforms → Sage crashes when inverting/dividing large number field elements 
Changing the summary and component of this ticket to reflect the fact that this is purely a bug in the number field code.
Is there a doctest to check that Newforms
now works correctly that takes less time? (It takes about 5 minutes on the system that I tested this on.) Alternatively, would you agree with just removing this doctest (given that the bug is in NumberFieldElement
)?
The current branch just fixes inversion; division is still susceptible to the same bug. I will upload a commit that fixes _div_()
in a similar way, and gets rid of the _invert_c_()
method.
comment:39 followup: 40 Changed 6 years ago by
Replying to pbruin:
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.
comment:40 Changed 6 years ago by
comment:41 Changed 6 years ago by
Authors:  Stephan Ehlen → Stephan Ehlen, Peter Bruin 

Branch:  u/ehlen/sage_crashes_when_computing_newforms → u/pbruin/20693sage_crashes_when_computing_newforms 
Commit:  a26dd5fd1a0e12c05e3b02f0715a11705f21e04c → 162035995bff5bb15d73bc13d2994655f24a8c23 
Pushing my commit and changing the branch does seem to work.
comment:42 Changed 6 years ago by
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.
comment:43 Changed 6 years ago by
@pbruin 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_()
?
comment:44 Changed 6 years ago by
Description:  modified (diff) 

comment:45 followup: 48 Changed 6 years ago by
@pbruin 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 ;)
comment:46 Changed 6 years ago by
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).
comment:47 Changed 6 years ago by
Commit:  162035995bff5bb15d73bc13d2994655f24a8c23 → eb3da684a613c81294dec6eff22972152950790f 

comment:48 Changed 6 years ago by
Authors:  Stephan Ehlen, Peter Bruin → Stephan Ehlen 

Reviewers:  → Peter Bruin 
Replying to ehlen:
@pbruin 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 nonzero. We call XGCD(e, h, t, g, f)
; then e is the resultant of g and f (a nonzero 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 righthand 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...
comment:49 followup: 50 Changed 6 years ago by
@pbruin: 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 doublecheck.
comment:50 Changed 6 years ago by
Replying to ehlen:
@pbruin: 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 doublecheck.
I agree. By the way, doctests still pass.
comment:51 Changed 6 years ago by
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: 20130830. 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?
comment:52 Changed 6 years ago by
Reviewers:  Peter Bruin → Peter Bruin, Fredrik Stromberg 

Status:  needs_review → positive_review 
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 wellformatted 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.
comment:53 Changed 6 years ago by
Branch:  u/pbruin/20693sage_crashes_when_computing_newforms → eb3da684a613c81294dec6eff22972152950790f 

Resolution:  → fixed 
Status:  positive_review → closed 
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.
And crash...
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.