Keith Vetter 2003-07-14 : How do you calculate the square root of a number? Well, now-a-days it's easy--just use sqrt. But suppose you don't have that function, how else can you do it?
Here are some various ways of doing it, from Newton's Method to the paper-and-pencil "long square root" method that I was taught in grade school which is akin to long division.
SQRT
Here's the best way but least fun. It's here for timing purposes:
proc SqrtA {num} {
# sqrt() -- for timing purposes
set ans [expr {sqrt($num)}]
return $ans
}Newton's Method
This is an iterative that takes an initial guess and successively makes it better using the recursion:
F(n+1) = (F(n) + num/F(n)) / 2
This converges quickly and is generally the method of choice. The tricky part is knowing when to stop (pick some epsilon) and determining the initial guess--any value will do but some converge faster. A lot has been written as to the best initial guess, the algorithm I use here is to use a number whose high-bit position is half that of the original number.
proc SqrtB {num} {
# Newton's method
set epsilon .0001
set ans [expr {1 << int((log($num)/log(2) / 2))}] ;# Initial guess
while {1} {
set last $ans
set ans [expr {(($num / $ans) + $ans) / 2.0}]
if {abs($last - $ans) < $epsilon} break
}
return $ans
}Successive Approximation
If your only dealing with integers then you can just do a brute force test of each possible bit.
proc SqrtC {num} {
# Successive integer approximation -- try adding each bit
set BitsPerInt [expr {4 * [string length [format %X -1]]}]
set bit [expr {1 << ($BitsPerInt/2)}]
set ans 0
for {set i 0} {$i < $BitsPerInt} {incr i 2} {
set bit [expr {$bit / 2}] ;# Next bit to possibly add
set tmp [expr {$ans + $bit}]
if {$tmp*$tmp <= $num} { ;# Is adding bit not too big???
set ans $tmp
}
}
return $ans
}
Long Square Root
For lack of a better name I'm calling this "Long Square Root" because it is a paper-and-pencil method very similar to long division. I learned this method in grade school but it seems to be a forgotten technique.
In assembly language, working in base 2, this is a blazingly fast and very compact algorithm involving just shifts, adds and comparisons.
As a refresher and to get a common terminology, in long division, you have the divisor, dividend and quotient. Each step involves bringing down the next digit from the dividend and adding it to the running dividend, then finding the maximum digit (D) which multiplied by the divisor is less than the running dividend. That digit (D) then becomes part of the quotient.
In long square root the method is almost identical except that you bring down pairs of digits and the divisor analog is a varying, more complex entity.
The square root analog to divisor, dividend and quotient terms I'll call cross term, square, root. Each step involves bringing down the next 2 digits from the square and adding it to the running square, then finding the maximum digit (D) which multiplied by the cross term (see below) is less than the running square. That digit (D) then becomes part of the root.
The cross term part changes at each step. It's value is: twice the current root with the maximum digit appended ==> 20*root + D
proc SqrtD {num} {
# "Long" square root - integer only
if {[string length $num] & 1} {set num " $num"} ;# Must be even length
set root 0
set v 0
foreach pair [regexp -inline -all {..} $num] {
scan $pair %d pair ;# Goddam octal nonsense
set v [expr {$v * 100 + $pair}] ;# Running square
# Find max D s.t. D * (2*(root*10) + D) < v
for {set D 1} {$D < 10} {incr D} {
if {(20 * $root + $D) * $D > $v} break
}
set v [expr {$v - (20 * $root + $D-1) * ($D-1)}]
set root [expr {10 * $root + $D-1}] ;# Add new digit to the root
if {$v < 0} {error "ERROR: overflow"}
}
return $root
}This is an algorithm that works much better in binary. "Bringing down" is just shifting bits in and finding the maximum digit is simply one comparison.
Long Square Root -- reals
The above method can easily handle real numbers--just like in long division, when you reach the decimal point just keep bringing down 0's (but for square root you bring down pair's of 0's). A tricky part is keeping track of where the decimal point should go.
This function will return as many decimal places as its input, to wit supplying 2.0000 will return 1.4142.
proc SqrtE {num} {
# "Long" square root - real numbers
foreach {int frac} [split $num "."] break ;# Get integer and fraction
if {[string length $int] & 1} {set int " $int"}
set flen [string length $frac]
if {$flen > 0} { append frac [string repeat "0" $flen]}
set root 0
set v 0
foreach pair [regexp -inline -all {..} "$int$frac"] {
scan $pair %d pair ;# Goddam octal nonsense
set v [expr {$v * 100 + $pair}] ;# Running square
# Find max D s.t. D * (2*(root*10) + D) < v
for {set D 1} {$D < 10} {incr D} {
if {(20 * $root + $D) * $D > $v} break
}
set v [expr {$v - (20 * $root + $D-1) * ($D-1)}]
set root [expr {10 * $root + $D-1}]
if {$v < 0} {error "ERROR: overflow"}
}
if {$flen > 0} { ;# Scale to proper decimal
set root [expr {$root / pow(10,$flen)}]
}
return $root
}Timing values to come