不使用sqrt函数如何求平方根？

Question

我正在研究一种不使用sqrt函数的算法来计算平方根，并尝试将其转化为编程代码。最终，我在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;
     } 

现在的问题是在声明中初始化迭代次数nCount（这里为50）。例如，要找出36的平方根需要22次迭代，因此没有问题，而找到15625的平方根需要超过50次迭代，因此它将在50次迭代后返回temp的值。请提供一个解决方案。

- Arun Pandey

---

:) root = pow(x,0.5); - jrok

通常情况下，您会检查每次迭代所产生的差异，当差异降至一定水平时，即可将结果视为足够准确。同时请注意，通常牛顿法比您正在使用的二分法快得多。 - Jerry Coffin

那么检查temp*temp与该数字之间的差异，如果足够接近，就停止循环怎么样？ - Martin Beckett

@jrok 真的很喜欢！说正经的，这是专业的方式：http://en.wikipedia.org/wiki/Newton's_method - Bathsheba

@MartinBeckett 是的，如果(temp*temp-num > 0.001)，它可以工作。那么0.001是检查足够接近的好数字吗？ - Arun Pandey

Answer 1

有一种更优秀的算法，最多只需要6次迭代即可收敛到双精度数的最高精度：

#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)
}

算法以1作为平方根的第一个近似值。然后，在每个步骤中，它通过取当前值y和x/y的平均值来改善下一个近似值。如果y = sqrt(x)，那么它将是相同的。如果y > sqrt(x)，那么x/y大约比sqrt(x)小相同的数量。换句话说，它会非常快速地收敛。

更新：为了加快对非常大或非常小的数字的收敛速度，将sqrt()函数更改为提取二进制指数，并从[1, 4)范围内的数字计算平方根。现在需要<math.h>中的frexp()获取二进制指数，但可以通过从IEEE-754数字格式中提取位而不使用frexp()来获得此指数。

- mvp

---

我猜停止规则的条件应该使用小于（或等于）符号：while (fabs(y-z) <= 0.00000000001) - NoChance

@EmmadKareem：不，规则很好。不过我会更改它以使用相对精度——它可能会对非常大的数字产生麻烦。另一个不错的改进是将回合数限制在20次左右。 - mvp

我进行了一些测试，对于数字1到4，只需要6轮就可以收敛到最大的双精度。然而，对于像1000000这样的大数，需要更长的时间。解决方案是将指数单独提取出来，并且仅从既不大也不小的数字中计算平方根 - 大约在1附近。这样可以将轮数固定为6轮。 - mvp

此外，我认为这是非常罕见的情况之一，即您可以使用双精度浮点数的精确相等性检查y == z，而不是检查fabs(y-z) < delta。这是因为该算法收敛非常快，一旦收敛，从此时起它将是完全相同的数字。 - mvp

@mvp 我知道 y = (y + x/y) / 2; 是来自于 y^2 = x，但是你的脑海中是什么启发你写下这个公式的？在数学中，这个技巧被称为什么？如果我们需要学习这种方法，应该搜索什么？ - codey modey

这被称为巴比伦方法。这个公式也可以从牛顿法中推导出来。 - mvp

Answer 2

为什么不尝试使用巴比伦算法求平方根。

以下是我的代码：

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;
}

祝你好运！

- Newskooler

---

使用 ( 和 ) 正确地，可以使代码难以阅读。我刚开始也有点困惑，为什么你要写 s-number，当 s = number。 - Abhinav Gauniyal

我不会使用固定的小数点。我会将行double error = 0.00001;替换为double error = number * 10e-8; - 这将使误差相对于输入值的大小。但是，这个简单明了的例子还是很好的，给你一个赞。 - Combinatix

很棒的答案！简单而高效。有趣的是，如果需要找到最近的整数解，可以将循环变成无限循环，并在s的最后两个值相等时退出。 - Marcus

Answer 3

我发现这个问题已经很老了，有很多答案，但是我有一个简单且非常有效的答案。

#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;
}

我希望这对未来的用户有所帮助。

- Himanshu

---

尝试使用值为10,000,000,000调用此函数，并查看获取答案需要多长时间。或者，使用0.000000001进行调用，并查看答案是否正确。 - mvp

Answer 4

完全删除您的nCount（因为有些根据此算法需要多次迭代）。

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;
 }

- Zac Howland

---

无法正常工作，即使是完美平方数也只显示最接近的值。 - Arun Pandey

哎呀，是的，你说得对（抱歉，打字太快了）。正如@mvp指出的那样，这是一种寻找平方根近似值非常缓慢的算法。还有其他的算法要快得多。 - Zac Howland

Answer 5

如果您需要在不使用sqrt()的情况下找到平方根，可以使用root=pow(x,0.5)。

其中x是您需要找到平方根的值。

- Mani Kant

Answer 6

//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;
}

- Kartik

Answer 7

这里有一种非常简单但不安全的方法来找到一个数字的平方根。 不安全是因为它只适用于自然数，其中您知道底数和指数都是自然数。我不得不在一个任务中使用它，那里既不允许我使用#include<cmath>库，也不允许我使用指针。

幂 = 底数 ^ 指数

// FUNCTION: square-root
int sqrt(int x)
{
    int quotient = 0;
    int i = 0;

    bool resultfound = false;
    while (resultfound == false) {
        if (i*i == x) {
          quotient = i;
          resultfound = true;
        }
        i++;
    }
    return quotient;
}

- PatrickSteiner

Answer 8

这是一种非常简单的递归方法。

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);
}

- Cameron Monks

Answer 9

这里有一段非常棒的代码，可以找到平方根，而且比原始的平方根函数还要快。

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;
}

- galaxyan

---

1）没有解释，就毫无价值。 2）以分歧告终，说明观点未被理解。 3）没有归属，就是剽窃。 - Pascal Cuoq

@PascalCuoq：作者不详。http://www.hackersdelight.org/hdcodetxt/rsqrt.c.txt - namar0x0309

实际上，这归功于Greg Walsh：http://en.wikipedia.org/wiki/Fast_inverse_square_root#History_and_investigation。如果你读了许多关于此函数的文章，你可以使用相同的技巧直接计算平方根，而不是计算倒数平方根并取倒数，这完全破坏了原始算法中的任何优雅之处：http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Approximations_that_depend_on_the_floating_point_representation - Pascal Cuoq

Answer 10

在查看之前的回复后，我希望这可以帮助解决任何不明确的问题。如果之前的解决方案和我的解决方案相似性不明显，或者这种求根方法不清楚，我还制作了一张图表，可以在这里找到。

这是一个工作的根函数，能够解决任何n次方根

（默认为平方根，为了这个问题）

#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)))));
}

解释：

点击此处查看图表

这是通过泰勒级数、对数性质和一些代数运算实现的。

例如：

log A = N
   x

*注意：对于平方根，N = 2；对于任何其他根，您只需要更改一个变量，即N。

1）更改底数，将以“x”为底的对数函数转换为自然对数。

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

2) 重新排列以分离ln(x)，并最终得出“x”。

ln(A)/N = ln(x)

3) 将两边都设置为'e'的指数，

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

4) 泰勒级数将“ln”表示为无限级数。

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] + . . .

*注意：为了提高精度，请继续进行该系列。为简洁起见，在我的函数中使用10^9来表示自然对数的级数收敛，其精度约为7位数字或百万分之十位数。

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

5) 现在，只需将自然对数方程式插入步骤3中得到的简化方程式中即可。

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

6) 这个实现可能看起来有些过度，但它的目的是展示如何在不必猜测和检查的情况下解决根问题。此外，它还可以让你用自己的pow函数替换cmath库中的pow函数：

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))));
}

- JReyn