Opened 8 years ago

Closed 7 years ago

#10281 closed defect (fixed)

Multimodular echelon form over cyclotomic fields fails

Reported by: mraum Owned by: jason, was
Priority: critical Milestone: sage-5.0
Component: linear algebra Keywords: cyclotomic, echelon
Cc: Merged in: sage-5.0.beta14
Authors: William Stein Reviewers: Martin Raum
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by was)

The following example comes from Eisenstein series. It leads to previous_prime being called with invalid arguments.

sage: K.<rho> = CyclotomicField(106)
sage: coeffs = [(18603/107*rho^51 - 11583/107*rho^50 - 19907/107*rho^49 - 13588/107*rho^48 - 8722/107*rho^47 + 2857/107*rho^46 - 19279/107*rho^45 - 16666/107*rho^44 - 11327/107*rho^43 + 3802/107*rho^42 + 18998/107*rho^41 - 10798/107*rho^40 + 16210/107*rho^39 - 13768/107*rho^38 + 15063/107*rho^37 - 14433/107*rho^36 - 19434/107*rho^35 - 12606/107*rho^34 + 3786/107*rho^33 - 17996/107*rho^32 + 12341/107*rho^31 - 15656/107*rho^30 - 19092/107*rho^29 + 8382/107*rho^28 - 18147/107*rho^27 + 14024/107*rho^26 + 18751/107*rho^25 - 8301/107*rho^24 - 20112/107*rho^23 - 14483/107*rho^22 + 4715/107*rho^21 + 20065/107*rho^20 + 15293/107*rho^19 + 10072/107*rho^18 + 4775/107*rho^17 - 953/107*rho^16 - 19782/107*rho^15 - 16020/107*rho^14 + 5633/107*rho^13 - 17618/107*rho^12 - 18187/107*rho^11 + 7492/107*rho^10 + 19165/107*rho^9 - 9988/107*rho^8 - 20042/107*rho^7 + 10109/107*rho^6 - 17677/107*rho^5 - 17723/107*rho^4 - 12489/107*rho^3 - 6321/107*rho^2 - 4082/107*rho - 1378/107, 1, 4*rho + 1), (0, 1, rho + 4)]
sage: m = matrix(2, coeffs)
sage: v = m.echelon_form()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)

/home/martin/<ipython console> in <module>()

/opt/sage/local/lib/python2.6/site-packages/sage/matrix/matrix_cyclo_dense.so in sage.matrix.matrix_cyclo_dense.Matrix_cyclo_dense.echelon_form (sage/matrix/matrix_cyclo_dense.cpp:13469)()

/opt/sage/local/lib/python2.6/site-packages/sage/matrix/matrix_cyclo_dense.so in sage.matrix.matrix_cyclo_dense.Matrix_cyclo_dense._echelon_form_multimodular (sage/matrix/matrix_cyclo_dense.cpp:14366)()

/opt/sage/local/lib/python2.6/site-packages/sage/rings/arith.pyc in previous_prime(n)
   1044     n = ZZ(n)-1
   1045     if n <= 1:
-> 1046         raise ValueError, "no previous prime"
   1047     if n <= 3:
   1048         return ZZ(n)

ValueError: no previous prime

Apply:

  1. trac_10281-sage-5.0-beta9.patch

Install:

  1. http://wstein.org/home/wstein/patches/linbox-1.1.6.p8.spkg

Attachments (1)

trac_10281-sage-5.0-beta9.patch (25.8 KB) - added by was 7 years ago.
This patch is for beta9. I've made the one change that robertwb suggests to this patch (the formal cython change)

Download all attachments as: .zip

Change History (17)

comment:1 Changed 7 years ago by was

This is not just an issue with cyclotomic linear algebra. I'm doing some computations of modular symbols spaces of weight 4 with trivial character, and just hit "no previous prime". We need real 64-bit mod-p linear algebra to fix this issue. We need more primes!

comment:2 Changed 7 years ago by was

This turned out to be easier to resolve than it might have otherwise been, since Burcin Erocal and Martin Albrecht had improved matrices modn so they can now handle a modulus up to 223. However, they hadn't adapted the chinese remainder theorem infrastructure to take this into account.

I've written code that makes it so for dense linear algebra we have all the primes up to 223 at our disposal, instead of just the primes up to 46341. There are more, bigger, and better primes up to 223:

sage: n = prod(prime_range(46431))
sage: len(str(n))                 
20025
sage: time n = prod(prime_range(2^23))
Time: CPU 0.90 s, Wall: 0.90 s
sage: len(str(n))
3641670

sage: prime_pi(2^23)
564163
sage: prime_pi(46431)
4797
sage: prime_pi(2^23)//prime_pi(46431)
117

The upshot is that in various contexts we can deal with answers with up to 3.5 million digits, instead of only 20,025 digits, which was pretty feeble. Moreover, since the linear algebra is just as fast or faster, and we start with the biggest prime we can, certain linear algebra computations should be much faster with this new code.

I did *not* fix anything for sparse matrices, except avoiding introducing a major bug (caught by doctests). See #12679 for sparse matrices, as this should be another ticket.

Also, obviously, this fix is only a partial fix, and we'll someday have to deal with running out of primes on huge problems. But it'll take around 150 times as long to ever reach such problems, for what it is worth.

Note that it would probably be trivial to change the *old* matrix_modn_dense code to use 64-bit ints, hence support a modulus up to 231. However the new Erocal-Albrecht code can't go beyond 223. So that code -- which they deprecated -- should probably be resurrected. That's a project for later.

comment:3 Changed 7 years ago by was

  • Status changed from new to needs_review

comment:4 Changed 7 years ago by was

Note - the patch is bigger than expected since I fixed some of the formatting of docstrings in a file, just because.

comment:5 follow-up: Changed 7 years ago by mraum

Apply trac_10281-sage-5.0-beta9.patch

(for the patchbot)

I want the latest version to be tested on a variety of machines.

In principle, the patch looks good. It would be create if someone had a 32-bit system to test this, though.

comment:6 Changed 7 years ago by robertwb

Looks good to me. The only comment I have is that the changes adding

cdef extern from "stdint.h": 
    ctypedef unsigned long uint_fast64_t 

cdef extern from "multi_modular.h": 
    ctypedef uint_fast64_t mod_int 

are unnecessary, declaring it as an unsigned long in Cython is sufficient.

Changed 7 years ago by was

This patch is for beta9. I've made the one change that robertwb suggests to this patch (the formal cython change)

comment:7 in reply to: ↑ 5 Changed 7 years ago by was

Replying to mraum:

Apply trac_10281-sage-5.0-beta9.patch

(for the patchbot)

I want the latest version to be tested on a variety of machines.

In principle, the patch looks good. It would be create if someone had a 32-bit system to test this, though.

I've started a 32-bit Ubuntu 10.4 virtual machine on boxen, and I'm running tests on it. I'll report back.

comment:8 Changed 7 years ago by was

I tested on 32-bit and it works, but *only* after making a 1-line change to the Linbox spkg (which isn't needed on 64-bit, evidently).

I've posted the new linbox spkg here:

http://wstein.org/home/wstein/patches/linbox-1.1.6.p8.spkg

So to apply this ticket:

  1. Build new linbox: http://wstein.org/home/wstein/patches/linbox-1.1.6.p8.spkg
  2. Apply the patch to the sage repo: trac_10281-sage-5.0-beta9.patch

William

comment:9 Changed 7 years ago by was

  • Priority changed from major to critical

comment:10 follow-up: Changed 7 years ago by mraum

  • Description modified (diff)

I have tested the spkg, and it is completely ok. The commit message could be a bit more informative, but since it refers to the ticket, it's ok also.

William, do you want to make the change suggested by Robert? Personally, I would keep the extern statements for clearness (they obviously don't do any harm).

If you don't intend to, then I would set this to positive review, so that everything can be run on the build farm.

comment:11 in reply to: ↑ 10 Changed 7 years ago by was

Replying to mraum:

I have tested the spkg, and it is completely ok. The commit message could be a bit more informative, but since it refers to the ticket, it's ok also.

William, do you want to make the change suggested by Robert? Personally, I would keep the extern statements for clearness (they obviously don't do any harm).

I did make the changes and updated the second patch. I didn't bother with the beta2 patch, since that one will never get applied by the release manager.

If you don't intend to, then I would set this to positive review, so that everything can be run on the build farm.

Please do. Thanks!

I have to point out that I was running some big modular symbols computations over cyclotomic fields and found some "not enough primes" issues even with this patch. That's because we only get to use split primes. So it will be very important to revive the old matrix_modn_dense code, but using a 64-bit int instead of 32-bit, so we get modulo up to 2^32 (instead of 2^23). (That is already #4968.) I plan to do that this week or next week. The current patch here is very good foundation for doing #4968 right.

comment:12 Changed 7 years ago by mraum

  • Status changed from needs_review to positive_review

Great! Let me know when you have made progress on #4968. I hope to have time to review it.

comment:13 Changed 7 years ago by mraum

  • Authors set to William Stein
  • Reviewers set to Martin Raum

comment:14 Changed 7 years ago by was

  • Description modified (diff)

comment:15 Changed 7 years ago by jdemeyer

  • Milestone set to sage-5.0

comment:16 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.0.beta14
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.