BigInteger to/from (uint[], bool)
用(uint[], bool)来表示有符号大数，其中uint[]是大数的绝对值，bool为false时是负数。
/// <summary>
/// (<see cref="uint"/>[], <see cref="bool"/>) to <see cref="BigInteger"/>
/// </summary>
private BigInteger ValueOf((uint[], bool) value)
{
var result = BigInteger.Zero;
var array = value.Item1.SkipWhile(num => num == 0u).ToArray();
foreach (var num in array)
{
result <<= 32;
result |= (num & 0xFFFF_FFFF);
}
return value.Item2 ? result : -result;
}
/// <summary>
/// <see cref="BigInteger"/> to (<see cref="uint"/>[], <see cref="bool"/>)
/// </summary>
private (uint[], bool) ToTuple(BigInteger value)
{
var positive = BigInteger.Abs(value);
var byteCount = positive.GetByteCount();
var len = (int)Math.Ceiling(byteCount / 4d);
var result = new uint[len];
for (var i = len - 1; i >= 0; --i)
{
result[i] = (uint)(positive & 0xFFFF_FFFF);
positive >>= 32;
}
return (result, value >= 0);
}Unsigned Knuth Division
/// <summary>
/// <paramref name="divident"/> / <paramref name="divisor"/> and
/// <paramref name="divident"/> % <paramref name="divisor"/> using Knuth division.
/// all integers must be positive. first is quotient, last is remainder.
/// </summary>
public static (uint[], bool)[] DivideAndRemainderKnuth((uint[], bool) dividend, (uint[], bool) divisor)
{
if (IsZero(divisor))
throw new ArithmeticException("BigInteger divide by zero");
// Dividend is zero
if (IsZero(dividend))
return new (uint[], bool)[] { ZERO, ZERO };
var cmp = Compare(dividend, divisor);
// Dividend less than divisor
if (cmp < 0)
return new (uint[], bool)[] { ZERO, dividend };
// Dividend equal to divisor
if (cmp == 0)
return new (uint[], bool)[] { ONE, ZERO };
// Special case one word divisor
if (divisor.Item1.Length == 1)
{
var r = DivideOneWord(dividend.Item1, divisor.Item1[0], out uint[] quotient);
return new (uint[], bool)[] { (quotient, true), (new uint[] { r }, true) };
}
// Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
if (dividend.Item1.Length >= KNUTH_POW2_THRESH_LEN)
{
var trailingZeroBits = Math.Min(GetLowestSetBit(dividend), GetLowestSetBit(divisor));
if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS * 32)
{
var aa = ShiftRight(dividend, trailingZeroBits);
var bb = ShiftRight(divisor, trailingZeroBits);
var r = DivideAndRemainderKnuth(aa, bb);
r[1] = ShiftLeft(r[1], trailingZeroBits);
return r;
}
}
return DivideAndRemainderKnuthAux(dividend, divisor);
}
/// <summary>
/// <paramref name="divident"/> / <paramref name="divisor"/> and
/// <paramref name="divident"/> % <paramref name="divisor"/> using Knuth division.
/// all integers must be positive. first tuple is quotient, last is remainder.
/// </summary>
private static (uint[], bool)[] DivideAndRemainderKnuthAux((uint[], bool) dividend, (uint[], bool) div)
{
// assert div.intLen > 1
// D1 normalize the divisor
var shift = div.Item1[0].NumberOfLeadingZeros();
// Copy divisor value to protect divisor
var dlen = div.Item1.Length;
uint[] divisor;
(uint[] value, int offset, int intLen) rem; // Remainder starts as dividend with space for a leading zero
if (shift > 0)
{
divisor = new uint[dlen];
CopyAndShift(div.Item1, 0, dlen, divisor, 0, shift);
if (dividend.Item1[0].NumberOfLeadingZeros() >= shift)
{
var remarr = new uint[dividend.Item1.Length + 1];
CopyAndShift(dividend.Item1, 0, dividend.Item1.Length, remarr, 1, shift);
rem.value = remarr;
rem.intLen = dividend.Item1.Length;
rem.offset = 1;
}
else
{
var remarr = new uint[dividend.Item1.Length + 2];
var c = 0u;
var n2 = 32 - shift;
for (int i = 1, rFrom = 0; i < dividend.Item1.Length + 1; i++, rFrom++)
{
var b = c;
c = dividend.Item1[rFrom];
remarr[i] = (b << shift) | (c >> n2);
}
remarr[dividend.Item1.Length + 1] = c << shift;
rem.value = remarr;
rem.intLen = dividend.Item1.Length + 1;
rem.offset = 1;
}
}
else
{
divisor = div.Item1;
rem.value = new uint[dividend.Item1.Length + 1];
Array.Copy(dividend.Item1, 0, rem.value, 1, dividend.Item1.Length);
rem.intLen = dividend.Item1.Length;
rem.offset = 1;
}
var nlen = rem.intLen;
// Set the quotient size
var limit = nlen - dlen + 1;
//(uint[] value, int intLen, int offset) quotient;
var quotient = new uint[limit];
// Must insert leading 0 in rem if its length did not change
rem.offset = 0;
rem.intLen++;
var dh = divisor[0];
var dhLong = dh & LONG_MASK;
var dl = divisor[1];
uint qhat, qrem, nh, nh2, nm; bool skipCorrection;
// D2 Initialize j
for (var j = 0; j < limit - 1; j++)
{
// D3 Calculate qhat
// estimate qhat
qhat = 0;
qrem = 0;
skipCorrection = false;
nh = rem.value[j + rem.offset];
nh2 = nh + 0x80000000;
nm = rem.value[j + 1 + rem.offset];
if (nh == dh)
{
qhat = uint.MaxValue;
qrem = nh + nm;
skipCorrection = (qrem + 0x80000000) < nh2;
}
else
{
var nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0)
{
qhat = (uint)(nChunk / dhLong);
qrem = (uint)(nChunk - (qhat * dhLong));
}
else
{
var tmp = DivWord(nChunk, dh);
qhat = (uint)(tmp & LONG_MASK);
qrem = (uint)((ulong)tmp >> 32);
}
}
if (qhat == 0)
continue;
if (!skipCorrection)
{ // Correct qhat
var nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
var rs = ((qrem & LONG_MASK) << 32) | nl;
var estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (UnsignedLongCompare(estProduct, rs))
{
qhat--;
qrem = (uint)((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong)
{
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (UnsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
// D4 Multiply and subtract
rem.value[j + rem.offset] = 0;
var borrow = MulSub(rem.value, divisor, qhat, dlen, j + rem.offset);
// D5 Test remainder
if (borrow + 0x80000000 > nh2)
{
// D6 Add back
DivAdd(divisor, rem.value, j + 1 + rem.offset);
qhat--;
}
// Store the quotient digit
quotient[j] = qhat;
} // D7 loop on j
// D3 Calculate qhat
// estimate qhat
qhat = 0;
qrem = 0;
skipCorrection = false;
nh = rem.value[limit - 1 + rem.offset];
nh2 = nh + 0x80000000;
nm = rem.value[limit + rem.offset];
if (nh == dh)
{
qhat = uint.MaxValue;
qrem = nh + nm;
skipCorrection = qrem + 0x80000000 < nh2;
}
else
{
var nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0)
{
qhat = (uint)(nChunk / dhLong);
qrem = (uint)(nChunk - (qhat * dhLong));
}
else
{
var tmp = DivWord(nChunk, dh);
qhat = (uint)(tmp & LONG_MASK);
qrem = (uint)((ulong)tmp >> 32);
}
}
if (qhat != 0)
{
if (!skipCorrection)
{ // Correct qhat
var nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
var rs = ((qrem & LONG_MASK) << 32) | nl;
var estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (UnsignedLongCompare(estProduct, rs))
{
qhat--;
qrem = (uint)((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong)
{
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (UnsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
// D4 Multiply and subtract
rem.value[limit - 1 + rem.offset] = 0;
var borrow = MulSub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
// D5 Test remainder
if (borrow + 0x80000000 > nh2)
{
// D6 Add back
DivAdd(divisor, rem.value, limit + rem.offset);
qhat--;
}
// Store the quotient digit
quotient[limit - 1] = qhat;
}
var remainder = rem.value;
// D8 Unnormalize
if (shift > 0)
remainder = ShiftRight((rem.value, true), shift).Item1;
remainder = StripLeadingZeroInts(remainder);
quotient = StripLeadingZeroInts(quotient);
return new (uint[], bool)[] { (quotient, true), (remainder, true) };
}
/// <summary>
/// This method is used for division of an n word dividend by a one word
/// divisor.The quotient is placed into quotient.The one word divisor is
/// specified by divisor.
/// </summary>
private static uint DivideOneWord(uint[] dividend, uint divisor, out uint[] quotient)
{
var divisorLong = divisor & LONG_MASK;
quotient = new uint[dividend.Length];
// Special case of one word dividend
if (dividend.Length == 1)
{
var dividendValue = dividend[0] & LONG_MASK;
var q = (uint)(dividendValue / divisorLong);
var r = (uint)(dividendValue - q * divisorLong);
quotient[0] = q;
return r;
}
// Normalize the divisor
var shift = divisor.NumberOfLeadingZeros();
var rem = dividend[0];
var remLong = rem & LONG_MASK;
if (remLong < divisorLong)
{
quotient[0] = 0;
}
else
{
quotient[0] = (uint)(remLong / divisorLong);
rem = (uint)(remLong - (quotient[0] * divisorLong));
remLong = rem & LONG_MASK;
}
var xlen = dividend.Length;
while (--xlen > 0)
{
var dividendEstimate = (remLong << 32) |
(dividend[dividend.Length - xlen] & LONG_MASK);
uint q;
if (dividendEstimate >= 0)
{
q = (uint)(dividendEstimate / divisorLong);
rem = (uint)(dividendEstimate - q * divisorLong);
}
else
{
var tmp = DivWord(dividendEstimate, divisor);
q = (uint)(tmp & LONG_MASK);
rem = (uint)((ulong)tmp >> 32);
}
quotient[dividend.Length - xlen] = q;
remLong = rem & LONG_MASK;
}
quotient = StripLeadingZeroInts(quotient);
// Unnormalize
if (shift > 0)
return rem % divisor;
else
return rem;
}
/// <summary>
/// This method divides a long quantity by an int to estimate
/// qhat for two multi precision numbers.It is used when
/// the signed value of n is less than zero.
/// Returns long value where high 32 bits contain remainder value and
/// low 32 bits contain quotient value.
/// </summary>
private static long DivWord(long dividend, uint divisor)
{
var dLong = divisor & LONG_MASK;
long r, q;
if (dLong == 1)
{
q = (uint)dividend;
r = 0;
return (r << 32) | (q & LONG_MASK);
}
// Approximate the quotient and remainder
q = (long)(((ulong)dividend >> 1) / ((ulong)dLong >> 1));
r = dividend - q * dLong;
// Correct the approximation
while (r < 0)
{
r += dLong;
q--;
}
while (r >= dLong)
{
r -= dLong;
q++;
}
// n - q*dlong == r && 0 <= r <dLong, hence we're done.
return (r << 32) | (q & LONG_MASK);
}
private static void CopyAndShift(uint[] src, int srcFrom, int srcLen, uint[] dst, int dstFrom, int shift)
{
var n2 = 32 - shift;
var c = src[srcFrom];
for (var i = 0; i < srcLen - 1; i++)
{
var b = c;
c = src[++srcFrom];
dst[dstFrom + i] = (b << shift) | (c >> n2);
}
dst[dstFrom + srcLen - 1] = c << shift;
}
private static int MulSub(uint[] q, uint[] a, uint x, int len, int offset)
{
var xLong = x & LONG_MASK;
var carry = 0L;
offset += len;
for (var j = len - 1; j >= 0; j--)
{
var product = (a[j] & LONG_MASK) * xLong + carry;
var difference = q[offset] - product;
q[offset--] = (uint)difference;
carry = (long)((ulong)product >> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1 : 0);
}
return (int)carry;
}
private static int DivAdd(uint[] a, uint[] result, int offset)
{
var carry = 0L;
for (var j = a.Length - 1; j >= 0; j--)
{
var sum = (a[j] & LONG_MASK) +
(result[j + offset] & LONG_MASK) + carry;
result[j + offset] = (uint)sum;
carry = (long)((ulong)sum >> 32);
}
return (int)carry;
}
private static bool UnsignedLongCompare(long left, long right)
=> (left + long.MinValue) > (right + long.MinValue);Unsigned Burnikel-Ziegler Division
/// <summary>
/// <paramref name="divident"/> / <paramref name="divisor"/> and
/// <paramref name="divident"/> % <paramref name="divisor"/> using Burnikel-Ziegler division.
/// all integers must be positive. first is quotient, last is remainder.
/// </summary>
public static (uint[], bool)[] DivideAndRemainderBurnikelZiegler((uint[], bool) divident, (uint[], bool) divisor)
{
var r = divident.Item1.Length;
var s = divisor.Item1.Length;
if (r < s)
{
return new (uint[], bool)[] { ZERO, divident };
}
else
{
// Unlike Knuth and Barrett division, we don't check for common powers of two here because
// BZ already runs faster if both numbers contain powers of two and cancelling them has no
// additional benefit.
// step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
var m = 1 << (32 - (s / BURNIKEL_ZIEGLER_THRESHOLD).NumberOfLeadingZeros());
var j = (s + m - 1) / m; // step 2a: j = ceil(s/m)
var n = j * m; // step 2b: block length in 32-bit units
var n32 = 32L * n; // block length in bits
var sigma = (int)Math.Max(0, n32 - BitLength(divisor)); // step 3: sigma = max{T | (2^T)*B < beta^n}
if (sigma > 0)
divisor = ShiftLeft(divisor, sigma); // step 4a: shift b so its length is a multiple of n
if (sigma > 0)
divident = ShiftLeft(divident, sigma); // step 4b: shift a by the same amount
// step 5: t is the number of blocks needed to accommodate a plus one additional bit
var t = (int)((BitLength(divident) + n32) / n32);
if (t < 2)
{
t = 2;
}
// step 6: conceptually split a into blocks a[t-1], ..., a[0]
var a1 = GetBlock(divident, t - 1, t, n); // the most significant block of a
// step 7: z[t-2] = [a[t-1], a[t-2]]
var z = GetBlock(divident, t - 2, t, n); // the second to most significant block
z.Item1 = Disjoint(a1.Item1, z.Item1); // z[t-2]
// do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
(uint[], bool) ri, qi, quotient = ZERO;
for (var i = t - 2; i > 0; i--)
{
//step 8a: compute(qi, ri) such that z = b * qi + ri
var tmp = divide2n1n(z, divisor);
qi = tmp[0];
ri = tmp[1];
//step 8b: z = [ri, a[i - 1]]
z = GetBlock(divident, i - 1, t, n); // a[i-1]
z.Item1 = Disjoint(ri.Item1, z.Item1);
qi = ShiftLeft(qi, i *32* n);
quotient = Add(quotient, qi); // update q (part of step 9)
}
// final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
var tmp2 = divide2n1n(z, divisor);
qi = tmp2[0];
ri = tmp2[1];
quotient = Add(quotient, qi);
ri = ShiftRight(ri, sigma); // step 9: a and b were shifted, so shift back
return new (uint[], bool)[] { quotient, ri };
}
}
/// <summary>
/// This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
/// It divides a 2n-digit number by a n-digit number.
/// </summary>
private static (uint[], bool)[] divide2n1n((uint[], bool) left, (uint[], bool) right)
{
var n = right.Item1.Length;
// step 1: base case
if ((n & 1) != 0 || n < BURNIKEL_ZIEGLER_THRESHOLD || 16*n >= left.Item1.Length)
return DivideAndRemainderKnuth(left, right);
// step 2: view left as [a1,a2,a3,a4] where each ai is n/2 ints or less
var aUpper = GetUpper(left, n >> 1);// aUpper = [a1,a2,a3]
left = GetLower(left, n >> 1); // left = a4
// step 3: r1[0]=aUpper/b, r1[1]=aUpper%b
var r1 = divide3n2n(aUpper, right);
// step 4: r2[0]=[r1[1],left]/b, r2[1]=[r1[1],left]%b
r1[1].Item1 = Disjoint(r1[1].Item1, left.Item1); // r1[1] = [r1[1],left]
var r2 = divide3n2n(r1[1], right);
// step 5: let quotient=[r1[0],r2[0]] and return (quotient, r2[1])
var quotient = Disjoint(r1[0].Item1, r2[0].Item1);
return new (uint[], bool)[] { (quotient, true), r2[1] };
}
/// <summary>
/// This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
/// It divides a 3n-digit number by a 2n-digit number.
/// </summary>
private static (uint[], bool)[] divide3n2n((uint[], bool) left, (uint[], bool) right)
{
var n = right.Item1.Length >> 1; // half the length of b in ints
// step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
var a12 = GetUpper(left, n);
var a3 = GetLower(left, n);
// step 2: view b as [b1,b2] where each bi is n ints or less
var b1 = GetUpper(right, n);
var b2 = GetLower(right, n);
(uint[], bool) r, d, quotient = (new uint[n], true);
if (CompareShifted(left, right, n) < 0)
{
// step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
var tmp = divide2n1n(a12, b1);
quotient = tmp[0];
r = tmp[1];
// step 4: d=quotient*b2
d = Multiply(quotient, b2);
}
else
{
// step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
Array.Fill(quotient.Item1, uint.MaxValue);
a12 = Add(a12, b1);
b1 = ShiftLeft(b1, 32*n);
a12 = Subtract(a12, b1);
r = a12;
// step 4: d=quotient*b2=(b2 << 32*n) - b2
d = ShiftLeft(b2, 32*n);
d = Subtract(d, b2);
}
// step 5: r = r*beta^n + a3 - d (paper says a4)
// However, don't subtract d until after the while loop so r doesn't become negative
r = ShiftLeft(r, 32*n);
r = Add(r, a3);
// step 6: add b until r>=d
while (Compare(r, d) < 0)
{
r = Add(r, right);
quotient = Subtract(quotient, ONE);
}
r = Subtract(r, d);
return new (uint[], bool)[] { quotient, r};
}
/// <summary>
/// Returns a <code>(uint[], bool)</code> containing <paramref name="blockLength"/>
/// <see cref="uint"/>s from <paramref name="value"/>
/// number, starting at <paramref name="index"/>*<paramref name="blockLength"/>.<br/>
/// Used by Burnikel-Ziegler division.
/// </summary>
private static (uint[], bool) GetBlock((uint[], bool) value, int index, int numBlocks, int blockLength)
{
var blockStart = index * blockLength;
if (blockStart >= value.Item1.Length)
{
return ZERO;
}
int blockEnd;
if (index == numBlocks - 1)
{
blockEnd = value.Item1.Length;
}
else
{
blockEnd = (index + 1) * blockLength;
}
if (blockEnd > value.Item1.Length)
{
return ZERO;
}
var newVal = new uint[blockEnd - blockStart];
Array.Copy(value.Item1, value.Item1.Length - blockEnd, newVal, 0, newVal.Length);
return (newVal, true);
}
private static uint[] Disjoint(uint[] left, uint[] right)
{
Array.Resize(ref left, left.Length + right.Length);
Array.Copy(right, 0, left, left.Length - right.Length, right.Length);
return left;
}
private static int CompareShifted((uint[], bool) left, (uint[], bool) right, int ints)
{
var blen = right.Item1.Length;
var alen = left.Item1.Length - ints;
if (alen < blen)
return -1;
if (alen > blen)
return 1;
// Add Integer.MIN_VALUE to make the comparison act as unsigned integer
// comparison.
for (int i = 0, j = 0; i < alen; i++, j++)
{
var b1 = left.Item1[i] + 0x80000000;
var b2 = right.Item1[j] + 0x80000000;
if (b1 < b2)
return -1;
if (b1 > b2)
return 1;
}
return 0;
}
/// <summary>
/// <paramref name="value"/> << <paramref name="k"/> (<paramref name="k"/> > 0)
/// </summary>
public static (uint[], bool) ShiftLeft((uint[], bool) value, int k)
{
if (IsZero(value)) return ZERO;
if (k == 0) return value;
if (k < 0) return ShiftRight(value, -k);
var nInts = k >> 5;
var nBits = k & 0x1F;
var len = value.Item1.Length;
uint[] result;
if (nBits == 0)//k = 0 (mod 32)
{
result = value.Item1.ToArray();
Array.Resize(ref result, len + nInts);
}
else
{
var i = 0;
var nBits2 = 32 - nBits;
var highBits = value.Item1[0] >> nBits2;
if (highBits != 0)
{
result = new uint[len + nInts + 1];
result[i++] = highBits;
}
else
{
result = new uint[len + nInts];
}
var j = 0;
while (j < len - 1)
result[i++] = (value.Item1[j++] << nBits) | (value.Item1[j] >> nBits2);
result[i] = value.Item1[j] << nBits;
}
return (result, value.Item2);
}
/// <summary>
/// <paramref name="value"/> >> <paramref name="k"/> (<paramref name="k"/> > 0)
/// </summary>
public static (uint[], bool) ShiftRight((uint[], bool) value, int k)
{
if (IsZero(value)) return ZERO;
if (k == 0) return value;
if (k < 0) return ShiftLeft(value, -k);
var nInts = k >> 5;
var nBits = k & 0x1F;
var len = value.Item1.Length;
uint[] result;
// Special case: entire contents shifted off the end
if (nInts >= len)
return (new uint[] { 0 }, true);
if (nBits == 0)// k = 0 (mod 32)
{
var newMagLen = len - nInts;
result = value.Item1.Take(newMagLen).ToArray();
}
else
{
var i = 0;
var highBits = value.Item1[0] >> nBits;
if (highBits != 0)
{
result = new uint[len - nInts];
result[i++] = highBits;
}
else
{
result = new uint[len - nInts - 1];
}
var nBits2 = 32 - nBits;
var j = 0;
while (j < len - nInts - 1)
result[i++] = (value.Item1[j++] << nBits2) | (value.Item1[j] >> nBits);
}
if (!value.Item2)
{
// Find out whether any one-bits were shifted off the end.
var onesLost = false;
for (int i = len - 1, j = len - nInts; i >= j && !onesLost; i--)
onesLost = (value.Item1[i] != 0);
if (!onesLost && nBits != 0)
onesLost = (value.Item1[len - nInts - 1] << (32 - nBits) != 0);
if (onesLost)
result = Increment(result);
}
return (result, value.Item2);
}
private static uint[] Increment(uint[] value)
{
var lastSum = 0u;
for (var i = value.Length - 1; i >= 0 && lastSum == 0; i--)
lastSum = (value[i] += 1);
if (lastSum == 0)
{
value = new uint[value.Length + 1];
value[0] = 1;
}
return value;
}
/// <summary>
/// <para>Returns the number of bits in the minimal two's-complement
/// representation of <paramref name="value"/>, excluding a sign bit.</para>
/// For positive BigIntegers, this is equivalent to the number of bits in
/// the ordinary binary representation.For zero this method returns 0
/// </summary>
public static int BitLength((uint[], bool) value)
{
var array = StripLeadingZeroInts(value.Item1);
return array.Length == 0
? 0
: ((array.Length - 1) << 5) + BitLength(array[0]);
}
/// <summary>
/// Returns the number of bits in the minimal two's-complement
/// representation of <paramref name="value"/>
/// </summary>
private static int BitLength(uint value)
=> 32 - value.NumberOfLeadingZeros();
/// <summary>
/// Returns the number of bits in the minimal two's-complement
/// representation of <paramref name="value"/>
/// </summary>
private static int BitLength(int value)
=> 32 - value.NumberOfLeadingZeros();
private static (uint[], bool) GetLower((uint[], bool) value, int n)
{
if (value.Item1.Length <= n)
{
return Abs(value);
}
var result = value.Item1.Skip(value.Item1.Length - n).Take(n).ToArray();
result = StripLeadingZeroInts(result);
return (result, value.Item2);
}
private static (uint[], bool) GetUpper((uint[], bool) value, int n)
{
if (value.Item1.Length <= n)
{
return ZERO;
}
var result = value.Item1.Take(value.Item1.Length - n).ToArray();
result = StripLeadingZeroInts(result);
return (result, value.Item2);
}Unsigned Barrett Division
/// <summary>
/// <paramref name="divident"/> / <paramref name="divisor"/> and
/// <paramref name="divident"/> % <paramref name="divisor"/> using Barrett division.
/// all integers must be positive. first tuple is quotient, last is remainder.
/// </summary>
public static (uint[], bool)[] DivideAndRemainderBarrett((uint[], bool) divident, (uint[], bool) divisor)
{
var m = BitLength(divident);
var n = BitLength(divisor);
if (m < n)
return new (uint[], bool)[] { ZERO, divident };
else
{
var lowestSetBit = Math.Min(GetLowestSetBit(divident), GetLowestSetBit(divisor));
// Cancel common powers of two if it will speed up division, which is
// the case iff it reduces bitLength() below the next lower power of two.
if (n.NumberOfLeadingZeros() < (n - lowestSetBit).NumberOfLeadingZeros())
{
var result = DivideAndRemainder(ShiftRight(divident, lowestSetBit), ShiftRight(divisor, lowestSetBit));
result[1] = ShiftLeft(result[1], lowestSetBit);
return result;
}
if (m <= 2 * n)
{
// this case is handled by Barrett directly
var mu = Inverse(divisor, m - n);
return BarrettBase(divident, divisor, mu);
}
else
{
// treat each n-bit piece of a as a digit and do long division by val
// (which is also n bits), reusing the inverse
var mu2n = Inverse(divisor, n);
var startBit = m / n * n; // the bit at which the current n-bit piece starts
var quotient = ZERO;
var remainder = ShiftRight(divident, startBit);
var mask = Subtract(ShiftLeft(ONE, n), ONE); // n ones
while (startBit > 0)
{
startBit -= n;
var ai = And(ShiftRight(divident, startBit), mask);
remainder = Add(ShiftLeft(remainder, n), ai);
var mu = ShiftRightRounded(mu2n, 2 * n - BitLength(remainder)); // mu = 2^(remainder.length-n)/val
var c = BarrettBase(remainder, divisor, mu);
quotient = Add(ShiftLeft(quotient, n), c[0]);
remainder = c[1];
}
return new (uint[], bool)[] { quotient, remainder };
}
}
}
/// <summary>
/// <paramref name="dividend"/> / <paramref name="divisor"/> and <paramref name="dividend"/> % <paramref name="divisor"/>.
/// The binary representation of <paramref name="divisor"/> must be at least half as
/// long, and no longer than, the binary representation of <paramref name="dividend"/>.
/// </summary>
private static (uint[], bool)[] BarrettBase((uint[], bool)dividend, (uint[], bool) divisor, (uint[], bool) mu)
{
var m = BitLength(dividend);
var n = BitLength(divisor);
var a1 = ShiftRight(dividend, n - 1);
var q = Multiply(a1, mu);
q = ShiftRight(q, m - n + 1);
var r = Subtract(dividend, Multiply(divisor, q));
while (!r.Item2 || Compare(r, divisor) >= 0)
{
if (!r.Item2)
{
r = Add(r, divisor);
q = Subtract(q, ONE);
}
else
{
r = Subtract(r, divisor);
q = Add(q, ONE);
}
}
return new (uint[], bool)[] { q, r };
}
/// <summary>
/// 2 ^ (bitLength(<paramref name="value"/>) + <paramref name="n"/>).
/// </summary>
private static (uint[], bool) Inverse((uint[], bool) value, int n)
{
var m = BitLength(value);
if (n <= NEWTON_THRESHOLD)
return Divide(ShiftLeft(ONE, n * 2), ShiftRightRounded(value, m - n));
// let numSteps = ceil(log2(n/NEWTON_THRESHOLD)) and initialize k
var numSteps = BitLength((n + NEWTON_THRESHOLD - 1) / NEWTON_THRESHOLD);
var k = new int[numSteps];
var ki = n;
for (var i = numSteps - 1; i >= 0; i--)
{
ki = (ki + 1) / 2;
k[i] = ki < NEWTON_THRESHOLD ? NEWTON_THRESHOLD : ki;
}
// calculate 1/this truncated to k0 fraction digits
var z = ShiftLeft(ONE, k[0] * 2);//.divide(shiftRightRounded(m - k[0])); // exp=k0 because exp(this)=m
z = Divide(z, ShiftRightRounded(value, m - k[0]));
for (var i = 0; i < numSteps; i++)
{
ki = k[i];
// the following BigIntegers represent numbers of the form a*2^(-exponent)
var s = Square(z); // exponent = 2ki
var t = ShiftRightRounded(value, m - 2 * ki - 3); // exponent = 2ki+3
var u = Multiply(t, s); // exponent = 4ki+3 > 2ki+1
var w = Add(z, z); // exponent = ki
w = ShiftLeft(w, 3 * ki + 3); // increase #fraction digits to 4ki+3 to match u
z = Subtract(w, u); // exponent = 4ki+3
if (i < numSteps - 1)
z = ShiftRightRounded(z, 4 * ki + 3 - k[i + 1]); // reduce #fraction digits to k[i+1]
else
z = ShiftRightRounded(z, 4 * ki + 3 - n); // final step: reduce #fraction digits to n
}
return z;
}
private static (uint[], bool) ShiftRightRounded((uint[], bool) value, int n)
{
var result = ShiftRight(value, n);
if (n > 0 && TestBit(value, n - 1))
result = Add(result, ONE);
return result;
}
/// <summary>
/// Returns <see cref="true"/> if and only if the designated bit is set.
/// (Computes ((value & (1 << n)) != 0).)
/// </summary>
public static bool TestBit((uint[], bool) value, int n)
{
if (n < 0)
throw new ArithmeticException("Negative bit address");
return (GetInt(value, n >> 5) & (1 << (n & 31))) != 0;
}
/// <summary>
/// Returns the specified int of the little-endian two's complement
/// representation(int 0 is the least significant). The int number can
/// be arbitrarily high(values are logically preceded by infinitely many
/// sign ints).
/// </summary>
private static int GetInt((uint[], bool) value, int n)
{
if (n < 0)
return 0;
if (n >= value.Item1.Length)
return value.Item2 ? 0 : -1;
var magInt = (int)value.Item1[value.Item1.Length - n - 1];
return (value.Item2 ? magInt :
(n <= FirstNonzeroIntNum(value) ? -magInt : ~magInt));
}
/// <summary>
/// Returns the index of the int that contains the first nonzero int in the
/// little-endian binary representation of the magnitude(int 0 is the
/// least significant). If the magnitude is zero, return value is undefined.
/// </summary>
private static int FirstNonzeroIntNum((uint[], bool) value)
{
int i;
var mlen = value.Item1.Length;
for (i = mlen - 1; i >= 0 && value.Item1[i] == 0; i--)
;
return mlen - i - 1;
}Signed Division
private static readonly int KNUTH_POW2_THRESH_ZEROS = 3;
private static readonly int KNUTH_POW2_THRESH_LEN = 6;
private static readonly int NEWTON_THRESHOLD = 100;
private static readonly int BURNIKEL_ZIEGLER_THRESHOLD = 80;
private static readonly int BURNIKEL_ZIEGLER_OFFSET = 40;
private static readonly (uint[], bool) ZERO = (new uint[] { 0 }, true);
private static readonly (uint[], bool) ONE = (new uint[] { 1 }, true);
private static readonly (uint[], bool) MINUSONE = (new uint[] { 1 }, false);
/// <summary>
/// stimates whether Barrett Division will be more efficient than Burnikel-Ziegler when
/// dividing two numbers of a given length in ints.
/// </summary>
private static bool ShouldDivideBarrett(int length)
{
if (Environment.Is64BitOperatingSystem)
{
// The following values were determined experimentally on a 64-bit JVM.
if (length <= 14863)
return false;
if (length <= 16383) // 2^14-1
return true;
if (length <= 23759)
return false;
if (length <= 32767) // 2^15-1
return true;
if (length <= 46063)
return false;
if (length <= 65535) // 2^16-1
return true;
if (length <= 73871)
return false;
if (length <= 131071) // 2^17-1
return true;
if (length <= 159887)
return false;
return true;
}
else
{
// The following values were determined experimentally on a 32-bit JVM.
if (length <= 16175)
return false;
if (length <= 16383) // 2^14-1
return true;
if (length <= 26991)
return false;
if (length <= 32767) // 2^15-1
return true;
if (length <= 51215)
return false;
if (length <= 65535) // 2^16-1
return true;
if (length <= 80367)
return false;
if (length <= 131071) // 2^17-1
return true;
if (length <= 155663)
return false;
return true;
}
}
private static (uint[], bool) Negate((uint[], bool) value)
=> (value.Item1, !value.Item2);
private static (uint[], bool) Abs((uint[], bool) value)
=> (value.Item1, true);
private static bool IsZero((uint[], bool) value)
=> StripLeadingZeroInts(value.Item1).Count() == 0;
private static bool IsOne((uint[], bool) value)
=> StripLeadingZeroInts(value) == ONE;
private static bool IsMinusOne((uint[], bool) value)
=> StripLeadingZeroInts(value) == MINUSONE;
private static uint[] StripLeadingZeroInts(uint[] value)
=> value.SkipWhile(num => num == 0u).ToArray();
private static (uint[], bool) StripLeadingZeroInts((uint[], bool) value)
=> (StripLeadingZeroInts(value.Item1), value.Item2);
/// <summary>
/// 大数除法，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) Divide((uint[], bool) divident, (uint[], bool) divisor)
{
if (IsZero(divident))
return ZERO;
if (IsZero(divisor))
throw new DivideByZeroException();
if (IsOne(divisor))
return divident;
if (IsMinusOne(divisor))
return Negate(divident);
if (IsOne(divident))
return ZERO;
if (IsMinusOne(divident))
return ZERO;
var quotientSign = divident.Item2 == divisor.Item2;
divident = Abs(divident);
divisor = Abs(divisor);
(uint[], bool) result;
if (divident.Item1.Length < BURNIKEL_ZIEGLER_THRESHOLD ||
divident.Item1.Length - divisor.Item1.Length < BURNIKEL_ZIEGLER_OFFSET)
{
result = DivideAndRemainderKnuth(divident, divisor)[0];
}
else if (!ShouldDivideBarrett(divident.Item1.Length) || !ShouldDivideBarrett(divisor.Item1.Length))
{
result = DivideAndRemainderBurnikelZiegler(divident, divisor)[0];
}
else
{
result = DivideAndRemainderBarrett(divident, divisor)[0];
}
result.Item1 = StripLeadingZeroInts(result.Item1);
result.Item2 = quotientSign;
return result;
}
/// <summary>
/// 大数除法，返回值第一个是系数，最后一个是余数，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool)[] DivideAndRemainder((uint[], bool) divident, (uint[], bool) divisor)
{
if (IsZero(divident))
return new (uint[], bool)[] { ZERO, ZERO };
if (IsZero(divisor))
throw new DivideByZeroException();
if (IsOne(divisor))
return new (uint[], bool)[] { divident, ZERO };
if (IsMinusOne(divisor))
return new (uint[], bool)[] { Negate(divident), ZERO };
if (IsOne(divident))
return new (uint[], bool)[] { ZERO, ONE };
if (IsMinusOne(divident))
return new (uint[], bool)[] { ZERO, MINUSONE };
var remainderSign = divident.Item2;
var quotientSign = divident.Item2 == divisor.Item2;
divident = Abs(divident);
divisor = Abs(divisor);
(uint[], bool)[] result;
if (divident.Item1.Length < BURNIKEL_ZIEGLER_THRESHOLD ||
divident.Item1.Length - divisor.Item1.Length < BURNIKEL_ZIEGLER_OFFSET)
{
result = DivideAndRemainderKnuth(divident, divisor);
}
else if (!ShouldDivideBarrett(divident.Item1.Length) || !ShouldDivideBarrett(divisor.Item1.Length))
{
result = DivideAndRemainderBurnikelZiegler(divident, divisor);
}
else
{
result = DivideAndRemainderBarrett(divident, divisor);
}
result[0].Item1 = StripLeadingZeroInts(result[0].Item1);
result[0].Item2 = quotientSign;
result[1].Item1 = StripLeadingZeroInts(result[1].Item1);
result[1].Item2 = remainderSign;
return result;
}
测试
[TestMethod]
public void DivideKnuthTest()
{
for (var i = 0; i < 100; ++i)
{
var dividend = BigInteger.Abs(RandomBigInteger());
var divisor = BigInteger.Abs(RandomBigInteger());
var test = DivideAndRemainderKnuth(ToTuple(dividend), ToTuple(divisor));
var expected = dividend / divisor;
Assert.AreEqual(expected, ValueOf(test[0]));
expected = BigInteger.Remainder(dividend, divisor);
Assert.AreEqual(expected, ValueOf(test[1]));
}
}
[TestMethod]
public void DivideAndRemainderBurnikelZieglerTest()
{
var ran = new Random();
for (var i = 0; i < 100; ++i)
{
var dividend = BigInteger.Abs(RandomBigInteger());
var divisor = BigInteger.Abs(RandomBigInteger());
var test = DivideAndRemainderBurnikelZiegler(ToTuple(dividend), ToTuple(divisor));
var expected = dividend / divisor;
Assert.AreEqual(expected, ValueOf(test[0]));
expected = BigInteger.Remainder(dividend, divisor);
Assert.AreEqual(expected, ValueOf(test[1]));
}
}
[TestMethod]
public void DivideAndRemainderBarrettTest()
{
var ran = new Random();
for (var i = 0; i < 100; ++i)
{
var dividend = BigInteger.Abs(RandomBigInteger());
var divisor = BigInteger.Abs(RandomBigInteger());
var test = DivideAndRemainderBarrett(ToTuple(dividend), ToTuple(divisor));
var expected = dividend / divisor;
Assert.AreEqual(expected, ValueOf(test[0]));
expected = BigInteger.Remainder(dividend, divisor);
Assert.AreEqual(expected, ValueOf(test[1]));
}
}
[TestMethod]
public void DivideAndRemainderTest()
{
var ran = new Random();
for (var i = 0; i < 100; ++i)
{
var dividend = RandomBigInteger();
var divisor = RandomBigInteger();
var test = DivideAndRemainder(ToTuple(dividend), ToTuple(divisor));
var expected = BigInteger.DivRem(dividend, divisor, out BigInteger rem);
Assert.AreEqual(expected, ValueOf(test[0]));
Assert.AreEqual(rem, ValueOf(test[1]));
}
}
[TestMethod]
public void DivideTest()
{
var ran = new Random();
for (var i = 0; i < 100; ++i)
{
var dividend = RandomBigInteger();
var divisor = RandomBigInteger();
var test = Divide(ToTuple(dividend), ToTuple(divisor));
var expected = BigInteger.DivRem(dividend, divisor, out BigInteger rem);
Assert.AreEqual(expected, ValueOf(test));
}
}
private BigInteger RandomBigInteger()
{
var ran = new Random(Guid.NewGuid().GetHashCode());
var bytes = new byte[ran.Next(100, 5000)];
ran.NextBytes(bytes);
return new BigInteger(bytes);
}
扫一扫
2764
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
