Big Square Roots

by Michael Gilleland

Introduction

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.

Algorithm

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:

1. 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.
2. Compute result = ((n/g) + g)/2. Let g be the result just computed.
3. 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.

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

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.

Source Code

//----------------------------------------------------------
// 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 ());
    }

  }

}

Sample Output

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

Links

• 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