尊重与诚信是软件发展的根本因素。
1、多年前 Skyivben 发布的大数计算的代码
本文的代码全部来自于 银河 的文章。
博客园——Skyivben
https://www.cnblogs.com/skyivben/
因为 博客园 的格式问题，大家用起来不方便，因而作者稍微整理一下，并修改了个别的 bug ，发在这里，直接复制保存即可使用。
如果您在发布相关成果的时候，请声明这是 银河 的成果。
再谈 BigInteger - 使用快速傅里叶变换https://www.cnblogs.com/skyivben/archive/2008/07/25/1251697.html
浅谈 BigIntegerhttps://www.cnblogs.com/skyivben/archive/2008/07/13/1241681.html
再谈 BigInteger - 优化https://www.cnblogs.com/skyivben/archive/2010/04/22/1717722.html浅谈 BigIntegerhttps://www.cnblogs.com/skyivben/archive/2008/07/13/1241681.html
BigArithmetic - 提供任意精度的算术运算的静态类https://www.cnblogs.com/skyivben/archive/2008/07/25/1250891.html
2、 Skyiv.Numeric.BigArithmetic 源代码
using System;
using System.Diagnostics;
namespace Skyiv.Numeric
{
/// <summary>
/// 提供任意精度的算术运算的静态类。使用快速傅里叶变换。
/// 本类对字节数组进行算术运算，字节数组以 100 为基。
/// 字节数组中第一个元素存储的数字是最高有效位。
/// https://www.cnblogs.com/skyivben/archive/2008/07/25/1250891.html
/// C语言数值算法程序大全(第二版)，ISBN 7-5053-2931-6 / TP 993
/// Numerical Recipes in C, The Art of Scientific Computing, Second Edition
/// Cambridge University Press 1988, 1992
/// [美] W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery 著
/// 傅祖芸 赵梅娜 丁岩 等译，傅祖芸 校，电子工业出版社，1995年10月第一版
/// </summary>
public static class BigArithmetic
{
// 字节数组的元素包含的十进制数字的个数
static readonly byte Len = 2;
// 字节数组的基
static readonly byte Base = (byte)Math.Pow(10, Len);
// 字节数组的元素的最大值
static readonly byte MaxValue = (byte)(Base - 1);
/// <summary>
/// 复函数的快速傅里叶变换
/// 变量 nn 是复数据点的个数，实型数组 data[1..2*nn] 的实际界长是两倍 nn，
/// 而每个复数值占据了两个相继的存储单元。nn 必须是 2 的整数幂
/// pp.431-432, four1, 12.2 快速傅里叶变换(FFT)
/// </summary>
/// <param name="data">实型数组 data[1..2*nn]。注意，下标从 1 开始</param>
/// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/nn</param>
public static void ComplexFFT(double[] data, bool isInverse)
{
// n 必须是 2 的正整数幂
int n = data.Length - 1;
// 变量 nn 是复数据点的个数
int nn = n >> 1;
// 这个循环实现位序颠倒
for (int i = 1, j = 1; i < n; i += 2)
{
if (j > i)
{
Utility.Swap(ref data[j], ref data[i]);
Utility.Swap(ref data[j + 1], ref data[i + 1]);
}
int m = nn;
for (; m >= 2 && j > m; m >>= 1) j -= m;
j += m;
}
// 执行 log2(nn) 次外循环
for (int mmax = 2, istep = 4; n > mmax; mmax = istep)
{
istep = mmax << 1;
// 下面是关于三角递归的初始赋值
double theta = (isInverse ? -2 : 2) * Math.PI / mmax;
double wtemp = Math.Sin(0.5 * theta);
double wpr = -2 * wtemp * wtemp;
double wpi = Math.Sin(theta);
double wr = 1;
double wi = 0;
for (int m = 1; m < mmax; m += 2)
{
for (int i = m; i <= n; i += istep)
{
int j = i + mmax;
// 下面是 Danielson-Lanczos 公式
double tempr = wr * data[j] - wi * data[j + 1];
double tempi = wr * data[j + 1] + wi * data[j];
data[j] = data[i] - tempr;
data[j + 1] = data[i + 1] - tempi;
data[i] += tempr;
data[i + 1] += tempi;
}
// 三角递归
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
}
}
/// <summary>
/// 单个实函数的快速傅里叶变换
/// 计算一组 n 个实值数据点的傅里叶变换。用复傅里叶变换的正半频率替换这些数据，
/// 它存储在数组 data[1..n] 中。复变换的第一个和最后一个分量的实数值分别返回
/// 单元 data[1] 和 data[2] 中。n 必须是 2 的幂次。这个程序也能计算复数据数组
/// 的逆变换，只要该数组是实值数据的变换(在这种情况下，其结果必须乘以 1/n)即可。
/// pp.436, realft, 12.3.2 单个实函数的 FFT
/// </summary>
/// <param name="data">实型数组 data[1..n]。注意，下标从 1 开始</param>
/// <param name="isInverse">是否逆变换。注意: 逆变换未乘上归一化因子 1/n</param>
public static void RealFFT(double[] data, bool isInverse)
{
// n 必须是 2 的整数幂
int n = data.Length - 1;
// 此处是正向变换
if (!isInverse) ComplexFFT(data, isInverse);
// 递归的初始赋值
double theta = (isInverse ? -2 : 2) * Math.PI / n;
double wtemp = Math.Sin(0.5 * theta);
double wpr = -2 * wtemp * wtemp;
double wpi = Math.Sin(theta);
double wr = 1 + wpr;
double wi = wpi;
double c1 = 0.5;
double c2 = isInverse ? 0.5 : -0.5;
int n3 = n + 3;
int n4 = n >> 2;
for (int i = 2; i <= n4; i++)
{
int i1 = i + i - 1, i2 = i1 + 1, i3 = n3 - i2, i4 = i3 + 1;
// 两个分离变换是从 data 中分离出来
double h1r = c1 * (data[i1] + data[i3]);
double h1i = c1 * (data[i2] - data[i4]);
double h2r = -c2 * (data[i2] + data[i4]);
double h2i = c2 * (data[i1] - data[i3]);
// 此处重新组合以形成原始实型数据的真实变换
data[i1] = h1r + wr * h2r - wi * h2i;
data[i2] = h1i + wr * h2i + wi * h2r;
data[i3] = h1r - wr * h2r + wi * h2i;
data[i4] = -h1i + wr * h2i + wi * h2r;
// 递归式
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
double tmp = data[1];
if (!isInverse)
{
// 同时挤压第一个和最后一个数据使它们都在原始数组中
data[1] = tmp + data[2];
data[2] = tmp - data[2];
}
else
{
data[1] = c1 * (tmp + data[2]);
data[2] = c1 * (tmp - data[2]);
// 此处是逆变换
ComplexFFT(data, isInverse);
}
}
/// <summary>
/// 比较 x[0..n-1] 和 y[0..n-1]
/// </summary>
/// <param name="x">第一操作数 x[0..n-1]</param>
/// <param name="y">第二操作数 y[0..n-1]</param>
/// <param name="n">两个操作数 x 和 y 的精度</param>
/// <returns>比较结果：-1:小于 1:大于 0:等于</returns>
public static int Compare(byte[] x, byte[] y, int n)
{
Debug.Assert(x.Length >= n && y.Length >= n);
for (int i = 0; i < n; i++)
{
if (x[i] != y[i])
{
return (x[i] < y[i]) ? -1 : 1;
}
}
return 0;
}
/// <summary>
/// 求补码。注意，被操作数被修改。
/// pp.775, mpneg, 20.6 任意精度的运算
/// </summary>
/// <param name="data">被操作数 data[0..n-1]</param>
/// <param name="n">被操作数 data 的精度</param>
/// <returns>被操作数的补码 data[0..n-1]</returns>
public static byte[] Negative(byte[] data, int n)
{
Debug.Assert(data.Length >= n);
for (int k = Base, i = n - 1; i >= 0; i--)
{
data[i] = (byte)((k = MaxValue + k / Base - data[i]) % Base);
}
return data;
}
/// <summary>
/// 减法。从 minuend[0..n-1] 中减去 subtrahend[0..n-1]，得到 difference[0..n-1]
/// pp.774, mpsub, 20.6 任意精度的运算
/// </summary>
/// <param name="difference">差 difference[0..n-1]</param>
/// <param name="minuend">被减数 minuend[0..n-1]</param>
/// <param name="subtrahend">减数 subtrahend[0..n-1]</param>
/// <param name="n">被减数 minuend 和减数 subtrahend 的精度</param>
/// <returns>差 difference[0..n-1]</returns>
public static byte[] Subtract(byte[] difference, byte[] minuend, byte[] subtrahend, int n)
{
Debug.Assert(minuend.Length >= n && subtrahend.Length >= n && difference.Length >= n);
for (int k = Base, i = n - 1; i >= 0; i--)
{
difference[i] = (byte)((k = MaxValue + k / Base + minuend[i] - subtrahend[i]) % Base);
}
return difference;
}
/// <summary>
/// 加法。augend[0..n-1] 与 addend[0..n-1] 相加，得到 sum[0..n]
/// pp.774, mpadd, 20.6 任意精度的运算
/// </summary>
/// <param name="sum">和 sum[0..n]</param>
/// <param name="augend">被加数 augend[0..n-1]</param>
/// <param name="addend">加数 addend[0..n-1]</param>
/// <param name="n">被加数 augend 和加数 addend 的精度</param>
/// <returns>和 sum[0..n]</returns>
public static byte[] Add(byte[] sum, byte[] augend, byte[] addend, int n)
{
Debug.Assert(augend.Length >= n && addend.Length >= n && sum.Length >= n + 1);
int k = 0;
for (int i = n - 1; i >= 0; i--)
{
sum[i + 1] = (byte)((k = k / Base + augend[i] + addend[i]) % Base);
}
sum[0] += (byte)(k / Base);
return sum;
}
/// <summary>
/// 捷加法。augend[0..n-1] 与整数 addend 相加，得到 sum[0..n]
/// pp.774, mpadd, 20.6 任意精度的运算
/// </summary>
/// <param name="sum">和 sum[0..n]</param>
/// <param name="augend">被加数 augend[0..n-1]</param>
/// <param name="n">被加数 augend 的精度</param>
/// <param name="addend">加数 addend</param>
/// <returns>和 sum[0..n]</returns>
public static byte[] Add(byte[] sum, byte[] augend, int n, byte addend)
{
Debug.Assert(augend.Length >= n && sum.Length >= n + 1);
int k = Base * addend;
for (int i = n - 1; i >= 0; i--)
{
sum[i + 1] = (byte)((k = k / Base + augend[i]) % Base);
}
sum[0] += (byte)(k / Base);
return sum;
}
/// <summary>
/// 捷除法。dividend[0..n-1] 除以整数 divisor，得到 quotient[0..n-1]
/// pp.775, mpsdv, 20.6 任意精度的运算
/// </summary>
/// <param name="quotient">商 quotient[0..n-1]</param>
/// <param name="dividend">被除数 dividend[0..n-1]</param>
/// <param name="n">被除数 dividend 的精度</param>
/// <param name="divisor">除数 divisor</param>
/// <returns>商 quotient[0..n-1]</returns>
public static byte[] Divide(byte[] quotient, byte[] dividend, int n, byte divisor)
{
Debug.Assert(quotient.Length >= n && dividend.Length >= n);
for (int r = 0, k = 0, i = 0; i < n; i++, r = k % divisor)
{
quotient[i] = (byte)((k = Base * r + dividend[i]) / divisor);
}
return quotient;
}
/// <summary>
/// 乘法。multiplicand[0..n-1] 与 multiplier[0..m-1] 相乘，得到 product[0..n+m-1]
/// pp.776-777, mpmul, 20.6 任意精度的运算
/// </summary>
/// <param name="product">积 product[0..n+m-1]</param>
/// <param name="multiplicand">被乘数 multiplicand[0..n-1]</param>
/// <param name="n">被乘数 multiplicand 的精度</param>
/// <param name="multiplier">乘数 multiplier[0..m-1]</param>
/// <param name="m">乘数 multiplier 的精度</param>
/// <returns>积 product[0..n+m-1]</returns>
public static byte[] Multiply(byte[] product, byte[] multiplicand, int n, byte[] multiplier, int m)
{
int mn = m + n, nn = 1;
Debug.Assert(product.Length >= mn && multiplicand.Length >= n && multiplier.Length >= m);
// 为变换找出最小可用的 2 的幂次
while (nn < mn) nn <<= 1;
double[] a = new double[nn + 1], b = new double[nn + 1];
for (int i = 0; i < n; i++) a[i + 1] = multiplicand[i];
for (int i = 0; i < m; i++) b[i + 1] = multiplier[i];
// 执行卷积，首先求出二个傅里叶变换
RealFFT(a, false);
RealFFT(b, false);
// 复数相乘的结果(实部和虚部)
b[1] *= a[1];
b[2] *= a[2];
for (int i = 3; i <= nn; i += 2)
{
double t = b[i];
b[i] = t * a[i] - b[i + 1] * a[i + 1];
b[i + 1] = t * a[i + 1] + b[i + 1] * a[i];
}
// 进行傅里叶逆变换
RealFFT(b, true);
byte[] bs = new byte[nn + 1];
// 执行最后完成所有进位的过程
long cy = 0;
for (int i = nn, n2 = nn / 2; i >= 1; i--)
{
long t = (long)(b[i] / n2 + cy + 0.5);
// 原书中这句使用循环，有严重的性能问题
bs[i] = (byte)(t % Base);
cy = t / Base;
}
if (cy >= Base) throw new OverflowException("FFT Multiply");
bs[0] = (byte)cy;
Array.Copy(bs, product, n + m);
return product;
}
/// <summary>
/// 除法。dividend[0..n-1] 除以 divisor[0..m-1]，m ≤ n，
/// 得到：商 quotient[0..n-m]，余数 remainder[0..m-1]。
/// pp.778, mpdiv, 20.6 任意精度的运算
/// </summary>
/// <param name="quotient">商 quotient[0..n-m]</param>
/// <param name="remainder">余数 remainder[0..m-1]</param>
/// <param name="dividend">被除数 dividend[0..n-1]</param>
/// <param name="n">被除数 dividend 的精度</param>
/// <param name="divisor">除数 divisor[0..m-1]</param>
/// <param name="m">除数 divisor 的精度</param>
/// <returns>商 quotient[0..n-m]</returns>
public static byte[] DivRem(byte[] quotient, byte[] remainder, byte[] dividend, int n, byte[] divisor, int m)
{
Debug.Assert(m <= n && dividend.Length >= n && divisor.Length >= m && quotient.Length >= n - m + 1 && remainder.Length >= m);
int MACC = 3;
byte[] s = new byte[n + MACC], t = new byte[n - m + MACC + n];
// s = 1 / divisor
Inverse(s, n - m + MACC, divisor, m);
// quotient = dividend / divisor
Array.Copy(Multiply(t, s, n - m + MACC, dividend, n), 1, quotient, 0, n - m + 1);
// s = quotient * divisor
Array.Copy(Multiply(t, quotient, n - m + 1, divisor, m), 1, s, 0, n);
// s = dividend - quotient * divisor
Subtract(s, dividend, s, n);
Array.Copy(s, n - m, remainder, 0, m);
// 调整商和余数
if (Compare(remainder, divisor, m) >= 0)
{
Subtract(remainder, remainder, divisor, m);
Add(s, quotient, n - m + 1, 1);
Array.Copy(s, 1, quotient, 0, n - m + 1);
}
return quotient;
}
/// <summary>
/// 求倒数。
/// pp.777 - 778, mpinv, 20.6 任意精度的运算
/// </summary>
/// <param name="inverse">倒数 inverse[0..n-1]，在 inverse[0] 后有基数的小数点</param>
/// <param name="n">倒数 inverse 的精度</param>
/// <param name="data">被操作数 data[0..m-1]，data[0] > 0，在 data[0] 后有基数的小数点</param>
/// <param name="m">被操作数 data 的精度</param>
/// <returns>倒数 inverse[0..n-1]，在 inverse[0] 后有基数的小数点</returns>
public static byte[] Inverse(byte[] inverse, int n, byte[] data, int m)
{
Debug.Assert(inverse.Length >= n && data.Length >= m);
InitialValue(inverse, n, data, m, false);
if (n == 1) return inverse;
byte[] s = new byte[n], t = new byte[n + n];
// 牛顿迭代法: inverse = inverse * ( 2 - data * inverse ) => inverse = 1 / data
for (; ; )
{
// s = data * inverse
Array.Copy(Multiply(t, inverse, n, data, m), 1, s, 0, n);
// s = -(data * inverse)
Negative(s, n);
// s = 2 - data * inverse
s[0] -= (byte)(Base - 2);
// inverse = inverse * s
Array.Copy(Multiply(t, s, n, inverse, n), 1, inverse, 0, n);
int i = 1;
// 判断 s 的小数部分是否为零
for (; i < n - 1 && s[i] == 0; i++) ;
// 若 s 收敛到 1 则返回 inverse = 1 / data
if (i == n - 1) return inverse;
}
}
/// <summary>
/// 求平方根 sqrt，以及平方根的倒数 invSqrt。invSqrt 也可设为 sqrt，此时，invSqrt 也是平方根。
/// pp.778-779, mpsqrt, 20.6 任意精度的运算
/// </summary>
/// <param name="sqrt">平方根 sqrt[0..n-1]，在 sqrt[0] 后有基数的小数点</param>
/// <param name="invSqrt">平方根的倒数 invSqrt[0..n-1]，在 invSqrt[0] 后有基数的小数点</param>
/// <param name="n">平方根的精度</param>
/// <param name="data">被操作数 data[0..m-1]，data[0] > 0，在 data[0] 后有基数的小数点</param>
/// <param name="m">被操作数 data 的精度</param>
/// <returns>平方根 sqrt[0..n-1]，在 sqrt[0] 后有基数的小数点</returns>
public static byte[] Sqrt(byte[] sqrt, byte[] invSqrt, int n, byte[] data, int m)
{
Debug.Assert(sqrt.Length >= n && invSqrt.Length >= n && data.Length >= m);
if (n <= 1) throw new ArgumentOutOfRangeException("n", "must greater than 1");
InitialValue(invSqrt, n, data, m, true);
byte[] s = new byte[n], t = new byte[n + Math.Max(m, n)];
// invSqrt = invSqrt * (3 - data * invSqrt * invSqrt) / 2 => invSqrt = 1 / sqrt(data)
for (; ; )
{
// s = invSqrt * invSqrt
Array.Copy(Multiply(t, invSqrt, n, invSqrt, n), 1, s, 0, n);
// s = data * invSqrt * invSqrt
Array.Copy(Multiply(t, s, n, data, m), 1, s, 0, n);
// s = -(data * invSqrt * invSqrt)
Negative(s, n);
// s = 3 - data * invSqrt * invSqrt
s[0] -= (byte)(Base - 3);
// s = (3 - data * invSqrt * invSqrt) / 2
Divide(s, s, n, 2);
// invSqrt = invSqrt * s
Array.Copy(Multiply(t, s, n, invSqrt, n), 1, invSqrt, 0, n);
int i = 1;
// 判断 s 的小数部分是否为零
for (; i < n - 1 && s[i] == 0; i++) ;
// 若 s 没有收敛到 1 则继续迭代
if (i < n - 1) continue;
// sqrt = invSqrt * data = sqrt(data)
Array.Copy(Multiply(t, invSqrt, n, data, m), 1, sqrt, 0, n);
return sqrt;
}
}
/// <summary>
/// 采用浮点运算以得到一个初始近似值 u[0..n-1]: u = 1 / data 或者 u = 1 / sqrt(data)
/// </summary>
/// <param name="u">初始近似值 u[0..n-1]</param>
/// <param name="n">所需的精度</param>
/// <param name="data">被操作数 data[0..m-1]</param>
/// <param name="m">被操作数 data 的精度</param>
/// <param name="isSqrt">是否求平方根</param>
/// <returns>初始近似值 u[0..n-1]</returns>
static byte[] InitialValue(byte[] u, int n, byte[] data, int m, bool isSqrt)
{
Debug.Assert(u.Length >= n && data.Length >= m);
// double 可达到 16 位有效数字
int scale = 16 / Len;
double fu = 0;
for (int i = Math.Min(scale, m) - 1; i >= 0; i--) fu = fu / Base + data[i];
fu = 1 / (isSqrt ? Math.Sqrt(fu) : fu);
for (int i = 0; i < Math.Min(scale + 1, n); i++)
{
int k = (int)fu;
u[i] = (byte)k;
fu = Base * (fu - k);
}
return u;
}
}
}
2、Skyiv.Numeric.BigInteger源代码
using System;
using System.Text;
namespace Skyiv.Numeric
{
/// <summary>
/// BigInterger类
/// https://www.cnblogs.com/skyivben/archive/2008/07/25/1251697.html
/// </summary>
public class BigInteger : IEquatable<BigInteger>, IComparable<BigInteger>
{
static readonly byte Len = 2;
static readonly byte Base = (byte)Math.Pow(10, Len);
// 符号，取值：-1, 0, 1。
sbyte sign;
// 字节数组以 100 为基，字节数组中第一个元素存储的数字是最高有效位。
byte[] data;
BigInteger()
{
}
BigInteger(long x)
{
sign = (sbyte)((x == 0) ? 0 : ((x > 0) ? 1 : -1));
data = new byte[10]; // long.MinValue = -9,223,372,036,854,775,808
ulong z = (x < 0) ? (ulong)-x : (ulong)x;
for (int i = data.Length - 1; z != 0; i--, z /= Base) data[i] = (byte)(z % Base);
Shrink();
}
BigInteger(BigInteger x)
{
sign = x.sign; // x != null
data = new byte[x.data.Length];
Array.Copy(x.data, data, data.Length);
}
public static implicit operator BigInteger(long x)
{
return new BigInteger(x);
}
public static BigInteger Parse(string s)
{
if (s == null) return null;
s = s.Trim().Replace(",", null);
if (s.Length == 0) return 0;
BigInteger z = new BigInteger();
z.sign = (sbyte)((s[0] == '-') ? -1 : 1);
if (s[0] == '-' || s[0] == '+') s = s.Substring(1);
int r = s.Length % Len;
z.data = new byte[s.Length / Len + ((r != 0) ? 1 : 0)];
int i = 0;
if (r != 0) z.data[i++] = byte.Parse(s.Substring(0, r));
for (; i < z.data.Length; i++, r += Len) z.data[i] = byte.Parse(s.Substring(r, Len));
z.Shrink();
return z;
}
public static void Swap(ref BigInteger x, ref BigInteger y)
{
BigInteger z = x;
x = y;
y = z;
}
public static ulong Abs(long x)
{
return (x < 0) ? (ulong)-x : (ulong)x;
}
public static BigInteger Abs(BigInteger x)
{
if (x == null) return null;
BigInteger z = new BigInteger(x);
z.sign = Math.Abs(x.sign);
return z;
}
public static BigInteger Pow(BigInteger x, int y)
{
if (x == null) return null;
BigInteger z = 1, n = x;
for (; y > 0; y >>= 1, n *= n) if ((y & 1) != 0) z *= n;
return z;
}
public static BigInteger operator +(BigInteger x)
{
if (x == null) return null;
return new BigInteger(x);
}
public static BigInteger operator -(BigInteger x)
{
if (x == null) return null;
BigInteger z = new BigInteger(x);
z.sign = (sbyte)-x.sign;
return z;
}
public static BigInteger operator ++(BigInteger x)
{
return x + 1;
}
public static BigInteger operator --(BigInteger x)
{
return x - 1;
}
public static BigInteger operator +(BigInteger x, BigInteger y)
{
if (x == null || y == null) return null;
if (x.AbsCompareTo(y) < 0) Utility.Swap(ref x, ref y);
BigInteger z = new BigInteger();
z.sign = x.sign;
byte[] bs = Utility.Expand(y.data, x.data.Length);
bool isAdd = x.sign * y.sign == 1;
z.data = new byte[x.data.Length + (isAdd ? 1 : 0)];
if (isAdd) BigArithmetic.Add(z.data, x.data, bs, bs.Length);
else BigArithmetic.Subtract(z.data, x.data, bs, bs.Length);
z.Shrink();
return z;
}
public static BigInteger operator -(BigInteger x, BigInteger y)
{
if (x == null || y == null) return null;
return x + (-y);
}
public static BigInteger operator *(BigInteger x, BigInteger y)
{
if (x == null || y == null) return null;
if (x.sign * y.sign == 0) return 0;
BigInteger z = new BigInteger();
z.sign = (sbyte)(x.sign * y.sign);
z.data = new byte[x.data.Length + y.data.Length];
BigArithmetic.Multiply(z.data, x.data, x.data.Length, y.data, y.data.Length);
z.Shrink();
return z;
}
public static BigInteger operator /(BigInteger dividend, BigInteger divisor)
{
BigInteger remainder;
return DivRem(dividend, divisor, out remainder);
}
public static BigInteger operator %(BigInteger dividend, BigInteger divisor)
{
BigInteger remainder;
DivRem(dividend, divisor, out remainder);
return remainder;
}
public static BigInteger DivRem(BigInteger dividend, BigInteger divisor, out BigInteger remainder)
{
remainder = null;
if (dividend == null || divisor == null) return null;
if (divisor.sign == 0) throw new DivideByZeroException();
if (dividend.AbsCompareTo(divisor) < 0)
{
remainder = new BigInteger(dividend);
return 0;
}
BigInteger quotient = new BigInteger();
remainder = new BigInteger();
quotient.data = new byte[dividend.data.Length - divisor.data.Length + 1];
remainder.data = new byte[divisor.data.Length];
BigArithmetic.DivRem(quotient.data, remainder.data, dividend.data,
dividend.data.Length, divisor.data, divisor.data.Length);
quotient.sign = (sbyte)(dividend.sign * divisor.sign);
remainder.sign = dividend.sign;
quotient.Shrink();
remainder.Shrink();
return quotient;
}
public static BigInteger Sqrt(BigInteger x)
{
if (x == null || x.sign < 0) return null;
if (x.sign == 0) return 0;
if (x.data.Length == 1) return new BigInteger((long)Math.Sqrt(x.data[0]));
BigInteger z = new BigInteger();
z.sign = 1;
z.data = new byte[x.data.Length / 2 + 3];
z.data = Adjust(BigArithmetic.Sqrt(z.data, z.data, z.data.Length, x.data, x.data.Length), x.data.Length);
z.Shrink();
// 平方根有可能比实际小 1。
BigInteger z1 = z + 1;
return (z1 * z1 <= x) ? z1 : z;
}
static byte[] Adjust(byte[] bs, int digits)
{
if (bs[0] >= 10) throw new OverflowException("sqrt adjust");
byte[] nbs = new byte[(digits + 1) / 2];
if (digits % 2 == 0)
for (int k = bs[0], i = 0; i < nbs.Length; i++, k = bs[i] % 10)
nbs[i] = (byte)(k * 10 + bs[i + 1] / 10);
else Array.Copy(bs, nbs, nbs.Length);
return nbs;
}
void Shrink()
{
int i;
for (i = 0; i < data.Length; i++) if (data[i] != 0) break;
if (i != 0)
{
byte[] bs = new byte[data.Length - i];
Array.Copy(data, i, bs, 0, bs.Length);
data = bs;
}
if (data.Length == 0) sign = 0;
}
public static bool operator ==(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return object.ReferenceEquals(y, null);
return x.Equals(y);
}
public static bool operator !=(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return !object.ReferenceEquals(y, null);
return !x.Equals(y);
}
public static bool operator <(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return !object.ReferenceEquals(y, null);
return x.CompareTo(y) < 0;
}
public static bool operator >(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return false;
return x.CompareTo(y) > 0;
}
public static bool operator <=(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return true;
return x.CompareTo(y) <= 0;
}
public static bool operator >=(BigInteger x, BigInteger y)
{
if (object.ReferenceEquals(x, null)) return object.ReferenceEquals(y, null);
return x.CompareTo(y) >= 0;
}
public bool Equals(BigInteger other)
{
return CompareTo(other) == 0;
}
public override bool Equals(object other)
{
if (other == null || GetType() != other.GetType()) return false;
return Equals(other as BigInteger);
}
public int CompareTo(BigInteger other)
{
if (object.ReferenceEquals(other, null)) return 1;
if (sign < other.sign) return -1;
if (sign > other.sign) return 1;
if (sign == 0) return 0;
return sign * AbsCompareTo(other);
}
public override string ToString()
{
StringBuilder sb = new StringBuilder();
if (sign < 0) sb.Append('-');
sb.Append((data.Length == 0) ? 0 : (int)data[0]);
for (int i = 1; i < data.Length; i++) sb.Append(data[i].ToString("D" + Len));
return sb.ToString();
}
public override int GetHashCode()
{
int hash = sign;
foreach (int n in data) hash ^= n;
return hash;
}
int AbsCompareTo(BigInteger other)
{
if (data.Length < other.data.Length) return -1;
if (data.Length > other.data.Length) return 1;
return BigArithmetic.Compare(data, other.data, data.Length);
}
}
public static class Utility
{
public static T[] Expand<T>(T[] x, int n)
{
T[] z = new T[n]; // assume n >= x.Length
Array.Copy(x, 0, z, n - x.Length, x.Length);
return z;
}
public static void Swap<T>(ref T x, ref T y)
{
T z = x;
x = y;
y = z;
}
}
}
![【LeetCode】从前序与中序遍历序列构造二叉树 [M]（二叉树重构）](https://img-blog.csdnimg.cn/img_convert/2bebde118125eb0db2c884ee7a64608d.jpeg)
