Computing sqrt(2) by hand

(This is Yan Sheng.)

Today's blog post will be a change of pace from the usual contest math content: one day I was wondering to myself how many digits of $\sqrt2$ I could calculate by hand (who hasn't?). Math is more fun when you try stuff out for yourself, so here's some warm-up:

Exercise: Come up with a method of calculating the decimal expansion of $\sqrt2$ by hand. Then spend 5 or 10 minutes (or until you get bored) actually executing it. Can you get the right answer to 6 digits? 8? 10? How long would it take (and how painful would it be) to get 20 correct digits using your method?

Note: throughout this post, when I refer to "decimal places" or "digits" I always mean after the decimal point, eg. $\pi\approx3.125$ is an approximation with 1 correct digit.

History of computation records

Here's a (probably incomplete) table of past records for number of digits of $\sqrt2$ calculated, partly to add some space to avoid spoilers for those who are trying the above exercise, and partly because it turned out to be surprisingly hard to compile such a list.

Source	Year	Digits
YBC 7289	1800-1600 BCE	5
Baudhayana	800-500 BCE	5
Zhu Zaiyu	1584	24
J. M. Boorman	1887	315 (486 claimed)
R. Coustal	1950	1,032
H. S. Uhler	1950	1,544
K. Takahashi, M. Sibuya	1966	14,000
M. Lal	1967	19,600
M. Lal	1968	39,000
M. Lal, W. F. Lunnon	1968	100,000
J. Dutka	1971	1,000,000
R. Nemiroff, J. Bonnell	1994	10,000,000
Y. Kanada, D. Takahashi	1997	137,438,953,444
S. Kondo	2010	1,000,000,000,000
A. Yee	2012	2,000,000,000,050
R. Watkins	2016	2,000,000,000,100
R. Watkins	2016	5,000,000,000,000
R. Watkins	2016	10,000,000,000,000
T. Hanselmann	2022	10,000,000,001,000
J. Ranous	2023	20,000,000,000,000 (unverified)

Some notes on the above table:

• Someone on Hacker News claims that Riemann computed 38 digits of $\sqrt2$ by hand, but I couldn't find any other source which corroborates this. If anyone wants to visit the Göttingen State and University Library to scan the relevant pages of Riemann's notes, please let me know.
• The introduction to the Uhler paper is pretty funny:

  The present investigation was begun with the single intention of testing a value for $\sqrt2$ as computed by J. Marcus Boorman about 65 years ago. This datum was chosen because the author had acquired the impression that Boorman seldom if ever produced a reliable result in his calculations and it seemed fair to check this particular approximation since it was claimed that he had verified many of his digits.

  Spoiler alert: not all of Boorman's digits were correct.

• Breaking news: the latest record was just announced yesterday (Nov 17, 2023)!

Back to calculation by hand: here are some of the things I tried, see how you compare!

Digit-by-digit

Let's start with some obvious facts: we have$$1^2=1<2<4=2^2,$$so $1<\sqrt2<2$, ie. the digit of $\sqrt2$ before the decimal point is $1$.

How should we find the next digit? Write $\sqrt2=1+c/10$, so we want to solve the equation$$\left(1+\frac c{10}\right)^2=2,$$which after expanding and clearing denominators gives$$(20+c)c=200-10^2=100.$$Now we check that $24\times4<100<25\times5$, so $4<c<5$, and thus we have $\sqrt2=1.4\ldots$.

We can repeat this process for the next digit: writing $c=4+c'/10$, we want to solve$$\left(24+\frac{c'}{10}\right)\left(4+\frac{c'}{10}\right)=100,$$or$$(280+c')c'=100(100-24\times4)=400.$$Check that $281\times1<400<282\times2$, so $1<c'<2$, and we have $\sqrt2=1.41\ldots$.

One more time, just to be sure: writing $c'=1+c''/10$, we want to solve$$\left(281+\frac{c''}{10}\right)\left(1+\frac{c''}{10}\right)=400,$$or$$(2820+c'')c''=100(400-281\times1)=11900.$$Check that $2824\times4<11900<2825\times5$, so $4<c''<5$, and we have $\sqrt2=1.414\ldots$.

Okay, so we can continue this process and just calculate the digits one by one. We don't want to write down a quadratic equation for every digit, so conventionally the above procedure is written in a way similar to long division, as follows: $$\begin{align*} 1\,\phantom{0}4\,\phantom{0}1\,\phantom{0}4\,\phantom{0}2\,\phantom{0}1\\ \sqrt{2\,00\,00\,00\,00\,00}\\ \underline{{}-1}\phantom{\,00\,00\,00\,00\,00}\\ 1\,00\phantom{\,00\,00\,00\,00}\\ \underline{{}-96}\phantom{\,00\,00\,00\,00}&&&(24\times4)\\ 4\,00\phantom{\,00\,00\,00}\\ \underline{{}-2\,81}\phantom{\,00\,00\,00}&&&(281\times1)\\ 1\,19\,00\phantom{\,00\,00}\\ \underline{{}-1\,12\,96}\phantom{\,00\,00}&&&(2824\times4)\\ 6\,04\,00\phantom{\,00}\\ \underline{{}-5\,65\,64}\phantom{\,00}&&&(28282\times2)\\ 38\,36\,00\\ \underline{{}-28\,28\,41}&&&(282841\times1)\\ 10\,07\,59 \end{align*}$$

Exercise: Figure out how the procedure we were performing with quadratic equations is equivalent to this "long division" method.

Exercise: Use this method to compute the square root of your favourite number between 1 and 100 (which is not a square) to 4 decimal places.

Of course, in principle this method gives us as many digits as we want, but in practice I get bored by the time the calculation gets as far as shown above. This suggests a rating system:

Digits Until I Got Bored: 5

Time to move on and brainstorm some other methods.

Heron's method

Suppose we have some approximation $x_0\approx\sqrt2$. If we had a procedure that uses $x_0$ to somehow generate a better approximation $x_1$, then by repeating same procedure we could get an even better approximation $x_2$, then $x_3$ and so on, giving us a sequence of numbers that converges to $\sqrt2$. Now the problem of calculating $\sqrt2$ to a given precision reduces to calculating a good enough approximation $x_n$.

Okay, so how do we generate a better approximation $x'$ from a given approximation $x$? Actually, how do we generate any other approximation from $x$? We probably need to think about some properties that characterise $\sqrt2$.

For example, we might consider the equation$$\sqrt2=\frac2{\sqrt2}.$$This suggests that if $x$ is a good approximation to $\sqrt2$, then so should $\dfrac2x$. In fact, if we write $x=\sqrt2+\varepsilon$, then we can expand the fraction $\dfrac2x$ as$$\frac2x=\sqrt2-\varepsilon+\frac{\varepsilon^2}{\sqrt2+\varepsilon}.$$This says that $\dfrac2x$ is nearly as good of an approximation as $x$, but in the opposite direction: the error changes sign. Hence the average $\dfrac12\left(x+\dfrac2x\right)$ should be an even better approximation of $\sqrt2$, with error on the order of $\varepsilon^2$ instead of $\varepsilon$.

Thus we expect that the sequence defined by$$x_{n+1}=\dfrac12\left(x_n+\dfrac2{x_n}\right)$$should converge to $\sqrt2$, provided that the initial approximation $x_0$ is good enough. The convergence turns out to be pretty fast, since the error is roughly squared with each iteration (this is known as quadratic convergence); equivalently, the number of correct digits doubles at every step!

Exercise: Show that in fact $x_n\to\sqrt2$ as long as $x_0>0$. What happens if $x_0<0$?

Exercise: Show that we also get the above recurrence by applying Newton's method to the function $f(x)=x^2-2$.

The above procedure is known as Heron's method (yes, same guy as Heron's formula), or the Babylonian method (even though we don't know if the Babylonians knew about it).

Let's see this in action. We'll start with $x_0=1$, and calculate $$\begin{align*} x_1=\frac12\left(1+\frac21\right)=\frac32&=1.5000000\ldots\\ x_2=\frac12\left(\frac32+\frac43\right)=\frac{17}{12}&=1.4166666\ldots\\ x_2=\frac12\left(\frac{17}{12}+\frac{24}{17}\right)=\frac{577}{408}&=1.4142156\ldots\\ \end{align*}$$ Here $x_2$ gives 5 correct digits of $\sqrt2$.

At this point I realise a key flaw of this method: I don't actually want to divide large integers by hand, and the next term in this sequence ($x_3=\dfrac{665857}{470832}$) already looks pretty bad to me, even though it probably has at least 10 correct digits. So maybe this is the end of the road, and thus

Digits Until I Got Bored: 5…

…unless?

Let's give the method one more shot, choosing a different starting value $x_0=1.4=\dfrac75$. Then$$x_1=\frac12\left(\frac75+\frac{10}7\right)=\frac{99}{70}=1.4142857\ldots,$$which is alright since dividing by 7 isn't too bad. Now the miracle is that dividing by 99 is also not too bad, so the next iteration is actually really easy: $$\begin{align*} x_1=\frac{99}{70}&=1.4142857142\ldots\\ \frac2{x_1}=\frac{140}{99}&=1.4141414141\ldots\\ \hline x_2&=1.4142135642\ldots \end{align*}$$ This approximation has 8 correct digits. But now our luck runs out, because $x_2=\dfrac{19801}{13860}$, and I'm not really up to dividing by 19801. Final verdict:

Digits Until I Got Bored: 8

Exercise: Could we have come up with a recurrence without using division by $x_n$? Consider the identity$$\sqrt2=\dfrac12{\sqrt2^3}.$$Hence if $x\approx\sqrt2$, then $\dfrac12x^3\approx\sqrt2$ as well. How good is the latter approximation compared to the former? (Note that here the two errors have the same sign.) Hence show that the sequence$$x_{n+1}=\frac32x_n-\frac14x_n^3$$converges quadratically to $\sqrt2$, if $x_0$ is close enough to $\sqrt2$. Which function is this the Newton's method iteration for?

Exercise: We can simplify the above recurrence even further for hand calculations. Let $b_n=\dfrac12x_n^2$; show that$$b_{n+1}=b_n\left(\frac{3-b_n}2\right)^2.$$Writing $Y_n=\dfrac{3-b_n}2$, we have $b_n=b_0Y_0^2Y_1^2\cdots Y_{n-1}^2$, and $b_n\to1$. Hence deduce the relation between the product $Y_0Y_1\cdots$ and $\sqrt2$. If $Y_n=1+\varepsilon_n$, express $\varepsilon_{n+1}$ as a function of $\varepsilon_n$. Now actually try to approximate $\sqrt2$ by calculating a few steps of this recurrence. (This method is known as Goldschmidt's algorithm.)

Rational approximations

If we want a good rational approximation $\dfrac ab\approx\sqrt2$, one possible requirement is for $\lvert a-b\sqrt2\rvert$ to be small. How do we generate numbers of the form $a-b\sqrt2$ (with $a,b$ nonzero integers) which are close to 0? Here's an idea: since $\lvert1-\sqrt2\rvert<1$, we could consider its powers$$(1-\sqrt2)^n=a_n-b_n\sqrt2,$$which converge to 0 exponentially fast. Then $\dfrac{a_n}{b_n}$ is a sequence of rational approximations converging to $\sqrt2$. The terms $a_n$ and $b_n$ can be computed with the recurrence $$\begin{align*} a_{n+1}-b_{n+1}\sqrt2&=(1-\sqrt2)(a_n-b_n\sqrt2)\\ &=(a_n+2b_n)-(a_n+b_n)\sqrt2, \end{align*}$$ ie. $$\begin{align*} a_{n+1}&=a_n+2b_n,\\ b_{n+1}&=a_n+b_n. \end{align*}$$ The first few values are given as follows: $$\begin{matrix} n&0&1&2&3&4&5\\ \hline a_n&1&1&3&7&17&41\\ b_n&0&1&2&5&12&29 \end{matrix}$$

How good is this approximation, ie. what is the error like? Note that the above recurrence remains unchanged when we replace $-\sqrt2$ with $\sqrt2$, so that we also have$$(1+\sqrt2)^n=a_n+b_n\sqrt2.$$But $(1-\sqrt2)(1+\sqrt2)=-1$, so $$\begin{align*} \left\lvert\frac{a_n}{b_n}-\sqrt2\right\rvert&=\frac1{b_n}\lvert a_n-b_n\sqrt2\rvert\\ &=\frac1{b_n(a_n+b_n\sqrt2)}, \end{align*}$$ which is close to $\dfrac1{2\sqrt2b_n^2}$ for large $n$.

Hey, here's an idea: why calculate all the powers of $1-\sqrt2$, when I could get to higher powers with less calculation by repeated squaring? Let's try it out: $$\begin{align*} (1-\sqrt2)^2&=3-2\sqrt2\\ (1-\sqrt2)^4=(3-2\sqrt2)^2&=17-12\sqrt2\\ (1-\sqrt2)^8=(17-12\sqrt2)^2&=577-408\sqrt2\\ (1-\sqrt2)^{16}=(577-408\sqrt2)^2&=665857-470832\sqrt2 \end{align*}$$ Hang on, don't these numbers look kinda familiar?

Exercise: Show that squaring $a-b\sqrt2$ is equivalent to performing one iteration of Heron's method on $\dfrac ab$.

There are several other approaches to getting good rational approximations of $\sqrt2$, which I'll leave as exercises:

Exercise: Notice that$$\sqrt2-1=\frac1{\sqrt2+1}=\frac1{2+(\sqrt2-1)}.$$Show that the sequence$$1,\quad1+\frac12,\quad1+\frac1{2+\dfrac12},\quad1+\frac1{2+\dfrac1{2+\dfrac12}},\quad\ldots$$converges to $\sqrt2$. Simplify each of these continued fractions; what do you see?

Exercise: Another possible indication that a rational approximation $\dfrac ab\approx\sqrt2$ is good is when $\lvert a^2-2b^2\rvert$ is small. Clearly $a^2-2b^2$ cannot be 0, unless $a=b=0$. Find some small integer solutions to $a^2-2b^2=\pm1$; what do you see? Guess what the complete set of solutions to $a^2-2b^2=\pm1$ are, then prove your conjecture. (In general, the Pell's equation $a^2-Nb^2=1$ is closely related to the continued fraction expansion of $\sqrt N$.)

Power series

Maybe we should generalise the problem, from computing $\sqrt2$ to eg. $\sqrt{1+x}$, where $x$ is small enough so that 1 is a good initial guess. What's a good first order correction to the initial guess? If $\sqrt{1+x}\approx1+a_1x,$ then we want$$1+x\approx(1+a_1x)^2=1+2a_1x+a_1x^2,$$so when $x$ is small our best choice is $a_1=\dfrac12$.

What's the second order correction? If $\sqrt{1+x}\approx1+\dfrac12x+a_2x^2$, then we want $$\begin{align*} 1+x&\approx\left(1+\frac12x+a_2x^2\right)^2\\ &=1+x+\left(2a_2+\frac14\right)x^2+O(x^3), \end{align*}$$ so we should choose $a_2=-\dfrac18$.

Of course, we could continue this procedure for higher order corrections, and obtain a power series expansion$$\sqrt{1+x}=1+a_1x+a_2x^2+a_3x^3+\cdots,$$by equating coefficients on both sides of$$1+x=(1+a_1x+a_2x^2+a_3x^3+\cdots)^2,$$and solving for the $a_n$ recursively.

Fortunately, we don't actually have to do that calculation ourselves, because Newton's binomial series does the hard work for us: for all $\lvert x\rvert<1$ and real $\alpha$, we have $$(1+x)^\alpha=1+\alpha x+\frac{\alpha(\alpha-1)}{2!}x^2+\frac{\alpha(\alpha-1)(\alpha-2)}{3!}x^3+\cdots.$$

Exercise: Derive the above series by computing the $n$-th derivative of $(1+x)^\alpha$ at $x=0$, and using Taylor's theorem.

Substituting $\alpha=\dfrac12$, we get$$\sqrt{1+x}=1+\frac12x-\frac18x^2+\frac1{16}x^3-\frac5{128}x^4+\frac7{256}x^5-\cdots.$$It is tempting at this point to just substitute $x=1$ to get $\sqrt2$, but that's a bad idea: the series has radius of convergence of 1, so even though the series does converge at $x=1$ (because it's alternating), it does so very slowly:

Exercise: Show that the $n$-th coefficient of the above power series, namely$$a_n=\frac{\frac12(-\frac12)(-\frac32)\cdots(-\frac{2n-3}2)}{n!},$$satisfies$$\dfrac1{4n^{3/2}}\leq\lvert a_n\rvert\leq\dfrac1{2n^{3/2}}.$$(Hence we need roughly $10^{2k/3}$ terms to approximate $\sqrt2$ to $k$ decimal places by using the series with $x=1$.)

A much better alternative is to substitute $x=-\dfrac12$ to get$$\sqrt{1-\frac12}=\sqrt{\frac12}=\frac12\sqrt2,$$so that the series converges geometrically, although it still takes roughly $3k$ terms to get $k$ decimal places of $\sqrt2$, and I don't think I'd be too happy to compute the decimal expansion of terms like $a_{3k}x^{3k}$.

Can we do better? Instead of using the series to calculate $\sqrt2$ directly, we would really want to calculate the square root of a number closer to 1, so that the series converges faster. Here's where the rational approximations of $\sqrt2$ that we found earlier can help us: for instance, since $\sqrt2\approx\dfrac32$, we could instead calculate$$\frac23\sqrt2=\sqrt{\frac89}=\sqrt{1-\frac19}$$by substituting $x=-\dfrac19$ into the series, and then multiply the sum by $\dfrac32$. Or we could calculate$$\frac34\sqrt2=\frac3{2\sqrt2}=\sqrt{\frac98}=\sqrt{1+\frac18},$$by substituting $x=\dfrac18$ into the series, and then multiply the sum by $\dfrac43$. In either case, we now only need roughly $k$ terms to get $k$ decimal places of $\sqrt2$, which is still unpleasant, but at least it's 3 times less unpleasant than before.

So now we go through each rational approximation of $\sqrt2$ from before to write down similar identities, and see if any catches our eye. Eventually we find a miracle:$$\frac7{10}\sqrt2=\frac7{5\sqrt2}=\sqrt{\frac{49}{50}}=\sqrt{1-\frac1{50}},$$so we should substitute $x=-\dfrac1{50}$ into the series. This is very good for us, since it is very easy to divide by powers of 50.

Let's actually do this. We'll rewrite the series a little, substituting $x\to-2x$ to get$$\sqrt{1-2x}=1-x-\frac1{2!}x^2-\frac{1\cdot3}{3!}x^3-\frac{1\cdot3\cdot5}{4!}x^4+\cdots.$$We now want to take $x=0.01$. It's easy to get each term in the series from the previous term by multiplying by $\dfrac{2n-3}nx$, so we compute: $$\begin{align*} 1-\tfrac7{10}\sqrt2&=0.01\\ &+\phantom{0.00}005\\ &+\phantom{0.0000}005\\ &+\phantom{0.000000}00625\\ &+\phantom{0.00000000}00875\\ &+\phantom{0.0000000000}013125\\ &+\phantom{0.000000000000}0206\\ &+\phantom{0.00000000000000}04\\ \hline 1-\tfrac7{10}\sqrt2&=0.0100505063388335\\ \tfrac7{10}\sqrt2&=0.9899494936611665\\ \sqrt2&=\phantom01.414213562373095 \end{align*}$$ That's 15 digits of $\sqrt2$ while barely breaking a sweat. Wow!

Digits Until I Got Bored: 15++

Exercise: If you're thinking to yourself, "hmm, I bet Euler could do something like that," you'd be correct! In fact, in his 1755 textbook on differential calculus, he did it even better, and avoids the division by 7 that we had to do at the end of our calculation. By using the binomial series with $\alpha=-\dfrac12$ and the identity$$\frac57\sqrt2=\left(1-\frac1{50}\right)^{-1/2},$$compute $\sqrt2$ to 12 (or 15, or 20) decimal places.

Exercise: Instead of using the binomial series, we could also custom-design our own power series expansion. For example, using the approximation $\sqrt2\approx\dfrac{99}{70}$, we can write$$\frac{99}{140}\sqrt2=\sqrt{1+\frac1{9800}}.$$Find the first few terms of the power series $1+\frac12x^2+\cdots$ whose square is$$1+\frac{x^2}{1-2x}=1+x^2+2x^3+4x^4+8x^5+\cdots$$By substituting the appropriate value for $x$, compute $\sqrt2$ to 12 (or 15, or 20) decimal places. (If you don't know the nice trick for dividing by 99, you'll figure it out by the end of this exercise. Hint: group the digits into pairs.)

More fun things to try

• Explain the pattern in the decimal expansion of$$\frac12-\frac15\sqrt6=0.0101020514\ldots$$Hence compute $\sqrt6$ to 12 (or 15, or 20) decimal places.
• (MIT Integration Bee 2023 Finals) Evaluate$$\left\lfloor10^{20}\int_2^\infty\frac{x^9}{x^{20}-48x^{10}+575}\,dx\right\rfloor.$$(Time limit: 4 mins)
• Generalise the methods here to computing arbitrary square roots $\sqrt N$.
  • In particular, there doesn't seem to be a great power series for $\sqrt3$, in the same way that we had miracles for $\sqrt2$ and $\sqrt6$. This comes from the fact that the expressions $\lvert2m^2-10^{2n}\rvert$ and $\lvert6m^2-10^{2n}\rvert$ can take surprisingly small values for $m,n$ positive integers. What is the smallest possible value of $\lvert3m^2-10^{2n}\rvert$? The second smallest?
• Of course, "Digits Until I Got Bored" isn't a very scientific metric. Come up with a more objective measurement for the efficiency the different methods we've seen (ideas: number of 1-digit operations required for $k$ decimal places, or maybe just multiplications).
  • (hard) Generalise to a different activity that you enjoy: come up with a metric that quantifies the amount of fun you'll have, and optimise for that metric. Was that a good idea?
• Suppose you're a handheld calculator (ie. suppose you're a programmer at Hewlett-Packard in the 1970s). What method would you implement for the square root function? Compare your proposal with CORDIC, the actual algorithm used for the first handheld calculators; why did HP use CORDIC instead of your proposal?
• Suppose you're feeling very inspired by this blog post, and you want to start a project of computing $\sqrt2$ by hand to 50 or 100 decimal places. Which method would you use? How would you check for mistakes in your intermediate calculations?