计算大整数(System.Numerics.BigInteger)的平方根

40

.NET 4.0 提供了 System.Numerics.BigInteger 类型来处理任意大的整数。我需要计算一个 BigInteger 的平方根(或者一个合理的近似值,例如整数平方根)。为了不必重复造轮子,请问是否有一个好的扩展方法可以做到这一点?


计算任意精度平方根的最简单可行方法可能是牛顿法 - mqp
抱歉,我刚开始思考这个背后的数学原理就感到脑袋疼 :-P。而且这些数字太大了,无法转换为长整型? - Alxandr
1
是的,我需要大约256位,可能是512位 - 所以不要使用无符号长整型进行欺骗。 - Anonym
8个回答

25

Check if BigInteger is not a perfect square提供了计算Java BigInteger的整数平方根的代码。 这里将其翻译为C#扩展方法。

    public static BigInteger Sqrt(this BigInteger n)
    {
        if (n == 0) return 0;
        if (n > 0)
        {
            int bitLength = Convert.ToInt32(Math.Ceiling(BigInteger.Log(n, 2)));
            BigInteger root = BigInteger.One << (bitLength / 2);

            while (!isSqrt(n, root))
            {
                root += n / root;
                root /= 2;
            }

            return root;
        }

        throw new ArithmeticException("NaN");
    }

    private static Boolean isSqrt(BigInteger n, BigInteger root)
    {
        BigInteger lowerBound = root*root;
        BigInteger upperBound = (root + 1)*(root + 1);

        return (n >= lowerBound && n < upperBound);
    }

非正式测试表明,对于小整数而言,这比Math.Sqrt慢了约75倍。VS分析器指出isSqrt中的乘法是热点。


7
BigInteger不优化除法运算符。使用右移一位操作符代替除以2可以提高性能(至少在我的情况下)。 - GeirGrusom
3
UpperBound的定义也可以重写为多项式展开形式,即BigInteger upperBound = lowerBound + root + root + 1,或者内联在返回语句中,如return n >= lowerBound && n <= lowerBound + root + root。请注意,翻译保持原文意思,尽量使语言通俗易懂,不包含任何解释或其他内容。 - Jesan Fafon

15
我不确定牛顿法是否是计算大数平方根的最佳方法,因为它涉及到除法运算,对于大数来说速度较慢。你可以使用CORDIC方法,它只使用加法和位移运算(这里展示的是无符号整数的情况)。
static uint isqrt(uint x)
{
    int b=15; // this is the next bit we try 
    uint r=0; // r will contain the result
    uint r2=0; // here we maintain r squared
    while(b>=0) 
    {
        uint sr2=r2;
        uint sr=r;
                    // compute (r+(1<<b))**2, we have r**2 already.
        r2+=(uint)((r<<(1+b))+(1<<(b+b)));      
        r+=(uint)(1<<b);
        if (r2>x) 
        {
            r=sr;
            r2=sr2;
        }
        b--;
    }
    return r;
}

有一种类似的方法只使用加法和移位,称为“迪杰斯特拉平方根”,例如在这里解释:http://lib.tkk.fi/Diss/2005/isbn9512275279/article3.pdf

这个程序计算一个整数的整数平方根。如果你需要小数,你可以预先缩放操作数。 - Nordic Mainframe
你可以通过继续循环负数 b 的值并将 -n 的左移转换为 n 的右移来计算任意精度。 - Chris Dodd
很容易适应64位长整型,这正是我所需要的。谢谢! - yoyo
那么如何进行扩展?将 b 更改为什么? - Ahmed Fwela
通过缩放,我指的是使用 BigInteger 到任意位数长度。 - Ahmed Fwela

9

好的,首先进行一些速度测试,测试一些在此发布的变量。 (我只考虑那些可以给出精确结果并且适用于 BigInteger 的方法):

+------------------------------+-------+------+------+-------+-------+--------+--------+--------+
| variant - 1000x times        |   2e5 | 2e10 | 2e15 |  2e25 |  2e50 |  2e100 |  2e250 |  2e500 |
+------------------------------+-------+------+------+-------+-------+--------+--------+--------+
| my version                   |  0.03 | 0.04 | 0.04 |  0.76 |  1.44 |   2.23 |   4.84 |  23.05 |
| RedGreenCode (bound opti.)   |  0.56 | 1.20 | 1.80 |  2.21 |  3.71 |   6.10 |  14.53 |  51.48 |
| RedGreenCode (newton method) |  0.80 | 1.21 | 2.12 |  2.79 |  5.23 |   8.09 |  19.90 |  65.36 |
| Nordic Mainframe (CORDIC)    |  2.38 | 5.52 | 9.65 | 19.80 | 46.69 |  90.16 | 262.76 | 637.82 |
| Sunsetquest (without divs)   |  2.37 | 5.48 | 9.11 | 24.50 | 56.83 | 145.52 | 839.08 | 4.62 s |
| Jeremy Kahan (js-port)       | 46.53 | #.## | #.## |  #.## |  #.## |   #.## |   #.## |   #.## |
+------------------------------+-------+------+------+-------+-------+--------+--------+--------+

+------------------------------+--------+--------+--------+---------+---------+--------+--------+
| variant - single             | 2e1000 | 2e2500 | 2e5000 | 2e10000 | 2e25000 |  2e50k | 2e100k |
+------------------------------+--------+--------+--------+---------+---------+--------+--------+
| my version                   |   0.10 |   0.77 |   3.46 |   14.97 |  105.19 | 455.68 | 1,98 s |
| RedGreenCode (bound opti.)   |   0.26 |   1.41 |   6.53 |   25.36 |  182.68 | 777.39 | 3,30 s |
| RedGreenCode (newton method) |   0.33 |   1.73 |   8.08 |   32.07 |  228.50 | 974.40 | 4,15 s |
| Nordic Mainframe (CORDIC)    |   1.83 |   7.73 |  26.86 |   94.55 |  561.03 | 2,25 s | 10.3 s |
| Sunsetquest (without divs)   |  31.84 | 450.80 | 3,48 s |  27.5 s |    #.## |   #.## |   #.## |
| Jeremy Kahan (js-port)       |   #.## |   #.## |   #.## |    #.## |    #.## |   #.## |   #.## |
+------------------------------+--------+--------+--------+---------+---------+--------+--------+

- value example: 2e10 = 20000000000 (result: 141421)
- times in milliseconds or with "s" in seconds
- #.##: need more than 5 minutes (timeout)

描述:

Jeremy Kahan (js-port)

Jeremy的简单算法有效,但由于简单的加减操作,计算工作量呈指数级快速增长... :)

Sunsetquest (without divs)

不分割的方法很好,但由于分而治之的变种,结果收敛相对较慢(特别是对于大的数字)。

Nordic Mainframe (CORDIC)

CORDIC算法已经非常强大,但不可变的BigIntegers逐位操作会产生许多开销。

我是这样计算所需位数的:int b = Convert.ToInt32(Math.Ceiling(BigInteger.Log(x, 2))) / 2 + 1;

RedGreenCode (newton method)

经过验证的牛顿法表明,旧事物并不一定缓慢。尤其是大量数据的快速收敛几乎无法被超越。

RedGreenCode (bound opti.)

Jesan Fafon提出的节省乘法的建议引起了很大的反响。

my version

首先,使用Math.Sqrt()计算小数字,一旦double类型的精度不足,就使用牛顿迭代法。但是,我会尝试使用Math.Sqrt()预先计算尽可能多的数字,这样可以使牛顿迭代法更快地收敛。

以下是源代码:

static readonly BigInteger FastSqrtSmallNumber = 4503599761588223UL; // as static readonly = reduce compare overhead

static BigInteger SqrtFast(BigInteger value)
{
  if (value <= FastSqrtSmallNumber) // small enough for Math.Sqrt() or negative?
  {
    if (value.Sign < 0) throw new ArgumentException("Negative argument.");
    return (ulong)Math.Sqrt((ulong)value);
  }

  BigInteger root; // now filled with an approximate value
  int byteLen = value.ToByteArray().Length;
  if (byteLen < 128) // small enough for direct double conversion?
  {
    root = (BigInteger)Math.Sqrt((double)value);
  }
  else // large: reduce with bitshifting, then convert to double (and back)
  {
    root = (BigInteger)Math.Sqrt((double)(value >> (byteLen - 127) * 8)) << (byteLen - 127) * 4;
  }

  for (; ; )
  {
    var root2 = value / root + root >> 1;
    if ((root2 == root || root2 == root + 1) && IsSqrt(value, root)) return root;
    root = value / root2 + root2 >> 1;
    if ((root == root2 || root == root2 + 1) && IsSqrt(value, root2)) return root2;
  }
}

static bool IsSqrt(BigInteger value, BigInteger root)
{
  var lowerBound = root * root;

  return value >= lowerBound && value <= lowerBound + (root << 1);
}

完整的基准源代码:

using System;
using System.Numerics;
using System.Diagnostics;

namespace MathTest
{
  class Program
  {
    static readonly BigInteger FastSqrtSmallNumber = 4503599761588223UL; // as static readonly = reduce compare overhead
    static BigInteger SqrtMax(BigInteger value)
    {
      if (value <= FastSqrtSmallNumber) // small enough for Math.Sqrt() or negative?
      {
        if (value.Sign < 0) throw new ArgumentException("Negative argument.");
        return (ulong)Math.Sqrt((ulong)value);
      }

      BigInteger root; // now filled with an approximate value
      int byteLen = value.ToByteArray().Length;
      if (byteLen < 128) // small enough for direct double conversion?
      {
        root = (BigInteger)Math.Sqrt((double)value);
      }
      else // large: reduce with bitshifting, then convert to double (and back)
      {
        root = (BigInteger)Math.Sqrt((double)(value >> (byteLen - 127) * 8)) << (byteLen - 127) * 4;
      }

      for (; ; )
      {
        var root2 = value / root + root >> 1;
        if ((root2 == root || root2 == root + 1) && IsSqrt(value, root)) return root;
        root = value / root2 + root2 >> 1;
        if ((root == root2 || root == root2 + 1) && IsSqrt(value, root2)) return root2;
      }
    }

    static bool IsSqrt(BigInteger value, BigInteger root)
    {
      var lowerBound = root * root;

      return value >= lowerBound && value <= lowerBound + (root << 1);
    }

    // newton method
    public static BigInteger SqrtRedGreenCode(BigInteger n)
    {
      if (n == 0) return 0;
      if (n > 0)
      {
        int bitLength = Convert.ToInt32(Math.Ceiling(BigInteger.Log(n, 2)));
        BigInteger root = BigInteger.One << (bitLength / 2);

        while (!isSqrtRedGreenCode(n, root))
        {
          root += n / root;
          root /= 2;
        }

        return root;
      }

      throw new ArithmeticException("NaN");
    }
    private static bool isSqrtRedGreenCode(BigInteger n, BigInteger root)
    {
      BigInteger lowerBound = root * root;
      //BigInteger upperBound = (root + 1) * (root + 1);

      return n >= lowerBound && n <= lowerBound + root + root;
      //return (n >= lowerBound && n < upperBound);
    }

    // without divisions
    public static BigInteger SqrtSunsetquest(BigInteger number)
    {
      if (number < 9)
      {
        if (number == 0)
          return 0;
        if (number < 4)
          return 1;
        else
          return 2;
      }

      BigInteger n = 0, p = 0;
      var high = number >> 1;
      var low = BigInteger.Zero;

      while (high > low + 1)
      {
        n = (high + low) >> 1;
        p = n * n;
        if (number < p)
        {
          high = n;
        }
        else if (number > p)
        {
          low = n;
        }
        else
        {
          break;
        }
      }
      return number == p ? n : low;
    }

    // javascript port
    public static BigInteger SqrtJeremyKahan(BigInteger n)
    {
      var oddNumber = BigInteger.One;
      var result = BigInteger.Zero;
      while (n >= oddNumber)
      {
        n -= oddNumber;
        oddNumber += 2;
        result++;
      }
      return result;
    }

    // CORDIC
    public static BigInteger SqrtNordicMainframe(BigInteger x)
    {
      int b = Convert.ToInt32(Math.Ceiling(BigInteger.Log(x, 2))) / 2 + 1;
      BigInteger r = 0; // r will contain the result
      BigInteger r2 = 0; // here we maintain r squared
      while (b >= 0)
      {
        var sr2 = r2;
        var sr = r;
        // compute (r+(1<<b))**2, we have r**2 already.
        r2 += (r << 1 + b) + (BigInteger.One << b + b);
        r += BigInteger.One << b;
        if (r2 > x)
        {
          r = sr;
          r2 = sr2;
        }
        b--;
      }
      return r;
    }

    static void Main(string[] args)
    {
      var t2 = BigInteger.Parse("2" + new string('0', 10000));

      //var q1 = SqrtRedGreenCode(t2);
      var q2 = SqrtSunsetquest(t2);
      //var q3 = SqrtJeremyKahan(t2);
      //var q4 = SqrtNordicMainframe(t2);
      var q5 = SqrtMax(t2);
      //if (q5 != q1) throw new Exception();
      if (q5 != q2) throw new Exception();
      //if (q5 != q3) throw new Exception();
      //if (q5 != q4) throw new Exception();

      for (int r = 0; r < 5; r++)
      {
        var mess = Stopwatch.StartNew();
        //for (int i = 0; i < 1000; i++)
        {
          //var q = SqrtRedGreenCode(t2);
          var q = SqrtSunsetquest(t2);
          //var q = SqrtJeremyKahan(t2);
          //var q = SqrtNordicMainframe(t2);
          //var q = SqrtMax(t2);
        }
        mess.Stop();
        Console.WriteLine((mess.ElapsedTicks * 1000.0 / Stopwatch.Frequency).ToString("N2") + " ms");
      }
    }
  }
}

你好,有点晚了,但我现在已经添加了基准测试代码。测试时间仍然在27-28秒之间。计算结果也是正确的。如果有性能问题,请告诉我...;-) - MaxKlaxx
感谢您的发布。我无法让ngMax.Zp.TickCount正常工作,所以我使用了内置的“Stopwatch timer = Stopwatch.StartNew(); ..... timer.Stop(); Console.WriteLine(timer.ElapsedMilliseconds);”,仍然得到3.5秒的结果...但是您的结果同时也只有11毫秒。看起来要么是“ngMax”的问题,要么就是我的电脑快了8倍(但我不认为是这样的!)。无论如何,基准比例似乎是准确的 - 这才是最重要的。在一个快速函数上做得好! - SunsetQuest
1
ngMax.Zp.TickCount返回一个(双精度)毫秒时间(带小数点)。有趣的是,这并不是那样:我将我的项目从.Net 4.5(vs2013)升级到了.Net Core 3.1(vs2019),现在我也能够达到3.8-3.9秒。BigInteger本身在过去几年中在性能方面得到了很大的改进。^^ - MaxKlaxx
我注意到一些舍入问题:Sqrt(4503599761588224)返回67108865,但应该是67108864; Sqrt(9511487747939116817576)返回97526856546,但应该是97526856547。 - SunsetQuest
1
好的,舍入误差现在应该已经修复了。再次感谢您的提示! - MaxKlaxx
显示剩余2条评论

3

简短回答:(但请注意,下面有更多详细信息)

Math.Pow(Math.E, BigInteger.Log(pd) / 2)

其中pd代表您想要执行平方根操作的BigInteger

详细回答和解释:

另一种理解这个问题的方法是理解平方根和对数的工作原理。

如果你有方程5^x = 25,要解出x,我们必须使用对数。在这个例子中,我将使用自然对数(其他底数的对数也是可能的,但自然对数是简单的方法)。

5^x = 25

重写后,我们有:
x(ln 5) = ln 25

为了隔离x,我们需要:
x = ln 25 / ln 5

我们看到这个结果是 x = 2。但既然我们已经知道了 x(在 5^2 中,x = 2),让我们改变我们不知道的部分并编写一个新方程来解决新的未知数。让 x 成为平方根操作的结果。这给我们带来了:
2 = ln 25 / ln x

重新组织以隔离变量x,我们得到:
ln x = (ln 25) / 2

为了移除日志,我们现在使用自然对数和特殊数字e的特殊身份。具体地说, e^ln x = x。现在重新编写方程,得到:
e^ln x = e^((ln 25) / 2)

简化左侧,我们有

x = e^((ln 25) / 2)

x将是25的平方根。您还可以将此想法扩展到任何根或数字,对于x的第y个根的一般公式成为e^((ln x) / y)

现在特别应用于C#,BigIntegers和这个问题,我们只需实现该公式。 警告:尽管数学是正确的,但存在有限限制。此方法仅会使您接近答案,且有一个大的未知范围(取决于您操作的数字大小)。也许这就是为什么Microsoft没有实现这样的方法的原因。

// A sample generated public key modulus
var pd = BigInteger.Parse("101017638707436133903821306341466727228541580658758890103412581005475252078199915929932968020619524277851873319243238741901729414629681623307196829081607677830881341203504364437688722228526603134919021724454060938836833023076773093013126674662502999661052433082827512395099052335602854935571690613335742455727");
var sqrt = Math.Pow(Math.E, BigInteger.Log(pd) / 2);

Console.WriteLine(sqrt);

注意: BigInteger.Log() 方法返回一个双精度浮点数,因此会有两个问题。1)这个数字是不精确的,2)Log() 对于 BigInteger 输入有一个上限。为了检查上限,我们可以查看自然对数的正常形式,即 ln x = y。换句话说,e^y = x。由于 doubleBigInteger.Log() 的返回类型,最大的 BigInteger 应该是 edouble.MaxValue 次幂。在我的电脑上,那将是 e^1.79769313486232E+308。不精确性没有得到解决。有人想实现 BigDecimal 并更新 BigInteger.Log() 吗?
消费者要谨慎,但它能让你接近答案,将结果平方确实会产生类似于原始输入的数字,只是精度不如 RedGreenCode 的答案。愉快地(平方)开根!;)

2
您可以将此转换为您选择的语言和变量类型。这是一个在JavaScript中截断的平方根(对我来说最新鲜),它利用了1 + 3 + 5... +第n个奇数 = n ^ 2。所有变量都是整数,只进行加法和减法运算。
var truncSqrt = function(n) {
  var oddNumber = 1;
  var result = 0;
  while (n >= oddNumber) {
    n -= oddNumber;
    oddNumber += 2;
    result++;
  }
  return result;
};

1
真的很好奇这种方法相对于其他方法的表现如何。 - Jeremy Kahan

2

更新:为了获得最佳性能,请使用Newton Plus版本。

那个版本比这个快数百倍。我将保留这个作为一种备选方式进行参考。

    // Source: http://mjs5.com/2016/01/20/c-biginteger-square-root-function/  Michael Steiner, Jan 2016
    // Slightly modified to correct error below 6. (thank you M Ktsis D) 
    public static BigInteger SteinerSqrt(BigInteger number)
    {
        if (number < 9)
        {
            if (number == 0)
                return 0;
            if (number < 4)
                return 1;
            else
                return 2;
        }

        BigInteger n = 0, p = 0;
        var high = number >> 1;
        var low = BigInteger.Zero;

        while (high > low + 1)
        {
            n = (high + low) >> 1;
            p = n * n;
            if (number < p)
            {
                high = n;
            }
            else if (number > p)
            {
                low = n;
            }
            else
            {
                break;
            }
        }
        return number == p ? n : low;
    }

更新: 感谢 M Ktsis D 发现了一个错误。已经通过守卫条款进行了更正。


1
添加针对数字1和4的守卫条款将使这个答案完美。 - Marlon Dumal-is
非常感谢!你是一个英雄。 - killereks

1

***世界上最快的Java/C#大整数平方根算法!!!***

详细介绍请查看:https://www.codeproject.com/Articles/5321399/NewtonPlus-A-Fast-Big-Number-Square-Root-Function

Github地址:https://github.com/SunsetQuest/NewtonPlus-Fast-BigInteger-and-BigFloat-Square-Root

public static BigInteger NewtonPlusSqrt(BigInteger x)
{
    if (x < 144838757784765629)    // 1.448e17 = ~1<<57
    {
        uint vInt = (uint)Math.Sqrt((ulong)x);
        if ((x >= 4503599761588224) && ((ulong)vInt * vInt > (ulong)x))  // 4.5e15 =  ~1<<52
        {
            vInt--;
        }
        return vInt;
    }

    double xAsDub = (double)x;
    if (xAsDub < 8.5e37)   //  long.max*long.max
    {
        ulong vInt = (ulong)Math.Sqrt(xAsDub);
        BigInteger v = (vInt + ((ulong)(x / vInt))) >> 1;
        return (v * v <= x) ? v : v - 1;
    }

    if (xAsDub < 4.3322e127)
    {
        BigInteger v = (BigInteger)Math.Sqrt(xAsDub);
        v = (v + (x / v)) >> 1;
        if (xAsDub > 2e63)
        {
            v = (v + (x / v)) >> 1;
        }
        return (v * v <= x) ? v : v - 1;
    }

    int xLen = (int)x.GetBitLength();
    int wantedPrecision = (xLen + 1) / 2;
    int xLenMod = xLen + (xLen & 1) + 1;

    //////// Do the first Sqrt on hardware ////////
    long tempX = (long)(x >> (xLenMod - 63));
    double tempSqrt1 = Math.Sqrt(tempX);
    ulong valLong = (ulong)BitConverter.DoubleToInt64Bits(tempSqrt1) & 0x1fffffffffffffL;
    if (valLong == 0)
    {
        valLong = 1UL << 53;
    }

    ////////  Classic Newton Iterations ////////
    BigInteger val = ((BigInteger)valLong << 52) + (x >> xLenMod - (3 * 53)) / valLong;
    int size = 106;
    for (; size < 256; size <<= 1)
    {
        val = (val << (size - 1)) + (x >> xLenMod - (3 * size)) / val;
    }

    if (xAsDub > 4e254) // 4e254 = 1<<845.76973610139
    {
        int numOfNewtonSteps = BitOperations.Log2((uint)(wantedPrecision / size)) + 2;

        //////  Apply Starting Size  ////////
        int wantedSize = (wantedPrecision >> numOfNewtonSteps) + 2;
        int needToShiftBy = size - wantedSize;
        val >>= needToShiftBy;
        size = wantedSize;
        do
        {
            ////////  Newton Plus Iterations  ////////
            int shiftX = xLenMod - (3 * size);
            BigInteger valSqrd = (val * val) << (size - 1);
            BigInteger valSU = (x >> shiftX) - valSqrd;
            val = (val << size) + (valSU / val);
            size *= 2;
        } while (size < wantedPrecision);
    }

    /////// There are a few extra digits here, lets save them ///////
    int oversidedBy = size - wantedPrecision;
    BigInteger saveDroppedDigitsBI = val & ((BigInteger.One << oversidedBy) - 1);
    int downby = (oversidedBy < 64) ? (oversidedBy >> 2) + 1 : (oversidedBy - 32);
    ulong saveDroppedDigits = (ulong)(saveDroppedDigitsBI >> downby);


    ////////  Shrink result to wanted Precision  ////////
    val >>= oversidedBy;


    ////////  Detect a round-ups  ////////
    if ((saveDroppedDigits == 0) && (val * val > x))
    {
        val--;
    }

    ////////// Error Detection ////////
    //// I believe the above has no errors but to guarantee the following can be added.
    //// If an error is found, please report it.
    //BigInteger tmp = val * val;
    //if (tmp > x)
    //{
    //    Console.WriteLine($"Missed  , {ToolsForOther.ToBinaryString(saveDroppedDigitsBI, oversidedBy)}, {oversidedBy}, {size}, {wantedPrecision}, {saveDroppedDigitsBI.GetBitLength()}");
    //    if (saveDroppedDigitsBI.GetBitLength() >= 6)
    //        Console.WriteLine($"val^2 ({tmp}) < x({x})  off%:{((double)(tmp)) / (double)x}");
    //    //throw new Exception("Sqrt function had internal error - value too high");
    //}
    //if ((tmp + 2 * val + 1) <= x)
    //{
    //    Console.WriteLine($"(val+1)^2({((val + 1) * (val + 1))}) >= x({x})");
    //    //throw new Exception("Sqrt function had internal error - value too low");
    //}

    return val;
}

以下是基于日志的图表。请注意,细微的差异在性能方面有很大的影响。除了为了比较而添加的GMP(C++/Asm)外,所有都使用C#编写,还添加了Java版本(转换为C#)。 enter image description here

1
我把你的函数移植到了Go(big.Int),由于更快的内部除法,你的算法也运行得更快了:1e39457: .Net: 6.37 毫秒 / Go: 261 微秒1e2525223: .Net: 23.1 秒 / Go: 173 毫秒 - undefined
为了完整性:1e2525223: Math.Gmp.Native(GMP 6.x的.Net端口):48毫秒 - undefined
@MaxKlaxx - 哇,我知道Go很快,但是25倍到125倍真的令人惊讶。这已经接近GMP(原生的C++,可能还有手动调整的128位宽度SSE汇编指令)了,只是性能差了3倍。感谢你将它移植到Go平台上。如果你在网上分享,请分享链接,我想看看。将它移植到C++也很酷,但是我已经有20年没有用那种语言编程了,而且也没有一个直接的BigInteger类。 - undefined
1
Math.Gmp.Native:是的,"wrapper" 是正确的术语。SQRT 和 Go:Go 本身的原生版本在内部使用了 PaulZimmermann 算法,比 .Net 中的 NewtonPlusSqrt 稍微快一点。然而,Go 中的 NewtonPlusSqrt 速度极快,接近于 GMP。至于发布方面,我还没有对该版本进行足够的测试(特别是对较小值的优化)。顺便说一下,我从 2003 年开始使用 .Net,使用了将近 20 年,去年四月份才自愿转向 Go,所以它值得一看;-)。 - undefined

1
以下两种方法使用巴比伦方法来计算所提供数字的平方根。Sqrt方法返回BigInteger类型,因此只会提供最后一个整数的答案(没有小数点)。
该方法将使用15次迭代,尽管在进行了几次测试后,我发现12-13次迭代就足以处理80位以上的数字,但我决定将其保持在15次以防万一。
由于巴比伦平方根逼近方法要求我们选择一个长度为要找到平方根的数字长度的一半的数字,因此RandomBigIntegerOfLength()方法提供了该数字。
RandomBigIntegerOfLength()以整数长度作为参数,并提供该长度的随机生成数字。该数字是使用Random类的Next()方法生成的,为了避免数字在开头有0(例如041657180501613764193159871),从而导致DivideByZeroException,Next()方法被调用两次。需要指出的重要一点是,最初数字是一个接一个地生成,然后从字符串转换为BigInteger类型。

Sqrt 方法使用 RandomBigIntegerOfLength 方法获取提供的参数 "number" 长度的一半的随机数,然后使用巴比伦方法进行 15 次迭代计算平方根。您可以更改迭代次数为更小或更大的值。由于巴比伦方法不能提供 0 的平方根,因为它需要除以 0,在提供 0 作为参数时,它将返回 0。

//Copy the two methods
public static BigInteger Sqrt(BigInteger number)
{
    BigInteger _x = RandomBigIntegerOfLength((number.ToString().ToCharArray().Length / 2));
    
    try
    {
        for (int i = 0; i < 15; i++)
        {
            _x = (_x + number / _x) / 2;
        }
        return _x;
    }
    catch (DivideByZeroException)
    {
        return 0;
    }
}

// Copy this method as well
private static BigInteger RandomBigIntegerOfLength(int length)
{
    Random rand = new Random();

    string _randomNumber = "";
    
    _randomNumber = String.Concat(_randomNumber, rand.Next(1, 10));

    for (int i = 0; i < length-1; i++)
    {
        _randomNumber = String.Concat(_randomNumber,rand.Next(10).ToString());
    }
    if (String.IsNullOrEmpty(_randomNumber) == false) return BigInteger.Parse(_randomNumber);
    else return 0;
}

请勿只发布代码答案。解释代码的作用及其实现方式。 - M-Chen-3

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接