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
}AMG: Or just use [::tcl::mathfunc::sqrt].
Newton's Method
This is an iterative process 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 you are 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 value.
The square root analog to the divisor, dividend and quotient terms I'll call cross term, square, and 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. Its 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.
AM It is a pity that this method can not easily be extended to arbitrary roots, such as the cubic root. Note that the method boils down to:
x^2 = x1^2 + 2 * x1 * d + d^2 d -> x2 x^2 = (x1+x2)^2 + 2 * (x1+x2) * d + d^2 d -> x3 etc.
with each d between 0 and 9*10^-n.
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 pairs 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.
Tcl 8.5 will have a proper "integer square root" function isqrt [L1 ] (given a number, find the integer part of its square root). Until 8.5 is out, here's a stopgap:
proc ::tcl::mathfunc::isqrt {x} {
set l [string length $x]
# if floating point will return an exact result, use it
if { $l <= 15 } {
return [expr {int(sqrt($x))}]
}
# get a first approximation in floating point, rescaling to avoid
# overflows. (Should use powers of 2 for this!)
set p10 [expr {($l - 34) / 2}]
if { $p10 > 0 } {
set x1 [expr {$x / 10**(2*$p10)}]
set s0 [expr {(10 ** $p10) * entier(sqrt($x1))}]
} else {
set s0 [expr {entier(sqrt($x))}]
}
# Improve the result with Newton's method
set s1 $s0
while {1} {
set s2 [expr {($s1 + $x/$s1) >> 1}]
if {$s2 == $s0 || $s2 == $s1} break
set s0 $s1
set s1 $s2
}
if {$s2 * $s2 > $x} {
return [expr {$s2 - 1}]
} else {
return $s2
}
}Sarnold -- Well, you could use the math::bigfloat package from tcllib, though I admit it makes use of some exp() and log()...
package require math::bigfloat namespace import ::math::bigfloat::* puts [tostr [sqrt [fromdouble $mydouble 15]]]
DKF: I was reading [L2 ] (Wayback Machine Web Archive) about some obscure code that calculates an approximation to the inverse square root of a number (i.e., it raises the number to the power -0.5) a value useful in much 3D code apparently. What particularly interested me was the method of calculation, which I reproduce here in Tcl code with some simplifications to make the code a bit more readable:
proc f2i x {binary scan [binary format f $x] i f; set f}
proc i2f x {binary scan [binary format i $x] f f; set f}
proc InvSqrt x {
set xhalf [expr {0.5 * $x}]
set i [f2i $x]
set i [expr {0x5f3759df - ($i>>1)}]
set y [i2f $i]
set y [expr {$y * (1.5 - $xhalf * $y * $y)}]
return $y
}Testing it out:
% InvSqrt 0.25 1.99661429917 % InvSqrt 0.01 9.98252138783
OK, that's an approximation, but it's deeply strange that it works at all...
Lars H: Not so strange, when one considers that the last line looks rather much like (and according to the text you link to indeed is) a Newton--Raphson iteration formula; for f(y) = 1/y**2 - x one gets f'(y) = -1/2y**3 and thus y - f(y)/f'(y) = 1.5*y - 0.5*x*y**3. Repeating the last calculation step will produce better approximations:
proc InvSqrt2 x {
set xhalf [expr {0.5 * $x}]
set i [f2i $x]
set i [expr {0x5f3759df - ($i>>1)}]
set y [i2f $i]
set y [expr {$y * (1.5 - $xhalf * $y * $y)}]
set y [expr {$y * (1.5 - $xhalf * $y * $y)}]
return $y
}
InvSqrt2 0.01 ; # Returns 9.99995420142It is impressive that it gets it that good in just one step, but N-R generally is pretty magical. There are for example tons of results that the asymptotic complexity of computing some function (e.g. the Square Root) is the same (up to a constant factor) as that of multiplying two numbers with the given number of digits, which are all based on the observation that the function can be computed by N-R iteration and since this converges so very fast it's only the very last iteration step that needs to be carried out at full precision...
As for the magical starting approximation (which is analysed in a paper [L3 ] linked to by the text linked to above), the idea is (if one simplifies heavily) to multiply the exponent of the floating point number by -0.5, by shifting the bits of its representation one step to the right and then adjusting for the way it is encoded. The trick relies on the fact that the float (i) is in IEEE-754 format, (ii) is 32 bits long (but changing the magic constant 0x5f3759df suffices for adapting the trick to other sizes), and (iii) has the same endianness as the integers. When $tcl_platform(byteOrder) is "bigEndian", this means one rather has to define
proc i2f x {binary scan [binary format I $x] f f; set f}
proc f2i x {binary scan [binary format f $x] I f; set f}FPX: Supposing that "sqrt" does not exist? I can't resist offering the following trivial alternative:
proc sqrt {n} {
return [expr {exp(log($n)/2)}]
}TP also likes using the B and C scales on A little slide-rule, which is FPX's sqrt proc in graphical form :-)
AM Here is a related page: What if you do not have exp()?