Dear sage-devel,
Something is wrong with the multi-modular matrix multiplication code for
matrices over ZZ. At random, and infrequently, it gives incorrect
results. For example, the following code chooses random 3x2 and 2x10
integer matrices and multiplies them together using the multi-modular
algorithm. It then does the same multiplication 100 more times, checks
that the answer is always the same, and if not raises an exception:
sage: for n in range(2000):
....: A = MatrixSpace(ZZ,3,2).random_element()
....: B = MatrixSpace(ZZ,2,10).random_element()
....: try_once = A._multiply_multi_modular(B)
....: for k in range(100):
....: try_again = A._multiply_multi_modular(B)
....: if try_once != try_again:
....: print "="*60
....: print "n = %s, k = %s"%(n,k)
....: print "A = "
....: print A
....: print "B ="
....: print B
....: print "first attempt = "
....: print try_once
....: print "k-th retry = "
....: print try_again
....: raise RuntimeError
....:
This fails with very high probability (on Sage 4.6.2 under OS X, built
from source, and on Sage 4.6.2 under Red Hat Enterprise Linux, binary
install) with output such as:
============================================================
n = 27, k = 43
A =
[-1 0]
[ 1 0]
[ 0 1]
B =
[ -2 1 -1 0 0 -1 -1 3 -1 -1]
[ 1 1 -116 3 0 -1 -1 0 -1 0]
first attempt =
[ 2 -1 1 0 0 1 1 -3 1 1]
[ -2 1 -1 0 0 -1 -1 3 -1 -1]
[ 1 1 -116 3 0 -1 -1 0 -1 0]
k-th retry =
[ 2 1102 1 0 0 1 1 1100 1 1]
[1101 1 1102 0 0 1102 1102 3 1102 1102]
[ 1 1 987 3 0 1102 1102 0 1102 0]
---------------------------------------------------------------------------
RuntimeError [...]
Note that the two candidates for the matrix product here are congruent
modulo 1103, which is prime. If you rerun the code with verbose
logging, using set_verbose(2), then every time it fails the two
candidates for the matrix product are congruent modulo the prime being
used in the multi-modular algorithm. Thus I suspect that the Chinese
Remainder Theorem code in sage/ext/multi_modular.pyx is not handling a
corner case properly.
The relevant routine (mpz_crt_vec_tail) involves both Cython and GMP, so
if someone with Cython and/or GMP experience could give me a hand
debugging it then I would be very grateful.
This is trac #11358.
Yours,
Tom
--
Tom Coates
Royal Society University Research Fellow
Reader in Pure Mathematics
Imperial College London
--
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
