【笔记】Big Integer - Modular Exponentiation

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;
}