Montgomery Reduce
/// <summary>
/// <paramref name="left"/> * <paramref name="right"/> * R ^ -1 (mod <paramref name="modulo"/>)
/// (R = 2 ^ bitLength(<paramref name="modulo"/>))
/// </summary>
public static BigInteger MonPro(BigInteger left, BigInteger right, BigInteger modulo, BigInteger inverseM, int length)
{
var bits = left.ToBinaryString().ToCharArray().Select(ch => ch - '0').Reverse().ToArray();
Array.Resize(ref bits, length);
var result = BigInteger.Zero;
for (var i = 0; i < length; ++i)
{
var tmp = (result.IsEven ? 0 : 1) + (bits[i] == 0 || right.IsEven ? 0 : 1);
tmp = tmp != 1 || inverseM.IsEven ? 0 : 1;
result = (result + (bits[i] == 0 ? 0 : right) + (tmp == 0 ? 0 : modulo)) >> 1;
}
if (result >= modulo) result -= modulo;
return result;
}
/// <summary>
/// Converts a <see cref="BigInteger"/> to a binary string.
/// </summary>
public static string ToBinaryString(this BigInteger value)
{
var bytes = value.ToByteArray();
var idx = bytes.Length - 1;
// Create a StringBuilder having appropriate capacity.
var base2 = new StringBuilder(bytes.Length << 3);
// Convert first byte to binary.
var binary = Convert.ToString(bytes[idx], 2);
// Ensure leading zero exists if value is positive.
if (binary[0] != '0' && value.Sign == 1)
{
base2.Append('0');
}
// Append binary string to StringBuilder.
base2.Append(binary);
// Convert remaining bytes adding leading zeros.
for (idx--; idx >= 0; idx--)
{
base2.Append(Convert.ToString(bytes[idx], 2).PadLeft(8, '0'));
}
return base2.ToString();
}
MonModPowOdd
/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod <paramref name="modulo"/>) <paramref name="modulo"/> is odd
/// </summary>
private static BigInteger MonModPowOdd(BigInteger value, BigInteger exponent, BigInteger modulo)
{
var len = modulo.BitLength();
var inverse = MonGCD(BigInteger.One << (len - 1), modulo);
var tmp = MonDomainIn(value, modulo, len);
var result = MonDomainIn(BigInteger.One, modulo, len);
var bits = exponent.ToBinaryString().ToCharArray();
for (var i = 0; i < bits.Length; ++i)
{
result = MonPro(result, result, modulo, inverse[1], len);
if (bits[i] == '1')
{
result = MonPro(result, tmp, modulo, inverse[1], len);
}
}
result = MonDomainOut(result, inverse[0], modulo);
return result;
}
/// <summary>
/// bit length of <see cref="BigInteger"/>
/// </summary>
public static int BitLength(this BigInteger value)
{
if (value.Sign < 0) value = -value;
var byteCount = value.GetByteCount();
var bitCount = (byteCount - 1) << 3;
value >>= bitCount;
return bitCount + 32 - ((uint)value).NumberOfLeadingZeros();
}
/// <summary>
/// 返回{R, K},R*(2left) - K*right = 1.
/// left = 2 ^ n and right is odd
/// </summary>
public static BigInteger[] MonGCD(BigInteger left, BigInteger right)
{
var result = new BigInteger[] { 1, 0 };
var s = left;
var t = right;
while (left.Sign > 0)
{
left >>= 1;
result[1] >>= 1;
if (result[0].IsEven)
{
result[0] >>= 1;
}
else
{
result[0] = (result[0] + t) >> 1;
result[1] += s;
}
}
return result;
}
/// <summary>
/// <paramref name="value"/> * 2 ^ <paramref name="k"/> (mod <paramref name="modulo"/>)
/// </summary>
private static BigInteger MonDomainIn(BigInteger value, BigInteger modulo, int k)
{
var result = value << k;
result %= modulo;
return result;
}
/// <summary>
/// <paramref name="value"/> * <paramref name="inverseR"/> (mod <paramref name="modulo"/>)
/// </summary>
private static BigInteger MonDomainOut(BigInteger value, BigInteger inverseR, BigInteger modulo)
{
var result = value;
result *= inverseR;
result %= modulo;
return result;
}
MonModPowEven
/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod <paramref name="modulo"/>) <paramref name="modulo"/> is even
/// </summary>
private static BigInteger MonModPowEven(BigInteger value, BigInteger exponent, BigInteger modulo)
{
//modulo = oddmodulo * evenmodulo
var k = modulo.NumberOfTrailingZeros();
var evenmodulo = BigInteger.One << k;
var oddmodulo = (modulo >> k);
var even = ModPowPowerOf2(value, exponent, k);
if (oddmodulo.IsOne) return even;
//china remainder theorem
var odd = MonModPowOdd(value, exponent, oddmodulo);
var oddinverse = ModInversePowerOf2(oddmodulo, k);
var eveninverse = ModInverse(evenmodulo, oddmodulo);
var result = evenmodulo * odd * eveninverse;
result += oddmodulo * even * oddinverse;
result %= modulo;
return result;
}
/// <summary>
/// <see cref="BigInteger"/>二进制尾数0的个数
/// </summary>
public static int NumberOfTrailingZeros(this BigInteger n)
{
n = ~n & (n - 1);
var c = 0;
while (!n.IsZero)
{
++c;
n >>= 1;
}
return c;
}
/// <summary>
/// <paramref name="value"/> ^ -1 (mod <paramref name="modulo"/>)
/// </summary>
public static BigInteger ModInverse(BigInteger value, BigInteger modulo)
{
if (modulo.Sign <= 0)
throw new ArithmeticException("modulo should be positive");
if (modulo.IsOne)
return BigInteger.Zero;
if (value.Sign < 0 || value >= modulo)
value = Mod(value, modulo);
if (value.IsOne)
return BigInteger.One;
var u = modulo; var v = value;
var r = BigInteger.Zero; var s = BigInteger.One;
var vLen = 0;
while ((vLen = v.BitLength()) > 1)
{
var shift = u.BitLength() - vLen;
if (u.Sign != v.Sign)
{
u += (v << shift);
r += (s << shift);
}
else
{
u -= (v << shift);
r -= (s << shift);
}
if (u.BitLength() < vLen)
{
Swap(ref u, ref v);
Swap(ref r, ref s);
}
}
if (v.IsZero) return BigInteger.Zero;//inverse does not exist
if (v.Sign < 0) s = -s;
if (s >= modulo) s -= modulo;
if (s.Sign < 0) s += modulo;
return s;
}
/// <summary>
/// <paramref name="value"/> ^ -1 (mod 2 ^ <paramref name="k"/>)
/// assumes <paramref name="value"/> is odd.
/// </summary>
public static BigInteger ModInversePowerOf2(BigInteger value, int k)
{
if (value.IsEven)
throw new ArithmeticException("value should be odd");
if (k == 0)
return BigInteger.Zero;
var modulo = BigInteger.One << k;
var mod = value & (modulo - 1);
if (mod.Sign < 0) mod += modulo;
if (mod.IsOne) return 1;
var bits = new int[k];
var tmp = BigInteger.One;
for (var i = 0; i < k; ++i)
{
bits[i] = tmp.IsEven ? 0 : 1;
tmp = (tmp - (bits[i] == 0 ? 0 : mod)) >> 1;
}
var bytes = new byte[(k + 8) / 8];
for (int i = 0, j = 0, t = 0; i < bits.Length; i += 8, ++j)
{
t = 0;
for (var m = i; m < i + 8 && m < bits.Length; ++m)
t += (bits[m] << (m - i));
bytes[j] = (byte)(t & 0xFF);
}
return new BigInteger(bytes);
}
/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod 2 ^ <paramref name="k"/>).
/// </summary>
public static BigInteger ModPowPowerOf2(BigInteger value, BigInteger exponent, int k)
{
var modulo = (BigInteger.One << k) - 1;
value &= modulo;
if (value.Sign < 0) value += modulo + 1;
var result = BigInteger.One;
var bits = exponent.ToBinaryString().ToCharArray();
for (var i = 0; i < bits.Length; ++i)
{
result = (result * result) & modulo;
if (bits[i] == '1')
{
result = (result * value) & modulo;
}
}
return result;
}
MonModPow
/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod 2 ^ <paramref name="k"/>).
/// </summary>
public static BigInteger MonModPow(BigInteger value, BigInteger exponent, BigInteger modulo)
{
if (modulo.Sign <= 0)
throw new ArgumentException("modulo should be positive");
if (exponent.Sign < 0)
throw new ArgumentException("exponent should be nonegative");
if (value.Sign < 0 || value >= modulo)
value = Mod(value, modulo);
if (modulo.IsOne)
return 0;
if (exponent.IsZero)
return 1;
if (value.IsOne)
return 1;
if (value.IsZero)
return 0;
return modulo.IsEven
? MonModPowEven(value, exponent, modulo)
: MonModPowOdd(value, exponent, modulo);
}
/// <summary>
/// <paramref name="value"/> (mod <paramref name="modulo"/>)
/// </summary>
public static BigInteger Mod(BigInteger value, BigInteger modulo)
{
if (modulo.Sign <= 0)
throw new ArithmeticException("modulo should be positive");
if (modulo.IsOne)
return BigInteger.Zero;
value = BigInteger.Remainder(value, modulo);
return value.Sign < 0 ? value + modulo : value;
}