sqrt()

ID:22445
 
Nov 10 2006, 9:32 pm
Keywords: computers, math
Daxsa asks how to compute a square root without using the C library sqrt() function. I started writing my comment in response, when soon it grew to a full blog post. Here is my solution:

It's pretty common to use the Newton-Raphson method to compute the square root of a number. Recently I had to use this method to solve a different equation at work, but it's also good for solving the square root problem.

Basically you frame the question in terms of a function that equals zero. If you want to find the square root X of number N, thenX2 = N. Expressed as a zero-value equation, that is

X2 - N = 0


The Newton-Raphson method starts with an initial guess, and then iterates through a loop, altering the guess until the current guess is really close to the previous guess. Each time through the loop, the new guess G is computed by subtracting F/F' from the old guess X. F and F' are the function above and its derivative. The derivative of this particular function is:

2*X


So each time through the loop, we compute the new guess as follows:

G = X - (X2 - N)/(2*X)
or
G = X - 0.5*X + 0.5*N/X
or
G = 0.5*(X + N/X)


So the only thing left is what to use as your initial guess. Ideally you'd want to pick something that's pretty close to the final answer because otherwise for some functions, you might get the wrong result. However for this particular application, the function has only one positive solution, so there's no real danger if you simply start with the input number N as your initial guess.

mob/verb/square_root(n as num)
    var/nval = abs(n)
    var/x = 0
    var/g = nval
    var/tolerance = 0.00001
    var/nloop = 0

    // The loop should converge quickly.  If it doesn't, maybe our
    // tolerance is too small.
    while (abs(g - x) > tolerance && nloop < 1000)
        x = g
        g = 0.5 * (x + nval/x)
        nloop++

    src << "The square root of [n] is roughly [g][n < 0 ? "i" : ""]"

#ifdef DEBUG
    src << "We were off by [sqrt(nval) - g] after [nloop] iterations"
#endif
Nov 10 2006, 9:35 pm
Uh, there is a power of function in BYOND. It's **.

If you already new this, then I am an idiot and should eat a pastrie of worms.
Nov 10 2006, 9:41 pm
Just out of curiosity, why are you defining this as a mob procedure?
Nov 10 2006, 10:03 pm
Eternal_Dread wrote:
Uh, there is a power of function in BYOND. It's **.

If you already new this, then I am an idiot and should eat a pastrie of worms.


There's also a squareroot function, he was just showing an example that could easily be adapted to other languages.

Nov 11 2006, 1:55 am
Popisfizzy wrote:
Just out of curiosity, why are you defining this as a mob procedure?

It's a verb, so presumably it's so that you can just chuck it in a DM file and test it by using the verb. =)
Nov 11 2006, 4:13 am
Out of interest, do you know how to compute log()?
Nov 11 2006, 7:22 am
Crispy wrote:
It's a verb

Ah, damn I thought it was a procedure. My bad.
Nov 11 2006, 7:37 am
Nadrew and Crispy are right. The code is just a demonstration verb, as DM already has a built in sqrt() function (which I even use in debug mode to show that my result is correct).

DeathAwaitsU wrote:
Out of interest, do you know how to compute log()?

Yep, just use the same method. The natural log, ln, is the logarithm of a number in base e. e has value 2.718281828459. So we want to solve for the natural log X of N, or ln(N) = X. To express this as a zero-value function that will work with the Newton-Raphson method, we take e to the power of both sides of that equation:

e^(ln(N)) = e^X

By the definition of natural log, e^(ln(N)) equals N, so:

N = e^X or
e^X - N = 0

There's your function of X. The derivative of that function is e^X. So now you can plug those two functions into the line of code that computes the guess G:

G = X - F/F' = X - ((e^X - N) / e^X)

or in the code, add this at the top:

var/e = 2.718281828459

and replace the g = ... line with:

g = x - ((e**x - nval)/e**x)

Easy peasy. That gets you the natural log. If you want the log in another base (such as base 10 or base 2), you simply divide by the natural log of the new base's number. So:

log10(N) = ln(N) / ln(10)
Nov 21 2006, 2:20 pm
thanks guys....

i was also thinking of having something like

Sqrt(x)= EXP^(1/2(x))

or some sort of infinite seris....
digit by digit calculation