FLINT3.7.0-dev
/

flint_mpn_sqrtrem

Synopsis#include <flint/mpn_extras.h>
mp_size_t flint_mpn_sqrtrem(mp_ptr s, mp_ptr r, mp_srcptr a, mp_size_t an)source

Sets (s, \lceil \mathit{an}/2 \rceil) to \lfloor \sqrt a \rfloor. If r is not NULL, the remainder a - s^2 \le 2s is written to (r, \lceil \mathit{an}/2 \rceil + 1), zero padded, and its number of limbs is returned; r must have room for \max(\mathit{an}, 2) limbs. If r is NULL, the return value is 0 if a is a perfect square and 1 otherwise (GMP's convention). Requires \mathit{an} \ge 1 and a_{an-1} \ne 0.

On 64-bit machines, inputs of up to four limbs have dedicated code. Two-limb inputs use the hardware double-precision square root for the initial approximation, followed for a \ge 2^{108} by one Newton step and a final adjustment. When the compiler reports a fast fused multiply-add (FP_FAST_FMA), the Newton step is evaluated in floating point: the residual a - D^2 of the integral approximation D = 2^{32} \sqrt{a_1} is computed exactly by an FMA from an exact splitting of a, so that no 128/64-bit division is needed. Three- and four-limb inputs use one fully unrolled step of the divide and conquer square root of Zimmermann [Zim1999], as GMP's mpn_dc_sqrtrem: the two-limb code gives the root of the top two limbs, one 2/1 division gives the remaining limb of the root (for three limbs the remaining 32 bits, using an unbalanced split instead of GMP's zero padding), and the remainder is corrected once without branching.

Larger inputs below FLINT_MPN_SQRTREM_NEWTON_CUTOFF limbs use _flint_mpn_sqrtrem_divconquer on 64-bit machines (GMP's mpn_sqrtrem on 32-bit machines): GMP's divide and conquer algorithm with the code above as basecases, FLINT's squaring and the divide and conquer division of _flint_mpn_divrem_preinv1. When no remainder is wanted, the last step uses an approximate quotient with a guard limb (flint_mpn_divapprox), and the remainder of the division is evaluated (from the low part of a product, flint_mpn_mulmid) only if that is inconclusive, in particular for perfect squares.

Above the cutoff, a is viewed as a fixed-point number in [B^{-2}, 1) and _mp_real_sqrt_newton is used with three extra fraction limbs, so that the truncated root is certified whenever the first fraction limb lies in [2, B-2]; otherwise it is corrected by O(1) steps. With a NULL remainder, the Newton path tests a = s^2 on the low two limbs before comparing in full.

References#

  • [Zim1999] P. Zimmermann. Karatsuba square root. INRIA Research Report RR-3805, 1999.
Cite this page

The FLINT team, “flint_mpn_sqrtrem (mpn_extras.h),” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/mpn_extras/flint_mpn_sqrtrem.html (generated 2026-10-20).