Editing Integer square root (section)

==Karatsuba square root algorithm==

The Karatsuba square root algorithm is a combination of two functions: a [[Access modifiers|public]] function, which returns the integer square root of the input, and a recursive [[Access modifiers|private]] function, which does the majority of the work.

The public function normalizes the actual input, passes the normalized input to the private function, denormalizes the result of the private function, and returns that.

The private function takes a normalized input, divides the input bits in half, passes the most-significant half of the input recursively to the private function, and performs some integer operations on the output of that recursive call and the least-significant half of the input to get the normalized output, which it returns.

For [[Arbitrary-precision arithmetic|big-integers]] of "50 to 1,000,000 digits", [[Burnikel-Ziegler division|Burnikel-Ziegler Karatsuba division]] and [[Karatsuba algorithm|Karatsuba multiplication]] are recommended by the algorithm's creator.<ref>{{cite web | last = Zimmermann | first = Paul | title = Karatsuba Square Root |date = 1999 <!-- This is the date the document was finished, not the date it was submitted to Inria --> | publisher = [[French Institute for Research in Computer Science and Automation|Inria]] | publication-date = 2006-05-24 | series = Research report #3805 | url = https://inria.hal.science/inria-00072854v1/file/RR-3805.pdf | archive-url = https://web.archive.org/web/20230511212802/https://inria.hal.science/inria-00072854v1/file/RR-3805.pdf | archive-date = 2023-05-11 }}</ref>

An example algorithm for 64-bit unsigned integers is below. The algorithm:

# Normalizes the input inside {{mono|u64_isqrt}}.
# Calls {{mono|u64_normalized_isqrt_rem}}, which requires a normalized input.
# Calls {{mono|u32_normalized_isqrt_rem}} with the most-significant half of the normalized input's bits, which will already be normalized as the most-significant bits remain the same.
# Continues on recursively until there's an algorithm that's faster when the number of bits is small enough.
# {{mono|u64_normalized_isqrt_rem}} then takes the returned integer square root and remainder to produce the correct results for the given normalized {{mono|u64}}.
# {{mono|u64_isqrt}} then denormalizes the result.

<syntaxhighlight lang="rust" line="1">
/// Performs a Karatsuba square root on a `u64`.
pub fn u64_isqrt(mut n: u64) -> u64 {
    if n <= u32::MAX as u64 {
        // If `n` fits in a `u32`, let the `u32` function handle it.
        return u32_isqrt(n as u32) as u64;
    } else {
        // The normalization shift satisfies the Karatsuba square root
        // algorithm precondition "a₃ ≥ b/4" where a₃ is the most
        // significant quarter of `n`'s bits and b is the number of
        // values that can be represented by that quarter of the bits.
        //
        // b/4 would then be all 0s except the second most significant
        // bit (010...0) in binary. Since a₃ must be at least b/4, a₃'s
        // most significant bit or its neighbor must be a 1. Since a₃'s
        // most significant bits are `n`'s most significant bits, the
        // same applies to `n`.
        //
        // The reason to shift by an even number of bits is because an
        // even number of bits produces the square root shifted to the
        // left by half of the normalization shift:
        //
        // sqrt(n << (2 * p))
        // sqrt(2.pow(2 * p) * n)
        // sqrt(2.pow(2 * p)) * sqrt(n)
        // 2.pow(p) * sqrt(n)
        // sqrt(n) << p
        //
        // Shifting by an odd number of bits leaves an ugly sqrt(2)
        // multiplied in.
        const EVEN_MAKING_BITMASK: u32 = !1;
        let normalization_shift = n.leading_zeros() & EVEN_MAKING_BITMASK;
        n <<= normalization_shift;

        let (s, _) = u64_normalized_isqrt_rem(n);

        let denormalization_shift = normalization_shift / 2;
        return s >> denormalization_shift;
    }
}

/// Performs a Karatsuba square root on a normalized `u64`, returning the square
/// root and remainder.
fn u64_normalized_isqrt_rem(n: u64) -> (u64, u64) {
    const HALF_BITS: u32 = u64::BITS >> 1;
    const QUARTER_BITS: u32 = u64::BITS >> 2;
    const LOWER_HALF_1_BITS: u64 = (1 << HALF_BITS) - 1;

    debug_assert!(
        n.leading_zeros() <= 1,
        "Input is not normalized: {n} has {} leading zero bits, instead of 0 or 1.",
        n.leading_zeros()
    );

    let hi = (n >> HALF_BITS) as u32;
    let lo = n & LOWER_HALF_1_BITS;

    let (s_prime, r_prime) = u32_normalized_isqrt_rem(hi);

    let numerator = ((r_prime as u64) << QUARTER_BITS) | (lo >> QUARTER_BITS);
    let denominator = (s_prime as u64) << 1;

    let q = numerator / denominator;
    let u = numerator % denominator;

    let mut s = (s_prime << QUARTER_BITS) as u64 + q;
    let mut r = (u << QUARTER_BITS) | (lo & ((1 << QUARTER_BITS) - 1));
    let q_squared = q * q;
    if r < q_squared {
        r += 2 * s - 1;
        s -= 1;
    }
    r -= q_squared;

    return (s, r);
}
</syntaxhighlight>