1 Figure 3. First and second bits of
1xxx)01111001 square root are
-1000 being calculated.
The first half
of Figure 3 shows
an example with
the first bit calculated. Notice in
the figure how
we don’t know
the divisor yet, so I’ve substituted “xxx.”
The binary division procedure requires
multiplying the guess by the divisor, so
“0” is used for each “x” for now.
The second half of Figure 3 shows
the partially-completed calculation of
the second bit. We make a guess of 1,
multiply the trial bit by the divisor
(which is 11xx now to stay equal to our
quotient), and subtract from the
remainder. We must also account for
what we should have subtracted when
we computed the previous bit. Back
then, we only subtracted 1000. Now
we have a divisor of 11xx, so we should
have subtracted 1100 before. I’ve put
this make-up value to be subtracted on
an additional line. Because we’re shifted right one bit, it’s written as 1000,
but you can follow the “1” bit up the
column and see it’s in the right place.
So at this point, we have two values
to subtract — 1100 and 1000 — whose
binary total is 10100. This total amount
to subtract is larger than our current
remainder of 1110, so we must conclude the trial was too big — instead of
1, it must be 0.
The first half
of Figure 4 shows
the second bit set
back to 0, and
the calculation of
the third bit.
done the same:
We subtract the
product of the
trial bit and the
divisor, and also
include a line
to subtract the
value of the trial bit had it been set
back when we computed the first bit.
This time, the total of the two values to
subtract (1010 and 1000) is 10010,
while the current remainder is 11100.
The remainder is big enough, so the
trial bit is correct.
The second half of Figure 4 shows
the calculation of the final bit. Note
how that extra value we subtract (to
correct for previously subtracting the
wrong value) has more “1” bits in it.
Each of the 1 bits corresponds to a
subtraction for a previous quotient bit
where we didn’t subtract the bit we’re
working on now. So, for each 1 bit in
the quotient (not counting the trial bit),
there will be a corresponding 1 bit in
the extra value to subtract to account
for not subtracting it originally.
Once more: For each 1 bit in the
quotient, there will be a 1 bit in the
extra value to subtract. This is the same
as saying the extra value to subtract
equals the quotient (and divisor) before
setting our trial bit. So, we subtract the
divisor with the trial bit, and we subtract
the divisor without the trial bit. Adding
that up, we subtract 2 * divisor + trial.
That’s it. Since the divisor and
quotient are always equal, I’m going to
start calling them the root. To calculate
each bit of the root, try to subtract two
times the root plus the trial bit from the
current remainder. If it fits, the trial bit
is part of the root; if not, that bit is zero.
Listing 1 represents this algorithm
as a function in C. You don’t need to
know C to understand it if you look in
the sidebar for the explanation of the
operators. This algorithm also includes
the extra step of rounding the result
based on whether the next bit would
be 0 or 1. This function has been
compiled and tested for both Windows
and the AVR microcontroller.
LISTING 1. Fast square root algorithm in C.
typedef unsigned long UINT;
#define INT_BITS ((sizeof(UINT)) * 8)
UINT SquareRoot(UINT arg)
trial = 1 << (INT_BITS - 1); // set up 100000... binary
root = 0; // Inside loop, really root * 2
trial >>= 1;
root |= trial;
// move trial to next position in root
// combine trial bit into root
if (arg < root)
root ^= trial;
// does trial root bit fit?
// no, remove trial bit
arg -= root;
root += trial;
// root fits, remove it from arg
// double the trial bit
root >>= 1;
trial >>= 1;
// move both root & trial to next
// position within arg
} while (trial != 0);
// Compute rounding
// See if next bit computed would be 1; round up if so.
if (arg > root)
Figure 4. Third
and fourth bits of
square root being
SERVO 01.2008 39