FFT space arithmetic

[fredrik-johansson]

This PR refactors fft_small for more easily working with transformed operands, and uses this in various places. This enables huge speedups for algorithms that can reuse the same operand for many multiplications, e.g. matrix multiplication -- the largest speedup I've measured so far is 15x for 64x64 fmpz_mat_mul with 100,000-digit entries.

• The previously monolithic fft_small multiplication routines simplify into function calls for creating a plan object, creating operands, transforming operands, doing pointwise operations, and transforming back
• We add gr rings for transformed nmod_poly, mpn_mod_poly and (signed) mpn integers. Others can easily be added in the future.
• Applications (so far): FFT based fmpz_mat_mul, nmod_poly_mat_mul, complex multiplication (fmpzi_mul), 2x2 matrix multiplication in gr_poly_gcd_hgcd. There is also truncating complex multiplication and matrix multiplication useful for acb_mul, arb_mat_mul, etc., but this is not yet tuned and integrated.

This was developed mainly using Claude Fable 5 over a couple of weeks; the code still needs some cleanup.

Important remarks

• Signed output conversions was missing and had to be implemented. The signed conversions are slower than than the unsigned ones.
• It is extremely important for performance that transforms are done fully in-place without allocating scratch space. The design consequence is that gr transform rings either need to maintain an internal scratch pool or use destructive output conversions instead of standard conversions. I opted for destructive conversions.
• The standard cost model assuming 1 input transform ~= 1 output transform is inaccurate: output conversion with CRT is generally much more expensive than input conversion, especially in the signed case. This means that the achievable speedups are often smaller than theory predict. For example, for 2x2 fmpz_mat multiplication transformed multiplication barely beats Strassen with unsigned coefficients and not at all with signed coefficients (currently FFT multiplication is only used by default here for n >= 3).
• Since FFT:ed operands consume a lot of memory, the generic transform rings decline when the expected memory usage for all operands exceeds the limit flint_fft_small_max_transformed_ring_size which currently defaults to 4 GB. Matrix multiplication still benefits for matrices much larger than this: the multiplication algorithm simply selects a smaller block size than the full matrix.
• Multithreading is not yet fully tuned. Hopefully existing code shouldn't slow down, but this may need checking.

Speedup for fmpz_mat_mul

Timings for integer matrices with uniform (fmpz_randbits) entries. FFT multiplication kicks in from around 7000-10000 bits for small n, the tuned cutoff increasing with larger n.

  n         bits       old        new   speedup

  3         8192   8.47e-05  8.46e-05   1.001
  3        16384   0.000258  0.000144   1.792
  3        32768    0.00073   0.00029   2.517
  3        65536    0.00133  0.000607   2.191
  3       131072    0.00262    0.0012   2.183
  3       262144    0.00552   0.00246   2.244
  3       524288     0.0112   0.00676   1.657
  3      1048576     0.0234    0.0122   1.918
  3      2097152     0.0482    0.0318   1.516

  4         8192   0.000175  0.000176   0.994
  4        16384   0.000527  0.000286   1.843
  4        32768    0.00155  0.000623   2.488
  4        65536    0.00275   0.00131   2.099
  4       131072    0.00532   0.00269   1.978
  4       262144     0.0112   0.00563   1.989
  4       524288     0.0226    0.0105   2.152
  4      1048576     0.0474    0.0274   1.730
  4      2097152     0.0983    0.0646   1.522

  8         8192    0.00126  0.000849   1.484
  8        16384    0.00363   0.00145   2.503
  8        32768     0.0109    0.0028   3.893
  8        65536     0.0191   0.00576   3.316
  8       131072     0.0369    0.0123   3.000
  8       262144      0.077    0.0282   2.730
  8       524288      0.157    0.0562   2.794
  8      1048576      0.331     0.112   2.955
  8      2097152      0.676     0.288   2.347

 16         8192    0.00912   0.00444   2.054
 16        16384     0.0264   0.00718   3.677
 16        32768     0.0792    0.0147   5.388
 16        65536       0.14    0.0318   4.403
 16       131072      0.272    0.0668   4.072
 16       262144      0.566     0.133   4.256
 16       524288      1.146     0.271   4.229
 16      1048576      2.403     0.563   4.268
 16      2097152      4.944     1.357   3.643

 32         8192     0.0652    0.0288   2.264
 32        16384      0.191    0.0488   3.914
 32        32768      0.573    0.0975   5.877
 32        65536      0.994     0.191   5.204
 32       131072      1.923     0.378   5.087
 32       262144      4.015     0.761   5.276
 32       524288      8.112     1.565   5.183
 32      1048576     16.972     3.238   5.242
 32      2097152     34.915     7.483   4.666

 64         8192       0.42     0.243   1.728
 64        16384      1.191     0.351   3.393
 64        32768      3.396     0.639   5.315
 64        65536      9.608     1.216   7.901
 64       131072     26.475     2.405   11.008
 64       262144     71.051     4.869   14.593
 64       524288    178.076    11.586   15.370
 64      1048576     118.14     21.37   5.528
 64      2097152    242.875    50.687   4.792

128         8192      1.812     1.815   0.998
128        16384      5.072     2.614   1.940
128        32768     14.157     4.591   3.084
128        65536     38.849     8.739   4.445
128       131072    108.001    17.541   6.157
128       262144     288.81    35.733   8.082
128       524288    719.603    92.443   7.784
128      1048576    1740.92   178.814   9.736

256         8192      8.073     8.067   1.001
256        16384     21.536    21.533   1.000
256        32768     59.463     34.65   1.716
256        65536    163.358    68.146   2.397

512         8192     38.683    38.687   1.000
512        16384      98.91    98.926   1.000
512        32768    262.065     261.9   1.001
512        65536    712.764   600.961   1.186

Note: matrix mulhigh is implemented but not yet used/profiled.

Speedup for nmod_poly_mat_mul

Example timings modulo p = nextprime(2^63), uniform entries.

  n          len        old        new  speedup

  2           64   2.39e-05  2.06e-05   1.160
  2          256   7.81e-05  4.32e-05   1.808
  2         1024   0.000401  0.000323   1.241
  2         4096    0.00178   0.00155   1.148
  2        16384    0.00779   0.00678   1.149
  2        65536     0.0342    0.0313   1.093
  2       262144       0.16     0.135   1.185

  3           64    8.5e-05  5.16e-05   1.647
  3          256   0.000311  0.000108   2.880
  3         1024    0.00141  0.000533   2.645
  3         4096    0.00616    0.0023   2.678
  3        16384     0.0268    0.0107   2.505
  3        65536      0.117    0.0695   1.683
  3       262144       0.54     0.311   1.736

  4           64   0.000202  9.71e-05   2.080
  4          256   0.000758  0.000218   3.477
  4         1024    0.00335     0.001   3.350
  4         4096     0.0145   0.00429   3.380
  4        16384     0.0635    0.0206   3.083
  4        65536      0.278     0.124   2.242
  4       262144      1.295     0.563   2.300

  8           64    0.00167  0.000516   3.236
  8          256    0.00621   0.00115   5.400
  8         1024     0.0272   0.00492   5.528
  8         4096      0.119     0.023   5.174
  8        16384      0.516     0.133   3.880
  8        65536      2.238     0.561   3.989
  8       262144      10.51      3.69   2.848

 16           64     0.0108   0.00297   3.636
 16          256     0.0511   0.00644   7.935
 16         1024      0.223    0.0294   7.585
 16         4096      0.963     0.145   6.641
 16        16384      4.181     0.664   6.297
 16        65536     18.082     2.811   6.433

 32           64     0.0766    0.0204   3.755
 32          256      0.414    0.0458   9.039
 32         1024      1.796     0.211   8.512
 32         4096       7.69     0.873   8.809
 32        16384     33.537     3.739   8.970

 64           64       0.51     0.512   0.996
 64          256      3.378     0.389   8.684
 64         1024     14.382     1.368   10.513

128           64      1.038     1.033   1.005
128          256       7.64     2.772   2.756
128         1024     52.733     9.307   5.666

Note: matrix mulmid is implemented but not yet used/profiled.

Polynomial HGCD

The generic gr_poly_gcd_hgcd gets an overload for its 2x2 matrix multiplication to allow using FFT instead of Strassen. The speedup for the GCD as a whole is small but measurable. Timings for nmod_poly, p being a 50-bit FFT prime:

     len        old       new   speedup

     100   1.38e-05  1.38e-05   1.000
     300   9.07e-05  9.08e-05   0.999
     900   0.000821  0.000822   0.999
    2700    0.00293   0.00289   1.014
    8100      0.012    0.0117   1.026
   24300     0.0379    0.0365   1.038
   72900      0.144     0.128   1.125
  218700      0.463     0.439   1.055
  656100      1.713     1.537   1.115
 1968300      5.511     5.197   1.060
 5904900     18.686    17.705   1.055

Timings for p = nextprime(2^63):

     len        old       new   speedup

     100   1.57e-05  1.56e-05   1.006
     300   0.000107  0.000108   0.991
     900   0.000967  0.000968   0.999
    2700    0.00582   0.00561   1.037
    8100     0.0237    0.0225   1.053
   24300     0.0978    0.0907   1.078
   72900      0.346     0.312   1.109
  218700      1.256     1.144   1.098
  656100      4.406     3.975   1.108
 1968300     15.085      13.7   1.101
 5904900     55.245    51.534   1.072

Timings for mpn_mod, p = nextprime(2^127):

     len        old       new   speedup

     100   9.06e-05  9.14e-05   0.991
     300    0.00107   0.00107   1.000
     900    0.00559   0.00559   1.000
    2700     0.0247    0.0236   1.047
    8100     0.0963    0.0911   1.057
   24300      0.366     0.334   1.096
   72900      1.346     1.204   1.118
  218700      4.719      4.23   1.116
  656100     16.905    15.366   1.100

Complex multiplication

Timings for fmpzi_mul and fmpzi_sqr, which become wrappers around the new flint_mpn_mul_complex and flint_mpn_sqr_complex:

      bits        old       new   speedup         old       new   speedup

        50   9.77e-09  9.37e-09   1.043      8.84e-09  8.19e-09   1.079
       150   5.31e-08  5.05e-08   1.051      4.62e-08  3.69e-08   1.252
       450   1.24e-07   1.3e-07   0.954         1e-07  9.02e-08   1.109
      1350   6.32e-07  6.08e-07   1.039      4.66e-07  3.91e-07   1.192
      4050   3.29e-06  3.28e-06   1.003      2.38e-06   2.1e-06   1.133
     12150   1.92e-05  1.91e-05   1.005      1.29e-05  1.25e-05   1.032
     36450   8.02e-05     8e-05   1.002      6.53e-05  5.31e-05   1.230
    109350   0.000249  0.000238   1.046      0.000189  0.000164   1.152
    328050   0.000884  0.000817   1.082      0.000681  0.000576   1.182
    984150    0.00275   0.00246   1.118       0.00209   0.00181   1.155
   2952450    0.00942   0.00773   1.219       0.00727   0.00625   1.163
   8857350     0.0283    0.0258   1.097        0.0205    0.0186   1.102
  26572050      0.099    0.0967   1.024        0.0725    0.0654   1.109
  79716150      0.342     0.303   1.129         0.251     0.218   1.151
 239148450      1.117      0.97   1.152         0.816     0.746   1.094
 717445350      3.857     3.232   1.193         2.815     2.584   1.089
2152336050     15.717    16.154   0.973        10.283    10.477   0.981
6457008150     51.441    52.725   0.976        34.021    34.203   0.995

Some nice 10-20% speedups. It should be noted that the speedup for squaring over most of the range just comes from a better Karatsuba variant (see #2781); FFT squaring is only used from 190M bits (vs 44K bits for multiplication).

At 2152336050 bits the 4 GB transform limit kicks in and we fall back to Strassen as before (the new code is 2% slower due to differences in memory allocation).

If I bump the 4 GB limit I get:

2152336050     16.275    12.275   1.326        10.335     8.045   1.285

but at 6457008150 bits my machine runs out of memory.

Complex mulhigh implemented but not yet used/profiled.