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

Question

.NET 4.0提供了System.Numerics.BigInteger任意大整数的类型.我需要计算a的平方根(或合理的近似值 - 例如,整数平方根)BigInteger.所以我没有重新实现轮子,有没有人有一个很好的扩展方法呢？

Answer 1

检查BigInteger是不是一个完美的正方形,有代码来计算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中的乘法作为热点.

Answer 2

我不确定Newton's Method是否是计算bignum平方根的最佳方法,因为它涉及对bignum来说缓慢的划分.您可以使用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;
}

有一种类似的方法只使用加法和移位,称为"Dijkstras Square Root",例如这里解释:

• http://lib.tkk.fi/Diss/2005/isbn9512275279/article3.pdf

Answer 3

最后更新：2023-09-18（添加 SunsetQuests“世界最快”NewtonPlusSqrt）

好的，首先对这里发布的一些变体进行一些速度测试。（我只考虑给出准确结果并且至少适合 BigInteger 的方法）：

+------------------------------+-------+--------+------+-------+-------+-------+--------+--------+
| variant - 1000x times        |   2e5 |   2e10 | 2e15 |  2e25 |  2e50 | 2e100 |  2e250 |  2e500 |
+------------------------------+-------+--------+------+-------+-------+-------+--------+--------+
| NewtonPlusSqrt (SunsetQuest) |  0.02 |   0.03 | 0.03 |  0.25 |  0.40 |  0.72 |   1.13 |   3.13 |
| SqrtMax (my version)         |  0.02 |   0.03 | 0.03 |  0.53 |  1.01 |  1.17 |   2.46 |   6.89 |
| RedGreenCode (newton method) |  0.24 |   0.62 | 0.82 |  1.30 |  1.56 |  2.35 |   5.09 |  12.06 |
| Nordic Mainframe (CORDIC)    |  1.48 |   3.52 | 5.58 | 10.84 | 22.72 | 44.02 | 124.65 | 345.51 |
| SteinerSqrt (without divs)   |  1.60 |   3.16 | 5.61 | 14.30 | 31.64 | 66.23 | 301.14 | 1.37 s |
| Jeremy Kahan (js-port)       | 22.06 | 8.12 s | #.## |  #.## |  #.## |  #.## |   #.## |   #.## |
+------------------------------+-------+--------+------+-------+-------+-------+--------+--------+

+------------------------------+--------+--------+--------+---------+---------+--------+--------+
| variant - single             | 2e1000 | 2e2500 | 2e5000 | 2e10000 | 2e25000 |  2e50k | 2e100k |
+------------------------------+--------+--------+--------+---------+---------+--------+--------+
| NewtonPlusSqrt (SunsetQuest) |   0.02 |   0.04 |   0.12 |    0.44 |    2.63 |   9.55 |  37.76 |
| SqrtMax (my version)         |   0.08 |   0.57 |   2.37 |   10.51 |   68.70 | 303.74 | 1.31 s |
| RedGreenCode (newton method) |   0.15 |   0.79 |   3.53 |   13.65 |   92.83 | 395.66 | 1.66 s |
| Nordic Mainframe (CORDIC)    |   1.21 |   5.36 |  19.54 |   61.62 |  379.39 | 1.39 s | 5.61 s |
| SteinerSqrt (without divs)   |   8.40 |  75.91 | 451.48 |  2.67 s | 28.63 s |  172 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)

杰里米的简单算法有效，但由于简单的加/减，计算工作量呈指数级快速增长......:)

SteinerSqrt (without divs)

不划分的方法很好，但由于分而治之的变体，结果收敛相对较慢（尤其是在数字很大时）

Nordic Mainframe (CORDIC)

尽管不可变 BigInteger 的逐位操作会产生很大的开销，但 CORDIC 算法已经相当强大了。

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

RedGreenCode (newton method)

经过验证的牛顿法表明，旧的东西不一定要慢。尤其是大数的快速收敛，很难被超越。

RedGreenCode (bound opti.)

Jesan Fafon 节省乘法的提议在这里带来了很多。

MaxSqrt (my version)

首先：首先使用 Math.Sqrt() 计算小数，一旦 double 的精度不再足够，则使用牛顿算法。但是，我尝试使用 Math.Sqrt() 预先计算尽可能多的数字，这使得牛顿算法收敛得更快。

NewtonPlusSqrt (SunsetQuest)

我不会想到如此高的速度增长仍然是可能的（尤其是在数量非常大的情况下）。因此，我用 SunsetQuest 中更快的 NewtonPlusSqrt 替换了我的版本。

注意：不幸的是，许多 .Net 版本不支持 BigInteger.GetBitLength()（例如，2022 年 8 月发布的 .Net 4.8.1 也不支持），因此我添加了一个可选的后备版本，它几乎同样快。

这里是 NewtonPlusSqrt + 可选的 GetBitLengthFallback() 的源代码：

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 xLen = GetBitLengthFallback(x);
  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, let's 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-up ////////
  if ((saveDroppedDigits == 0) && (val * val > x))
  {
    val--;
  }

  return val;
}

public static int GetBitLengthFallback(BigInteger x) // only support x > 0
{
  byte[] bytes = x.ToByteArray();
  byte msb = bytes[bytes.Length - 1];
  int msbBits = 0;
  while (msb != 0)
  {
    msb >>= 1;
    msbBits++;
  }
  return (bytes.Length - 1) * 8 + msbBits;
}

完整基准-来源：

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

namespace BigSqrt
{
  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 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;
    }

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

    public static int GetBitLengthFallback(BigInteger x) // only support x > 0
    {
      byte[] bytes = x.ToByteArray();
      byte msb = bytes[bytes.Length - 1];
      int msbBits = 0;
      while (msb != 0)
      {
        msb >>= 1;
        msbBits++;
      }
      return (bytes.Length - 1) * 8 + msbBits;
    }

    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 xLen = GetBitLengthFallback(x);
      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, let's 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-up ////////
      if ((saveDroppedDigits == 0) && (val * val > x))
      {
        val--;
      }

      return val;
    }

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

      //var q1 = SqrtRedGreenCode(t2);
      //var q2 = SteinerSqrt(t2);
      //var q3 = SqrtJeremyKahan(t2);
      //var q4 = SqrtNordicMainframe(t2);
      //var q5 = SqrtMax(t2);
      //var q6 = NewtonPlusSqrt(t2);
      //if (q6 != q1) throw new Exception();
      //if (q6 != q2) throw new Exception();
      //if (q6 != q3) throw new Exception();
      //if (q6 != q4) throw new Exception();
      //if (q6 != q5) 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 = SteinerSqrt(t2);
          //var q = SqrtJeremyKahan(t2);
          //var q = SqrtNordicMainframe(t2);
          //var q = SqrtMax(t2);
          var q = NewtonPlusSqrt(t2);
        }
        mess.Stop();
        Console.WriteLine((mess.ElapsedTicks * 1000.0 / Stopwatch.Frequency).ToString("N2") + " ms");
      }
    }
  }
}

Answer 4

计算平方根到任意精度的最简单可行方法可能是牛顿方法.