Two-prime and cyclic fft_small-based mpn multiplication

[fredrik-johansson]

The fft_small based integer multiplication previously only supported using 3 to 8 FFT primes. We add a specialized 2-prime version which has much less overhead for CRT and conversions.

This speeds up fft_small-based flint_mpn_mul and friends at every size up to about 12,000 limbs, usually by about 10-20%. The peak speedup I've measured is more than 50% occuring at 689 limbs (13,000 digits). The GMP Toom -> fft_small crossover for full multiplication drops from 400 to 240 limbs (7700 -> 4600 digits) on Zen 3 and to 300 limbs for mulhigh. Plot of the speedup for flint_mpn_mul:

Speedup for flint_mpn_mul, flint_mpn_sqr, flint_mpn_mulhigh_n, flint_mpn_sqrhigh, n in limbs:

     n    mul    sqr  mulhi   sqrhi
   220  1.010  0.998  0.995  1.000
   231  1.000  1.004  1.005  1.000
   242  1.031  1.003  0.996  1.000
   254  1.073  1.003  1.000  1.007
   266  1.117  1.001  0.999  0.998
   279  0.991  1.001  1.001  1.000
   292  1.034  1.001  0.999  1.004
   306  1.083  0.998  1.009  1.002
   321  1.172  1.000  1.075  1.001
   337  1.260  0.995  1.148  1.003
   353  1.294  1.000  1.220  1.001
   370  1.112  1.000  1.000  0.998
   388  1.141  1.000  1.076  0.998
   407  1.175  1.056  1.151  1.001
   427  1.190  1.125  1.291  1.000
   448  1.078  1.070  1.118  0.991
   470  1.104  1.123  1.228  0.984
   493  1.204  1.181  1.308  1.080
   517  1.220  1.230  1.385  1.142
   542  1.068  1.124  1.130  1.040
   569  1.093  1.190  1.199  1.126
   597  1.402  1.253  1.256  1.209
   626  1.346  1.256  1.222  1.182
   657  1.370  1.322  1.250  1.275
   689  1.541  1.374  1.362  1.349
   723  1.362  1.257  1.203  1.223
   759  1.401  1.307  1.290  1.314
   796  1.116  1.323  1.265  1.281
   835  1.135  1.152  1.292  1.354
   876  1.161  1.160  1.361  1.438
   919  1.095  1.099  1.245  1.290
   964  1.373  1.443  1.275  1.320
  1012  1.151  1.143  1.270  1.303
  1062  1.088  1.097  1.254  1.249
  1115  1.111  1.124  1.279  1.309
  1170  1.302  1.349  1.242  1.315
  1228  1.337  1.384  1.289  1.369
  1289  1.309  1.365  1.265  1.338
  1353  1.362  1.441  1.309  1.333
  1420  1.201  1.267  1.141  1.188
  1491  1.201  1.280  1.141  1.223
  1565  1.083  1.072  1.075  1.060
  1643  1.063  1.066  1.050  1.048
  1725  1.037  1.043  1.020  1.022
  1811  1.067  1.050  1.050  1.026
  1901  1.090  1.081  1.082  1.070
  1996  1.086  1.084  1.076  1.066
  2095  1.047  1.054  1.050  1.031
  2199  1.052  1.052  1.043  1.028
  2308  1.233  1.270  1.186  1.209
  2423  1.211  1.218  1.149  1.152
  2544  1.223  1.240  1.170  1.176
  2671  1.212  1.240  1.163  1.175
  2804  1.148  1.178  1.074  1.085
  2944  1.142  1.200  1.071  1.100
  3091  1.083  1.046  1.062  1.039
  3245  1.060  1.062  1.044  1.025
  3407  1.076  1.068  1.057  1.040
  3577  1.048  1.040  1.040  1.000
  3755  1.109  1.100  1.090  1.081
  3942  1.181  1.221  1.130  1.147
  4139  1.061  1.061  1.040  1.024
  4345  1.058  1.045  1.033  1.015
  4562  1.153  1.188  1.099  1.101
  4790  1.179  1.192  1.121  1.120
  5029  1.157  1.194  1.100  1.113
  5280  1.156  1.198  1.110  1.110
  5544  1.120  1.155  1.036  1.056
  5821  1.100  1.149  1.019  1.027
  6112  1.102  1.155  1.015  1.041
  6417  1.073  1.064  1.028  1.000
  6737  1.072  1.072  1.020  0.995
  7073  1.080  1.073  1.019  1.000
  7426  1.112  1.114  1.061  1.043
  7797  1.108  1.117  1.053  1.041
  8186  1.083  1.088  0.992  0.953
  8595  1.081  1.091  1.026  1.004
  9024  1.169  1.210  1.074  1.087
  9475  1.165  1.202  1.071  1.076
  9948  1.200  1.241  1.097  1.105
 10445  1.189  1.241  1.084  1.095
 10967  1.110  1.166  1.008  1.024
 11515  1.109  1.165  1.002  1.021

The PR also adds cyclic mpn products (multiplying mod 2^N-1) for future use and makes some tweaks to the negacyclic code.

Developed using Claude Fable 5.

The updated tuning parameters for Zen 3 and Zen 4 are from my Zen 3 machine; the Skylake and later parameters are Claude's. Someone may want to double check these and do some tuning for arm64 as well.

Something to note is that there are three versions of the code for 2-prime input conversion. The one chosen for x86-64 uses byte-unaligned word loads, which in practice performed better than a SIMD version. The SIMD version is used on arm64 out of caution that unaligned loads may be much slower there, but neither version has been timed.

[fredrik-johansson]

Two comments:

The speedup actually extends to 100,000 limbs or so, maybe half that for high products. Not sure if I'll bump the cutoff here or in a followup PR with further tuning.

This speeds up transform-based complex multiplication, but 2x2 matrix multiplication slows down at some sizes. The reason seems to be that complex multiplication uses the cached mpn_ctx buffer for its scratch space while matrix multiplication mallocs its operands and suffers extra page faults with the larger allocations. I think the solution is to redesign memory allocation for transformed rings to allow using the mpn_ctx buffer everywhere.

[fredrik-johansson]

Fixing some loose ends here forced refactoring the transformed mpn ring memory management for the better (simpler interface, faster).