Iterative Methods for Computing the Square Root
Computing the Square Root using Halley’s Method
Explanation, Python Implementation and Numerical Analysis
While Edmond Halley is widely known for Halley’s Comet and his contributions to the Inverse Square Law in physics, his improvement to Netwon’s method for computing the square root (aka the Babylonian / Heron’s method) remains a lesser known side achievement.
In this article we explore Halley’s iterative method to calculate the square root of a positive number.
Also in this series:
How does it work?
It is supposed that Halley used Newton’s method to calculate square roots whenever he need it in the calculation of orbits. Later on, he is thought to have come up with his amazing idea to speed up that method, now known as Halley’s method.
Factorising out Newton’s correction term of -f(xₙ)/f’(xₙ) demonstrates how Newton’s method was improved upon by Halley:
Thus, we can see Newton’s correction multiplied by Halley’s correction:
How did Edmond Halley come up with this formula?
Edmond Halley most certainly used Taylor’s expansion for the infinitely differentiable f(x) =x² — S at any point a ∈ ℝ. We will not consider complex numbers a ∈ ℂ, as it is beyond the scope of this article.
Taylor’s expansion looks as follows for f(x) at any point a ∈ ℝ:
Evaluate the factorials to get the following:
Lastly, simplify to obtain the expansion of f(x):
In order to use Taylor’s expansion to approximate the roots of the continuous function i.e., when f(x) = 0, we need to replace the x in blue with xₙ₊₁ and the pink a with xₙ i.e., the approximations at steps n+1 and n.
The goal is for f(xₙ₊₁) = 0 i.e., xₙ₊₁ is a root of f(x). Thus, we replace the left-hand side with zero:
We can identify the following from the right-hand side:
- The first two terms give us Netwon’s method.
- The first three terms are the basis for Halley’s method.
- The more terms one uses becomes Householder’s method.
How do the first two terms i.e., expansion to 1ˢᵗ order, give us Netwon’s method? Simple. All you do is rearrange the terms to obtain the famous root finding method:
Fits like a glove! 🧤
Why would Halley chose a Taylor’s expansion to 2ⁿᵈ order instead of fitting a parabola like f(x) = x² — S with a boring straigh? The answer is simple again. Taylor’s expansion fits the parabola perfectly. How, I hear you asking? Here is how:
We can also use desmos to demonstrate that Taylor’s expansion to second order fits our function like a glove.
A bit of Mathemagic
We have Taylor’s expansion to second order. How would Edmond Halley have derived his fantastic approximation from this expansion? 🤔
One might be tempted to solve this a quaratic by letting δ = xₙ₊₁ — xₙ be the difference term.
We soon realize that we get stuck in a pickle. Good luck trying to solve it…
I’m guessing that Halley realized this, so he probably used a bit of maths magic or “mathemagic” as I like to call it, to derive his formula. What do I mean by “magic”? In mathematics, we sometimes have to make some assumptions and lateral thinking where conventional methods don’t work.
I’m literally exhausted as I finish writing this… 🥵
Fantastic, we have derived Halley’s method:
- xₙ₊₁ = xₙ — [f(xₙ)f’(xₙ)] / [[f’(xₙ)]² — ½f(xₙ)f’’(xₙ)].
![xₙ₊₁ = xₙ — [f(xₙ)f’(xₙ)] / [[f’(xₙ)]² — ½f(xₙ)f’’(xₙ)].](https://miro.medium.com/v2/resize:fit:640/format:webp/1*Tpfp3gVQZ6CvUTV4aKHnig.png%20640w)
In my article on the Babylonian / Heron’s Method, I show that if we substitute f(x) = x²— S and f’(x) = 2x into xₙ₊₁ = xₙ — f(xₙ) / f’(xₙ) i.e., Newton’s method, we obtain the simplified formula xₙ₊₁ = ½(xₙ + S / xₙ) aka Heron’s method.
For computational purpose, we are also require a simplified formula for Halley’s method.
Thus, the computational form of Halley’s method is:
- xₙ₊₁ = xₙ · (xₙ² + 3S) / (3xₙ² + S)
[…] we are using exhaustive enumeration. Notice that this method for finding the square root has nothing in common with the way of finding square roots using a pencil that you might have learned in middle school. It is often the case that the best way to solve a problem with a computer is quite different from how one would approach the problem by hand. (Guttag 2016, p. 31)
In numerical methods, we have to use the simplest starting point. It is impossible to store ‘all’ of the square numbers in this universe and, it is even more computationally costly to search for the nearest square numbers to ours. Thus, we always choose the initial guess as max(1, S) / 2.
# Initial guess:
x = max(1, S) / 2;Then we refine our guess at each step with the simplified form:
# Refinded guess:
x = x * (x**2 + 3 * S) / (3 * x**2 + S)With this compact formula, we can easily compute the first few guesses in desmos. Notice how rapidly the guesses gravitate towards the true answer for √2. Brilliant!
When do we reach a ‘good-enough’ approximation?
That depends on the threshold value tau for the approximation errors abs(x[n] - x[n-1]) and abs(x[n]² - S). Therefore, we loop until either condition is met:
abs(x[n]² - S) < tau- or
abs(x[n] - x[n-1]) < tau - or
abs(x[n] - x[n-2]) < tausince Halley’s method includes a second derivative. - In finite difference methods, this would mean that you also need to consider the difference two steps back.
Since most modern computers work with 64-bit double precision, we should not set the threshold lower than 1e-15 i.e, 15 decimal places (d.p.). For consistency, I have chosen 12 d.p. i.e., tau = 10^(-12) in all tests.
Get Aurel Nicolae’s stories in your inbox
Join Medium for free to get updates from this writer.
Does it converge?
Yes, Halley’s method is cubically convergent. The convergence error term eₙ = |xₙ — α| that is, the distance of our approximation xₙ from the true value α = √S, depends on the previous error term in a cubic polynomial fashion:
- eₙ₊₁ ≈ |eₙ|³/(4S)
- The error at step n+1 is proportional to the previous error cubed.
- Cubing the error at step n means that it is getting closer to the true value faster than the quadraticaly convergent Babylonian method with error convergence eₙ₊₁ ≈ eₙ²/(2√S), and even faster than the linearly convergent Bisection method with convergence error eₙ₊₁ ≈ eₙ/2.
- For example, if eₙ = ⅟₂ = 0.5 then eₙ² = ¼ = 0.25 and eₙ³ = ⅟₈ = 0.125. So you see, the error cubed is the smallest. Thus the distance between our guess with Halley’s method and the true value √S of gets smaller faster with every step.
- We call the power of 3 in this proportional relationship eₙ₊₁ ∼ eₙ³ the order of convergence, notated with a ‘q’. The order q = 3 implies that Halley’s method converges cubically.
- The coefficient of 1/(4S) is called the rate of convergence, notated with a ‘μ’. This rate combined with the cubed error, dictates how fast the iterative method approaches the true value.
A sequence of real numbers (xₙ) is said to converge to x ∈ ℝ (and we write lim_{n → ∞} xₙ = x) if, for any ε > 0, there exists an integer N ∈ ℕ such that, for all n > N, |xₙ — x| < ε. (Chongchitnan 2026, p. 10)
Proof of Cubic Convergence
We have already derived the compact rational form of Halley’s method as:
- xₙ₊₁ = xₙ · (xₙ² + 3S) / (3xₙ² + S)
Thus, we have the convergence error eₙ₊₁ ≈ |eₙ|³/(4S) that is, the next error is proportional to the previous error cubed and the rate of convergence μ is 1/(4S).
Comparing the rate of convergence alone to the Babylonian / Heron’s method which has a rate of convergence μ = 1/(2√S), we can see that Halley’s rate is smaller. Obvioulsy, we cannot compare rate in isolation to their order of convergence. If the Babylonian’s and Halley’s methods had the exact same order, then Halley’s would have a faster rate of convergence.
Need for Speed
As we were discussing the rate of convergence, the impact of order of convergence can be much greater. Cubic convergence i.e., order 3, implies that we approach the true value α = √S faster than the Babylonian’s method of order 2 and much faster than the Bisection method of order 1.
For example, to approximate √10 with 12 d.p. precision, the Halley’s method takes 3 iterations to converge, while the Babylonian’s method takes 5 iterations and the Bisection method converges at a snail pace of 41.
Caution: Halley’s method can be ‘unstable’ or volatile. Due the limitations of 64-bit architecture and because we use the second derivative of f(x), the guesses can oscillate between two values forever. That is why I am checking for the difference between guesses 2 steps appart |xₙ — xₙ₋₂| as well.
To quantify the number of iterations n in terms of the desired precision d i.e. the number of decimal places, the estimated number of iterations as the error gets smaller for that desired precision is:
- n ≥ log₃ |log₁₀(ρ×10⁻ᵈ) / log₁₀(ρe₀)|
- where rho ρ = 1/(2√S) i.e. μ = ρ², and e₀ is the initial error |x₀ — α|.
Derivation of the estimated number of iterations next:
Nota Bene: This estimation applies when the guess is very close to the true value α of the square root of S.
Example for S = 2 with 12 d.p. and ρ = 1/(2√2), the estimated number of iterations is 3:
In [1]: import math
In [2]: S = 2
In [3]: d = 12
In [4]: x0 = max(1, S) / 2
In [5]: e0 = abs(x0 - math.sqrt(S))
In [6]: rho = 1 / (2 * math.sqrt(S))
In [7]: math.log(9, 3)
Out[7]: 2.0
In [8]: n = math.log(abs( math.log10(rho * math.pow(10, -d)) / math.log10(rho * e0) ), 3)
In [9]: print(n)
2.4603603859954797
In [10]: math.ceil(n)
Out[10]: 3This means that it takes Halley’s method roughly 3 iterations for the approximation error of √2 to get below τ = 1e-12, importanlty as we get closer to the true answer.
Decimal Precision
How does the decimal precision improve with each iteration? The precision improves by thrice the previous precision plus a bonus or penalty of log₁₀(4S) of a decimal digit. Tripling the precision relates directly to the ‘speed’ that we approach the true value of √S with a precision d (d.p.).
Here is the derivation of the precision at every step:
The derived formula dₙ₊₁ ≈ 3dₙ + log₁₀(4S) demonstrates the tripling of the precision from step to step using Halley’s method.
Let’s see this in action in ipython:
In [1]: S = 2
In [2]: x0 = max(1, S) / 2
In [3]: print(x0)
1.0
In [4]: x1 = x0 * (x0**2 + 3*S) / (3 * x0**2 + S)
In [5]: print(x1)
1.4
In [6]: x2 = x1 * (x1**2 + 3*S) / (3 * x1**2 + S)
In [7]: print(x2)
1.4142131979695434
In [8]: x3 = x2 * (x2**2 + 3*S) / (3 * x2**2 + S)
In [9]: print(x3)
1.414213562373095
In [10]: x4 = x3 * (x3**2 + 3*S) / (3 * x3**2 + S)
In [11]: print(x4)
1.4142135623730951
In [12]: import numpy as np
In [13]: guesses = np.array([x0, x1, x2, x3, x4])
In [14]: guesses[1:] - guesses[:-1]
Out[14]: array([4.00000000e-01, 1.42131980e-02, 3.64403552e-07, 2.22044605e-16])Recall that we estimated 3 iterations. We analysie the differences:
- from 4.00000000e-01 to 1.42131980e-02, precision doubles,
- from 1.42131980e-02 to 3.64403552e-07, precision more than triples,
- and from 3.64403552e-07 to 2.22044605e-16, even though we hit the IEEE double precision limitation, the precision more than doubles.
That’s good enough for me. If you want to see precision tripling properly, you should use the multiprecision math library.
In [1]: import mpmath as mpm
In [2]: mpm.mp.dps = 60
In [4]: S = mpm.mpmathify(2)
In [5]: x0 = mpm.fdiv(max(1, S), 2)
In [6]: print(x0)
1.0
In [8]: x1 = mpm.fmul(x0, mpm.fdiv(mpm.fadd(mpm.power(x0, 2), mpm.fmul(3, S)), mpm.fadd(mpm.fmul(3, mpm.power(x0, 2)), S)))
In [9]: print(x1)
1.4
In [10]: x2 = mpm.fmul(x1, mpm.fdiv(mpm.fadd(mpm.power(x1, 2), mpm.fmul(3, S)), mpm.fadd(mpm.fmul(3, mpm.power(x1, 2)), S)))
In [11]: print(x2)
1.4142131979695431472081218274111675126903553299492385786802
In [12]: x3 = mpm.fmul(x2, mpm.fdiv(mpm.fadd(mpm.power(x2, 2), mpm.fmul(3, S)), mpm.fadd(mpm.fmul(3, mpm.power(x2, 2)), S)))
In [13]: print(x3)
1.4142135623730950487956400807542599463542382401452485875064
In [14]: x4 = mpm.fmul(x3, mpm.fdiv(mpm.fadd(mpm.power(x3, 2), mpm.fmul(3, S)), mpm.fadd(mpm.fmul(3, mpm.power(x3, 2)), S)))
In [15]: print(x4)
1.41421356237309504880168872420969807856967187537694807317668
In [16]: import numpy as np
In [17]: vfsub = np.vectorize(mpm.fsub)
In [18]: guesses = np.array([x0, x1, x2, x3, x4])
In [19]: vfsub(guesses[1:], guesses[:-1])
Out[19]:
array([mpf('0.400000000000000000000000000000000000000000000000000000000000062'),
mpf('0.0142131979695431472081218274111675126903553299492385786802029042'),
mpf('0.000000364403551901587518253343092433663882910196010008826198356773409'),
mpf('6.04864345543813221543363523169948567027845529940506684702000768e-21')],
dtype=object)
In [20]: Notice that the difference between the 3rd and 4th guesses is triple the difference before of 3.64403552e-07 i.e.,
- 6.048643455438132215433635231699485670278455299405066847(…)e-21.
Python implementation
The following is my implementation of xₙ₊₁ = xₙ · (xₙ² + 3S) / (3xₙ² + S). The stopping criteria was challenging because of the IEEE double precision limitation. The algorithm arrives at 64-bit limit very quickly.
Thus, if we carry on the algorithm from above, note that we oscillate between two values infinitely:
In [1]: S = 2
In [2]: x0 = max(1, S) / 2
In [3]: x1 = x0 * (x0**2 + 3*S) / (3 * x0**2 + S)
In [4]: x2 = x1 * (x1**2 + 3*S) / (3 * x1**2 + S)
In [5]: x3 = x2 * (x2**2 + 3*S) / (3 * x2**2 + S)
In [6]: print(x3)
1.414213562373095
In [7]: x4 = x3 * (x3**2 + 3*S) / (3 * x3**2 + S)
In [8]: print(x4)
1.4142135623730951
In [9]: x5 = x4 * (x4**2 + 3*S) / (3 * x4**2 + S)
In [10]: print(x5)
1.414213562373095
In [11]: x6 = x5 * (x5**2 + 3*S) / (3 * x5**2 + S)
In [12]: print(x6)
1.4142135623730951Thus, I had to carefully check for duplicates in the list of guesses to avoid an infinite loop i.e., look two steps behind abs(x[n] - x[n-2]) < tau.
I am also using a custom utility function to print the results to a file. You can remove that functionality if you do not want it.
Disclaimer: My implementation might not be the best and most efficient. I encourage you to rewrite it in your own way.
Analysis
I ran tests using Halley’s method for various numbers with the same threshold 1e-12 or 12 d.p. and the following is a summary of my tests:
╭─────────────────────────────────┬───────────────────────────────────┬──────────────────┬─────────────────────╮
│ Number │ Approximation │ Actual Guesses │ Estimated Guesses │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 1 │ 0.9999999999997379 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 2 │ 1.4142135623730949 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 3 │ 1.7320508075688772 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 4 │ 2.0000000000000000 │ 0 │ 0 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 5 │ 2.2360679774997902 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 7 │ 2.6457513110645903 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 8 │ 2.8284271247461898 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 9 │ 3.0000000000000000 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 10 │ 3.1622776601683791 │ 3 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 25 │ 5.0000000000000000 │ 4 │ 5 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 37 │ 6.0827625302982211 │ 4 │ 7 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 49 │ 7.0000000000000000 │ 4 │ 5 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 73 │ 8.5440037453175304 │ 4 │ 4 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 81 │ 9.0000000000000018 │ 4 │ 4 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 100 │ 10.0000000000000000 │ 5 │ 4 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 121 │ 11.0000000000000000 │ 5 │ 4 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 1000 │ 31.6227766016837926 │ 6 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 10000 │ 99.9999999999999716 │ 8 │ 3 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 100000 │ 316.2277660168379043 │ 9 │ 2 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 987654 │ 993.8078285060950066 │ 10 │ 2 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 987543201 │ 31425.2000948283530306 │ 15 │ 2 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 673003460632639326586069439443 │ 820367881278051.7500000000000000 │ 37 │ 1 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 1000000000000000000000000000000 │ 1000000000000000.0000000000000000 │ 34 │ 1 │
├─────────────────────────────────┼───────────────────────────────────┼──────────────────┼─────────────────────┤
│ 1367999732000000000000000000071 │ 1169615206809487.5000000000000000 │ 36 │ 1 │
╰─────────────────────────────────┴───────────────────────────────────┴──────────────────┴─────────────────────╯For the number S under 100 and in the lower one hundreds, the actual number of iterations matches quite well the number of estimated guesses. Looking at the square root of 2 once again, we can see how the guess play out:
(.venv) C:\Users\...\Math_Square_Root>python c3Halley\halley_sqrt.py
Enter positive number: 2
Enter precision between 0 and 15 (Default 15 d.p.): 12
iteration: 0 guess: 1.0 guess²: 1.0
iteration: 1 guess: 1.4 guess²: 1.9599999999999997
iteration: 2 guess: 1.4142131979695434 guess²: 1.9999989693112428
Number of guesses: 3, Threshold: 1e-12, Estimate iterations: 3.
Approximation for the square root of 2.0 is: 1.414213562373095
Compared to math.sqrt(2.0): 1.4142135623730951This is because the initial error term e₀ = |x₀ — α| is quite small for these ‘low’ numbers. However, as we test numbers much bigger than one hundred, that is when we can see what is ment by:
This estimation applies when the guess is very close to the true value α of the square root of S.
For S = 100,000 as an example, the estimated number of guesses is 2 but in practice it takes us 9 guesses:
(.venv) C:\Users\...\Math_Square_Root>python c3Halley\halley_sqrt.py
Enter positive number: 100_000
Enter precision between 0 and 15 (Default 15 d.p.): 12
iteration: 0 guess: 50000.0 guess²: 2500000000.0
iteration: 1 guess: 16668.444420741056 guess²: 277837039.4073336
iteration: 2 guess: 5561.480265041226 guess²: 30930062.73844302
iteration: 3 guess: 1869.7925037381722 guess²: 3496124.0070354627
iteration: 4 guess: 670.3546302931796 guess²: 449375.3303555056
iteration: 5 guess: 346.8946939909042 guess²: 120335.92871904308
iteration: 6 guess: 316.29032685507946 guess²: 100039.570862093
iteration: 7 guess: 316.2277660174499 guess²: 100000.00000038702
iteration: 8 guess: 316.22776601683796 guess²: 100000.00000000001
Number of guesses: 9, Threshold: 1e-12, Estimate iterations: 2.
Approximation for the square root of 100000.0 is: 316.2277660168379
Compared to math.sqrt(100000.0): 316.22776601683796This means that the last 2 iterations are very close as we can see above in the output.
Conclusion
Edmond Halley may be famous for calculating the orbits of comets, but his method of calculating square roots is the true “shooting star” in my eyes (excuse the corny pun 🍿). Improving on Newton’s method by tripling the precision is an amazing achievement.
Download
- Here is the link to my GitHub repo.
Appendix: Multiprecision
Here is my multiple precision implementation of the Babylonian aka Heron’s method using the mpmath library accompanied by my disclaimer that there are probably better and more efficient implementations out there.
Also in this series:
