Serious bug in cyclotomic matrix multiplication

[daveloeffler]

I've been investigating trac #8541, where a user reported an
unexpected error when calculating some modular forms. This turns out
to be due to an extremely evil bug in matrix multiplication over
cyclotomic fields:

sage: K.<zeta4> = CyclotomicField(4)
sage: m = matrix(K, [186])
sage: n = matrix(K, [125])
sage: m * n
[-23087]

Investigating with random integer matrices coerced into K, the answer
seems to be wrong in a fairly substantial fraction of cases. I'm
posting here because I think this rates as a sufficiently serious bug
to rate as a blocker for the next release, but I didn't feel I had the
authority to mark it as a blocker myself -- do people agree with my
assessment of this bug, and more importantly does anyone know how to
fix it?

David

[Tom Boothby]

+1 that's a blocker if I've ever seen one.

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org
>
> To unsubscribe, reply using "remove me" as the subject.
>

[Dima Pasechnik]

+10^6 for this being a blocker!

[William Stein]

On Friday, April 9, 2010, Dima Pasechnik <dim...@gmail.com> wrote:
> +10^6 for this being a blocker!

And I'll fix this this weekend if nobody else does. It is probably a
problem doing a multi modular matrix multiply and not having enough
primes. I designed and mostly implemented the cycle linalg
algorithms in sage, so this is my fault, probably.

It is very fast. It gets the wrong answer faster than magma.

William

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org
>
> To unsubscribe, reply using "remove me" as the subject.
>

--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

[Craig Citro]

> And I'll fix this this weekend if nobody else does.  It is probably a
> problem doing a multi modular matrix multiply and not having enough
> primes.   I designed and mostly implemented the cycle linalg
> algorithms in sage, so this is my fault, probably.
>

I helped with the cyclo linalg code, so I took a quick look -- it
looks like this is an honest problem with the code for matrices over
ZZ:

sage: from sage.matrix.matrix_integer_dense import _lift_crt
sage: _lift_crt(matrix(ZZ,2,1), [matrix(Integers(46337),2,1,[23250,
0])])
[-23087]
[ 0]

Obviously there's just an issue with which residue is getting
returned, because 23087 + 23250 = 46337. Probably a mishandled corner
case in _lift_crt -- I glanced at mpz_crt_vec in sage/ext/
multi_modular.pyx, and lo, and behold! There's a comment that's
telling us what's going on -- around line 589:

# normalize to be between -prod/2 and prod/2.
if mpz_cmp(z[j], self.half_product) > 0:
mpz_sub(z[j], z[j], self.product)

Indeed, that's where the residue is getting changed. So we just need
to change which normalization we return -- or document what _lift_crt
does better, and change the way we use it in cyclotomic matrix
multiply.

> It is very fast.   It gets the wrong answer faster than magma.
>

:)

-cc

[Craig Citro]

> And I'll fix this this weekend if nobody else does.  It is probably a
> problem doing a multi modular matrix multiply and not having enough
> primes.   I designed and mostly implemented the cycle linalg
> algorithms in sage, so this is my fault, probably.
>

... and a fix is up. Robert Bradshaw pointed out that the bound is
just wrong -- we need an extra factor of 2 to deal with the sign. I've
created a ticket and posted a patch:

http://trac.sagemath.org/sage_trac/ticket/8666

Probably wouldn't hurt if someone wanted to do some speed tests; I'd
expect the factor of 2 wouldn't hurt us, but those are famous last
words ...

-cc