sqrt_nm1 (n, A): Input: A = {a_{n-1}, ..., a_0}, n >= 2 Output: X = {x_{n-2}, ..., x_0} \approx sqrt(B^{n-2} A) if (n == 2) x_0 <-- floor (sqrt (A)) return X = {x_0} if (n == 3) // Needs normalization c <-- count_leading_zeros (a_2) & ~1 A' <-- 2^c A x_1 = floor (sqrt ({a'_2, a'_1})) E = A' - B x_1^2 return X = (B x_1 + floor (E / 2 x_1)) 2^{-c/2} k <-- floor (n/2) H <-- sqrt_nm1 (n - (k-1), {a_{n-1},...,a_{k-1}}) // n - (k-1) limbs if (n odd) // n = 2k+1 // We have a (k+1)-limb H, and produce k-1 low limbs. E <-- B A - H^2 X <-- B^{k-1} H + floor (B^{k-1} E / 2H) // XXX else // n = 2k // We have a k-limb H, and produce k-1 low limbs E <-- A - H^2 X <-- B^{k-1} H + floor (B^{k-1} E / 2H) return X sqrt_n (n, A): Input: A = {a_{n-1}, ..., a_0}, n >= 1 Output: X = {x_{n-1}, ..., x_0} \approx sqrt(B^{n-1} A) if (n == 1) x_0 <-- floor (sqrt (a_0)) return X = {x_0} if (n == 2) h <-- floor (sqrt(A)) E <-- A - h*h X <-- sqrt(B) h + floor (sqrt(B) E / 2h) k <-- floor (n/2) H <-- sqrt_n(k+1, {a_{n-1},...,a_{n-1-k}}) // k+1 limbs if (n odd) // n = 2k+1 // We have (k+1)-limb H, and produce k low limbs E <-- A - H^2 X <-- B^k H + floor (B^k E / 2 H) else // n = 2k // We have (k+1)-limb H, and produce k-1 low limbs E <-- B A - H^2 X <-- B^{k-1} H + floor (B^{k-1} E / 2H) return X Size of remainder: Assume X = floor(sqrt(A)) = sqrt(A) - e, where A is n limbs. Define the remainder R = A - X^2. Then R = A - (sqrt(A) - e)^2 = 2 e sqrt(A) - e^2 = e (2 sqrt(A) - e) < 2 sqrt(A) If n is even, n = 2k, then A - X^2 < B^{k+1} If n is odd, n = 2k+1, then A - X^2 < B^{k+1} So in all cases, R < B^{floor(n/2) + 1}, and X < B^{ceil(n/2)}, and with almost half a limb margin. Correctness test: 0 > A - (X+1)^2 = A - X^2 - 2X - 1 = R - 2X - 1 R < 2X + 1 R <= 2X Small cases. Consider sqrtrem, size n. n = 1: sqrt_1 n = 2: sqrt_2 n = 3: sqrt_n(2) -> sqrt_2 + special update n = 4: sqrt_nm1(3) -> sqrt_nm1(2) -> sqrt_2 n = 5: sqrt_n (3) -> sqrt_n(2) -> sqrt_2 + special update n = 6: sqrt_nm1 (4) -> sqrt_nm1 (3) -> sqrt_nm1(2) -> sqrt_2