Algorithm

Source Code

Sample Output

Links

The base Java class java.lang.Math contains a useful grab bag of mathematical goodies:

- Trigonometric functions, such as sine, cosine, tangent, arcsine, arccosine, arctangent
- Constants, such as PI and E
- Other miscellaneous functions such as square root, ceiling, floor, etc.

In the absence of such a class, I needed a simple way to find square roots of big numbers in Java, and this paper presents a straightforward implementation of what turns out to be a very old technique.

There has to be a simpler way, I thought, and indeed there is. You can find it at the Argonne National Laboratory's web site, in the archives of the "Ask a Scientist" service for K-12 educators and students. Someone asked a scientist "How do you find square roots?", and R.C. Winther answered with a lucid description, which I'll paraphrase here:

- Start with an initial guess (g) of what the square root of some number (call it n) might be. The initial guess doesn't even have to be close. For simplicity's sake, let's always choose g = 1.
- Compute result = ((n/g) + g)/2. Let g be the result just computed.
- Repeat step 2 until the last two results obtained are the same.

Let's look at some examples. First, let's find the square root of 9 by this method. In our calculations, we'll round off everything to three decimal places.

n | g | n/g | (n/g) + g | ((n/g) + g)/2 |
---|---|---|---|---|

9 | 1 | 9 | 10 | 5 |

9 | 5 | 1.8 | 6.8 | 3.4 |

9 | 3.4 | 2.65 | 6.05 | 3.025 |

9 | 3.025 | 2.975 | 6 | 3 |

9 | 3 | 3 | 6 | 3 |

Next, let's find the square root of 30 by this method. This time, we'll round off everything to eight decimal places.

n | g | n/g | (n/g) + g | ((n/g) + g)/2 |
---|---|---|---|---|

30 | 1 | 30 | 31 | 15.5 |

30 | 15.5 | 1.93548371 | 17.43548387 | 8.71774194 |

30 | 8.71774194 | 3.44125809 | 12.15900003 | 6.07950002 |

30 | 6.07950002 | 4.93461632 | 11.01411634 | 5.50705817 |

30 | 5.50705817 | 5.44755459 | 10.95461276 | 5.47730638 |

30 | 5.47730638 | 5.47714477 | 10.95445115 | 5.47722558 |

30 | 5.47722558 | 5.47722557 | 10.95445115 | 5.47722558 |

Professional mathematicians are of course familiar with this algorithm,
which is known as Newton's method, but its pedigree is in fact much
older than Newton.
In the first century AD, the Greek mathematician
Heron of Alexandria used essentially the
same method in book 1, chapter 8, of his *Metrica* to determine
the square root of 720. See Jean-Luc Chambert et al., *op. cit.*,
pp. 202-203, and Thomas Heath, *A History of Greek Mathematics, Volume II:
From Aristarchus to Diophantus* (NY: Dover, 1981), pp. 323-326.

There is evidence that this method goes even further back, as
far as the Babylonians.
In the Yale University Babylonian Collection, there is a cuneiform
tablet (dated ca. 2000-1700 BC and catalogued as YBC 7289) that seems
to use this method to compute the square root of 2.
See H.L. Resnikoff and R.O. Wells, Jr.,
*Mathematics in Civilization* (NY: Dover, 1984), pp. 76-79,
and also the following links:

- Bill Casselman, YBC 7289, with splendid full-size color photographs
- Franz Gnaedinger, Babylonian Mathematics (YBC 7289, Plimpton 322)
- Duncan J. Melville, YBC 7289
- J.J. O'Connor and E.F. Robertson, Babylonian Pythagoras, a very clear explanation

There's one obvious refinement we can make, and that is in the choice of our initial guess. We've used 1 as the initial guess so far in the examples, but let's try to do a little better. We know that the square root of 100 is 10, so a better initial guess for the square root of any 3 or 4 digit number is 10. Likewise, we know that the square root of 10,000 is 100, so a better initial guess for the square root of any 5 or 6 digit number is 100. And so forth.

//---------------------------------------------------------- // Compute square root of large numbers using Heron's method //---------------------------------------------------------- import java.math.*; public class BigSquareRoot { private static BigDecimal ZERO = new BigDecimal ("0"); private static BigDecimal ONE = new BigDecimal ("1"); private static BigDecimal TWO = new BigDecimal ("2"); public static final int DEFAULT_MAX_ITERATIONS = 50; public static final int DEFAULT_SCALE = 10; private BigDecimal error; private int iterations; private boolean traceFlag; private int scale = DEFAULT_SCALE; private int maxIterations = DEFAULT_MAX_ITERATIONS; //--------------------------------------- // The error is the original number minus // (sqrt * sqrt). If the original number // was a perfect square, the error is 0. //--------------------------------------- public BigDecimal getError () { return error; } //------------------------------------------------------------- // Number of iterations performed when square root was computed //------------------------------------------------------------- public int getIterations () { return iterations; } //----------- // Trace flag //----------- public boolean getTraceFlag () { return traceFlag; } public void setTraceFlag (boolean flag) { traceFlag = flag; } //------ // Scale //------ public int getScale () { return scale; } public void setScale (int scale) { this.scale = scale; } //------------------- // Maximum iterations //------------------- public int getMaxIterations () { return maxIterations; } public void setMaxIterations (int maxIterations) { this.maxIterations = maxIterations; } //-------------------------- // Get initial approximation //-------------------------- private static BigDecimal getInitialApproximation (BigDecimal n) { BigInteger integerPart = n.toBigInteger (); int length = integerPart.toString ().length (); if ((length % 2) == 0) { length--; } length /= 2; BigDecimal guess = ONE.movePointRight (length); return guess; } //---------------- // Get square root //---------------- public BigDecimal get (BigInteger n) { return get (new BigDecimal (n)); } public BigDecimal get (BigDecimal n) { // Make sure n is a positive number if (n.compareTo (ZERO) <= 0) { throw new IllegalArgumentException (); } BigDecimal initialGuess = getInitialApproximation (n); trace ("Initial guess " + initialGuess.toString ()); BigDecimal lastGuess = ZERO; BigDecimal guess = new BigDecimal (initialGuess.toString ()); // Iterate iterations = 0; boolean more = true; while (more) { lastGuess = guess; guess = n.divide(guess, scale, BigDecimal.ROUND_HALF_UP); guess = guess.add(lastGuess); guess = guess.divide (TWO, scale, BigDecimal.ROUND_HALF_UP); trace ("Next guess " + guess.toString ()); error = n.subtract (guess.multiply (guess)); if (++iterations >= maxIterations) { more = false; } else if (lastGuess.equals (guess)) { more = error.abs ().compareTo (ONE) >= 0; } } return guess; } //------ // Trace //------ private void trace (String s) { if (traceFlag) { System.out.println (s); } } //---------------------- // Get random BigInteger //---------------------- public static BigInteger getRandomBigInteger (int nDigits) { StringBuffer sb = new StringBuffer (); java.util.Random r = new java.util.Random (); for (int i = 0; i < nDigits; i++) { sb.append (r.nextInt (10)); } return new BigInteger (sb.toString ()); } //----- // Test //----- public static void main (String[] args) { BigInteger n; BigDecimal sqrt; BigSquareRoot app = new BigSquareRoot (); app.setTraceFlag (true); // Generate a random big integer with a hundred digits n = BigSquareRoot.getRandomBigInteger (100); // Build an array of test numbers String testNums[] = {"9", "30", "720", "1024", n.toString ()}; for (int i = 0; i < testNums.length; i++) { n = new BigInteger (testNums[i]); if (i > 0) { System.out.println ("----------------------------"); } System.out.println ("Computing the square root of"); System.out.println (n.toString ()); int length = n.toString ().length (); if (length > 20) { app.setScale (length / 2); } sqrt = app.get (n); System.out.println ("Iterations " + app.getIterations ()); System.out.println ("Sqrt " + sqrt.toString ()); System.out.println (sqrt.multiply (sqrt).toString ()); System.out.println (n.toString ()); System.out.println ("Error " + app.getError ().toString ()); } } }

Computing the square root of 9 Initial guess 1 Next guess 5.0000000000 Next guess 3.4000000000 Next guess 3.0235294118 Next guess 3.0000915542 Next guess 3.0000000014 Next guess 3.0000000000 Next guess 3.0000000000 Iterations 7 Sqrt 3.0000000000 9.00000000000000000000 9 Error 0.00000000000000000000 ---------------------------- Computing the square root of 30 Initial guess 1 Next guess 15.5000000000 Next guess 8.7177419355 Next guess 6.0795000150 Next guess 5.5070581682 Next guess 5.4773063790 Next guess 5.4772255757 Next guess 5.4772255751 Next guess 5.4772255751 Iterations 8 Sqrt 5.4772255751 30.00000000052952574001 30 Error -0.00000000052952574001 ---------------------------- Computing the square root of 720 Initial guess 10 Next guess 41.0000000000 Next guess 29.2804878049 Next guess 26.9351210370 Next guess 26.8330100187 Next guess 26.8328157307 Next guess 26.8328157300 Next guess 26.8328157300 Iterations 7 Sqrt 26.8328157300 720.00000000013543290000 720 Error -0.00000000013543290000 ---------------------------- Computing the square root of 1024 Initial guess 10 Next guess 56.2000000000 Next guess 37.2103202847 Next guess 32.3647837114 Next guess 32.0020557400 Next guess 32.0000000661 Next guess 32.0000000000 Next guess 32.0000000000 Iterations 7 Sqrt 32.0000000000 1024.00000000000000000000 1024 Error 0.00000000000000000000 ---------------------------- Computing the square root of 4091003901585987357290452999329796564377935868371155936881162018216821980804517176141145508100318254 Initial guess 10000000000000000000000000000000000000000000000000 Next guess 209550195079299367864522649966489828218896793418557.79684405810091084109904022585880705727540501591270 Next guess 114536491223959647614245503583340019583427857317171.39481623120113435750889678042997539315029228108131 Next guess 75127199810195208191946961814969421616682217230416.51246197944058126421950875584398655254556449079256 Next guess 64790782017046982096534183561372152130639527179695.34728948529268210319975000563917742728801615162872 Next guess 63966270184743427189126780349337659174048381406572.22692340512880828760541344926488086150358076940314 Next guess 63960956292283814083342094802488714859091297948561.14711574274512116433752843287967792579630701229591 Next guess 63960956071544047910243450536455793538763694049291.28632379657071976628102300495529141197728862876571 Next guess 63960956071544047529338853709828068953366630815793.34341832696373266833349284427958104709380671225364 Next guess 63960956071544047529338853709828067819172266042961.42582721313752717361185018697161569618604233390746 Next guess 63960956071544047529338853709828067819172266042961.42582721313752717360179407668111645754167252074208 Next guess 63960956071544047529338853709828067819172266042961.42582721313752717360179407668111645754167252074208 Iterations 11 Sqrt 63960956071544047529338853709828067819172266042961.42582721313752717360179407668111645754167252074208 4091003901585987357290452999329796564377935868371155936881162018216821980804517176141145508100318253.9941843294265790729454364307536302861737946912377682142782519700429387214937782005481990878338827264 4091003901585987357290452999329796564377935868371155936881162018216821980804517176141145508100318254 Error 0.0058156705734209270545635692463697138262053087622317857217480299570612785062217994518009121661172736

- Harold Abelson and Gerald Jay Sussman,
Square Roots by Newton's Method,
from their book
*Structure and Interpretation of Computer Programs* - Kevin Brown, Ancient Square Roots
- Kevin Brown, Archimedes and the Square Root of 3
- Paul Hsieh, How to Calculate Square Roots
- Mohamed A. Khamsi and Helmut Knaust, The Newton-Raphson Method
- Romuauld Ireneus 'Scibor-Marchocki, Square Roots -- How to Find Them
- Ching-Kuang Shene, Computing the Square Root of a Positive Number