GFA BASIC 32‎ > ‎How To‎ > ‎

04 Calculating Square Roots

posted Jul 18, 2015, 3:15 PM by Troy Cheek   [ updated Jul 18, 2015, 3:25 PM ]
Back in college, I took a Computer Sciences course called something like Numerical Methods.  The long and the short of it was that this class taught you how to use a computer program to guess at answers to questions you weren't smart enough to answer on your own.  Since you have to thoroughly understand how to solve a problem in order to program a computer to solve it for you, this knocked the difficulty down a notch or three.  You only had to understand how to solve the problem in the most general of terms, then let the computer make multiple guesses, fine tuning the answer until it was "close enough."  I grew to love "close enough."

y = 10x^3 - 111x^2 + 14x - 72  Given that equation, let's say you wanted to find some value for x such that y=42.  Once upon a time, I could do that calculation.  Heck, at one time, I could probably do it in my head.  It's a simple matter of algebra.  Or maybe calculus.  Differential equations?  Sorry.  Long time, no math.

What Numerical Methods taught us was to set up a program that tries a guess for the value of x and sees what the resultant value of y comes out.  Using that result, we refine the guess for x and try again.  Assuming the results are what we call "convergent" (meaning they're getting closer to the answer instead of farther away or "divergent") we will eventually get an answer which, if not exact, is close enough for jazz.  Either the solution we get is arbitrarily close to y or the last couple guesses for x were close enough to each other that there really isn't much point continuing.

The solution for the above equation is left as an exercise for the reader, mainly because I have no idea how to do it.  I'm going to focus on calculating square roots.  A square root is a number which, when multiplied by itself, equals a target number you're aiming for.

Now, GFA BASIC 32 has a perfectly good function for calculating square roots.  It's called Sqr(), and has an identical sister called Sqrt().  This is different from some languages where Sqr() actually gives you the square of a number instead of the square root.  In GB32, that function is called Square().  The main point of this article is that we're going to ignore all these built-in functions and program one for ourselves.

Now, there are a couple of "approved" ways to calculate a square root.  These involve fancy formulas or methods that nobody remembers and you have to look up every time you want to use them.  If you're going to do that, you might as well look up the square root itself or at least look up the syntax for the Sqrt() command in your language of choice.  Instead, let's focus on a method that uses our basic understanding of what a square root is, which is a number that when multiplied by itself equals our target number.

Two common methods for doing this are called "bracketing" and "successive approximation."  Bracketing is when you know the answer is in some range, say higher than one number but definitely lower than another, and you keep refining the low and high values, bracketing them in, until you get the right answer or your lower and upper brackets are so close as to be effectively the same number, which is your answer.

If you're looking for the square root of 25, well that's easy.  It's easy to remember that 5 time 5 equals 25.  And 6 times 6 is 36.  The square root of 30 is tougher, but since we know the square roots of 25 and 36, it's easy to see that the square root of 30 is somewhere between 5 and 6.  To get our initial low and high bracket numbers for square roots, we can do something like this:

Dim lo, hi, med, tar, tmp As Double
Dim i%
tar = 30
For i% = 0 To tar + 1
  If i% * i% < tar Then lo = i%
  If i% * i% > tar Then hi = i% : Exit For
Next i%
Print "Calculating the square root of "; tar; " somewhere between "; lo; " ("; lo * lo; ") and "; hi; " ("; hi * hi; ")"

In this example, we dimension a lot of variables as floating point double precision.  I normally work solely with integers when I can, since GB32 works faster with integers, but square roots are often fractional by nature.  In this case, we start at 0 and work our way up to just larger than the target value, checking each number to see if the square of the number is less than or larger than our target.  We set the lo and hi ends of our brackets appropriately, and exit out when we're done.  Now, if you know anything about square roots, you might think that going larger or even anywhere near as large as the target value is a bit extreme, since the square root is always much smaller than the square.  However, if you start looking for the square roots of numbers in the range 0>x>1, also know as fractions or decimal numbers, you'll find that the square roots are actually larger than the squares.  Sqrt(0.25) -> 0.5, for example.  Now, to do the actual bracketing.

For i% = 1 To 40
  med = (lo + hi) / 2
  tmp = med * med
  Print "Step "; i%, lo, med, hi, tmp
  If tmp > tar Then
    hi = med
    lo = med
  If lo NEAR hi Then Exit For
  If tmp NEAR tar Then Exit For
  Delay .25
Next i%

I've limited this to 40 iterations, but you'll probably find your answer quicker than that.  I picked 40 because it pretty much fills the window I was using for output, and if you haven't found an answer in the first 40 or 50 iterations, you're probably doing something wrong.

Knowing that lo is too low and that hi is too high, we pick what we hope is a better guess by taking the average of the two values, which we'll call med.  Our first run through, we know that 5 is too low and 6 is too high, so we'll guess 5.5, which is actually pretty close at 30.25.  It's still a little high, so we set the high end of our bracket to this value and try again.  Looking between 5 and 5.5, we guess 5.25, which squares to 27.56, which is a little low, so we look between 5.25 and 5.5, and continue.

Notice the use of the NEAR operation.  In GB32, it's used to compare floating point variables.  Because GB32 uses a temporary 64 bit floating point number for internal calculations, which is many more bits than those used to store single or double precision variables, it's sometimes hard to get floating point math to generate exact matches.  NEAR checks floating point variables or calculations "only" out to 10 or 12 significant digits.  This is generally "close enough."

This method is fairly easy to understand, but is also fairly inefficient.  It takes about 35 steps to calculate a NEAR value, and each step is not always closer to the answer than the previous one, especially in the first few steps.  In math terms, the method is convergent, but not quickly convergent.

Successive approximation is a technique where you guess at an answer, check it, and use how close it is to fine tune your next guess.  How do we turn our knowledge of square roots into a successive approximation?  Well, we know that a root times a root equals the square.  If we divide both sides of that equation by the root, it turns out that a root equals the square divided by the root.  If we were again looking for the square root of 30 and guessed 5, then x=30/5, x=6.  5 and 6 aren't the same, but we know the root we're trying to calculate is somewhere in between, so we can guess (5+6)/2=5.5 for our new approximation.  x=30/5.5, x=5.45, which refines to 5.47 for our next guess.  Very quickly, we converge on an answer.

For i% = 0 To tar + 1
  If i% * i% >= tar Then old_app = i% : Exit For
Next i%
Print "Calculating the square root of "; tar; " somewhere near "; old_app; " ("; old_app * old_app; ")"
For i% = 1 To 40
  new_app = (tar / old_app)
  new_app = (new_app + old_app) / 2
  tmp = new_app * new_app
  Print "Step "; i%, new_app, tmp
  If old_app NEAR new_app Then Exit For
  old_app = new_app
  Delay .25
Next i%

This might be a little harder to understand than the bracketing method, but it's the same general principle.  We're using the results of our last guess to refine the next guess.  Or, as known in math circles, approximations.  In this case, we're dividing the square by an approximated root.  If our approximation was exact, the answer we'd get would be the exact root.  We assume it isn't, and that the real answer is somewhere between our approximation and our calculated root, so we average these two numbers and use the average as our next approximation.

This sequence converges rapidly and generally takes less than 10 steps.  Only 5 in this example.  Again, we use NEAR to compare our old approximation to our new one.  If they're NEAR enough, we know we've found our answer out to 10 or 12 significant digits, which is close enough for government work.

Useless Info:  This writer independently derived the bracket method shown above in a 7th grade science class when a problem called for a square root and I had left my calculator at home that day.  The science teacher, who was also my great aunt, was not impressed and demanded that I use the method taught in math class which vaguely resembles long division.  I protested that the long division method wasn't taught in math class and luckily the rest of the class backed me up.  The teacher left the room for a few minutes, probably checking this with the math teacher, and returned ranting about how they'd changed the math curriculum without telling her.  She spent the rest of the class ignoring the science assignment and teaching us how to "correctly" calculate square roots without a calculator.

This writer independently derived the division method shown above in a college level numerical methods class when the instructor decided that the bracket method didn't count as successive approximation.  The division method received a passing grade, but included a notation that I'd used Newton's method wrong.

Newton's method, which it turns out Newton neither invented nor probably ever used, is also known as the Babylonian method, and uses the same division and averaging of approximations, just rewritten slightly, but mathematically identical.  The difference being that I can quickly derive the division method from the definition of square roots, while I always have to look up the Newton method.

The reader is encouraged to download examples of all three methods and to try them out with other target values.  Suggested values include very large (12345678), very small (0.0005), your birth year (1967), the current year (2015), and important dates from history (1492, 1776, 2525).  The reader is also encouraged to un-comment the lines which set standard default starting values instead of the calculated ones to see how many extra iterations this requires.
Troy Cheek,
Jul 18, 2015, 3:15 PM