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

Multimodular echelon form over cyclotomic fields fails #10281

Closed
sagetrac-mraum mannequin opened this issue Nov 16, 2010 · 16 comments
Closed

Multimodular echelon form over cyclotomic fields fails #10281

sagetrac-mraum mannequin opened this issue Nov 16, 2010 · 16 comments

Comments

@sagetrac-mraum
Copy link
Mannequin

sagetrac-mraum mannequin commented Nov 16, 2010

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. attachment: trac_10281-sage-5.0-beta9.patch​

Install:

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

Component: linear algebra

Keywords: cyclotomic, echelon

Author: William Stein

Reviewer: Martin Raum

Merged: sage-5.0.beta14

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

@williamstein
Copy link
Contributor

comment:1

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!

@williamstein
Copy link
Contributor

comment:2

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 2^23. 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 2^23 at our disposal, instead of just the primes up to 46341. There are more, bigger, and better primes up to 2^23:

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 2^31. However the new Erocal-Albrecht code can't go beyond 2^23. So that code -- which they deprecated -- should probably be resurrected. That's a project for later.

@williamstein
Copy link
Contributor

comment:4

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

@sagetrac-mraum
Copy link
Mannequin Author

sagetrac-mraum mannequin commented Mar 17, 2012

comment:5

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.

@robertwb
Copy link
Contributor

comment:6

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.

@williamstein
Copy link
Contributor

Attachment: trac_10281-sage-5.0-beta9.patch.gz

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

@williamstein
Copy link
Contributor

comment:7

Replying to @sagetrac-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.

@williamstein
Copy link
Contributor

comment:8

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

@sagetrac-mraum
Copy link
Mannequin Author

sagetrac-mraum mannequin commented Mar 20, 2012

comment:10

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.

@sagetrac-mraum

This comment has been minimized.

@williamstein
Copy link
Contributor

comment:11

Replying to @sagetrac-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.

@sagetrac-mraum
Copy link
Mannequin Author

sagetrac-mraum mannequin commented Mar 20, 2012

comment:12

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

@sagetrac-mraum
Copy link
Mannequin Author

sagetrac-mraum mannequin commented Mar 20, 2012

Author: William Stein

@sagetrac-mraum
Copy link
Mannequin Author

sagetrac-mraum mannequin commented Mar 20, 2012

Reviewer: Martin Raum

@williamstein

This comment has been minimized.

@jdemeyer jdemeyer added this to the sage-5.0 milestone Mar 28, 2012
@jdemeyer
Copy link

Merged: sage-5.0.beta14

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

4 participants