如何改进这个求平方根的方法？

Question

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

到目前为止，我有以下代码：

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个刻度）。我知道这种特定的方法永远不会比本地方法更快，但我想知道是否有任何优化可以进行。

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

在你说“为什么不使用Math.Sqrt()呢？”之前……我正在进行这个学习练习，并不打算在任何生产代码中实际使用这个方法。

- David Brown

Answer 1

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

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

正如你提供的维基百科文章所述，使用牛顿迭代法不能高效地计算平方根，相反，应该使用对数。

- Ankur Goel

---

笔误：将“牛顿迭代法”更正为“巴比伦迭代法”-- 尽管牛顿迭代法在收敛速度方面表现良好（但是否收敛存在一些限制）。 - Jason S

如果根大于2^52*eps，则r可能会在根周围振荡，并且Math.Abs(r-last)永远不会小于eps。因此，虽然这个提案比原始程序稍微好一点，但仍可能导致比必要的迭代次数多得多。 - Accipitridae

这实际上可以节省约100个时钟周期，这似乎很奇怪，因为它执行了一个额外的方法以及一个比较。但是，我不会抱怨。;) - David Brown

@David：我相信Math.Abs会被JIT内联。 - BlueRaja - Danny Pflughoeft

Answer 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/

- Chris H

---

真是疯狂，但我认为在C#中这甚至不可能。 - BlueRaja - Danny Pflughoeft

你可以通过使用 StructLayoutAttribute 和 LayoutKind.Explicit 来创建一个联合体，从而解决指针问题。 - John Leidegren

Answer 3

你在这里做的是执行牛顿迭代法寻找根。因此，你可以使用一些更有效的寻根算法。你可以从这里开始搜索。

- Jasiu

---

算法优化通常比微小的优化更加重要，例如用位移替代除法。 - Richard

“我不认为‘使用另一个算法’是如何更快地执行该算法的好答案？他已经解释过了，他只是这样做是因为他想这样做，所以告诉他完全做另一件事情并没有什么用。” - Chad Birch

牛顿法收敛速度快，根本不是问题所在。真正的问题在于第一次近似值。C#似乎不允许进行C/C++中可能存在的位操作。 - Accipitridae

Answer 4

通过您的方法，每次迭代都会使正确位数加倍。

使用表格获取初始的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 次迭代以获取 mantissa 的平方根 r。结果为 r*2^(exponent/2)，使用类似于 ldexp 的函数构造。

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

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

- Eric Bainville

---

在C#中是否存在frexp和ldexp函数？我不知道这些函数或任何可以替代它们的东西。 OP的解决方案问题在于，在C#中很难找到初始近似值。 - Accipitridae

使用Jon Carmack的魔数进行近似：http://www.codemaestro.com/reviews/9 - Ankur Goel

我在Google Code Search上找到了frexp和ldexp的C#版本，但是这个例子实际上比我的原始代码要慢得多。当然，这也可能是frexp和ldexp实现的问题。不过，我觉得这段代码非常有趣。感谢您的发布！ - David Brown

Answer 5

不需要打破循环然后返回r，你可以直接返回r。这样做可能不会提供明显的性能增益。

- AlbertoPL

---

位移对于int等类型完全有效 - 但是它对于double类型有效吗？ 它甚至似乎没有定义... - Marc Gravell

“break/return”与“return”的区别非常微小；我认为你永远不会在这里发现差异。 - Marc Gravell

他试图节省时间，所以我建议甚至最琐碎的事情。 - AlbertoPL

+1！我希望我能够标记两个答案。使用“return r”而不是“break”肯定会使速度有所提高（尽管非常小，但正如您所说，我在这里处理的是时刻）。 - David Brown

Answer 6

将除以2的操作改为位移操作可能不会有太大的差别；考虑到这是一个常量除法，我希望编译器足够聪明以自动进行优化，但您可以尝试一下。

退出循环时，通过存储新的r值并与旧的r值进行比较，或者将x/r存储在一个变量中并在执行加法和除法之前与r进行比较，这样更有可能获得提高。

- Richard

Answer 7

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

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

应该非常准确，而且希望速度快。

这里是描述此处的初始估计的代码。它看起来相当不错。使用此代码，然后您还应该迭代，直到值在差异的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;
    }

- Unknown

---

它可以工作，但计算速度比简单地使用x / 2要慢三倍。 - David Brown

你的意思是它比x/2慢三倍，还是整个程序都要慢三倍？因为这样会比x/2给出更好的估计。 - Unknown

Answer 8

我也一直在研究这个，以便学习。你可能会对我尝试的两个修改感兴趣。

第一个修改是在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;
    }
}

原始方法更快，但是可能会有趣 :)

- ccook

Answer 9

定义一个公差，并在后续迭代在该公差范围内时提前返回。

- Sinan Ünür

Answer 10

将“/ 2”替换为“* 0.5”可以使我的机器运行速度快约1.5倍，但当然不如本地实现快。

- Kristof Verbiest