Square Roots: Babylonian Style!!

[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]

This is cool :-)

[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]

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]

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]

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]

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]

Oh, very good then.

Yours,
Anselme

[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]

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]

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]

Oops, I forgot about this. I've moved on to linear algebra now. =P

[rolliebollocks]

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]

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]

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