Square Roots: Babylonian Style!!

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

Square Roots: Babylonian Style!!

Postby KristopherWindsor » Dec 02, 2007 9:26

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

Postby anonymous1337 » Dec 02, 2007 9:52

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

Postby Mysoft » Dec 02, 2007 11:39

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

Postby integer » Dec 02, 2007 12:42

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

Postby Bunuel66 » Dec 02, 2007 19:22

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

Postby parakeet » Dec 06, 2007 17:02

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:

Postby KristopherWindsor » Dec 06, 2007 18:56

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

Postby parakeet » Dec 06, 2007 20:24

Oh, very good then.

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

Postby h4tt3n » Dec 08, 2007 15:07

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

Postby rolliebollocks » Oct 12, 2008 13:27

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

Postby rolliebollocks » Oct 12, 2008 13:28

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:

Postby KristopherWindsor » Oct 12, 2008 20:51

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

Postby rolliebollocks » Oct 13, 2008 17:13

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

Postby rolliebollocks » Aug 29, 2010 16:08

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:

Postby KristopherWindsor » Aug 29, 2010 20:16

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

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 3 guests