Radix integer division

[fredrik-johansson]

Implements

• radix_divrem
• radix_divrem_preinv
• radix_divexact
• radix_div (checked exact division)
• radix_inv_approx

and various algorithm variants. Plus wrappers for radix_integer_t:

• radix_integer_div
• radix_integer_divexact
• radix_integer_tdiv_q
• radix_integer_tdiv_r
• radix_integer_tdiv_qr
• radix_integer_fdiv_q
• radix_integer_fdiv_r
• radix_integer_fdiv_qr
• radix_integer_cdiv_q
• radix_integer_cdiv_r
• radix_integer_cdiv_qr

Asymptotically, division is implemented using an approximate Newton reciprocal (radix_inv_approx) with $O(M(n))$ complexity. The code uses fixed-point arithmetic with short and middle products to get a good practical constant factor. Middle product isn't optimized to use the FFT trick yet, but that's a missing optimization in radix_mulmid which is independent of the division code.

Shortcomings of this code:

• Division doesn't use the Karp-Markstein trick (possible ~10-20% speedup).
• The Newton reciprocal doesn't have an error bound, so the result must be validated by computing a remainder. This is fine for division with remainder, but means that approximate, quotient-only or exact division is a full multiplication slower than necessary.
• There is no radix-native basecase division (except for unbalanced division). Instead the basecase algorithm is to convert to binary, divide using mpn_tdiv_qr, and convert back. This solution isn't awful, but it's not necessarily optimal either.
• Other minor optimizations (tightest possible memory management, etc.).

This can be regarded as something of a prototype for implementing similar division code for mpn, a priori error bounds for the Newton reciprocal being the main missing ingredient. Empirically, the relative error is smaller $3 B^{-n}$. It should be possible to prove such a bound either for the current algorithm or for a slightly tweaked algorithm using similar techniques as in https://members.loria.fr/PZimmermann/papers/invert.pdf (note that radix_inv_approx is more general as it doesn't assume an even radix or normalized input). Ideally an error bound should be done for the full Karp-Markstein division, though.

Profile of division with remainder versus our mpn division code (there is actually no flint_mpn_divrem yet; I mocked one up for this profile program using mpn_tdiv_qr for small input and fmpz_tdiv_qr for large input):

$ build/radix/profile/p-divrem 
   decimal      mpn  decimal        time        time    relative
    digits    limbs    limbs flint_mpn_divrem radix_divrem  time

        19        1        1    1.19e-08    5.44e-09    0.457x
        38        2        2    1.71e-08    7.58e-08    4.433x
        57        3        3    4.03e-08    1.23e-07    3.052x
        76        4        4    4.55e-08    1.68e-07    3.692x
       114        6        6    6.49e-08    2.57e-07    3.960x
       171        9        9    1.15e-07     4.3e-07    3.739x
       247       13       13    1.85e-07    7.47e-07    4.038x
       361       19       19    3.15e-07    1.38e-06    4.381x
       532       28       28    5.49e-07    2.81e-06    5.118x
       798       42       42    1.06e-06    5.77e-06    5.443x
      1197       63       63    2.15e-06    7.06e-06    3.284x
      1786       93       94    4.02e-06    1.57e-05    3.905x
      2679      140      141    7.85e-06    2.84e-05    3.618x
      4009      209      211    1.55e-05    4.36e-05    2.813x
      6004      312      316    2.91e-05    6.58e-05    2.261x
      9006      468      474     5.5e-05    9.68e-05    1.760x
     13509      702      711    0.000107    0.000158    1.477x
     20254     1052     1066    0.000189    0.000233    1.233x
     30381     1577     1599    0.000324    0.000347    1.071x
     45562     2365     2398    0.000488    0.000531    1.088x
     68343     3548     3597    0.000765    0.000805    1.052x
    102505     5321     5395     0.00125     0.00122    0.976x
    153748     7981     8092     0.00187     0.00181    0.968x
    230622    11971    12138     0.00277     0.00288    1.040x
    345933    17956    18207     0.00426     0.00446    1.047x
    518890    26934    27310     0.00642      0.0067    1.044x
    778335    40400    40965      0.0101      0.0103    1.020x
   1167493    60599    61447      0.0154      0.0158    1.026x
   1751230    90898    92170       0.024       0.025    1.042x
   2626845   136347   138255      0.0357      0.0405    1.134x
   3940258   204520   207382      0.0581      0.0622    1.071x
   5910387   306780   311073      0.0905      0.0994    1.098x
   8865571   460169   466609       0.143       0.152    1.063x
  13298347   690253   699913       0.222       0.261    1.176x
  19947511  1035379  1049869       0.343       0.431    1.257x
  29921257  1553067  1574803       0.555       0.643    1.159x
  44881876  2329600  2362204       0.878       1.227    1.397x
  67322814  3494400  3543306       1.367       1.777    1.300x
 100984221  5241599  5314959       2.229       2.985    1.339x
 151476322  7862398  7972438       3.535       4.392    1.242x
 227214483 11793597 11958657       5.606       7.384    1.317x
 340821715 17690395 17937985       8.979      11.866    1.322x
 511232563 26535591 26906977      15.772      17.688    1.121x
 766848835 39803386 40360465      24.719      27.809    1.125x
1150273243 59705079 60540697      39.525      42.158    1.067x
1725409855 89557617 90811045      60.889      70.002    1.150x