bug in cyclotomic matrix multiplication

[dimpase]

As reported on ​https://ask.sagemath.org/question/64194/determinants-over-cyclotomic-fields-are-broken/ and on ​https://groups.google.com/d/msgid/sage-devel/44f0ff3c-71db-4555-a494-80b2ae222c22n%40googlegroups.com

K.<z> = CyclotomicField(16)
OK = K.ring_of_integers()
L = [[-575*z^7 - 439*z^6 - 237*z^5 + 237*z^3 + 439*z^2 + 575*z + 623, 0],
    [0,     -114*z^7 - 88*z^6 - 48*z^5 + 48*z^3 + 88*z^2 + 114*z + 123]]
U = [[-1926*z^7 - 1474*z^6 - 798*z^5 + 798*z^3 + 1474*z^2 + 1926*z + 2085, 0],
    [0,   -1014*z^7 - 777*z^6 - 421*z^5 + 421*z^3 + 777*z^2 + 1014*z + 1097]]
L, U = matrix(K,L), matrix(K,U)
LU = matrix( [ [L[i].inner_product(U.transpose()[j]) for j in range(2)] for i in range(2)] )
-LU[0][0]+(L*U)[0][0]

does not report 0, but 8388593*z^7 - 8388593*z - 8388593

I checked using another cyclotomics implementation, in GAP, that basic arithmetic over K, (i.e. computation of LU) is correct.

[dimpase]

Sage's matrix multiplication is clearly wrong here:

sage: libgap(L)*libgap(U)-libgap(LU)
[ [ 0, 0 ], [ 0, 0 ] ]
sage: libgap(L)*libgap(U)-libgap(L*U)
[ [ 8388593+8388593*E(16)-8388593*E(16)^7, 0 ], [ 0, 0 ] ]

[dimpase]

note that 8388593 is a prime, so there is some over-optimisation going on, I guess.

[dimpase]

The implementation is in src/sage/matrix/matrix_cyclo_dense.pyx

cdef _matrix_times_matrix_(self, baseMatrix right):
        """
        Return the product of two cyclotomic dense matrices.

        INPUT:
            self, right -- cyclotomic dense matrices with compatible
                           parents (same base ring, and compatible
                           dimensions for matrix multiplication).

        OUTPUT:
            cyclotomic dense matrix

        ALGORITHM:
            Use a multimodular algorithm that involves multiplying the
            two matrices modulo split primes.

so this explains that appearance of a largish prime in the error is not a fluke.

[dimpase]

the prime that pops up here is quite special:

sage: from sage.matrix.matrix_modn_dense_double import MAX_MODULUS as MAX_MODULUS_modn_dense_double
....: from sage.arith.multi_modular import MAX_MODULUS as MAX_MODULUS_multi_modular
....: MAX_MODULUS = min(MAX_MODULUS_modn_dense_double, MAX_MODULUS_multi_modular)
....: 
sage: MAX_MODULUS
8388608
sage: previous_prime(MAX_MODULUS)
8388593

Could it be one actually need "more primes" in this case, whatever this exactly means?

[dimpase]

If I take half MAX_MODULUS (which is 2²³ on my machine)

src/sage/matrix/matrix_cyclo_dense.pyx

@@ -652 +652 @@
         bound = 1 + 2 * A.height() * B.height() * self._ncols
 
         n = self._base_ring._n()
-        p = previous_prime(MAX_MODULUS)
+        p = previous_prime(2**22)
         prod = 1

then this particular error goes away, while all the tests in src/sage/matrix/matrix_cyclo_dense.pyx still pass.

[dimpase]

here is a fix, but why it works needs to be understood.

---

New commits:

​f4f1bf2

roughly halve maximum prime to start with, cf #34597

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​20f69aa

moved test in the right place, reduce MAX_MODULUS globally