Square roots, n-th roots and perfect powers

fmpz_sqrt and fmpz_sqrtrem compute the integer square root s = \lfloor \sqrt g \rfloor and the remainder g - s^2 \in [0, 2s]. When only squares are of interest, fmpz_perfect_sqrt is the function to use: it decides whether a is a square and returns the root if it is, rejecting almost all non-squares without computing a root. fmpz_is_square only answers the question. fmpz_root computes n-th roots rounded towards zero, and fmpz_is_perfect_power decides whether an integer is a nontrivial power.

Square roots#

Inline values use n_sqrt (a floating-point square root with a correction). Large values use FLINT's own square root code (_flint_mpn_sqrtrem) on 64-bit machines:

• one to four limbs: the hardware square root of a double, a Newton correction and one step of Zimmermann's divide-and-conquer square root, specialised to a two-limb root;
• up to FLINT_MPN_SQRTREM_NEWTON_CUTOFF limbs (7423 on recent x86-64): the divide-and-conquer ("Karatsuba") square root of Zimmermann [Zim1999], as in GMP, which computes the root of the top half recursively and the remaining limbs by one division, at a cost of a small multiple of a multiplication;
• above: Newton iteration for 1/\sqrt{\alpha} in fixed-point arithmetic with the Karp–Markstein trick [KM1997], followed by a check of the last limb that certifies the integer part.

The algorithms are described in Square roots. fmpz_sqrt does not form the remainder, which saves a squaring.

Recognising squares#

fmpz_is_square and fmpz_perfect_sqrt first reject a large input unless it is a square modulo 256 (low limb), then modulo the factors of 2^{48} - 1 (4095 = 9 \cdot 5 \cdot 7 \cdot 13, 17, 257, 97 and 241, from a single reduction modulo 2^{48} - 1), which lets about one random non-square in 900 through. A surviving candidate is decided by computing the root: for inputs below 2^{100}, rounding the floating-point root is exact for squares; otherwise the divide-and-conquer root is computed with an early exit as soon as its low bits show that the input cannot be a square. fmpz_perfect_sqrt returns this root, so for a square it costs one root computation, about half the cost of fmpz_is_square followed by fmpz_sqrt.

n-th roots and perfect powers#

fmpz_root uses n_sqrtrem, n_cbrtrem and n_rootrem for inline values and GMP's mpz_root (Newton iteration) otherwise. The root of a negative number is negative for odd n and an error for even n.

fmpz_is_perfect_power returns an exponent k \ge 2 with f = r^k if there is one. For one-limb values it uses n_is_perfect_power repeatedly, which finds the largest exponent. For larger values it follows the algorithm of GMP's mpz_perfect_power_p: the multiplicities of 2 and of small primes in f restrict the possible exponents to the divisors of their GCD, and the remaining candidates are tested with mpz_root. The exponent found is then not necessarily the largest: 2^{120} is reported as the square of 2^{60}. Apply the function again to the root to obtain a larger exponent.

Example#

Fifty digits of \sqrt 2, checked square roots, a cube root of a negative number and a perfect power:

#include <flint/flint.h>
#include <flint/fmpz.h>

int main(void)
{
    fmpz_t a, s, r;
    int k;

    fmpz_init(a); fmpz_init(s); fmpz_init(r);

    /* 50 digits of sqrt(2): floor(sqrt(2 * 10^98)) */
    fmpz_ui_pow_ui(a, 10, 98);
    fmpz_mul_ui(a, a, 2);
    fmpz_sqrtrem(s, r, a);
    flint_printf("isqrt(2*10^98) = "); fmpz_print(s); flint_printf("\n");
    flint_printf("remainder has %wu bits\n", fmpz_bits(r));

    /* checked square roots */
    fmpz_mul(a, s, s);
    flint_printf("perfect_sqrt(s^2) = %d, ", fmpz_perfect_sqrt(r, a));
    flint_printf("equal to s: %d\n", fmpz_equal(r, s));
    fmpz_add_ui(a, a, 1);
    flint_printf("is_square(s^2 + 1) = %d\n", fmpz_is_square(a));

    /* n-th roots round towards zero */
    fmpz_set_si(a, -1000001);
    k = fmpz_root(r, a, 3);
    flint_printf("root(-1000001, 3) = "); fmpz_print(r);
    flint_printf(", exact: %d\n", k);

    /* perfect powers */
    fmpz_ui_pow_ui(a, 7, 45);
    k = fmpz_is_perfect_power(r, a);
    flint_printf("7^45 = r^%d with r = ", k); fmpz_print(r); flint_printf("\n");

    fmpz_clear(a); fmpz_clear(s); fmpz_clear(r);
    return 0;
}

isqrt(2*10^98) = 14142135623730950488016887242096980785696718753769
remainder has 164 bits
perfect_sqrt(s^2) = 1, equal to s: 1
is_square(s^2 + 1) = 0
root(-1000001, 3) = -100, exact: 0
7^45 = r^45 with r = 7

The remainder 2 \cdot 10^{98} - s^2 is at most 2s; here both s and the remainder have 164 bits.