Opened 4 months ago

# bug in cyclotomic matrix multiplication

Reported by: Owned by: dimpase critical linear algebra was Dima Pasechnik N/A u/dimpase/matrix/cyclot_matrix_mult_fix 20f69aab69ebb4abc96d5f249a2cbc887a1c5fd2

### Description

```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.

### comment:1 Changed 4 months ago by 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 ] ]
```

### comment:2 Changed 4 months ago by dimpase

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

### comment:3 Changed 4 months ago by 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.

### comment:4 Changed 4 months ago by 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?

### comment:5 Changed 4 months ago by dimpase

If I take half `MAX_MODULUS` (which is 223 on my machine)

• ## src/sage/matrix/matrix_cyclo_dense.pyx

 a cdef class Matrix_cyclo_dense(Matrix_dense): 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.

### comment:6 Changed 4 months ago by dimpase

Authors: → Dima Pasechnik → u/dimpase/matrix/cyclot_matrix_mult_fix → f4f1bf23a3bc843e62e2137a2b82c312f5180665 new → needs_info

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

New commits:

 ​f4f1bf2 `roughly halve maximum prime to start with, cf #34597`

### comment:7 Changed 4 months ago by git

Commit: f4f1bf23a3bc843e62e2137a2b82c312f5180665 → 20f69aab69ebb4abc96d5f249a2cbc887a1c5fd2

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

 ​20f69aa `moved test in the right place, reduce MAX_MODULUS globally`

### comment:8 Changed 10 days ago by mkoeppe

Milestone: sage-9.8
Note: See TracTickets for help on using tickets.