尊重与诚信是软件发展的根本因素。

1、多年前 Skyivben 发布的大数计算的代码

本文的代码全部来自于 银河 的文章。

博客园——Skyivbenhttps://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;
        }
    }
}