Calculating a SQRT digit-by-digit? 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.
 A: 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.
A: 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.
A: As explained by two other answers, there is a method that yields the decimal places one by one. Unfortunately for you, the method still requires you to store and manipulate an approximation and the associated increment. Your problem is that you can only store them to 28 decimal places, so you will still not be able to go very far.
You can't have "arbitrary precision"; you will need to decide how many decimal places you are going to have. Suppose you decide on 100 dp. Since 100/28 is a little less than 4, you will need at least four 128 bit numbers to represent a single number. You will also need to create routines for operating on such numbers that are combinations of other numbers. As one of the operations you need to perform is multiplying, you may have difficulties when multiplying two 28 dp numbers. For that reason you could simplify some stages by using your (potentially) 28 dp numbers to hold number sthat never exceed 14 dp, except temporarily in the middle of a calculation. 
A: According to https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Mathematical_constants, square root of 2 can be computed as fast as the optimal multiplication algorithm. According to https://en.wikipedia.org/wiki/F%C3%BCrer%27s_algorithm, multiplication can be done in time n log n 2^(log * n), where log * n is the super-logarithm of n, the inverse tetration of n to the base 2.
A: Very good software implementations of sqrt for single/double/quadruple-precision floating-point numbers can be found in:


*

*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).


*jdk7u-jdk: 


*

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


*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.
