## Square Roots: Babylonian Style!!

KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:

### Square Roots: Babylonian Style!!

I did some searching on this forum and found that a faster method exists for calculating square roots without using SQR(), but I will post this anyway since this method still perplexes me. According to my calculus book, this method was used by the ancient Babylonians (except, without the computer :-P).

Tomorrow I will ask my teacher how this is supposedly derivable from Newton's method. ^_^

Code: Select all

`' Square Roots!' by Kristopher WindsorFunction squareroot (Byref number As Double) As Double  'BTW, just use SQR() instead! :-P  Static As Double r1, r2  r1 = 1  Do    r2 = r1    r1 = (r1 + number / r1) * .5  Loop Until Abs(r1 - r2) < 1E-14  Return r1End FunctionSub showroot (Byref number As Double)  Print "The square root of " & Right("0" & number, 2) & " is " & squareroot(number)End SubFor a As Double = 1 To 16  showroot aNext aSleep`
anonymous1337
Posts: 5494
Joined: Sep 12, 2005 20:06
Location: California
This is cool :-)
Mysoft
Posts: 779
Joined: Jul 28, 2005 13:56
Location: Brazil, Santa Catarina, Indaial (ouch!)
Contact:
after take a quick look this code. i guess that i looks like an "inverse" newtwon method... as far i remeber ;P (Bad memory haeueahuae)
integer
Posts: 391
Joined: Feb 01, 2007 16:54
Location: usa
That is Newton's Method.

if the SquareRoot of n = x then x^2 = n

x^2 - n = 0

x*x/x - n/x = 0

x - n/x =0

x = n/x

Since we do not know what x is,
Newton makes a guess and takes the
average of the guess and the calculation

first guess that the square root of n = 1

guess2 = (Guess1 + n/guess1)/2

thus,

x2 = ( x1 + n/x2)/2

or

x2 = ( x1^2 + n)/(2*x1)

Then the 2nd guess replaces the first
and the routine is iterated until an acceptable
error tolerance is reached.

-----
Incidentally, the method that used to be taught
in middle school, the paper & pencil routine
is based on:

(x+e)^2 = n

x^2 +2xe + e^2 = n

The intermediate values:
iterated with this
(2x*(base)+e)*e

if you are in base 10 then
(10*2*x + e )*e = (20x+e)*e

with x representing all of the previous digits (estimates)

The series of e's are the sucessive digits of the
square root of the number, until you reach the desired
tolerance wanted.
Bunuel66
Posts: 76
Joined: May 19, 2006 19:56

### Generalized Newton Method

In France there is a generic algorithm named Generalized Newton Method:

Let x_n an estimated of X such that f(X)=a. If x_n is close enough we can consider that a small h added to x_n will reach X, then, at first order:
f(x_n+h)=f(x_n)+g(x_n).h with g(x) the first derivative of f: d/dx[f(x)]

Then h=[a-f(x_n)]/g(x_n) and x_n+1=x_n+[a-f(x_n)]/g(x_n)

This sequence is converging (generally...) to X=f-1(a) f-1: inverse function of f.

Using f(x)=x^2 g(x)=2.x and applying the above algorithm provides the exact expression used for the so called babylonian algorithm.

This can be applied to a large number of functions for generating sequences allowing to compute by iterations the inverses.

Best regards.
parakeet
Posts: 48
Joined: Mar 30, 2006 15:46
Location: Lille, France

### by far this is the fastest

Hi,
If you aim at computing distances, then dividing the points' x and y coords by those distances (typically for 3D rendering), then the fastest algo (used by J. Carmack for Quake, but reportedly invented by... someone else) is :

float InvSqrt(float x){
float xhalf = 0.5f * x;
int i = *(int*)&x; // store floating-point bits in integer
i = 0x5f3759d5 - (i >> 1); // initial guess for Newton's method
x = *(float*)&i; // convert new bits into float
x = x*(1.5f - xhalf*x*x); // One round of Newton's method
return x;
}

and is explained everywhere.

Any volunteer translate it to FB ? :)

Yours,
Anselme
KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:
Parakeet, I think that is the same method I mentioned before.
I "found that a faster method exists for calculating square roots without using SQR()" (said in original post).

http://www.freebasic.net/forum/viewtopic.php?t=6597

;-)
parakeet
Posts: 48
Joined: Mar 30, 2006 15:46
Location: Lille, France

### fast inv square

Oh, very good then.

Yours,
Anselme
h4tt3n
Posts: 694
Joined: Oct 22, 2005 21:12
Location: Denmark
It's a nice function, but the accuracy is too low for it to be useful for anything else than making a fair estimate. In any form of math / physics implementation demanding accurate calulations the function seems pretty much worthless.

But isn't there a way to increase accuracy? As far as I understand it takes a good initial guess and then performs one iteration with the Newton-Rapson approximation method.

Having not yet checked up on it myself, does anyone here know how to make it perform more iterations, thus increasing the functions accuracy?

EDIT:

Ah, it's really as simple as just repeating the " x *= (1.5 - xhalf * x * x) " step any desired nuber of times. Anyway, in order to get a decent accuracy for single values you need to do three iterations, making the algo more than ten times slower than the native sqr() fuction. Too bad.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york
You can increase the accuracy of the result by repeating the process over and over and over and over and over again...

If you were to do this for an infinite amount of time, the number derived would approach perfect accuracy.

How this is distinct from any other methodology is beyond my comprehension, since to some extent, they all involve some degree of probability and approximation.

Math doesn't ever get to be reality itself, merely a model for it, and always a reductive model with certain variables considered negligible.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york
Kristopher,

Does it help that the babylonians were using a base of 60 instead of 10?

rb
KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

I would assume it's similar to Newton's implementation of the limit, since it approaches a discrepancy of 0 through each (recursive) manipulation.

I would assume it's related to the calculus of infinitesimals.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york
Hey Kris,

Did you realize this is faster than the native SQR function, and produces identical output?

<edit>

Optimized:

Code: Select all

`' Square Roots!' by Kristopher WindsorFunction squareroot (Byval number As Double ) As Double  dim As Double r1=1, r2=any  Do    r2 = r1    r1 = (r1 + number / r1) * .5  Loop Until Abs(r1 - r2) <  .001  Return r1End Functiondim as double tt=timer? sqr(666666669911)? timer-tt=timer? squareroot (666666669911)? timer-tsleep`
KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:
I'm quite sure it's not that fast; I think your test is flawed. :)