Faster matrix multiplication mod p <= 255

[fredrik-johansson]

Adds dedicated matrix multiplication code modulo $p \le 255$:

• nmod_mat_mul_u8, called from nmod_mat_mul when applicable
• _nmod_mat_mul_u8 acting directly on byte arrays, called internally (with conversion) from nmod_mat_mul_u8, and directly from gr_mat_mul with the nmod8 context

There are three internal algorithms, all of which beat nmod_mat_mul_blas; small GF(p) now reach teraop/second-equivalent speeds on my humble Zen3 laptop.

• GF(2) and GF(3) use bit-packing with M4R and Strassen
• $4 \le p \le 16$ use AVX2/AVX512/NEON vector byte-shuffling instructions (with a plain-C uint16 fallback) + Strassen with uint8 modular additions
• Larger moduli use flint_sgemm + block accumulation + Strassen with uint8 modular additions

I considered bit-packing for GF(5) and GF(7) too but preliminary tests didn't show this to be competitive with the byte-shuffling approach, so I've dropped this idea for now.

Single-threaded performance, n is the matrix dimension:

                 nmod_mat_mul                gr_mat_mul/nmod8
  p      n       old       new  speedup       old       new  speedup

  2     10  6.53e-07  3.96e-07 (  1.6)   9.57e-07  2.67e-07 (  3.6)     7.5 Gop/s
  2     30  7.73e-06  1.64e-06 (  4.7)   1.28e-05  1.11e-06 ( 11.5)    48.6 Gop/s
  2    100  0.000125  2.23e-05 (  5.6)    0.00026  1.72e-05 ( 15.1)   116.3 Gop/s
  2    300  0.000961  0.000121 (  7.9)    0.00507  7.23e-05 ( 70.1)   746.9 Gop/s
  2   1000    0.0242   0.00164 ( 14.8)      0.158  0.000591 (267.3)  3384.1 Gop/s
  2   3000     0.585    0.0208 ( 28.1)       3.98   0.01006 (395.6)  5367.8 Gop/s
  2  10000    19.527     0.446 ( 43.8)    108.562     0.323 (336.1)  6192.0 Gop/s

  3     10  6.74e-07  4.04e-07 (  1.7)   9.88e-07  2.68e-07 (  3.7)     7.5 Gop/s
  3     30  3.43e-06  1.71e-06 (  2.0)   1.31e-05  1.12e-06 ( 11.7)    48.2 Gop/s
  3    100  6.31e-05  1.83e-05 (  3.4)   0.000266   1.4e-05 ( 19.0)   142.9 Gop/s
  3    300  0.000687   0.00035 (  2.0)    0.00509  0.000266 ( 19.1)   203.0 Gop/s
  3   1000    0.0209    0.0038 (  5.5)      0.161   0.00249 ( 64.7)   803.2 Gop/s
  3   3000     0.589    0.0502 ( 11.7)      4.007    0.0381 (105.2)  1417.3 Gop/s
  3  10000    19.804      2.04 (  9.7)    109.389     1.978 ( 55.3)  1011.1 Gop/s

 13     10  6.81e-07   4.2e-07 (  1.6)   9.96e-07  2.82e-07 (  3.5)     7.1 Gop/s
 13     30  4.62e-06  1.75e-06 (  2.6)   1.32e-05  1.19e-06 ( 11.1)    45.4 Gop/s
 13    100  8.92e-05  1.93e-05 (  4.6)   0.000269  1.48e-05 ( 18.2)   135.1 Gop/s
 13    300  0.000694  0.000384 (  1.8)    0.00511  0.000358 ( 14.3)   150.8 Gop/s
 13   1000    0.0214   0.00992 (  2.2)       0.16   0.00845 ( 18.9)   236.7 Gop/s
 13   3000     0.594     0.233 (  2.5)      4.023     0.203 ( 19.8)   266.0 Gop/s
 13  10000    19.905     7.746 (  2.6)    109.195     7.409 ( 14.7)   269.9 Gop/s

251     10  6.78e-07  5.68e-07 (  1.2)   9.89e-07  4.31e-07 (  2.3)     4.6 Gop/s
251     30  5.19e-06  2.78e-06 (  1.9)   1.32e-05  2.17e-06 (  6.1)    24.9 Gop/s
251    100   0.00016  3.26e-05 (  4.9)   0.000268  2.77e-05 (  9.7)    72.2 Gop/s
251    300  0.000693  0.000599 (  1.2)    0.00517  0.000559 (  9.2)    96.6 Gop/s
251   1000    0.0212    0.0202 (  1.0)       0.16    0.0187 (  8.6)   107.0 Gop/s
251   3000     1.134     0.528 (  2.1)      4.084     0.499 (  8.2)   108.2 Gop/s
251  10000    38.623    14.269 (  2.7)    110.587     13.98 (  7.9)   143.1 Gop/s

Multithreaded performance (8 threads):

                 nmod_mat_mul                gr_mat_mul/nmod8
  p      n       old       new  speedup       old       new  speedup

  2     10  5.99e-07  3.96e-07 (  1.5)   9.26e-07  2.65e-07 (  3.5)     7.5 Gop/s
  2     30  1.72e-05  1.65e-06 ( 10.4)   1.26e-05  1.12e-06 ( 11.2)    48.2 Gop/s
  2    100   3.6e-05  2.21e-05 (  1.6)   0.000258  1.72e-05 ( 15.0)   116.3 Gop/s
  2    300   0.00025   0.00012 (  2.1)    0.00508  7.11e-05 ( 71.4)   759.5 Gop/s
  2   1000    0.0059   0.00152 (  3.9)      0.159  0.000497 (319.9)  4024.1 Gop/s
  2   3000      0.13    0.0126 ( 10.3)      3.969     0.005 (793.8) 10800.0 Gop/s
  2  10000     3.722     0.293 ( 12.7)    107.873     0.182 (592.7) 10989.0 Gop/s

  3     10   6.3e-07  3.98e-07 (  1.6)   9.73e-07  2.68e-07 (  3.6)     7.5 Gop/s
  3     30  1.91e-05  1.69e-06 ( 11.3)    1.3e-05  1.11e-06 ( 11.7)    48.6 Gop/s
  3    100  0.000136  1.82e-05 (  7.5)   0.000262  1.38e-05 ( 19.0)   144.9 Gop/s
  3    300   0.00022   0.00035 (  0.6)    0.00512  0.000267 ( 19.2)   202.2 Gop/s
  3   1000    0.0056   0.00344 (  1.6)       0.16   0.00215 ( 74.4)   930.2 Gop/s
  3   3000     0.132    0.0256 (  5.2)      3.978    0.0183 (217.4)  2950.8 Gop/s
  3  10000     3.758     1.131 (  3.3)    108.308     1.058 (102.4)  1890.4 Gop/s

 13     10  6.33e-07  4.16e-07 (  1.5)   9.84e-07   2.8e-07 (  3.5)     7.1 Gop/s
 13     30   2.6e-05  1.79e-06 ( 14.5)   1.31e-05  1.18e-06 ( 11.1)    45.8 Gop/s
 13    100   4.2e-05  1.91e-05 (  2.2)   0.000265  1.45e-05 ( 18.3)   137.9 Gop/s
 13    300    0.0002  0.000383 (  0.5)    0.00517  0.000349 ( 14.8)   154.7 Gop/s
 13   1000    0.0057    0.0044 (  1.3)       0.16    0.0026 ( 61.5)   769.2 Gop/s
 13   3000     0.131     0.062 (  2.1)      4.004     0.045 ( 89.0)  1200.0 Gop/s
 13  10000     3.776     1.823 (  2.1)    108.626     1.672 ( 65.0)  1196.2 Gop/s

251     10   6.4e-07  5.74e-07 (  1.1)   9.96e-07  4.31e-07 (  2.3)     4.6 Gop/s
251     30     3e-05  2.79e-06 ( 10.8)   1.33e-05   2.2e-06 (  6.0)    24.5 Gop/s
251    100   6.6e-05  3.08e-05 (  2.1)   0.000287  2.53e-05 ( 11.3)    79.1 Gop/s
251    300   0.00024   0.00028 (  0.9)    0.00551   0.00022 ( 25.0)   245.5 Gop/s
251   1000     0.006     0.008 (  0.8)      0.198    0.0056 ( 35.4)   357.1 Gop/s
251   3000     0.252     0.206 (  1.2)      4.155     0.197 ( 21.1)   274.1 Gop/s
251  10000     8.067     5.788 (  1.4)    111.191     5.523 ( 20.1)   362.1 Gop/s

Note that the nmod_mat_mul "old" timings are with flint_sgemm/flint_dgemm from #2810 enabled, so the speedups are essentially versus nmod_mat_mul_blas. The much higher speedup ratios for nmod8 are indicative of the improvement of the new nmod_mat_mul over the old nmod_mat_mul without a BLAS backend.

Note that gr_mat_mul with nmod8 coefficients now always beats nmod_mat_mul, as it should be. For GF(2), the cost of streaming ulong coefficients in and out of memory is virtually as high as that of actually multiplying in packed representation.

Developed using Claude Fable 5 and Opus 5.