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