Finding square root without using sqrt function?

Question

I was finding out the algorithm for finding out the square root without using sqrt function and then tried to put into programming. I end up with this working code in C++

    #include <iostream>
    using namespace std;

    double SqrtNumber(double num)
    {
             double lower_bound=0; 
             double upper_bound=num;
             double temp=0;                    /* ek edited this line */

             int nCount = 50;

        while(nCount != 0)
        {
               temp=(lower_bound+upper_bound)/2;
               if(temp*temp==num) 
               {
                       return temp;
               }
               else if(temp*temp > num)

               {
                       upper_bound = temp;
               }
               else
               {
                       lower_bound = temp;
               }
        nCount--;
     }
        return temp;
     }

     int main()
     {
     double num;
     cout<<"Enter the number\n";
     cin>>num;

     if(num < 0)
     {
     cout<<"Error: Negative number!";
     return 0;
     }

     cout<<"Square roots are: +"<<sqrtnum(num) and <<" and -"<<sqrtnum(num);
     return 0;
     } 

Now the problem is initializing the number of iterations nCount in the declaratione ( here it is 50). For example to find out square root of 36 it takes 22 iterations, so no problem whereas finding the square root of 15625 takes more than 50 iterations, So it would return the value of temp after 50 iterations. Please give a solution for this.

Selfstarter

:) root = pow(x,0.5); – Gunfight

Usually you check the difference each iteration is making, and when it drops below a certain level, you treat the result as accurate enough. Also note that Newton's method is usually quite a bit faster than bisection like you're using. – Urina

what about checking the difference between temp*temp and the number and stopping the loop if is close enough ? – Fizzy

@Gunfight love it! On a serious note, this is the professional way: en.wikipedia.org/wiki/Newton's_method – Vestry

@MartinBeckett yeah it worked if(temp*temp-num > 0.001) {upper_bound = temp;}else { lower_bound = temp; } Is 0.001 good to check close enough? – Selfstarter

@Vestry Yes, it was a tongue-in-cheek comment. – Gunfight

@ArunPandey: How accurate do you want your result? You typically want to get N significant digits, in which case you're looking for something like input / 10eN. – Urina

or root = exp(0.5 * log(x)) – Metronymic

@Vestry moreover, n-th root of x === x^(1/n) – Milomilon

Answer 1

There is a better algorithm, which needs at most 6 iterations to converge to maximum precision for double numbers:

#include <math.h>

double sqrt(double x) {
    if (x <= 0)
        return 0;       // if negative number throw an exception?
    int exp = 0;
    x = frexp(x, &exp); // extract binary exponent from x
    if (exp & 1) {      // we want exponent to be even
        exp--;
        x *= 2;
    }
    double y = (1+x)/2; // first approximation
    double z = 0;
    while (y != z) {    // yes, we CAN compare doubles here!
        z = y;
        y = (y + x/y) / 2;
    }
    return ldexp(y, exp/2); // multiply answer by 2^(exp/2)
}

Algorithm starts with 1 as first approximation for square root value. Then, on each step, it improves next approximation by taking average between current value y and x/y. If y = sqrt(x), it will be the same. If y > sqrt(x), then x/y < sqrt(x) by about the same amount. In other words, it will converge very fast.

UPDATE: To speed up convergence on very large or very small numbers, changed sqrt() function to extract binary exponent and compute square root from number in [1, 4) range. It now needs frexp() from <math.h> to get binary exponent, but it is possible to get this exponent by extracting bits from IEEE-754 number format without using frexp().

Emmer

I guess the condition for the stopping rule should use less than (or equal): while (fabs(y-z) <= 0.00000000001 – Endosteum

@EmmadKareem: nope, rule is fine. I would change it to use relative precision though - it may have trouble for very large numbers. Another nice improvement is to limit number of rounds to say 20 – Emmer

I did some tests, and for numbers 1 to 4 it only takes 6 rounds to converge to maximum double precision. However, for large numbers like 1000000 it will take longer. Solution is to take exponent separately, and compute sqrt only from number that neither large nor small - somewhere around 1. This way number of rounds can be even fixed to be say 6. – Emmer

Also, I think this is one of the very rare cases where instead of checking fabs(y-z) < delta you actually can use exact equality check for doubles: y == z. This is because algorithm converges very fast, and once converged, it will be exact same number from this point on – Emmer

@Emmer I know y = (y + x/y) / 2; comes from y^2 = x, but what clicked in your head to write this, in Mathematics what is this trick called as? As if we need to study this method in general what should we search upon? – Thulium

This is called Babylonian method. This formula can be also derived from Newton's method. – Emmer

Answer 2

Why not try to use the Babylonian method for finding a square root.

Here is my code for it:

double sqrt(double number)
{
    double error = 0.00001; //define the precision of your result
    double s = number;

    while ((s - number / s) > error) //loop until precision satisfied 
    {
        s = (s + number / s) / 2;
    }
    return s;
}

Good luck!

Avitzur

use ( & ) properly, makes reading code hard.. I was confused for a second why are you writing s-number when s = number – Faceless

I would not use fixed decimal points. I would replace the line double error = 0.00001; by double error = number * 10e-8; - that would make the error relative to size of input value. But +1 for nice and simple example. – Peterpeterborough

Great answer! Simple and efficient enough. Interestingly, if one needs to find the nearest integer solution, one can turn the loop into an infinite one, and exit when the last two values of s are equal. – Landaulet

Answer 3

Remove your nCount altogether (as there are some roots that this algorithm will take many iterations for).

double SqrtNumber(double num)
{
    double lower_bound=0; 
    double upper_bound=num;
    double temp=0;

    while(fabs(num - (temp * temp)) > SOME_SMALL_VALUE)
    {
           temp = (lower_bound+upper_bound)/2;
           if (temp*temp >= num)
           {
                   upper_bound = temp;
           }
           else
           {
                   lower_bound = temp;
           }
    }
    return temp;
 }

Pinna

not working it just shows nearest values even for the perfect square – Selfstarter

Oops, yes, you are correct (sorry, typing too fast). As @Emmer pointed out, this is a very slow algorithm for finding square root approximations. There are others that are much faster. – Pinna

Answer 4

As I found this question is old and have many answers but I have an answer which is simple and working great..

#define EPSILON 0.0000001 // least minimum value for comparison
double SquareRoot(double _val) {
    double low = 0; 
    double high = _val;
    double mid = 0; 

    while (high - low > EPSILON) {
            mid = low + (high - low) / 2; // finding mid value
            if (mid*mid > _val) {
                high = mid;
            } else {
                low = mid;
            }    
    }   
    return mid;
}

I hope it will be helpful for future users.

Hilbert

Try to call this function with value of say 10,000,000,000, and see how long does it take to get an answer. Or, call with 0.000000001, and see if answer is correct. – Emmer

Answer 5

if you need to find square root without using sqrt(),use root=pow(x,0.5).

Where x is value whose square root you need to find.

Mudd

Answer 6

Here is a very awesome code to find sqrt and even faster than original sqrt function.

float InvSqrt (float x) 
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f375a86 - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    x = x*(1.5f - xhalf*x*x);
    x = x*(1.5f - xhalf*x*x);
    x=1/x;
    return x;
}

Tacheometer

1) without explanation, it is worthless. 2) it ends in a division, showing that the idea is not understood 3) Without attribution, it is plagiarism. – Lombardy

@PascalCuoq: Author is unknown. hackersdelight.org/hdcodetxt/rsqrt.c.txt – Discerning

Actually it is attributed to Greg Walsh: en.wikipedia.org/wiki/… . If you read the many essays that have been written on the subject of this function, you could use the same tricks to computed sqrt directly, instead of computing the inverse square root and taking the inverse, which completely undoes any elegance there was in the original algorithm: en.wikipedia.org/wiki/… – Lombardy

Answer 7

//long division method.
#include<iostream>
using namespace std;
int main() {
int n, i = 1, divisor, dividend, j = 1, digit;
cin >> n;
while (i * i < n) {
    i = i + 1;
}
i = i - 1;
cout << i << '.';

divisor  = 2 * i;
dividend = n - (i * i );
while( j <= 5) {
    dividend = dividend * 100;
    digit = 0;
    while ((divisor * 10 + digit) * digit < dividend) {
        digit = digit + 1;
    }
    digit = digit  - 1;
    cout << digit;
    dividend = dividend - ((divisor * 10 + digit) * digit);
    divisor = divisor * 10 + 2*digit;
    j = j + 1;
}
cout << endl;
return 0;
}

Arcadia

Answer 8

This a very simple recursive approach.

double mySqrt(double v, double test) {
    if (abs(test * test - v) < 0.0001) {
        return test;
    }

    double highOrLow = v / test;
    return mySqrt(v, (test + highOrLow) / 2.0);
}
double mySqrt(double v) {
    return mySqrt(v, v/2.0);
}

Shena

Answer 9

After looking at the previous responses, I hope this will help resolve any ambiguities. In case the similarities in the previous solutions and my solution are illusive, or this method of solving for roots is unclear, I've also made a graph which can be found here.

This is a working root function capable of solving for any nth-root

(default is square root for the sake of this question)

#include <cmath> 
// for "pow" function

double sqrt(double A, double root = 2) {
    const double e = 2.71828182846;
    return pow(e,(pow(10.0,9.0)/root)*(1.0-(pow(A,-pow(10.0,-9.0)))));
}

Explanation:

click here for graph

This works via Taylor series, logarithmic properties, and a bit of algebra.

Take, for example:

log A = N
   x

*Note: for square-root, N = 2; for any other root you only need to change the one variable, N.

1) Change the base, convert the base 'x' log function to natural log,

log A   =>   ln(A)/ln(x) = N
   x

2) Rearrange to isolate ln(x), and eventually just 'x',

ln(A)/N = ln(x)

3) Set both sides as exponents of 'e',

e^(ln(A)/N) = e^(ln(x))  >~{ e^ln(x) == x }~>  e^(ln(A)/N) = x

4) Taylor series represents "ln" as an infinite series,

ln(x) = (k=1)Sigma: (1/k)(-1^(k+1))(k-1)^n

           <~~~ expanded ~~~>

[(x-1)] - [(1/2)(x-1)^2] + [(1/3)(x-1)^3] - [(1/4)(x-1)^4] + . . .

*Note: Continue the series for increased accuracy. For brevity, 10^9 is used in my function which expresses the series convergence for the natural log with about 7 digits, or the 10-millionths place, for precision,

ln(x) = 10^9(1-x^(-10^(-9)))

5) Now, just plug in this equation for natural log into the simplified equation obtained in step 3.

e^[((10^9)/N)(1-A^(-10^-9)] = nth-root of (A)

6) This implementation might seem like overkill; however, its purpose is to demonstrate how you can solve for roots without having to guess and check. Also, it would enable you to replace the pow function from the cmath library with your own pow function:

double power(double base, double exponent) {
    if (exponent == 0) return 1;
    int wholeInt = (int)exponent;
    double decimal = exponent - (double)wholeInt;
    if (decimal) {
        int powerInv = 1/decimal;
        if (!wholeInt) return root(base,powerInv);
        else return power(root(base,powerInv),wholeInt,true);
    }
    return power(base, exponent, true);
}

double power(double base, int exponent, bool flag) {
    if (exponent < 0) return 1/power(base,-exponent,true);
    if (exponent > 0) return base * power(base,exponent-1,true);
    else return 1;
}

int root(int A, int root) {
    return power(E,(1000000000000/root)*(1-(power(A,-0.000000000001))));
}

Multiphase