Calculating a SQRT digit-by-digit?

Question

I need to calculate the SQRT of $x$ to $y$ decimal places. I'm dealing with $128$-bit precision, giving a limit of $28$ decimal places. Obviously, if $\,y > 28$, the Babylonian method, which I'm currently using, becomes futile, as it simply doesn't offer $y$ decimal places.

My question to you is, can I calculate the Square Root digit-by-digit? For example, calculating the $n$-th digit, assuming I already have digit $\,n - 1$. I'm well aware this may be impossible, as I've not found any such method through Google yet, just checking here before I give up in search of arbitrary precision methods.

Answer 1

It used to be taught in school. You can see the method in Wikipedia. If you have high precision available, you can use Newton's method, which (once you get close) doubles the number of correct digits each time around.

Answer 2

This method is perhaps not very practical, and I don't know how to properly typeset (or even explain) this method, but I will show in an example how to compute $\sqrt{130}$ digit by digit. It is very similar to ordinary long division. (If someone has ideas on how to improve the typesetting, please show me!)

Step 1:

Insert delimiters, grouping the digits two by two from the right: $$\sqrt{1/30}$$

Step 2:

Start from with the leftmost digit (-group). What is the square root of $1$? It is $1$, so the first digit is $1$. Put $1$ in a memory column (to the left in this example). Subtract $1^2=1$, and move down the next digit group, $$ \begin{array}{rcl} 1 & \qquad & \sqrt{1/30}=1\ldots \\ +1 & & -1 \\ \overline{\phantom{+}2} & & \overline{\phantom{-}030} \end{array} $$

Step 3

Add a symbol $x$ to (two places in) the memory column: $$ \begin{array}{rcl} 1\phantom{1} & \qquad & \sqrt{1/30}=1\ldots \\ +1\phantom{1} & & -1 \\ \overline{\phantom{+}2x} & & \overline{\phantom{-0}30} \\ \phantom{}x \end{array} $$ We want to find a digit $x$ such that $x\cdot 2x$ is as large as possible, but below $30$ (our current remainder). This $x$ will be the next digit in the result. In this case, we get $x=1$ ($x=3$ would for example give $3\cdot23=69$, which is too much), so we replace $x$ with $1$ in the memory column and put a $1$ in the result. Finish the step by subtracting $1\cdot 21=21$ from the remainder, and moving down the next digit group (which is $00$, since all the decimals are zero in our case) $$ \begin{array}{rcl} 1\phantom{1} & \qquad & \sqrt{1/30}=11\ldots \\ +1\phantom{1} & & -1 \\ \overline{\phantom{+}21} & & \overline{\phantom{-0}30} \\ \phantom{+2}1 & & \phantom{}-21 \\ \overline{\phantom{+}22} & & \overline{\phantom{-00}900} \end{array} $$ As we have come to moving down decimals, we should also add a decimal point to the result.

Step 4

Add a symbol $x$ to (two places in) the memory column: $$ \begin{array}{rcl} 1\phantom{1x} & \qquad & \sqrt{1/30}=11.\ldots \\ +1\phantom{1x} & & -1 \\ \overline{\phantom{+}21}\phantom{x} & & \overline{\phantom{-0}30} \\ \phantom{+2}1\phantom{x} & & \phantom{}-21 \\ \overline{\phantom{+}22x} & & \overline{\phantom{-00}900} \\ \phantom{+22}x & & \end{array} $$ Which digit $x$ now makes $x\cdot 22x$ as large as possible, but less than $900$? The answer is $x=4$, which is the next digit in the result. $$ \begin{array}{rcl} 1\phantom{1x} & \qquad & \sqrt{1/30}=11.4\ldots \\ +1\phantom{1x} & & -1 \\ \overline{\phantom{+}21}\phantom{x} & & \overline{\phantom{-0}30} \\ \phantom{+2}1\phantom{x} & & \phantom{}-21 \\ \overline{\phantom{+}224} & & \overline{\phantom{-00}900} \\ \phantom{22}+4 & & \phantom{0}-896 \\ \overline{\phantom{+}228} & & \overline{\phantom{-0000}400} \end{array} $$ Subtract, move down the next digit group, add the memory column, ...

Step n

Imitate what we did in step 4.

Answer 3

Very good software implementations of sqrt for single/double/quadruple-precision floating-point numbers can be found in:

1. glibc:

   • e_sqrtf.c: Generates the square root of a single-precision floating-point number bit by bit.
   • e_sqrt.c: Calculates the square root of a double-precision floating-point number using some method I couldn't understand (If someone can lead us to the method used here. That would be much appreciated).

2. jdk7u-jdk:

   • e_sqrt.c: Generates the square root of a double-precision floating-point number bit by bit.

3. gnulib:

   • sqrtl.c: Calculates the square root of a long double floating-point number using Newton-Raphson method.

You can dive into these implementations written by amazingly smart people and try to extend them to fit your needs. If your system supports long double as 128-bit precision, the implementation in gnulib is exactly what you need.