如何改进这种平方根方法？答案

【问题标题】：How can I improve this square root method?如何改进这种平方根方法？
【发布时间】：2009-05-13 05:34:45
【问题描述】：

我知道这听起来像是一项家庭作业，但事实并非如此。最近我对用于执行某些数学运算的算法感兴趣，例如正弦、平方根等。目前，我正在尝试用 C# 编写计算平方根的Babylonian method。

到目前为止，我有这个：

public static double SquareRoot(double x) {
    if (x == 0) return 0;

    double r = x / 2; // this is inefficient, but I can't find a better way
                      // to get a close estimate for the starting value of r
    double last = 0;
    int maxIters = 100;

    for (int i = 0; i < maxIters; i++) {
        r = (r + x / r) / 2;
        if (r == last)
            break;
        last = r;
    }

    return r;
}

它工作得很好，每次都产生与 .NET Framework 的 Math.Sqrt() 方法完全相同的答案。不过，正如您可能猜到的那样，它比本机方法慢（大约 800 滴答声）。我知道这种特殊方法永远不会比本机方法快，但我只是想知道是否可以进行任何优化。

我立即看到的唯一优化是计算将运行 100 次，即使在已经确定答案之后（此时，r 始终是相同的值）。因此，我添加了一个快速检查，以查看新计算的值是否与先前计算的值相同并跳出循环。不幸的是，它在速度上并没有太大的区别，但似乎是正确的做法。

在你说“为什么不直接使用 Math.Sqrt() 来代替？”之前...我将其作为一个学习练习，并不打算在任何生产代码中实际使用此方法。

【问题讨论】：

标签： c# algorithm optimization math performance

【解决方案1】：

首先，您应该检查收敛性，而不是检查是否相等（r == last），其中 r 接近最后一个，其中 close 由任意 epsilon 定义：

eps = 1e-10  // pick any small number
if (Math.Abs(r-last) < eps) break;

正如您链接到的维基百科文章所提到的那样 - 您不能使用牛顿方法有效地计算平方根 - 相反，您使用对数。

【讨论】：

typo: s/"Newton's method"/"Babylonian method" -- 牛顿法适用于收敛速度（对它是否收敛有一些警告）
如果根大于 2^52*eps，那么 r 很可能围绕根振荡，并且 Math.Abs(r-last) 永远不会小于 eps。因此，虽然这个提议比原始程序好一点，但它仍然可能导致不必要的迭代。
这实际上减少了大约 100 个刻度，这看起来很奇怪，因为它正在执行一个附加方法以及一个比较。但是，我不会抱怨。 ;)
@David：我确定 Math.Abs 被 JIT 内联

【解决方案2】：

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

这是我最喜欢的快速平方根。实际上它是平方根的倒数，但是如果你愿意的话，你可以在之后把它取反...... .
http://www.beyond3d.com/content/articles/8/

【讨论】：

发疯了，虽然我认为这在 C# 中是不可能的
您可以创建一个联合来解决指针问题，只需使用StructLayoutAttribute 和LayoutKind.Explicit。

【解决方案3】：

你在这里做的是你执行Newton's method of finding a root。所以你可以使用一些更有效的寻根算法。你可以开始搜索它here。

【讨论】：

+1，算法改进通常胜过微优化，例如用位移位替换除法。
我看不出“使用不同的算法”对于“如何更快地执行此算法？”是一个很好的答案。他已经解释说他这样做只是因为他想这样做，所以告诉他完全做其他事情并不是一个有用的建议。
牛顿法收敛速度很快，根本不是问题。真正的问题是第一个近似值。 C# 似乎不允许 C/C++ 中可能的位摆弄。

【解决方案4】：

用位移代替除以 2 不太可能有那么大的区别；考虑到除法是一个常数，我希望编译器足够聪明，可以为你做到这一点，但你不妨试试看。

通过提前退出循环，您更有可能获得改进，因此要么将新 r 存储在变量中并与旧 r 进行比较，要么将 x/r 存储在变量中并在执行之前将其与 r 进行比较加法和除法。

【讨论】：

【解决方案5】：

您可以只返回 r，而不是中断循环然后返回 r。可能不会显着提高性能。

【讨论】：

位移位适用于 int （等） - 但它适用于 double 吗？它甚至似乎没有被定义......
"break/return" vs "return" so 最小；我不认为你会发现这里的区别
他正在努力节省蜱虫，所以我建议即使是最琐碎的事情。
+1！我希望我可以标记两个答案。使用“return r”而不是“break”肯定会产生速度差异（虽然非常小，但正如你所说，我在这里工作）。

【解决方案6】：

使用您的方法，每次迭代都会使正确位数翻倍。

使用表格获取最初的 4 位（例如），第一次迭代后您将有 8 位，第二次后有 16 位，第四次迭代后您需要的所有位（因为 double存储 52+1 位尾数）。

对于表查找，您可以提取 [0.5,1[ 中的尾数并从输入中提取指数（使用类似 frexp 的函数），然后使用适当的 2 幂乘以规范 [64,256[ 中的尾数。

mantissa *= 2^K
exponent -= K

在此之后，您输入的号码仍然是mantissa*2^exponent。 K 必须是 7 或 8，以获得偶数指数。您可以从包含尾数整数部分的所有平方根的表中获取迭代的初始值。执行 4 次迭代以获得尾数的平方根 r。结果是r*2^(exponent/2)，使用类似ldexp 的函数构造。

编辑。我在下面放了一些 C++ 代码来说明这一点。 OP 的函数 sr1 经过改进的测试需要 2.78s 来计算 2^24 平方根；我的函数 sr2 耗时 1.42s，硬件 sqrt 耗时 0.12s。

#include <math.h>
#include <stdio.h>

double sr1(double x)
{
  double last = 0;
  double r = x * 0.5;
  int maxIters = 100;
  for (int i = 0; i < maxIters; i++) {
    r = (r + x / r) / 2;
    if ( fabs(r - last) < 1.0e-10 )
      break;
    last = r;
  }
  return r;
}

double sr2(double x)
{
  // Square roots of values in 0..256 (rounded to nearest integer)
  static const int ROOTS256[] = {
    0,1,1,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9,
    9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,11,11,11,11,11,
    11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,12,12,12,
    12,12,12,12,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,
    13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,
    14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,
    15,15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16 };

  // Normalize input
  int exponent;
  double mantissa = frexp(x,&exponent); // MANTISSA in [0.5,1[ unless X is 0
  if (mantissa == 0) return 0; // X is 0
  if (exponent & 1) { mantissa *= 128; exponent -= 7; } // odd exponent
  else { mantissa *= 256; exponent -= 8; } // even exponent
  // Here MANTISSA is in [64,256[

  // Initial value on 4 bits
  double root = ROOTS256[(int)floor(mantissa)];

  // Iterate
  for (int it=0;it<4;it++)
    {
      root = 0.5 * (root + mantissa / root);
    }

  // Restore exponent in result
  return ldexp(root,exponent>>1);
}

int main()
{
  // Used to generate the table
  // for (int i=0;i<=256;i++) printf(",%.0f",sqrt(i));

  double s = 0;
  int mx = 1<<24;
  // for (int i=0;i<mx;i++) s += sqrt(i); // 0.120s
  // for (int i=0;i<mx;i++) s += sr1(i);  // 2.780s
  for (int i=0;i<mx;i++) s += sr2(i);  // 1.420s
}

【讨论】：

C# 中是否存在 frexp 和 ldexp？我不知道这些功能或任何可以替代它们的东西。 OP的解决方案的问题是在C#中很难找到一个初始近似值。
使用 Jon Carmack 的幻数进行近似：codemaestro.com/reviews/9
我在 Google 代码搜索上找到了 frexp 和 ldexp 的 C# 版本，但这个示例实际上比我的原始代码慢得多。当然，这也可能是 frexp 和 ldexp 实现的问题。然而，我发现这段代码真的很有趣。感谢发布！

【解决方案7】：

定义一个容差，并在后续迭代落在该容差范围内时尽早返回。

【讨论】：

【解决方案8】：

既然你说下面的代码不够快，试试这个：

    static double guess(double n)
    {
        return Math.Pow(10, Math.Log10(n) / 2);
    }

它应该非常准确，希望速度很快。

这是here 描述的初始估计代码。看起来还不错。使用此代码，然后您还应该进行迭代，直到值收敛在一个 epsilon 的差异内。

    public static double digits(double x)
    {
        double n = Math.Floor(x);
        double d;

        if (d >= 1.0)
        {
            for (d = 1; n >= 1.0; ++d)
            {
                n = n / 10;
            }
        }
        else
        {
            for (d = 1; n < 1.0; ++d)
            {
                n = n * 10;
            }
        }


        return d;
    }

    public static double guess(double x)
    {
        double output;
        double d = Program.digits(x);

        if (d % 2 == 0)
        {
            output = 6*Math.Pow(10, (d - 2) / 2);
        }
        else
        {
            output = 2*Math.Pow(10, (d - 1) / 2);
        }

        return output;
    }

【讨论】：

它有效，但计算时间比简单地使用 x / 2 长 3 倍。
你的意思是比 x/2 长 3 倍，还是整个程序？因为这应该比 x/2 给出更好的估计。

【解决方案9】：

出于学习目的，我也一直在研究这个。您可能对我尝试的两个修改感兴趣。

第一个是在 x0 中使用一阶泰勒级数逼近：

    Func<double, double> fNewton = (b) =>
    {
        // Use first order taylor expansion for initial guess
        // http://www27.wolframalpha.com/input/?i=series+expansion+x^.5
        double x0 = 1 + (b - 1) / 2;
        double xn = x0;
        do
        {
            x0 = xn;
            xn = (x0 + b / x0) / 2;
        } while (Math.Abs(xn - x0) > Double.Epsilon);
        return xn;
    };

第二个是尝试第三个订单（更贵），迭代

    Func<double, double> fNewtonThird = (b) =>
    {
        double x0 = b/2;
        double xn = x0;
        do
        {
            x0 = xn;
            xn = (x0*(x0*x0+3*b))/(3*x0*x0+b);
        } while (Math.Abs(xn - x0) > Double.Epsilon);
        return xn;
    };

我创建了一个辅助方法来为函数计时

public static class Helper
{
    public static long Time(
        this Func<double, double> f,
        double testValue)
    {
        int imax = 120000;
        double avg = 0.0;
        Stopwatch st = new Stopwatch();
        for (int i = 0; i < imax; i++)
        {
            // note the timing is strictly on the function
            st.Start();
            var t = f(testValue);
            st.Stop();
            avg = (avg * i + t) / (i + 1);
        }
        Console.WriteLine("Average Val: {0}",avg);
        return st.ElapsedTicks/imax;
    }
}

原来的方法更快，但同样，可能很有趣:)

【讨论】：

【解决方案10】：

将“/ 2”替换为“* 0.5”使我的机器上的速度提高了约 1.5 倍，但当然不如本机实现快。

【讨论】：

【解决方案11】：

好吧，原生 Sqrt() 函数可能没有用 C# 实现，它很可能用低级语言完成，而且肯定会使用更有效的算法。所以试图跟上它的速度可能是徒劳的。

但是，关于仅尝试优化 heckuvit 的功能，您链接的 Wikipedia 页面建议“起始猜测”为 2^floor(D/2)，其中 D 表示二进制数字的数量数字。您可以尝试一下，我没有看到您的代码中可以显着优化的其他内容。

【讨论】：

【解决方案12】：

你可以试试 r = x >> 1;

而不是 / 2（也可以在您设备 2 的其他位置）。它可能会给你一点优势。我还将100 移动到循环中。可能没什么，但我们在这里讨论的是蜱虫。

现在检查一下。

编辑：将 > 修复为 >>，但它不适用于双打，所以没关系。 100 的内联没有给我速度提升。

【讨论】：

我认为这行不通，因为 x>1 将是“真”或“假”，它应该是 >>。