Thursday, April 15, 2010

Computing square / cube roots

We are very much used to using the square root functions provided in the library. But many of us don't have a clue how it can be done ourselves ( at least some of my friends don't have clue because we were discussing about it sometime back). So how can it be implemented if you need one? I am obviously at advantage because I had previlege to take Numerical Methods course during my undergrad. But even those of you who didn't take Numerical methods course know a searching technique which can be used to compute square roots. What is the technique? Of Course, binary search. How can you use binary search to compute square root of say 10. The idea is simple: consider the interval 0 - 10. Find the midpoint which is 5. Check if 5*5 is less than the given number. If it is less then the root greater than 5 else it is less than 5. Keep on proceeding like this until u find the actual number. This would translate to code as follows:



double find_root(double number) {
double lo = 0;
double hi =number;
double mid;
double epsilon = 0.000001; //precision
while ( lo <>
mid = (lo + hi)/2;
if (abs((mid * mid) - number) <= epsilon)
break;
else if (mid * mid < number)
lo = mid;
else
hi = mid;

}
return mid;
}


This works but takes time in converging. For those of you familiar with numerical algorithms you know that there is a better way of doing this. One most widely used method for root finding is Newton Raphson method which exploits concept of derivative. In this method a root is guessed and then the guess is converged to the actual root. For example your initial guess is X0 then the next value closer to the root would be:
X1 = X0 - f(X0) / f`(X0)
Similarly X2 = X1 - f(X1) / f`(X1) and so on until you get the precision that you want. Generally, it can be written as:

X = X - f(X) / f`(X)

The idea is that at each iteration it takes you closer to the root. We know that root is a poin in XY plane where the curve meets X- axis. The above formula is taking you close to that X-axis whenever you apply it. Now lets express finding square root as a root finding problem. Say you want to find square root of c. Let the square root be x then we can write the equation as:

x^2 = c or x^2 - c = 0

We require computing derivative for this at each new value of x that we obtain iteratively. This cane be done by referring to the post: http://codinggeeks.blogspot.com/2010/04/computing-derivatives.html


double newton_raphson(double number) {
double fx , fdashx;
double x = 3.00; //make an smart guess here
double epsilon = 0.00001; //choose your precision
do {
fx = computer fx at x;
fdashx = compute_derivative(x);
x -= fx / fdashx;

} while ( abs (x * x - c) > epsilon);

return x;
}


You can choose your own precision. Note that this method can be uesd to compute any root though I illustrated code for square root computation. For computing any root replace the while condition with your function.

No comments:

Post a Comment