# Square root estimation

I happened to read a something online about estimating a square root in a very simple way, known as Heron's method or the Babylonian method for square roots. The following (python?) code was included in the post:

def heron_sqrt(n, precision=0.000000000001):
if n == 0:
return 0
guess = n / 3.0
while True:
lastguess = guess
guess = (guess + (n / guess)) / 2.0
if abs(lastguess - guess) <= precision:
break
return guess


What is a good approach to implement this very simple algorithm in a Mathematica function?

It's fairly simple using FixedPoint.

--- edited for safety, per recommendation by @dionys ---

heronSqrt[x_/;Element[x,Reals]&&x>=0] := FixedPoint[(# + x/#)/2. &, x/3.]

In[268]:= heronSqrt[2]

(* Out[268]= 1.41421356237 *)


Alternatively, seed the initial point with a complex value.

heronSqrtForComplexes[x_] :=
FixedPoint[(# + x/#)/2. &, x/3. + x*I/100.]

heronSqrtForComplexes[-2]

(* Out[201]= 0. - 1.41421356237 I *)

• Nice. This still has difficulties with negative numbers though ... and just keeps running, since it doesn't converge. – dionys Oct 13 '14 at 22:34
• @dionys I'm not sure why you would expect it to work with negative numbers though. There's nothing in the algorithm that will give you a complex number. – wxffles Oct 13 '14 at 23:22
• @wxffles It's not a serious question, really. There are better ways to get square roots and no one needs this. It's more like code golf. Negative numbers are just an obvious failing. Maybe there's a better way to do it that can showcase a cool idea that might actually be useful in something else. Anyway, I thought it was fun to puzzle about, which is why I posted it (and my answer below) in the first place. – dionys Oct 13 '14 at 23:28
• Tougher, meaner version now also in place. Not intended for children or people with heart conditions. Ask your doctor if heronSqrtForComplexes is suitable for your healthy usage. – Daniel Lichtblau Oct 14 '14 at 0:55
• I just wanted it to be a smallish offset. Would probably be better to use a random small value, in case there remains a "bad ray" as there was with the negative reals in the earlier version. – Daniel Lichtblau Oct 14 '14 at 14:51

Using NestList seems like a decent approach:

HeronSqrt[x_, n_: 5] := Module[{f},
f[num_, est_] := N@Mean[{est, num/est}];
NestList[f[x, #] &, x/3., n]]

HeronSqrt[36]

{12., 7.5,6.15, 6.00183, 6., 6.}


But I'm sure there are nicer ways to implement this. For instance, the above doesn't do so well with negative numbers:

HeronSqrt[-4, 100] // ListPlot