How to Manually Compute Square Root

09 Jul 2018

Despite the popularity of computers and calculators, I would always consider how we can possibly compute irrational numeric values, such as square root, trigonometry functions and so on. This topic may be worthless in modern age, but it might help in extreme conditions.

Using Taylor Series

The Taylor series method is probably the most universal method in computing transcendental functions. It stats that a function can be approximated by its derivatives, that is

\[\begin{align*} f(x) = \sum_{n = 0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^n. \end{align*}\]

Consider function \((x+1)^{\frac{1}{2}}\), whose Maclaurin series is

\[\begin{align*} (x+1)^{\frac{1}{2}} &= \sum_{n = 0}^{\infty} \frac{\left.\frac{d}{da^n}(1+a)^{\frac{1}{2}}\right|_{a = 0}}{n!} \cdot x^n \\ &= \sum_{n = 0}^{\infty} \binom{\frac{1}{2}}{n}x^n, \end{align*}\]

where

\[\begin{equation*} \binom{\alpha}{n} =  \begin{cases} 0 &, n = 0\\ \displaystyle \prod_{k = 1}^n \frac{\alpha - k + 1}{k} &, n \geq 1 . \end{cases} \end{equation*}\]

That is to say, 

\[\begin{align*} (1+x)^{\frac{1}{2}} = 1 + \frac{1}{2}x - \frac{1}{8}x^2 + \frac{1}{16}x^3 - \frac{5}{128}x^4 + \frac{7}{256}x^5 - \cdots \ . \end{align*}\]

Using Newton’s Method

Newton’s method can be used to find the approximation of roots very efficiently. To find \(x\) so that \(f(x) = 0\), the Newton’s method starts with an initial estimation \(x_0\), and repeats the iteration

\[\begin{align*} x_{n+1} = x_n - \frac{f(x_n)}{f^{\prime}(x_n)} \end{align*}\]

until convergence. Then the lastest \(x_{n+1}\) is a valid approximation of \(x\). It is worth noticing that the Newton’s method only finds one root, and which root it finds depends on the initial value \(x_0\). Finding the value of \(\sqrt{a}\) is equivalent to solving \(f(x) = x^2 - a = 0\). It can be seen that now

\[\begin{align*} x_{n+1} &= x_n - \frac{x^2 - a}{2x} \\ &= x_n -  \frac{1}{2} \left( x - \frac{a}{x} \right). \end{align*}\]

To guarantee that the the Newton’s Method finds the positive root of \(f(x)\), we can let \(x_0 > \sqrt{a}\), which means we can simply let \(x_0 = a + 1\). It practice, Newton’s method gives higher precision after same number of iterations. We can get very accurate estimation of square root after 7-8 iterations.

When Computing By Hand

The methods above are more suitable for computers. Human beings are better at handling integers rather than decimals, and there is indeed a square root algorithm that is based on integer arithmetic. Because it is a more indirect method, I am going to spend more time describing and analyzing its methodology. This method is originally documented in The Nine Chapters on the Mathematical Art, an ancient Chinese book on math. It is really amazing to know that Chinese people are able to formulate such a ingenious method at such a long time ago. The method is based on a trail-and-error approach, but it is much more efficient than bisection, and it is based on integer arithmetic.

Assume we are to square root \(n\), where \(n \geq 0, n \in \mathbb{Z}\). The basic procedure of the method is as follows:

1. Divide \(n\) into 2-digit parts from right to left. When the digit of \(n\) is odd, assume that there is a zero in front of the highest digit. Usually a vertical line is used to separate each part.

   

2. Based on the number in the leftmost part, find out the highest digit of the result. If we denote the number in the leftmost part by \(a_1\) and the highest digit of the result by \(b_1\), we are to find the largest \(b_1\) so that \(b_1^2 \leq a_1\).

3. Subtract the number in the leftmost part by the highest digit of the result squared, and then combine the difference with the next part to form the first remainder.

4. Divide the remainder by derived digits of the result multiplied by 20. The integral part of the quotient is called the hypothetical quotien. (Note that the hypothetical quotient is greater than or equal to 0 and less than or equal to 9)

5. Multiply the derived digits of the result by 20, and subsequently add it by the hypothetical quotient. Then multiply the sum by the hypothetical quotient. If the final product is less than or equal to the remainder, then the hypothetical quotient is the next digit of the result. Otherwise, repetitively subtract the hypothetical quotient by 1 and find the largest hypothetical quotient whose corresponding final product is less than or equal to the remainder. The new remainder is then the last remainder subtracted by the final product.

6. Repeat step 4 and 5 to compute the remaining digits of the result.

It can be seen that the amount of calculation grows almost quadratically. After 4-5 digits, it becomes quite difficult to compute by hand.

How Does It Work?

A very important task is to figure out why this method can accurately compute square root. For simplicity, we use the notation \(\mathbb{Z}(p,q)\) to denote the set \(\mathbb{Z} \cap [p,q]\), that is, the integers in range \([p,q]\).

Assume that \(n\) has \(N\) digits, which is denoted by \(n = o_1o_2\cdots o_N\). It can be seen that

\[n = \sum_{i = 1}^{N}o_i \times 10^{N - i}, \]

where \(o_1 \in \mathbb{Z}(1,9), o_i (i\geq 2) \in  \mathbb{Z}(0,9)\).

After the first step of the algorithm, the number is divided into several 2-digit segments. Let \(K = \lceil \frac{N}{2} \rceil\), then it can be denoted by \(n = a_1a_2\cdots a_{K}\). Similarly, we have

\[n = \sum_{i = 1}^{K} a_i \times 10^{2(K - i)}, \]

where \(a_i \in \mathbb{Z}(1,99), a_i (i\geq 2) \in \mathbb{Z}(0,99)\).

Given that \(n\) is a \(N\) digit number, it can be seen that \(10^N \leq n < 10^{N+1}\). Therefore \(10^{\lfloor \frac{N}{2}  \rfloor} \leq \sqrt{n} < 10^{\lceil \frac{N+1}{2} \rceil}\). Discuss the situation based on whether \(N\) is odd or even, it can be easily concluded that the integral part of \(\sqrt{n}\) has \(K\) digits. Denote the integral part of \(\sqrt{n}\)  by \(m =  b_1b_2\cdots b_K\). As mentioned above, 

\[m = \sum_{i = 1}^{K}b_i \times 10^{K - i},\]

where \(b_1 \in \mathbb{Z}(1,9), b_i (i\geq 2) \in  \mathbb{Z}(0,9)\).

Now \(b_1\) is selected to be the largest integer that satifies \(b_1^2 \leq a_1\). After step 3, \(n\) is converted into the following form:

\[\begin{align}\label{one} n = b_1^2 \times 10^{2(K-1)} + (a_1 - b_1^2) \times 10^{2(K-1)} + a_2 \times 10^{2(K-2)} + \sum_{i = 3}^{K} a_i \times 10^{2(K - i)}. \end{align}\]

Now \(\sqrt{\left[b_1 \times 10^{(K-1)}\right]^2}\) is our first approximation for \(\sqrt{n}\). We want to push the precision furthur by one digit, that is, to approximate \(n\) with \(\sqrt{\left[b_1 \times 10^{(K-1)} + b_2 \times 10^{(K-2)}\right]^2}\). Insert this into the equation \ref{one}, we have

\[\begin{multline} n = \left[ b_1 \times 10^{(K-1)} + b_2 \times 10^{(K-2)} \right] ^2  \\ - 2 \cdot b_1^2 \times 10^{(K-1)} \cdot b_2^2 \times 10^{(K-2)} - b_2^2 \times 10^{2(K-2)} \\ + (a_1 - b_1^2) \times 10^{2(K-1)} + a_2 \times 10^{2(K-2)} \\+ \sum_{i = 3}^{K} a_i \times 10^{2(K - i)}. \end{multline}\]

After a bit of arranging, it can be seen that 

\[\begin{multline}\label{two} n = \left[ b_1 \times 10^{(K-1)} + b_2 \times 10^{(K-2)} \right] ^2  \\ + 10^{2(K-2)}\left[  10^2 \cdot (a_1 - b_1^2) + a_2 - b_2 (20b_1 + b_2) \right] \\+ \sum_{i = 3}^{K} a_i \times 10^{2(K - i)}. \end{multline}\]

It can be seen that \(10^2 \cdot (a_1 - b_1^2) + a_2\) is the new remainder, and \(b_2 (20b_1 + b_2)\) is just identical to the procedure in step 5. Obviously there is no negative number throughout the entire computation, therefore, we must guarantee that

\[\begin{align}\label{three} 10^2 \cdot (a_1 - b_1^2) + a_2 - b_2 (20b_1 + b_2) \geq 0. \end{align}\]

The left side is a quadratic funtion of \(b_2\), which monotonically decreases on its domain.Therefore, we are to find the largest \(b_2\) that satifies \(\ref{three}\). Rearrange \(\ref{three}\) to make it the following inequality:

\[\begin{align}\label{four} 10^2 \cdot (a_1 - b_1^2) + a_2 - 20b_1b_2  - b_2^2 \geq 0. \end{align}\]

A necessary condition for this inequality is that

\[\begin{align} 10^2 \cdot (a_1 - b_1^2) + a_2 - 20b_1b_2 \geq 0, \end{align}\]

which gives

\[\begin{align} b_2 \leq \frac{10^2 \cdot (a_1 - b_1^2) + a_2}{20b_1}. \end{align}\]

Just as it is done in step 5, this inequality yields an upper bound for \(b_2\) and facilitates calculation. In extreme cases, the right side can be greater than or equal to 10, then we might need to start experimenting hypothetical quotient from 9. From now on it can be seen that the method is based on trail-and-error, and it tries to find the closest estimation at each step, which is the reason why we need \(b_i\) values to be as great as possible. Failure in finding such values might not cause contradiction mathematically, but it will result in bad estimation. It is not difficult to tell that the estimation will always be less than the actual value.

For the estimation at each step to be accurate, we need to assert that \(b_i  < 10 (i \geq 1)\) holds. In fact, solving \(\ref{four}\) with regards to \(b_2\) directly yields

\[\begin{align} \frac{-20b_1 - \sqrt{400a_1 + 4a_2}}{2} \leq  0 \leq  b_2 \leq \frac{-20b_1 + \sqrt{400a_1 + 4a_2}}{2} . \end{align}\]

Thus, we need to prove that

\[\begin{align} \frac{-20b_1 + \sqrt{400a_1 + 4a_2}}{2} < 10, \end{align}\]

which is equivalent to proving

\[\begin{align}\label{five} 400(a_1 - (b_1^2 + 2b_1) - 1) + 4a_2 < 0. \end{align}\]

Fortunately, according to the property of \(b_1\), we know that \(a_1 \leq (b_1 + 1)^2 - 1 = b_1^2 + 2b_1\). Therefore 

\[\begin{align} 400(a_1 - (b_1^2 + 2b_1) - 1) + 4a_2 < -400 + 4a_2. \end{align}\]

Because \(a_2 \in \mathbb{Z}(0,99)\), it can be seen that \(\ref{five}\) always holds. That is to say this algorithm is guaranteed to be able to find the best estimation at each step.

In fact, if we extend the definition of \(a_i\) and \(b_i\), the proof can be easily extended to the following steps of this algorithm. I am not going to cover too much detail on that.

Epilogue

I write this article because I cannot describe how stunned I am when I first read about this method on the Internet. Ancient China is known to be primitive and undeveloped. But the width and depth of our ancient civilization had come to a point where no other civilization can compare. There is no doubt that Chinese people have succeeded in demonstrating their wisdom and capability to the entire world. We used to be the past, we possess the potential now, and we will be in the future. I am proud to be born Chinese.