This is the results of yet another adventure that started with me reading something and wanting to know more about it. July 30, 2003 John N. Power (jpower@bcpl.net) posted a message to the PICList which mentioned the ENIAC (Electronic Numerical Integrator and Computer) and it's square root algorithm.
John published a paper on the algorithm. The original website hosting the paper appears to be gone, but wayback machine archive has a copy of the article.
Additional searches revealed Brian Shelburne's well written How the ENIAC took a Square Root.
Both papers make it clear how the algorithm works, the latter even provides some improvements. Something just made me want to code it up and play with it.
I knew the ENIAC algorithm was not the most efficient algorithm, and wanted to compare it to two other methods that I knew about: bisection and successive approximation. So I ended up writing a library providing simple code for comparison testing of the algorithms.
Click here for a link to my square root library source. The rest of this page discusses the techniques for finding a square root and my implementation.
My general laziness and the fact that Brian Shelburne's paper did such a good job of explaining the ENIAC algorithm are preventing me from describing the algorithm. My source implements his algorithm with only one deviation. Brian Shelburne noted that the sum of odd multiples of 100 is also a square. This generalizes to the sum of odd multiples of a square being a square.
Most modern computers perform binary math, which is not well suited to
multiplying and dividing by 100. In an attempt to ease computation, I
experimented by trying odd multiples of 64 and 256, allowing multiplications
and divisions to be converted to shifts. As it turns out, there are only
three multiplications and divisions. The time saved by changing the square
used is less significant than the time saved by allowing quicker convergence
to an answer. The smaller the square multiple used, the faster the
algorithm converges (64 converges faster than 100). However, the smaller
the multiple, the less accurate the result. If you play with the source,
you might want to tweak the multiple (radix) defined in
sqrt.h.
Successive approximation is an iterative algorithm for approximating the square root of a number. The algorithm is based on the following observations:
If x is the square root of N then,
(1) (N/x) = x
so
(2) (x + (N/x)) / 2
= x.
It can be shown that if y is an approximation of the square
root of N then,
(3) (y + (N/y)) / 2
is a better approximation.
Using equation (3) repeatedly allows an approximation
for the square root of N to be continuously refined. The
refinement of the approximation may be stopped once the difference between
successive approximations is less than some threshold.
Bisection is another iterative process used to approximate the square
root of number. The algorithm works by repeatedly reducing the range of
values which the square root of a value is known to exist in. Given an
initial value N and an approximate square root guess of
x, the bisection algorithm is implemented using the following
steps:
Step 1. Set the square root
upper bound to N.
Step 2. Set the square root lower bound to 1.
*
Step 3.If (x)2 is
close enough to N, x is the square root.
Exit.
Step 4.
x)2 < N, set the
lower bound to x, and let x be half way
between the upper bound and the old x.
x)2 > N, set the
upper bound to x, and let x be half way
between the lower bound and the old x.
Step 5.Go to Step 3.
*NOTE: For N between 0 and 1, set the
lower bound to 0 and the upper bound to 1.
I am releasing my implementation of the ENIAC square root algorithm and all of the other approximation methods under the under the GNU LGPL. The source code repository is available on GitHub. I recommend that you checkout the latest revision of the master branch, unless you're looking for something specific.
| Repository Location | https://github.com/michaeldipperstein/sqrt |
|---|
All the source code that I have provided is written in strict ANSI C. I would expect it to build correctly on any machine with an ANSI C compiler. I have tested the code compiled with gcc on Linux and mingw on Windows 98.
If you have any further questions or comments, you may contact me by e-mail. My e-mail address is: mdipperstein@gmail.com
HomeLast updated on January 6, 2019