Notation

Let ${\displaystyle B}$ be the base of the number system you are using, and ${\displaystyle n}$ be the degree of the root to be extracted. Let ${\displaystyle x}$ be the radicand processed thus far, ${\displaystyle y}$ be the root extracted thus far, and ${\displaystyle r}$ be the remainder. Let ${\displaystyle \alpha }$ be the next ${\displaystyle n}$ digits of the radicand, and ${\displaystyle \beta }$ be the next digit of the root. Let ${\displaystyle x'}$ be the new value of ${\displaystyle x}$ for the next iteration, ${\displaystyle y'}$ be the new value of ${\displaystyle y}$ for the next iteration, and ${\displaystyle r'}$ be the new value of ${\displaystyle r}$ for the next iteration. These are all integers.

Invariants

At each iteration, the invariant ${\displaystyle y^{n}+r=x}$ will hold. The invariant ${\displaystyle (y+1)^{n}>x}$ will hold. Thus ${\displaystyle y}$ is the largest integer less than or equal to the ${\displaystyle n}$th root of ${\displaystyle x}$, and ${\displaystyle r}$ is the remainder.

Initialization

The initial values of ${\displaystyle x,y}$, and ${\displaystyle r}$ should be 0. The value of ${\displaystyle \alpha }$ for the first iteration should be the most significant aligned block of ${\displaystyle n}$ digits of the radicand. An aligned block of ${\displaystyle n}$ digits means a block of digits aligned so that the decimal point falls between blocks. For example, in 123.4 the most significant aligned block of two digits is 01, the next most significant is 23, and the third most significant is 40.

Main loop

On each iteration we shift in ${\displaystyle n}$ digits of the radicand, so we have ${\displaystyle x'=B^{n}x+\alpha }$ and we produce one digit of the root, so we have ${\displaystyle y'=By+\beta }$. The first invariant implies that ${\displaystyle r'=x'-y'^{n}}$. We want to choose ${\displaystyle \beta }$ so that the invariants described above hold. It turns out that there is always exactly one such choice, as will be proved below.

Proof of existence and uniqueness of ${\displaystyle \beta }$

By definition of a digit, ${\displaystyle 0\leq \beta <B}$, and by definition of a block of digits, ${\displaystyle 0\leq \alpha <B^{n}}$

The first invariant says that:

${\displaystyle x'=y'^{n}+r'}$

or

${\displaystyle B^{n}x+\alpha =(By+\beta )^{n}+r'.}$

So, pick the largest integer ${\displaystyle \beta }$ such that

${\displaystyle (By+\beta )^{n}\leq B^{n}x+\alpha .}$

Such a ${\displaystyle \beta }$ always exists, since ${\displaystyle 0\leq \beta <B}$ and if ${\displaystyle \beta =0}$ then ${\displaystyle B^{n}y^{n}\leq B^{n}x+\alpha }$, but since ${\displaystyle y^{n}\leq x}$, this is always true for ${\displaystyle \beta =0}$. Thus, there will always be a ${\displaystyle \beta }$ that satisfies the first invariant

Now consider the second invariant. It says:

${\displaystyle (y'+1)^{n}>x'}$

or

${\displaystyle (By+\beta +1)^{n}>B^{n}x+\alpha .}$

Now, if ${\displaystyle \beta }$ is not the largest admissible ${\displaystyle \beta }$ for the first invariant as described above, then ${\displaystyle \beta +1}$ is also admissible, and we have

${\displaystyle (By+\beta +1)^{n}\leq B^{n}x+\alpha .}$

This violates the second invariant, so to satisfy both invariants we must pick the largest ${\displaystyle \beta }$ allowed by the first invariant. Thus we have proven the existence and uniqueness of ${\displaystyle \beta }$.

To summarize, on each iteration:

1. Let ${\displaystyle \alpha }$ be the next aligned block of digits from the radicand
2. Let ${\displaystyle x'=B^{n}x+\alpha }$
3. Let ${\displaystyle \beta }$ be the largest ${\displaystyle \beta }$ such that ${\displaystyle (By+\beta )^{n}\leq B^{n}x+\alpha }$
4. Let ${\displaystyle y'=By+\beta }$
5. Let ${\displaystyle r'=x'-y'^{n}}$

Now, note that ${\displaystyle x=y^{n}+r}$, so the condition

${\displaystyle (By+\beta )^{n}\leq B^{n}x+\alpha }$

is equivalent to

${\displaystyle (By+\beta )^{n}-B^{n}y^{n}\leq B^{n}r+\alpha }$

and

${\displaystyle r'=x'-y'^{n}=B^{n}x+\alpha -(By+\beta )^{n}}$

is equivalent to

${\displaystyle r'=B^{n}r+\alpha -((By+\beta )^{n}-B^{n}y^{n}).}$

Thus, we do not actually need ${\displaystyle x}$, and since ${\displaystyle r=x-y^{n}}$ and ${\displaystyle x<(y+1)^{n}}$, ${\displaystyle r<(y+1)^{n}-y^{n}}$ or ${\displaystyle r<ny^{n-1}+O(y^{n-2})}$, or ${\displaystyle r<nx^{{n-1} \over n}+O(x^{{n-2} \over n})}$, so by using ${\displaystyle r}$ instead of ${\displaystyle x}$ we save time and space by a factor of 1/${\displaystyle n}$. Also, the ${\displaystyle B^{n}y^{n}}$ we subtract in the new test cancels the one in ${\displaystyle (By+\beta )^{n}}$, so now the highest power of ${\displaystyle y}$ we have to evaluate is ${\displaystyle y^{n-1}}$ rather than ${\displaystyle y^{n}}$.

Summary

1. Initialize ${\displaystyle r}$ and ${\displaystyle y}$ to 0.
2. Repeat until desired precision is obtained:
   1. Let ${\displaystyle \alpha }$ be the next aligned block of digits from the radicand.
   2. Let ${\displaystyle \beta }$ be the largest ${\displaystyle \beta }$ such that ${\displaystyle (By+\beta )^{n}-B^{n}y^{n}\leq B^{n}r+\alpha .}$
   3. Let ${\displaystyle y'=By+\beta }$.
   4. Let ${\displaystyle r'=B^{n}r+\alpha -((By+\beta )^{n}-B^{n}y^{n}).}$
   5. Assign ${\displaystyle y\leftarrow y'}$ and ${\displaystyle r\leftarrow r'.}$
3. ${\displaystyle y}$ is the largest integer such that ${\displaystyle y^{n}<xB^{k}}$, and ${\displaystyle y^{n}+r=xB^{k}}$, where ${\displaystyle k}$ is the number of digits of the radicand after the decimal point that have been consumed (a negative number if the algorithm has not reached the decimal point yet).

Square root of 2 in binary

      1. 0  1  1  0  1
    ------------------
_  / 10.00 00 00 00 00     1
 \/   1                  + 1
     -----               ----
      1 00                100
         0               +  0
     --------            -----
      1 00 00             1001
        10 01            +   1
     -----------         ------
         1 11 00          10101
         1 01 01         +    1
         ----------      -------
            1 11 00       101100
                  0      +     0
            ----------   --------
            1 11 00 00    1011001
            1 01 10 01          1
            ----------
               1 01 11 remainder

Square root of 3

     1. 7  3  2  0  5
    ----------------------
_  / 3.00 00 00 00 00
 \/  1 = 20×0×1+1^2
     -
     2 00
     1 89 = 20×1×7+7^2 (27 x 7)
     ----
       11 00
       10 29 = 20×17×3+3^2  (343 x 3)
       -----
          71 00
          69 24 = 20×173×2+2^2 (3462 x 2)
          -----
           1 76 00
                 0 = 20×1732×0+0^2 (34640 x 0)
           -------
           1 76 00 00
           1 73 20 25 = 20×17320×5+5^2 (346405 x 5)
           ----------
              2 79 75

Cube root of 5

     1.  7   0   9   9   7
    ----------------------
_ 3/ 5. 000 000 000 000 000
 \/  1 = 300×(0^2)×1+30×0×(1^2)+1^3
     -
     4 000
     3 913 = 300×(1^2)×7+30×1×(7^2)+7^3
     -----
        87 000
             0 = 300×(17^2)×0+30×17×(0^2)+0^3
       -------
        87 000 000
        78 443 829 = 300×(170^2)×9+30×170×(9^2)+9^3
        ----------
         8 556 171 000
         7 889 992 299 = 300×(1709^2)×9+30×1709×(9^2)+9^3
         -------------
           666 178 701 000
           614 014 317 973 = 300×(17099^2)×7+30×17099×(7^2)+7^3
           ---------------
            52 164 383 027

Fourth root of 7

     1.   6    2    6    5    7
    ---------------------------
_ 4/ 7.0000 0000 0000 0000 0000
 \/  1 = 4000×(0^3)×1+600×(0^2)×(1^2)+40×0×(1^3)+1^4
     -
     6 0000
     5 5536 = 4000×(1^3)×6+600×(1^2)×(6^2)+40×1×(6^3)+6^4
     ------
       4464 0000
       3338 7536 = 4000×(16^3)×2+600×(16^2)×(2^2)+40×16×(2^3)+2^4
       ---------
       1125 2464 0000
       1026 0494 3376 = 4000×(162^3)×6+600×(162^2)×(6^2)+40×162×(6^3)+6^4
       --------------
         99 1969 6624 0000
         86 0185 1379 0625 = 4000×(1626^3)×5+600×(1626^2)×(5^2)+
         -----------------   40×1626×(5^3)+5^4
         13 1784 5244 9375 0000
         12 0489 2414 6927 3201 = 4000×(16265^3)×7+600×(16265^2)×(7^2)+
         ----------------------   40×16265×(7^3)+7^4
          1 1295 2830 2447 6799