Calculate square root of a BigInteger (System.Numerics.BigInteger)
.NET 4.0 provides the System.Numerics.BigInteger type for arbitrarily-large integers. I need to compute the square root (or a reasonable approximation -- e.g., integer square root) of a BigInteger. So that I don't have to reimplement the wheel, does anyone have a nice extension method for this?
Check if BigInteger is not a perfect square has code to compute the integer square root of a Java BigInteger. Here it is translated into C#, as an extension method.
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);
}
Informal testing indicates that this is about 75X slower than Math.Sqrt, for small integers. The VS profiler points to the multiplications in isSqrt as the hotspots.
I am not sure if Newton's Method is the best way to compute bignum square roots, because it involves divisions which are slow for bignums. You can use a CORDIC method, which uses only addition and shifts (shown here for unsigned ints)
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;
}
There's a similar method which uses only addition and shifts, called 'Dijkstras Square Root', explained for example here:
- http://lib.tkk.fi/Diss/2005/isbn9512275279/article3.pdf
Ok, first a few speed tests of some variants posted here. (I only considered methods which give exact results and are at least suitable for 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)
Descriptions:
Jeremy Kahan (js-port)
Jeremy's simple algorithm works, but the computational effort increases exponentially very fast due to the simple adding/subtracting... :)
Sunsetquest (without divs)
The approach without dividing is good, but due to the divide and conquer variant the results converges relatively slowly (especially with large numbers)
Nordic Mainframe (CORDIC)
The CORDIC algorithm is already quite powerful, although the bit-by-bit operation of the imuttable BigIntegers generates much overhead.
I have calculated the required bits this way: int b = Convert.ToInt32(Math.Ceiling(BigInteger.Log(x, 2))) / 2 + 1;
RedGreenCode (newton method)
The proven newton method shows that something old does not have to be slow. Especially the fast convergence of large numbers can hardly be topped.
RedGreenCode (bound opti.)
The proposal of Jesan Fafon to save a multiplication has brought a lot here.
my version
First: calculate small numbers at the beginning with Math.Sqrt() and as soon as the accuracy of double is no longer sufficient, then use the newton algorithm. However, I try to pre-calculate as many numbers as possible with Math.Sqrt(), which makes the newton algorithm converge much faster.
Here the source:
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);
}
full Benchmark-Source:
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");
}
}
}
}