.NET 4.0 提供了 System.Numerics.BigInteger
类型来处理任意大的整数。我需要计算一个 BigInteger
的平方根(或者一个合理的近似值,例如整数平方根)。为了不必重复造轮子,请问是否有一个好的扩展方法可以做到这一点?
.NET 4.0 提供了 System.Numerics.BigInteger
类型来处理任意大的整数。我需要计算一个 BigInteger
的平方根(或者一个合理的近似值,例如整数平方根)。为了不必重复造轮子,请问是否有一个好的扩展方法可以做到这一点?
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中的乘法是热点。
BigInteger upperBound = lowerBound + root + root + 1
,或者内联在返回语句中,如return n >= lowerBound && n <= lowerBound + root + root
。请注意,翻译保持原文意思,尽量使语言通俗易懂,不包含任何解释或其他内容。 - Jesan Fafonstatic 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;
}
好的,首先进行一些速度测试,测试一些在此发布的变量。 (我只考虑那些可以给出精确结果并且适用于 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");
}
}
}
}
简短回答:(但请注意,下面有更多详细信息)
Math.Pow(Math.E, BigInteger.Log(pd) / 2)
其中pd
代表您想要执行平方根操作的BigInteger
。
详细回答和解释:
另一种理解这个问题的方法是理解平方根和对数的工作原理。
如果你有方程5^x = 25
,要解出x
,我们必须使用对数。在这个例子中,我将使用自然对数(其他底数的对数也是可能的,但自然对数是简单的方法)。
5^x = 25
x(ln 5) = ln 25
x = ln 25 / ln 5
x = 2
。但既然我们已经知道了 x(在 5^2 中,x = 2),让我们改变我们不知道的部分并编写一个新方程来解决新的未知数。让 x 成为平方根操作的结果。这给我们带来了:2 = ln 25 / ln x
ln x = (ln 25) / 2
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
。由于 double
是 BigInteger.Log()
的返回类型,最大的 BigInteger
应该是 e 的 double.MaxValue
次幂。在我的电脑上,那将是 e^1.79769313486232E+308
。不精确性没有得到解决。有人想实现 BigDecimal
并更新 BigInteger.Log()
吗?var truncSqrt = function(n) {
var oddNumber = 1;
var result = 0;
while (n >= oddNumber) {
n -= oddNumber;
oddNumber += 2;
result++;
}
return result;
};
更新:为了获得最佳性能,请使用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 发现了一个错误。已经通过守卫条款进行了更正。
***世界上最快的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;
}
1e39457: .Net: 6.37 毫秒 / Go: 261 微秒
和 1e2525223: .Net: 23.1 秒 / Go: 173 毫秒
- undefined1e2525223: Math.Gmp.Native(GMP 6.x的.Net端口):48毫秒
- undefinedSqrt 方法使用 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;
}