Overhaul solving code for fmpz_mat and fmpq_mat

[fredrik-johansson]

Near-completely rewrite of the fmpz_mat and fmpq_mat linear solving code, done using Claude Fable 5.

This is generally on par with the old code for small systems with small bit sizes while giving order-of-magnitude speedups in large dimensions and/or large bit sizes.

• Termination uses heuristic probing (checking stabilization for a random linear combination of output coefficients) which runs frequently, instead of verifying AX = B every 1.4^k bits, so there is less overshoot. We also use continued fractions to identify unbalanced solutions, so there is even less overshoot than before e.g. when the solution has small denominator.
• Certification is based on a cheap bound certificate obtained from the reconstructed modular solution (essentially the same strategy as IML) instead of explicitly verifying AX = B
• At high precision, the Dixon solver switches from p-adic to p^k-adic lifting with an adaptive k
• The multimodular solver uses batched moduli and fast multimodular reduction/CRT
• Dixon starts with a small (30-bit) prime and a mod p LU factorization and only switches to a large (60-bit) prime and a mod p inverse after a few iterations when this pays off
• The can_solve algorithm is simplified to use the nonsingular solver as a kernel; fixes the broken logic noted in Bug in fmpq_mat_can_solve_multi_mod #2705 and generally improves performance
• Certifying that the system is singular is done by checking a single vector, not by effectively computing the RREF to certify the rank; this is faster for systems of low rank. Also adds an optimized fmpz_mat_is_singular (much faster than computing the rank or determinant) as a new function.
• Certifying that the can_solve system is inconsistent is also done by checking a single vector (the Fredholm alternative), not trying primes up to the Hadamard bound
• The fmpq_mat functions become wrappers around fmpz_mat algorithms and not vice versa
• Adds the global flag flint_fmpz_mat_force_small_primes to allow testing all modular fmpz_mat functions with starting prime p=2 instead of the default p=nextprime(2^59) or so (catches bugs like Bug in fmpq_mat_can_solve_multi_mod #2705)
• Some more loops are multithreaded

Odds and ends:

• I think the multimodular/Dixon solvers now guarantee a minimal positive denominator. In the future, we may want all solvers (including the default one) to guarantee this; need to profile first whether extra GCD computations are a performance issue for small n.
• The Dixon code uses a lot of logic to set its internal parameters by estimating the cost of various subroutines for the given input sizes. This has been calibrated by testing lots of balanced and unbalanced systems and fixing any regressions against the old code. I'm sure this logic could be simplified, but this is not easy.

Timings

n x n system; n x cols solution and right hand side

bits(A) = bits of the system matrix A
bits(B) = bits of the RHS
bits(X) = bits of solution numerator
bits(d) = bits of solution denominator

old fmpz_mat_solve vs new, 1T = single thread; 8T = multithreaded

n x 1, generic A and B:

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20   1 |       10      10     220     218 |  0.00103  0.00064 (  1.622x)   0.00041  0.00056 (  0.735x)
  40   1 |       10      10     460     459 |  0.00103  0.00083 (  1.248x)   0.00103  0.00171 (  0.598x)
  80   1 |       10      10     967     965 |  0.00449  0.00573 (  0.783x)   0.00461  0.00469 (  0.983x)
 160   1 |       10      10    2006    2005 |  0.02803  0.01280 (  2.190x)   0.01271  0.01104 (  1.151x)
 320   1 |       10      10    4180    4179 |  0.10650  0.07382 (  1.443x)   0.08711  0.06132 (  1.421x)
 640   1 |       10      10    8676    8673 |  0.82775  0.52926 (  1.564x)   0.65145  0.40112 (  1.624x)
1280   1 |       10      10   17985   17983 |  7.03068  4.61658 (  1.523x)   5.50908  2.69299 (  2.046x)
2560   1 |       10      10   37248   37246 | 54.73336 34.15226 (  1.603x)  43.69410 20.35498 (  2.147x)
  20   1 |      100     100    2020    2017 |  0.00144  0.00085 (  1.683x)   0.00128  0.00079 (  1.629x)
  40   1 |      100     100    4066    4065 |  0.00506  0.00294 (  1.719x)   0.00510  0.00293 (  1.738x)
  80   1 |      100     100    8168    8167 |  0.02864  0.01358 (  2.109x)   0.02644  0.01354 (  1.953x)
 160   1 |      100     100   16405   16404 |  0.23493  0.07378 (  3.184x)   0.19578  0.08267 (  2.368x)
 320   1 |      100     100   32979   32976 |  1.75300  0.46862 (  3.741x)   1.24135  0.44186 (  2.809x)
 640   1 |      100     100   66279   66277 | 14.54457  3.04329 (  4.779x)  10.92958  3.04483 (  3.590x)
  20   1 |     1000    1000   20023   20020 |  0.01846  0.01227 (  1.504x)   0.01823  0.01213 (  1.503x)
  40   1 |     1000    1000   40067   40066 |  0.11781  0.07261 (  1.623x)   0.11889  0.07272 (  1.635x)
  80   1 |     1000    1000   80169   80165 |  0.72406  0.45147 (  1.604x)   0.72059  0.43735 (  1.648x)
 160   1 |     1000    1000  160406  160403 |  5.06811  2.94102 (  1.723x)   5.02011  2.98940 (  1.679x)
 320   1 |     1000    1000  320978  320975 | 47.91358 23.12362 (  2.072x)  48.79206 23.61065 (  2.067x)
  20   1 |    10000   10000  200020  200018 |  0.90343  0.41354 (  2.185x)   0.90537  0.39631 (  2.285x)
  40   1 |    10000   10000  400062  400059 |  5.09594  1.93877 (  2.628x)   5.25380  1.60430 (  3.275x)
  80   1 |    10000   10000  800165  800163 | 47.13387  9.97633 (  4.725x)  46.74485  5.81098 (  8.044x)
  20   1 |   100000  100000 2000021 2000019 | 65.76680  5.76757 ( 11.403x)  63.72588  5.04621 ( 12.628x)

n x n RHS, generic A and B:

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20  20 |       10      10     224     223 |  0.00277  0.00114 (  2.432x)   0.00243  0.00156 (  1.556x)
  40  40 |       10      10     465     463 |  0.01058  0.00562 (  1.883x)   0.00934  0.00521 (  1.791x)
  80  80 |       10      10     968     965 |  0.09065  0.05389 (  1.682x)   0.07784  0.02339 (  3.327x)
 160 160 |       10      10    2013    2009 |  0.91407  0.53771 (  1.700x)   0.80475  0.21515 (  3.740x)
 320 320 |       10      10    4181    4176 | 12.31040  5.91530 (  2.081x)  10.34404  1.77520 (  5.827x)
  20  20 |      100     100    2023    2021 |  0.01139  0.00715 (  1.593x)   0.01234  0.00621 (  1.988x)
  40  40 |      100     100    4065    4061 |  0.12648  0.04971 (  2.544x)   0.12982  0.03453 (  3.759x)
  80  80 |      100     100    8168    8164 |  1.57117  0.48838 (  3.217x)   1.42706  0.25554 (  5.585x)
 160 160 |      100     100   16411   16408 | 18.97795  5.54495 (  3.423x)  18.73801  2.07305 (  9.039x)
  20  20 |     1000    1000   20024   20020 |  0.26848  0.10715 (  2.506x)   0.27815  0.08308 (  3.348x)
  40  40 |     1000    1000   40067   40065 |  4.24528  0.88322 (  4.807x)   4.08469  0.47626 (  8.577x)
  80  80 |     1000    1000   80163   80156 | 65.37644  8.22718 (  7.946x)  66.73743  3.72099 ( 17.935x)
  20  20 |    10000   10000  200024  200022 | 16.64655  2.68531 (  6.199x)  17.22796  1.55916 ( 11.049x)

n x 1, small determinant:

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20   1 |       96     100     912      10 |  0.00049  0.00044 (  1.126x)   0.00028  0.00036 (  0.769x)
  40   1 |       99     100    1830      10 |  0.00149  0.00127 (  1.177x)   0.00135  0.00122 (  1.110x)
  80   1 |      101     100    3528       9 |  0.00970  0.00432 (  2.243x)   0.00872  0.00432 (  2.022x)
 160   1 |      102     100    6826      10 |  0.08223  0.02149 (  3.827x)   0.06166  0.01971 (  3.128x)
 320   1 |      104     100   13701       5 |  0.62074  0.11417 (  5.437x)   0.43828  0.10274 (  4.266x)
 640   1 |      105     100   27544       9 |  7.75616  0.91027 (  8.521x)   5.51922  0.77342 (  7.136x)
1280   1 |      106     100   56921      10 | 68.53736  7.82182 (  8.762x)  47.43033  4.26698 ( 11.116x)
  20   1 |      968    1000   11229      95 |  0.00971  0.00386 (  2.516x)   0.00940  0.00386 (  2.434x)
  40   1 |      996    1000   22084     100 |  0.05177  0.02444 (  2.118x)   0.05087  0.02401 (  2.119x)
  80   1 |      997    1000   41596      99 |  0.28447  0.14169 (  2.008x)   0.28478  0.14008 (  2.033x)
 160   1 |     1001    1000   84753     100 |  3.09494  0.95828 (  3.230x)   3.11044  0.95597 (  3.254x)
 320   1 |     1001    1000  168213     100 | 31.27370  6.48781 (  4.820x)  31.36860  6.65025 (  4.717x)
  20   1 |     9850   10000  109103     997 |  0.33029  0.32717 (  1.010x)   0.31754  0.31705 (  1.002x)
  40   1 |     9930   10000  211720    1000 |  2.26403  1.34591 (  1.682x)   2.29024  1.21394 (  1.887x)
  80   1 |     9943   10000  437680    1000 | 15.44322  6.58073 (  2.347x)  16.36928  3.46972 (  4.718x)
  20   1 |    97441  100000 1227500    9998 | 25.09662  3.74995 (  6.693x)  26.43390  3.52797 (  7.493x)

n x n, small determinant:

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20  20 |       96     100     912      10 |  0.00777  0.00486 (  1.600x)   0.00508  0.00359 (  1.416x)
  40  40 |       99     100    1734      10 |  0.04502  0.01126 (  4.000x)   0.02280  0.00860 (  2.653x)
  80  80 |      101     100    3526      10 |  0.32361  0.06034 (  5.363x)   0.27439  0.02649 ( 10.356x)
 160 160 |      102     100    7196      10 |  4.85117  0.83238 (  5.828x)   4.33364  0.25891 ( 16.738x)
 320 320 |      103     100   13754      10 | 59.18597  8.66243 (  6.832x)  48.96883  1.73348 ( 28.249x)
  20  20 |      970    1000   11596     100 |  0.11373  0.02900 (  3.922x)   0.11974  0.02010 (  5.956x)
  40  40 |      986    1000   20693     100 |  0.94684  0.12392 (  7.640x)   0.93912  0.05447 ( 17.241x)
  80  80 |      997    1000   42608     100 | 28.93984  1.25332 ( 23.091x)  30.11140  0.28372 (106.129x)
  20  20 |     9839   10000  112506    1000 |  4.84765  0.50858 (  9.532x)   5.33179  0.26248 ( 20.313x)
  40  40 |     9907   10000  222181    1000 | 89.43461  3.65090 ( 24.497x)  90.06951  1.30027 ( 69.270x)
  20  20 |   100000  100000 2000025 2000022 |1452.32931 95.23054 ( 15.251x) 1522.16664 38.46118 ( 39.577x)

n x 1, generic A, generic X, B = A*X:

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20   1 |       10      23      10       1 |  0.00082  0.00016 (  5.199x)   0.00014  0.00007 (  1.949x)
  40   1 |       10      24      10       1 |  0.00056  0.00018 (  3.031x)   0.00018  0.00009 (  1.883x)
  80   1 |       10      24      10       1 |  0.00146  0.00085 (  1.717x)   0.00137  0.00038 (  3.641x)
 160   1 |       10      25      10       1 |  0.00894  0.00350 (  2.552x)   0.00455  0.00227 (  2.002x)
 320   1 |       10      25      10       1 |  0.03810  0.00778 (  4.898x)   0.01592  0.00675 (  2.358x)
 640   1 |       10      26      10       1 |  0.20787  0.03977 (  5.227x)   0.08127  0.03031 (  2.681x)
1280   1 |       10      27      10       1 |  1.41989  0.23946 (  5.930x)   0.43965  0.13331 (  3.298x)
2560   1 |       10      27      10       1 |  9.42601  1.69375 (  5.565x)   2.91836  0.80334 (  3.633x)
5120   1 |       10      28      10       1 | 63.44343 11.45415 (  5.539x)  18.00183  4.42136 (  4.072x)
  20   1 |      100     203     100       1 |  0.00020  0.00013 (  1.499x)   0.00010  0.00009 (  1.080x)
  40   1 |      100     203     100       1 |  0.00033  0.00025 (  1.292x)   0.00027  0.00024 (  1.138x)
  80   1 |      100     204     100       1 |  0.00134  0.00080 (  1.674x)   0.00147  0.00083 (  1.773x)
 160   1 |      100     205     100       1 |  0.00750  0.00421 (  1.783x)   0.00577  0.00367 (  1.571x)
 320   1 |      100     206     100       1 |  0.04390  0.01556 (  2.822x)   0.02406  0.01434 (  1.678x)
 640   1 |      100     206     100       1 |  0.28657  0.09678 (  2.961x)   0.14388  0.08673 (  1.659x)
1280   1 |      100     207     100       1 |  1.87550  0.63311 (  2.962x)   0.80944  0.47886 (  1.690x)
2560   1 |      100     207     100       1 | 11.63600  3.40288 (  3.419x)   4.67591  2.44869 (  1.910x)
  20   1 |     1000    2003    1000       1 |  0.00106  0.00074 (  1.442x)   0.00076  0.00066 (  1.158x)
  40   1 |     1000    2004    1000       1 |  0.00269  0.00177 (  1.518x)   0.00269  0.00187 (  1.445x)
  80   1 |     1000    2004    1000       1 |  0.01099  0.00677 (  1.623x)   0.01075  0.00643 (  1.673x)
 160   1 |     1000    2005    1000       1 |  0.04678  0.02384 (  1.962x)   0.04532  0.02471 (  1.834x)
 320   1 |     1000    2006    1000       1 |  0.27438  0.10276 (  2.670x)   0.23718  0.10372 (  2.287x)
 640   1 |     1000    2006    1000       1 |  1.29518  0.43254 (  2.994x)   1.15453  0.41624 (  2.774x)
1280   1 |     1000    2007    1000       1 |  6.56882  1.83656 (  3.577x)   5.44730  1.66148 (  3.279x)
2560   1 |     1000    2007    1000       1 | 30.18662  7.58505 (  3.980x)  23.84664  6.79392 (  3.510x)
  20   1 |    10000   20003   10000       1 |  0.02795  0.01766 (  1.582x)   0.02713  0.01754 (  1.547x)
  40   1 |    10000   20004   10000       1 |  0.09568  0.04203 (  2.276x)   0.09529  0.04251 (  2.241x)
  80   1 |    10000   20004   10000       1 |  0.35571  0.13817 (  2.574x)   0.34446  0.14097 (  2.443x)
 160   1 |    10000   20005   10000       1 |  2.48280  0.71611 (  3.467x)   2.49610  0.71937 (  3.470x)
 320   1 |    10000   20006   10000       1 | 10.81882  3.34673 (  3.233x)  11.16567  3.48687 (  3.202x)
  20   1 |   100000  200003  100000       1 |  2.16060  2.65309 (  0.814x)   2.13498  2.25072 (  0.949x)
  40   1 |   100000  200004  100000       1 |  9.39496  3.28451 (  2.860x)   9.34413  3.01681 (  3.097x)
  80   1 |   100000  200004  100000       1 | 35.77884 11.59471 (  3.086x)  36.02857 11.56898 (  3.114x)

rank n / 2 (detecting singular system)

  n cols |  bits(A) bits(B) bits(X) bits(d) |  old (1T) new(1T) ( speedup)  old (8T) new (8T) ( speedup)
--------------------------------------------------------------------------------------------------------
  20   1 |       13      10       0       0 |  0.00042  0.00004 ( 11.447x)   0.00009  0.00002 (  3.632x)
  40   1 |       13      10       0       0 |  0.00068  0.00020 (  3.325x)   0.00057  0.00019 (  3.036x)
  80   1 |       14      10       0       0 |  0.00475  0.00079 (  6.041x)   0.00447  0.00081 (  5.538x)
 160   1 |       15      10       0       0 |  0.04090  0.00335 ( 12.198x)   0.01582  0.00317 (  4.994x)
 320   1 |       15      10       0       0 |  0.46259  0.01593 ( 29.034x)   0.18946  0.01841 ( 10.294x)
 640   1 |       16      10       0       0 |  5.30634  0.09378 ( 56.581x)   1.46231  0.07434 ( 19.671x)
1280   1 |       17      10       0       0 | 65.11364  0.66078 ( 98.541x)  17.87687  0.45692 ( 39.124x)
  20   1 |      103     100       0       0 |  0.00093  0.00025 (  3.730x)   0.00043  0.00010 (  4.477x)
  40   1 |      104     100       0       0 |  0.00435  0.00073 (  5.942x)   0.00569  0.00189 (  3.012x)
  80   1 |      104     100       0       0 |  0.05263  0.00300 ( 17.545x)   0.02294  0.00311 (  7.381x)
 160   1 |      105     100       0       0 |  0.37787  0.01713 ( 22.056x)   0.14739  0.02225 (  6.624x)
 320   1 |      106     100       0       0 |  4.28961  0.07300 ( 58.764x)   1.51953  0.07170 ( 21.192x)
 640   1 |      106     100       0       0 | 48.09464  0.38822 (123.885x)  14.83424  0.40610 ( 36.529x)
  20   1 |     1003    1000       0       0 |  0.00737  0.00141 (  5.243x)   0.00759  0.00148 (  5.128x)
  40   1 |     1004    1000       0       0 |  0.06924  0.01491 (  4.643x)   0.05406  0.01503 (  3.595x)
  80   1 |     1004    1000       0       0 |  0.63615  0.07165 (  8.879x)   0.43693  0.07355 (  5.941x)
 160   1 |     1005    1000       0       0 |  6.47013  0.42994 ( 15.049x)   2.57732  0.43489 (  5.926x)

fmpz_mat_can_solve: A = 2k x k with rank k / 2, B = A * random (consistent)

    m    k    r cols | bits(A) bits(X) bits(d) s |  old (1T) new (1T) ( speedup)  old (8T) new (8T) ( speedup)
------------------------------------------------------------------------------------------------------------
   40   20   10    1 |      12      18       1 1 |  0.00035  0.00004 (  7.689x)   0.00006  0.00004 (  1.579x)
   80   40   20    1 |      14      32       1 1 |  0.00019  0.00017 (  1.093x)   0.00016  0.00012 (  1.269x)
  160   80   40    1 |      16      59       1 1 |  0.00112  0.00062 (  1.816x)   0.00114  0.00059 (  1.953x)
  320  160   80    1 |      17     103       1 1 |  0.00896  0.00339 (  2.645x)   0.00624  0.00265 (  2.357x)
  640  320  160    1 |      19     173       1 1 |  0.08318  0.01871 (  4.446x)   0.04844  0.01250 (  3.874x)
 1280  640  320    1 |      20     327       1 1 |  1.19549  0.11441 ( 10.450x)   0.48687  0.05994 (  8.123x)
   40   20   10    1 |      92     161       1 1 |  0.00021  0.00010 (  2.018x)   0.00016  0.00008 (  2.102x)
   80   40   20    1 |      98     207       1 1 |  0.00086  0.00048 (  1.799x)   0.00082  0.00039 (  2.130x)
  160   80   40    1 |     100     485       1 1 |  0.00781  0.00168 (  4.642x)   0.00771  0.00158 (  4.897x)
  320  160   80    1 |     101     904       1 1 |  0.08581  0.00683 ( 12.557x)   0.06522  0.00569 ( 11.452x)
  640  320  160    1 |     103    1725       1 1 |  1.10340  0.03470 ( 31.797x)   0.72505  0.02741 ( 26.450x)
   40   20   10    1 |     981    1398       1 1 |  0.00218  0.00047 (  4.633x)   0.00209  0.00047 (  4.462x)
   80   40   20    1 |     984    2760       1 1 |  0.01522  0.00155 (  9.800x)   0.01490  0.00142 ( 10.516x)
  160   80   40    1 |     996    5731       1 1 |  0.19993  0.00847 ( 23.598x)   0.19992  0.00830 ( 24.085x)
  320  160   80    1 |    1000   10937       1 1 |  1.42628  0.04619 ( 30.877x)   1.27982  0.04995 ( 25.620x)
   40   20   10    1 |    9447   13862       1 1 |  0.06240  0.01192 (  5.233x)   0.06186  0.01194 (  5.179x)
   80   40   20    1 |    9790   28414       1 1 |  0.48931  0.03782 ( 12.937x)   0.48923  0.03762 ( 13.006x)

fmpz_mat_can_solve, A = 2k x k with rank k / 2, B = random (certifying inconsistent) system

    m    k    r cols | bits(A) bits(X) bits(d) s |  old (1T) new (1T) ( speedup)  old (8T) new (8T) ( speedup)
------------------------------------------------------------------------------------------------------------
   40   20   10    1 |      12       0       0 0 |  0.00080  0.00021 (  3.843x)   0.00021  0.00011 (  1.855x)
   80   40   20    1 |      15       0       0 0 |  0.00149  0.00069 (  2.148x)   0.00166  0.00053 (  3.158x)
  160   80   40    1 |      16       0       0 0 |  0.01522  0.00236 (  6.437x)   0.01687  0.00090 ( 18.842x)
  320  160   80    1 |      17       0       0 0 |  0.10072  0.00426 ( 23.636x)   0.07023  0.00360 ( 19.527x)
  640  320  160    1 |      19       0       0 0 |  1.46347  0.03076 ( 47.580x)   0.77032  0.02352 ( 32.752x)
   40   20   10    1 |      95       0       0 0 |  0.00049  0.00014 (  3.434x)   0.00048  0.00012 (  3.972x)
   80   40   20    1 |      99       0       0 0 |  0.00478  0.00096 (  4.988x)   0.00483  0.00086 (  5.630x)
  160   80   40    1 |     100       0       0 0 |  0.05314  0.00288 ( 18.444x)   0.05248  0.00267 ( 19.640x)
  320  160   80    1 |     101       0       0 0 |  0.64543  0.01241 ( 52.010x)   0.48111  0.01164 ( 41.327x)
   40   20   10    1 |     987       0       0 0 |  0.01028  0.00116 (  8.834x)   0.01083  0.00105 ( 10.275x)
   80   40   20    1 |     979       0       0 0 |  0.08982  0.00674 ( 13.322x)   0.08694  0.00632 ( 13.759x)
  160   80   40    1 |     998       0       0 0 |  0.88347  0.03250 ( 27.181x)   0.88358  0.03403 ( 25.965x)
   40   20   10    1 |    9803       0       0 0 |  0.29096  0.03107 (  9.366x)   0.29041  0.03106 (  9.350x)
   80   40   20    1 |    9838       0       0 0 |  2.76369  0.34701 (  7.964x)   2.78548  0.33526 (  8.309x)