Square Roots: Babylonian Style!!

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:

Square Roots: Babylonian Style!!

Post by KristopherWindsor »

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 Windsor

Function 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 r1
End Function

Sub showroot (Byref number As Double)
  Print "The square root of " & Right("0" & number, 2) & " is " & squareroot(number)
End Sub

For a As Double = 1 To 16
  showroot a
Next a
Sleep
anonymous1337
Posts: 5494
Joined: Sep 12, 2005 20:06
Location: California

Post by anonymous1337 »

This is cool :-)
Mysoft
Posts: 836
Joined: Jul 28, 2005 13:56
Location: Brazil, Santa Catarina, Indaial (ouch!)
Contact:

Post by Mysoft »

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: 408
Joined: Feb 01, 2007 16:54
Location: usa

Post by integer »

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

Post by Bunuel66 »

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

Post by parakeet »

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:

Post by KristopherWindsor »

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

Post by parakeet »

Oh, very good then.

Yours,
Anselme
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

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

Post by rolliebollocks »

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

Post by rolliebollocks »

Kristopher,


Did you ever get an answer from your prof?

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:

Post by KristopherWindsor »

Oops, I forgot about this. I've moved on to linear algebra now. =P
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

KristopherWindsor wrote:Oops, I forgot about this. I've moved on to linear algebra now. =P
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

Post by rolliebollocks »

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 Windsor

Function 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 r1
End Function

dim as double t

t=timer
? sqr(666666669911)
? timer-t
t=timer
? squareroot (666666669911)
? timer-t
sleep
KristopherWindsor
Posts: 2428
Joined: Jul 19, 2006 19:17
Location: Sunnyvale, CA
Contact:

Post by KristopherWindsor »

I'm quite sure it's not that fast; I think your test is flawed. :)
Post Reply