flint_mpn division, square roots and powers

[fredrik-johansson]

Adds flint_mpn_* functions for division, square root, etc. and uses these in fmpz, arf and elsewhere. GMP is used for small sizes, Newton iteration based on flint_mpn_mul* for huge input (replacing old inefficient arf-based Newton code).

Implemented mostly using Claude Fable 5.1.

@albinahlback The interfaces are now in place for inserting faster basecase division code :-)

In passing this fixes #2317: fmpz_ndiv_qr deliberately changes tiebreaking to even. This breaks a Nemo test, which has to be updated downstream.

New public functions

mpn_extras — Hensel (2-adic) arithmetic

Function	Purpose
flint_mpn_binv	x^-1 mod B^n for odd x
flint_mpn_bdiv_qr	Hensel division, n quotient limbs, optional remainder
flint_mpn_bdiv_q	inline wrapper, quotient only
flint_mpn_bdiv_qr_1	single-limb divisor
flint_mpn_bdiv_qr_classical	block algorithm (short divisors)
flint_mpn_bdiv_qr_karp_markstein	Karp–Markstein algorithm
flint_mpn_brsqrt	2-adic reciprocal square root
flint_mpn_bsqrt	2-adic square root
_flint_mpn_mulhigh_known_low	high product given known low limbs
_flint_mpn_bdiv_qr_classical_preinv	block algorithm with a precomputed inverse

mpn_extras — Euclidean division

Function	Purpose
flint_mpn_tdiv_qr, flint_mpn_tdiv_q, flint_mpn_tdiv_r	truncating division (inline dispatch)
flint_mpn_cdiv_qr, flint_mpn_cdiv_q, flint_mpn_cdiv_r	ceiling division
flint_mpn_ndiv_qr, flint_mpn_ndiv_q, flint_mpn_ndiv_r	round to nearest, ties to even
flint_mpn_div	checked exact division
flint_mpn_divexact	exact division (assumed exact)
flint_mpn_divisible	divisibility test
flint_mpn_inv	correctly truncated floor(B^n / x)
_flint_mpn_tdiv_qr and _newton, _unbalanced, _preinv, _preinvn, _gmp	individual algorithms
_flint_mpn_divexact, _flint_mpn_divexact_hensel, _flint_mpn_divisible	out-of-line backends
flint_mpn_divexact_preinv_init / _clear / flint_mpn_divexact_preinv	exact division by a fixed divisor (+ type flint_mpn_divexact_preinv_t)

mpn_extras — square root, powering, misc

Function	Purpose
flint_mpn_sqrtrem	root and optional remainder; with r == NULL returns 0 iff a perfect square
flint_mpn_sqrt	checked exact square root
flint_mpn_is_square	perfect square test
_flint_mpn_sqrtrem, _flint_mpn_sqrtrem_newton, _flint_mpn_sqrtrem_gmp, _flint_mpn_sqrt, _flint_mpn_is_square	backends
flint_mpn_pow, flint_mpn_pow_bound_limbs	binary powering and its size bound
flint_mpn_mod_2exp48m1	value congruent mod 2^48 - 1 (AVX2 kernel; 64-bit only)

mpn_extras — mpz-like interface

Intermediate wrappers used by the fmpz layer; if there is interest we could promote these to an mpz_extras module in the future.

flint_mpz_tdiv_qr, _tdiv_q, _tdiv_r, flint_mpz_fdiv_qr, _fdiv_q, _fdiv_r,
flint_mpz_cdiv_qr, _cdiv_q, _cdiv_r, flint_mpz_mod, flint_mpz_divexact,
flint_mpz_sqrt, flint_mpz_sqrtrem (inline dispatch; _flint_mpz_* backends).

fmpz

Function	Purpose
fmpz_div, fmpz_div_ui, fmpz_div_si	checked exact division (fmpz_divides is now an alias of fmpz_div)
fmpz_perfect_sqrt	checked exact square root
fmpz_invmod_2exp	inverse modulo 2^N
fmpz_sqrtmod_2exp, fmpz_rsqrtmod_2exp	2-adic square root and reciprocal square root
fmpz_divmod_2exp	Hensel division modulo 2^N

arf

Function	Purpose
ARF_RND_FAST, ARF_RND_ACCURATE	relaxed rounding modes (≤ 1 ulp, ≤ 0.51 ulp)
arf_rnd_relaxed_to_strict	mode mapping helper
_arf_inv_newton, _arf_div_newton, _arf_sqrt_newton, _arf_rsqrt_newton	Newton backends
_arf_want_newton_inv/div/sqrt/rsqrt (inline) and _large variants	dispatch predicates

Changed behaviour of existing functions

Function	Change
fmpz_ndiv_qr	ties now round to even (was: towards zero); documented, test updated
fmpz_divides	now an inline alias of fmpz_div
fmpz_tdiv_q/qr, fdiv_q/qr/r, cdiv_q/qr, fmpz_mod, fmpz_divexact	Newton/Hensel for large operands; 1×1 and 2×1 hardware paths
fmpz_tdiv_q_si/ui, fdiv_q_si/ui, cdiv_q_si/ui, divexact_si/ui, tdiv_ui, fdiv_ui, cdiv_ui	1×1 / 2×1 paths without mpz promotion
fmpz_sqrt, fmpz_sqrtrem	Newton above the cutoff; 1- and 2-limb paths
fmpz_is_square	residue screening plus Newton certification
fmpz_pow_ui	uses flint_mpn_pow
arf_div, arf_ui_div, arf_sqrt, arf_rsqrt (and arf_div_ui/si/fmpz, arf_fmpz_div_fmpz, …)	accept the relaxed rounding modes and dispatch to Newton
arf_set_round / _arf_set_round_mpn	in-place rounding supported; no temporary for aliased arguments
arf_div	exact quotient via flint_mpn_tdiv_q instead of __gmpn_div_q
arb_div, arb_div_arf, arb_div_fmpz, arb_sqrt, arb_sqrt_arf, arb_rsqrt, arb_rsqrt_arf	use ARF_RND_FAST; all Newton machinery removed
flint_mpn_preinvn	uses flint_mpn_inv
flint_mpn_divides	uses _flint_mpn_tdiv_qr
mpn_tdiv_q (compat shim)	calls GMP's mpn_div_q
_fmpz_vec_scalar_divexact_fmpz, fmpz_mat_scalar_divexact_fmpz	precomputed 2-adic inverse for multi-limb divisors
_padic_inv, _padic_sqrt	p = 2 cases use fmpz_invmod_2exp / fmpz_sqrtmod_2exp
fmpz_poly_sqrt_KS, fmpz_mpoly_sqrt_heap, mpn_extras/multi_mod, multi_crt, multi_crt_once	routed through the new mpn functions
gr fmpz ring	sqrt via fmpz_perfect_sqrt, div via fmpz_div, new div_ui / div_si methods
radix_sqrtrem_newton_karp_markstein	drops an unnecessary zeroing

Removed

arb_div_newton, arb_div_arf_newton, _arb_fmpz_divapprox_newton,
arb_fmpz_divapprox, arb_sqrt_newton, arb_sqrt_arf_newton,
arb_rsqrt_arf_newton, the _fmpz_*_newton hooks and
MPZ_WANT_FLINT_DIVISION (plus tests t-div_newton.c, t-sqrt_newton.c
in arb, t-div_newton.c in fmpz).

Build system

• configure.ac probes __gmpn_divexact, __gmpn_mod_34lsub1,
  __gmpn_divisible_p (FLINT_HAVE_NATIVE_mpn_*).
• All 12 flint-mparam.h files define FLINT_MPN_TDIV_QR_NEWTON_CUTOFF,
  FLINT_MPN_DIVEXACT_NEWTON_CUTOFF, FLINT_MPN_SQRTREM_NEWTON_CUTOFF.

Performance

This mostly speeds up arithmetic with huge numbers.

I have tried very hard to ensure that there are no major regressions for small or unbalanced inputs. Unfortunately, due to the number of functions touched, it is possible that I've missed some spots. There are also some improvements for small input. For example:

• flint_mpn_sqrtrem has 1- and 2-limb code which beats GMP
• The fmpz division functions for 2-limb by 1-limb division and square root functions for 2-limb input have been optimized: they now compute the result in registers and write a small result without an unnecessary promote-demote cycle

Known regressions:

• arf_div/arb_div and arf_sqrt/arb_sqrt are 3-4 nanoseconds slower than before, which is significant at 64-256 bits of precision (~5-10% overhead). I couldn't figure out how to fix this easily. However, the code is suboptimal here to begin with here, and I think a future rewrite of the division for low precision will give a speedup that is much bigger than the current regression.

• Ceiling division is slightly slower than before at high precision as it doesn't use the quotient-only path

Benchmarks: fmpz vs mpz

Input: a random with 2n bits, b random with n bits, b2 = b^2. Each entry is the ratio time(mpz) / time(fmpz), so higher is better and values below 1.00 mean the fmpz function is slower than its GMP counterpart. The speedup columns are new/old, i.e. the improvement this PR gives over the current main branch.

fmpz_tdiv_q — truncating quotient

tdiv_q(q, a, b) for the inexact case, tdiv_q(q, b2, b) for the exact case (remainder 0).

bits	old, inexact	new, inexact	speedup	old, exact	new, exact	speedup
30	5.14	5.18	1.01	5.18	5.19	1.00
32	0.56	4.46	7.96	0.54	4.34	8.04
60	0.57	3.27	5.74	0.56	3.25	5.80
64	0.81	0.80	0.99	0.79	0.85	1.08
100	0.88	0.87	0.99	0.88	0.87	0.99
240	0.91	0.92	1.01	0.92	0.92	1.00
1,000	0.98	0.97	0.99	0.98	0.98	1.00
2,000	0.99	0.99	1.00	0.99	0.99	1.00
4,000	1.00	1.00	1.00	1.00	1.00	1.00
10,000	1.00	1.00	1.00	1.00	0.99	0.99
20,000	1.00	1.00	1.00	1.00	1.00	1.00
40,000	1.00	1.20	1.20	1.01	1.26	1.25
100,000	1.11	1.63	1.47	1.29	1.71	1.33
200,000	1.47	1.90	1.29	1.72	2.04	1.19
400,000	1.84	2.08	1.13	2.03	2.32	1.14
1,000,000	2.09	2.45	1.17	2.11	2.35	1.11
10,000,000	2.23	2.78	1.25	2.35	2.72	1.16

fmpz_cdiv_q — ceiling quotient

bits	old, inexact	new, inexact	speedup	old, exact	new, exact	speedup
30	8.64	8.60	1.00	7.27	7.27	1.00
32	0.86	6.32	7.35	0.77	5.78	7.51
60	0.86	4.95	5.76	0.79	4.55	5.76
64	0.83	0.82	0.99	0.88	0.87	0.99
100	0.88	0.86	0.98	0.91	0.90	0.99
240	0.93	0.91	0.98	0.98	0.94	0.96
1,000	0.98	0.98	1.00	0.99	0.98	0.99
2,000	0.99	0.99	1.00	1.00	0.99	0.99
4,000	1.00	1.00	1.00	1.00	1.00	1.00
10,000	1.00	1.00	1.00	0.99	1.00	1.01
20,000	1.00	1.00	1.00	1.00	1.00	1.00
40,000	1.00	1.11	1.11	1.00	1.05	1.05
100,000	1.43	1.47	1.03	1.08	1.42	1.31
200,000	1.74	1.60	0.92	1.30	1.54	1.18
400,000	2.10	1.84	0.88	1.56	1.79	1.15
1,000,000	2.38	2.02	0.85	1.77	1.97	1.11
10,000,000	2.66	2.35	0.88	1.93	2.21	1.15

fmpz_tdiv_qr — truncating quotient and remainder

bits	old, inexact	new, inexact	speedup	old, exact	new, exact	speedup
30	5.02	4.61	0.92	4.98	4.67	0.94
32	0.34	3.79	11.15	0.34	3.82	11.24
60	0.36	3.10	8.61	0.36	3.11	8.64
64	0.43	0.44	1.02	0.43	0.44	1.02
100	0.79	0.78	0.99	0.54	0.56	1.04
240	0.90	0.90	1.00	0.74	0.75	1.01
1,000	0.98	0.97	0.99	0.92	0.93	1.01
2,000	0.97	0.99	1.02	0.95	0.96	1.01
4,000	1.00	1.00	1.00	0.99	0.99	1.00
10,000	1.01	1.01	1.00	0.99	0.99	1.00
20,000	1.00	1.00	1.00	0.99	1.00	1.01
40,000	1.00	1.11	1.11	1.00	1.05	1.05
100,000	1.08	1.48	1.37	1.08	1.41	1.31
200,000	1.30	1.60	1.23	1.31	1.55	1.18
400,000	1.55	1.85	1.19	1.57	1.79	1.14
1,000,000	1.76	2.04	1.16	1.77	1.97	1.11
10,000,000	1.91	2.34	1.23	1.93	2.21	1.15

fmpz_divexact — exact division

Only the exact case exists: divexact(q, b2, b).

bits	old	new	speedup
30	2.89	2.90	1.00
32	0.37	2.46	6.65
60	0.42	2.07	4.93
64	0.73	0.72	0.99
100	0.79	0.84	1.06
240	0.89	0.86	0.97
1,000	0.96	0.97	1.01
2,000	1.00	1.00	1.00
4,000	0.99	1.00	1.01
10,000	1.01	1.00	0.99
20,000	1.00	1.00	1.00
40,000	1.00	1.00	1.00
100,000	1.05	1.57	1.50
200,000	1.34	1.71	1.28
400,000	1.64	1.97	1.20
1,000,000	1.96	2.27	1.16
10,000,000	2.06	2.68	1.30

fmpz_divisible — divisibility test

divisible(a, b) for the negative case, divisible(b2, b) for the positive case.

bits	old, not divisible	new, not divisible	speedup	old, divisible	new, divisible	speedup
30	2.27	2.23	0.98	2.58	2.47	0.96
32	1.06	1.03	0.97	1.13	1.03	0.91
60	1.18	1.13	0.96	1.11	1.04	0.94
64	0.88	0.82	0.93	0.88	0.84	0.95
100	0.94	0.89	0.95	0.93	0.92	0.99
240	0.97	0.93	0.96	0.95	0.94	0.99
1,000	0.99	1.59	1.61	0.99	0.82	0.83
2,000	1.00	1.68	1.68	1.00	0.92	0.92
4,000	1.00	2.11	2.11	1.01	0.97	0.96
10,000	1.00	2.04	2.04	1.00	0.99	0.99
20,000	1.00	2.16	2.16	1.00	0.99	0.99
40,000	1.00	1.70	1.70	1.00	1.00	1.00
100,000	1.00	1.95	1.95	1.00	1.00	1.00
200,000	1.00	2.92	2.92	1.00	1.54	1.54
400,000	1.04	3.51	3.37	1.01	1.77	1.75
1,000,000	1.00	3.33	3.33	1.00	2.01	2.01
10,000,000	1.28	4.88	3.81	1.00	2.08	2.08

fmpz_sqrtrem — square root with remainder

sqrtrem(s, r, a) for a generic (non-square) input, sqrtrem(s, r, b2) for a perfect square.

bits	old, non-square	new, non-square	speedup	old, square	new, square	speedup
30	1.60	1.40	0.87	1.64	1.55	0.95
32	0.20	0.78	3.90	0.20	0.79	3.95
60	0.30	0.96	3.20	0.30	0.99	3.30
64	0.68	0.60	0.88	0.37	0.72	1.95
100	0.83	0.85	1.02	0.72	0.70	0.97
240	0.91	0.87	0.96	0.80	0.77	0.96
1,000	0.94	0.96	1.02	0.91	0.91	1.00
2,000	0.97	0.97	1.00	0.96	0.96	1.00
4,000	1.00	1.00	1.00	0.96	0.95	0.99
10,000	1.00	1.00	1.00	0.98	0.98	1.00
20,000	1.00	1.00	1.00	1.00	0.99	0.99
40,000	1.00	1.00	1.00	1.00	1.00	1.00
100,000	1.00	1.11	1.11	1.00	1.06	1.06
200,000	1.00	1.39	1.39	1.00	1.33	1.33
400,000	1.00	1.64	1.64	0.99	1.61	1.63
1,000,000	1.00	2.07	2.07	1.00	1.99	1.99
10,000,000	1.00	2.33	2.33	1.00	2.18	2.18

fmpz_is_square — perfect square test

is_square(a) for a generic (non-square) input, is_square(b2) for a perfect square.

bits	old, non-square	new, non-square	speedup	old, square	new, square	speedup
30	0.91	0.94	1.03	3.42	1.82	0.53
32	0.76	0.83	1.09	0.96	1.47	1.53
60	0.78	0.75	0.96	0.97	1.47	1.52
64	0.79	0.69	0.87	0.94	1.36	1.45
100	0.81	0.88	1.09	0.94	0.96	1.02
240	0.81	0.83	1.02	1.00	0.97	0.97
1,000	0.83	0.83	1.00	0.99	0.97	0.98
2,000	0.93	1.55	1.67	1.00	0.98	0.98
4,000	0.91	1.01	1.11	1.01	0.99	0.98
10,000	1.01	7.75	7.67	1.00	1.00	1.00
20,000	0.90	1.30	1.44	1.00	1.00	1.00
40,000	0.90	1.35	1.50	1.00	1.00	1.00
100,000	0.98	1.42	1.45	1.00	0.92	0.92
200,000	0.95	1.30	1.37	1.00	1.12	1.20
400,000	1.02	1.12	1.10	1.00	1.45	1.45
1,000,000	1.00	78.00	78.00	1.00	1.79	1.79
10,000,000	1.02	1.12	1.10	1.00	1.95	1.95

Note: the speedup here depends on the distribution of inputs. The lone 78x speedup in the table is for a set of input where the FLINT's residue screen happened to reject all random inputs while GMP's didn't.