How Does a Calculator Find a Square Root?

Press Start again, then Next step slowly, and count the green digits at each step: 1, 1, 3, 6, 12, 24. After the first steps, each step about doubles them. Six steps, and all 25 digits shown are right.

Now give it a bad first guess. Choose 5,476 (one ropani, the land measure, is 5,476 square feet) and First guess: 1.

So the method has two speeds. Far from the answer it is slow, about halving the guess each step. Close to the answer it is astonishingly fast. A good first guess skips the slow part: try First guess: 100 for 5,476.

The average of two sides

If your guess g for √2 is too big, then 2 ÷ g is too small, because the two multiply to make 2. For g = 1.5, 2 ÷ 1.5 = 1.333…, and indeed 1.333… < 1.414… < 1.5. The true root is somewhere between them, so a sensible next guess is halfway between them, their average:

$$\text{new guess} = \frac{g + \dfrac{n}{g}}{2}$$

That is the whole method: guess, average the guess with the number divided by the guess, repeat. Each repeat is an iteration, and a method that improves an answer by repeating a step is an iterative method. It is often called the Babylonian method or Heron’s method, after Heron of Alexandria in Egypt, who described it in his Metrica about 2,000 years ago.

Why the digits double

Suppose the guess is too big by a small amount, which we can call e: the guess is √n + e. Working the average out (it takes a few lines of algebra) gives a new guess that is too big by about

$$\frac{e^2}{2\sqrt{n}}$$

The new error is the old error squared (then divided by about 2√n, a number that doesn’t change from step to step). Squaring a small number makes it far smaller: if the error is 0.001 (2 correct decimals), its square is 0.000001 (5 or 6 correct decimals). An error of 10⁻⁶ becomes about 10⁻¹². The number of zeros after the point doubles, and so do the correct digits. A method whose error is squared each step is said to converge quadratically.

Far from the root it is a different story. If the guess is much too big, n ÷ g is tiny, and the average of a big number and a tiny one is about half the big one. That is the halving you saw starting from 1 for 5,476.

You will build it in four steps, each one running here on the page.

Step 1: the method in a loop

n = 2guess = 1for step in range(1, 7):    guess = (guess + n / guess) / 2    print(step, guess)

1 1.52 1.41666666666666653 1.41421568627450974 1.41421356237468995 1.4142135623730956 1.414213562373095

After five steps the guess stops changing. That isn’t because the method has finished: a Python float keeps only about 16 significant digits (Grade 11, Floating point), and all 16 are already as right as a float can be. To watch the doubling go on, the program needs more digits.

Step 2: watch the digits double

Python’s decimal module does arithmetic with as many digits as you ask for. This program works with 60, and measures each guess against the true root (from the decimal module’s own square root, used here only as a ruler). gap.adjusted() gives the place of the gap’s first digit as a power of ten: −6 for a gap of 0.0000021, whose first digit is in the millionths. A guess that far off is right to 5 decimal places, so the program counts −(−6) − 1 = 5.

from decimal import Decimal, getcontextgetcontext().prec = 60          # work with 60 digits, not 16 n = Decimal(2)true_root = n.sqrt()            # decimal's own square root, as a rulerguess = Decimal(1)for step in range(1, 8):    guess = (guess + n / guess) / 2    gap = abs(guess - true_root)    correct = -gap.adjusted() - 1 if gap else 60    print(f'Step {step}: {str(guess)[:32]}  about {correct} correct decimals')

Step 1: 1.5  about 1 correct decimalsStep 2: 1.416666666666666666666666666666  about 2 correct decimalsStep 3: 1.414215686274509803921568627450  about 5 correct decimalsStep 4: 1.414213562374689910626295578890  about 11 correct decimalsStep 5: 1.414213562373095048801689623502  about 24 correct decimalsStep 6: 1.414213562373095048801688724209  about 48 correct decimalsStep 7: 1.414213562373095048801688724209  about 60 correct decimals

There it is, made by your program: 2, 5, 11, 24, 48. Each step doubles the correct decimals (and a little more). Six steps of adding and dividing give √2 to 48 decimal places; at step 7 it runs out of the 60 digits it was given.

Step 3: a function that knows when to stop

A real square-root function can’t be told how many steps to take: it must stop by itself. Here is a neat rule. If the first guess is at or above the root, every guess after it is too (an average of two sides is never below the square’s side), and each one is smaller than the last. With floats, the guesses eventually stop getting smaller, because there are no more digits to improve. So: keep going while the guess shrinks. Starting from n itself (or 1, for a number less than 1) puts the first guess above the root.

import math def my_sqrt(n):    guess = max(n, 1)           # never below the root    steps = 0    while True:        better = (guess + n / guess) / 2        if better >= guess:     # it has stopped getting smaller            return guess, steps        guess = better        steps = steps + 1 for n in [2, 10, 5476, 0.25, 10 ** 10, 10 ** 100]:    root, steps = my_sqrt(n)    print(f'√{n:g}: {root!r:20} {steps:3} steps, math.sqrt: {math.sqrt(n)!r}')

√2: 1.414213562373095      5 steps, math.sqrt: 1.4142135623730951√10: 3.162277660168379      6 steps, math.sqrt: 3.1622776601683795√5476: 74.0                  11 steps, math.sqrt: 74.0√0.25: 0.5                    6 steps, math.sqrt: 0.5√1e+10: 100000.0              21 steps, math.sqrt: 100000.0√1e+100: 1e+50                171 steps, math.sqrt: 1e+50

Two things to notice.

• For √2 and √10 your answer differs from Python’s math.sqrt in the last digit. Each float division and average is rounded, so your last guess can end one tiny step (in the 16th digit) below the closest float. math.sqrt is built to give the closest float there is, every time.
• The step counts: 5 or 6 for small numbers, but 171 for 10¹⁰⁰. From a first guess of 10¹⁰⁰, the guesses halve about 165 times before they get near 10⁵⁰. The fast part of the method is only ever the last 5 or 6 steps.

Step 4: a good first guess

So the secret of a fast square root is the first guess. A computer can make a very good one almost for free. A float is stored as a number between 0.5 and 1 times a power of 2 (Grade 11’s Floating point); math.frexp(n) gives the two parts. The square root of 2 to some power is 2 to half that power, so halving the power (rounded up) gives a power of 2 at or just above √n: never more than twice too big.

import math def first_guess(n):    m, e = math.frexp(n)        # n = m × 2**e, with m from 0.5 to 1    return 2.0 ** ((e + 1) // 2)  # a power of 2 at or above √n def my_sqrt(n):    guess = first_guess(n)    steps = 0    while True:        better = (guess + n / guess) / 2        if better >= guess:            return guess, steps        guess = better        steps = steps + 1 for n in [2, 10, 5476, 0.25, 10 ** 10, 10 ** 100]:    root, steps = my_sqrt(n)    print(f'√{n:g}: first guess {first_guess(n):g}, {steps} steps, {root!r}')

√2: first guess 2, 5 steps, 1.414213562373095√10: first guess 4, 5 steps, 3.162277660168379√5476: first guess 128, 5 steps, 74.0√0.25: first guess 1, 6 steps, 0.5√1e+10: first guess 131072, 5 steps, 100000.0√1e+100: first guess 1.87072e+50, 5 steps, 1e+50

From 171 steps to 5, for any number at all. The better the first guess, the fewer steps: start with half the digits right and a single step makes them all right.

• a.Each step doubles the correct digits (Step 2): 2, 5, 11, 24, 48. A method that adds one digit a step would need 48 steps; this one needs 6.
• b.The slow part is getting close (Step 3): from a bad guess the method only halves its error, 165 times for 10¹⁰⁰.
• c.A good first guess removes the slow part (Step 4): 5 or 6 steps for any number.

Now look back 3,800 years. The method’s first three steps from 1, done with fractions instead of decimals, give 3/2, then 17/12, then 577/408. This program works them out exactly with Python’s fractions module, writes 577/408 in base 60 (the Babylonians’ number system: one place for sixtieths, the next for sixtieths of sixtieths, as minutes and seconds still are), and compares it with two of the oldest values of √2 we know:

from fractions import Fraction guess = Fraction(1)for step in range(1, 4):    guess = (guess + 2 / guess) / 2    print(step, guess, float(guess)) # Write 577/408 in base 60, as a Babylonian scribe would.x = guessplaces = []for i in range(5):    whole = int(x)    places.append(whole)    x = (x - whole) * 60print('In base 60:', places)tablet = 1 + Fraction(24, 60) + Fraction(51, 60 ** 2) + Fraction(10, 60 ** 3)print('The tablet, 1;24,51,10:', float(tablet)) # The Sulba Sutra's rule: add a third, a fourth of that, take away a 34th of that.sulba = 1 + Fraction(1, 3) + Fraction(1, 3 * 4) - Fraction(1, 3 * 4 * 34)print('The Sulba Sutra:', sulba)

1 3/2 1.52 17/12 1.41666666666666673 577/408 1.4142156862745099In base 60: [1, 24, 51, 10, 35]The tablet, 1;24,51,10: 1.414212962962963The Sulba Sutra: 577/408

Your program has found something remarkable. The value of √2 in an Indian Sulba Sutra is exactly the method’s third step, 577/408. And a Babylonian clay tablet’s value, 1;24,51,10 in base 60, is the same 577/408 written in base 60 and cut off after three places (the next place, 35, is dropped). Neither text tells us how its value was found, so this doesn’t prove that either one used the averaging method. It does show that the oldest good values of √2 sit exactly on the path the method takes.

Cube roots

A cube-shaped water tank that holds 5,000 liters holds 5 cubic meters, so each side is ∛5 meters. The same kind of improvement works: if g is a guess, n ÷ g² is another, and a weighted average of them, (2g + n ÷ g²) ÷ 3, is better than either.

n = 5           # 5,000 liters is 5 cubic metersguess = 2for step in range(1, 7):    guess = (2 * guess + n / (guess * guess)) / 3    print(step, guess, '  cubed:', guess ** 3)

1 1.75   cubed: 5.3593752 1.7108843537414966   cubed: 5.00797282909646853 1.7099764289169748   cubed: 5.0000042302385724 1.709975946676833   cubed: 5.0000000000011935 1.709975946676697   cubed: 5.0000000000000016 1.709975946676697   cubed: 5.000000000000001

Each side is 1.71 meters. Look at the “cubed” column: the error goes 0.36, 0.008, 0.000004, 0.000000000001, doubling its zeros each time, just like the square root. For any root, the kth root, the rule is ((k − 1)g + n ÷ g^(k−1)) ÷ k: that is Challenge ★.

Newton’s method: any equation at all

Where did those rules come from? Finding √2 means finding the x where x² − 2 = 0: where the curve y = x² − 2 crosses the axis. Isaac Newton’s idea (1669, published later; Joseph Raphson described a simpler version in 1690) works for any curve: at your guess, follow the straight line that touches the curve there, and see where that line crosses the axis. That is your next guess.

The line’s slope at x is called the curve’s derivative, from calculus, which you meet in Grade 11 and 12 math. For x² − 2 it is 2x; for x⁵ − x − 1 it is 5x⁴ − 1 (the slope of xᵏ is k times x to the power k − 1). A line through the point (g, f(g)) with slope s crosses the axis at g − f(g) ÷ s, so the program needs the curve and its slope, each as a function, and passes them to newton like any other values (a function’s name without brackets is the function itself):

def newton(f, slope, guess, steps):    for step in range(1, steps + 1):        guess = guess - f(guess) / slope(guess)        print(step, guess)    return guess def square(x):    return x * x - 2 def square_slope(x):    return 2 * x def fifth(x):    return x ** 5 - x - 1 def fifth_slope(x):    return 5 * x ** 4 - 1 print('x * x - 2 = 0')newton(square, square_slope, 1, 5)print('x**5 - x - 1 = 0')x = newton(fifth, fifth_slope, 1, 6)print('check:', x ** 5 - x - 1)

x * x - 2 = 01 1.52 1.41666666666666673 1.41421568627450994 1.41421356237468995 1.4142135623730951x**5 - x - 1 = 01 1.252 1.17845939351690483 1.167537389396114 1.16730408282300835 1.16730397826143966 1.1673039782614187check: 6.661338147750939e-16

For x² − 2, Newton’s method gives exactly the averaging guesses, 1.5, 1.41666…: g − (g² − 2) ÷ 2g works out to (g + 2 ÷ g) ÷ 2. The ancient method is Newton’s method for square roots.

The second equation is the interesting one. No formula exists for the solution of x⁵ − x − 1 = 0: there is a formula for any equation with x², another (long) one for x³ and x⁴, but Paolo Ruffini (1799, with a gap in his proof) and Niels Henrik Abel (1824) proved that no formula of that kind can solve every equation with x⁵, and this one is a known example that no such formula solves. Yet six steps of your program give x = 1.1673039782614187, so close that x⁵ − x − 1 is 0.0000000000000007. A formula isn’t needed when you can iterate.

An open question: are √2’s digits fair?

Python’s math.isqrt(n) gives the whole-number square root of a whole number, exactly, however big. The whole-number square root of 2 × 10²⁰⁰⁰ is √2 × 10¹⁰⁰⁰ without its fraction: the first 1,001 digits of √2.

import mathdigits = str(math.isqrt(2 * 10 ** 2000))   # √2 × 10**1000, as a whole numberprint(digits[:40], '...')for d in '0123456789':    print(d, digits.count(d))

1414213562373095048801688724209698078569 ...0 1081 992 1093 824 1005 1046 907 1048 1139 92

Each digit comes about 100 times in 1,001, as if they had been drawn at random. Mathematicians think that in √2 every digit, every pair of digits and every longer string of them comes equally often in the long run (a normal number), and every test of its digits so far agrees. Nobody has proved it. It is not even known whether the digit 7 keeps appearing for ever.

Three challenges, each checked on the page: write your function in the program above the checker (press Edit), then press Check my program. The examples show what your function should give; there are more tests you can’t see.

★ Any root

Write nth_root(n, k): the kth root of a positive number n (the number that, multiplied by itself k times, makes n), by the method, rounded to 6 decimal places. nth_root(5, 3) is ∛5. (n ** (1 / k) would pass too, but the point is to build it.)

def nth_root(n, k):    # Give back the kth root of n, rounded to 6 decimal places.    return 0 print(nth_root(5, 3))

★★ The exact whole-number root

Write isqrt(n): the biggest whole number whose square is at most n, for any whole number n from 0 up. isqrt(99) is 9, isqrt(100) is 10. Don’t use floats (or math.isqrt): for big numbers a float’s 16 digits aren’t enough, and one of the tests shows it.

def isqrt(n):    # Give back the biggest k with k * k <= n.    return 0 print(isqrt(99), isqrt(100))

★★★ What rate is the bank charging?

A bank in Nepal works out a loan’s monthly payment (its EMI, equated monthly installment) with the formula in emi() below: borrow 10 lakh rupees at 12% a year for 60 months and you pay Rs. 22,244.45 a month. But if a lender tells you only the payment, what yearly rate are you paying? There’s no formula to turn emi() round. Write loan_rate(amount, payment, months): the yearly rate in percent, rounded to 2 decimal places.

def emi(amount, rate, months):    r = rate / 1200                 # the yearly percent, as a monthly fraction    grow = (1 + r) ** months    return amount * r * grow / (grow - 1) def loan_rate(amount, payment, months):    # Give back the yearly rate in percent that makes emi() equal payment.    return 0 print(round(emi(1000000, 12, 60), 2))print(loan_rate(1000000, 22244.45, 60))

• Grade 8, Decisions and Loops in Python: while loops, and loops that never end.
• Grade 11, From Code to Circuits: why a float keeps only about 16 digits, and how it is stored as a fraction times a power of 2.
• How Did a Rounding Error Stop a Missile Defense?: when the last digits of a calculation really matter.
• Why Can’t Anyone Forecast the Weather a Month Ahead?: another program that repeats one step, where tiny differences grow instead of shrinking.