1. 程式人生 > >【日記】Montgomery Moduluar Exponentiation

【日記】Montgomery Moduluar Exponentiation

對正整數x, y, m, R, m',如果存在如下約束:

0 <= x, y < R*m
R = 2 ^ n, R > m
gcd(m, 2) = 1
R^-1*R - m'*m = 1

則MonPro(x, y, m, m') = (xy)R^-1 (mod m)。

function REDC(x, m, m')
    a <- (x (mod R))m' (mod R) [so 0 < m < R]
    t <- (x + am)/R
    if t > m
        then return t - m
    else
        return t 

如果資料型別為uint, 因為m最小值為1,為滿足條件1和條件2,則R >= 2 ^ 32,取R = 2 ^ 32。

條件3表明m必須為奇數。

已知m和R,可以根據擴充套件歐幾里得方法求得R^-1和m'。因為R = 2 ^ 32溢位,所以第一個引數將被*2。

/// <summary>
/// 返回{R, K},R*(2left) - K*right = 1
/// left = 2 ^ n and right is odd
/// </summary>
private static uint[] MonGCD(uint left, uint right)
{
    var result = new uint[] { 1u, 0u };
    var s = left;
    var t = right;

    while (left > 0)
    {
        left >>= 1;
        result[1] >>= 1;

        if (IsEven(result[0]))
        {
            result[0] >>= 1;
        }
        else
        {
            result[0] = ((result[0] ^ t) >> 1) + (result[0] & t);//(result[0] + t) / 2
            result[1] += s;
        }
    }

    return result;
}

private static bool IsEven(uint n) => (n & 1) == 0;

MonPro(x, y, m, m') c# 程式碼如下:

/// <summary>
/// <paramref name="left"/> * <paramref name="right"/> * R ^ -1  (mod <paramref name="modulo"/>) (R = 2 ^ 32)
/// </summary>
private static uint MonPro(uint left, uint right, uint modulo, uint inverseM)
{
    ulong result = left;
    result *= right;
    var temp = result;

    result &= 0xFFFF_FFFF;//x (mod R)
    result *= inverseM;//(x (mod R)) * m'
    result &= 0xFFFF_FFFF;//(x (mod R)) * m' (mod R)
    result *= modulo; //((x (mod R)) * m' (mod R)) * m
    result += temp;//((x (mod R)) * m' (mod R)) * m + x

    var condition = result < temp;//overflow?
    result >>= 32;//(((x (mod R)) * m' (mod R)) * m + x) / R
    
    if (condition || (result >= modulo))
    {
        result -= modulo;
    }

    return (uint)result;
}

MonPro第二個版本

/// <summary>
/// <paramref name="left"/> * <paramref name="right"/> * R ^ -1  (mod <paramref name="modulo"/>) (R = 2 ^ 32)
/// </summary>
private static uint MonPro2(uint left, uint right, uint modulo, uint inverseM)
{
    var bits = Convert.ToString(left, 2).ToCharArray().Select(ch => (uint)(ch - '0')).Reverse().ToArray();
    Array.Resize(ref bits, 32);

    var result = 0UL;
    for (var i = 0; i < 32; ++i)
    {
        var tmp = (((result & 1) + bits[i]*(right & 1)) * inverseM) & 1;
        result = (result + bits[i]*right + tmp*modulo) >> 1;
    }

    if (result >= modulo) result -= modulo;

    return (uint)result;
}

模指運算中我們要計算xy (mod m),而MonPro(x, y, m, m') = (xy) R^-1 (mod m)。根據條件4,R^-1R = 1 (mod m)。

當同時對輸入引數乘以R時:z = MonPro(xR, yR, m, m')  = (xR)(yR) R^-1 (mod m) = (xy)R (mod m)。

對z再zR^-1 (mod m) = (xyR)R^-1 (mod m) = xy (mod m)。

所以在連乘時我們要對第一個輸入引數和最終結果都加以調整。程式碼如下:

/// <summary>
/// <paramref name="value"/> * 2 ^ 32 (mod <paramref name="modulo"/>)
/// </summary>
private static uint MonDomainIn(uint value, uint modulo)
{
    ulong result = value;
    result <<= 32;
    result %= modulo;

    return (uint)result;
}

/// <summary>
/// <paramref name="value"/> * <paramref name="inverseR"/> (mod <paramref name="modulo"/>)
/// </summary>
private static uint MonDomainOut(uint value, uint inverseR, uint modulo)
{
    ulong result = value;
    result *= inverseR;
    result %= modulo;
    
    return (uint)result;
}

指數運算可以參考

 

這樣就可以得到奇數模指運算的程式碼

/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod <paramref name="modulo"/>) <paramref name="modulo"/> is odd
/// </summary>
private static uint MonModPowOdd(uint value, uint exponent, uint modulo)
{
    var inverse = MonGCD(1u << 31, modulo);

    var tmp = MonDomainIn(value % modulo, modulo);
    var result = MonDomainIn(1, modulo);
    var bits = Convert.ToString(exponent, 2).ToCharArray();
    for (var i = 0; i < bits.Length; ++i)
    {
        result = MonPro(result, result, modulo, inverse[1]);
        if (bits[i] == '1')
        {
            result = MonPro(result, tmp, modulo, inverse[1]);
        }
    }

    result = MonDomainOut(result, inverse[0], modulo);

    return result;
}

當模為偶數時,我們可以把m分解為2^k * oddM,只需要計算m二進位制字串字尾為0的個數就可以得到k。

當整數n=0時,其後綴為0的個數是32,其它時候其二進位制至少包含一個1,所以0字尾最多隻有31個。然後通過二分法,每次把二進位制字串長度減半,測試其是否=0,直至最後測試完整個字串共需要log2(32) = 5步。

計算整數尾數0的程式碼:

/// <summary>
/// <paramref name="n"/>二進位制尾數0的個數
/// </summary>
public static int NumberOfTrailingZeros(this uint n)
{
    if (n == 0) return 32;

    var c = 31;
    var t = n << 16; if (t != 0u) { c -= 16; n = t; }
    t = n << 8; if (t != 0u) { c -= 8; n = t; }
    t = n << 4; if (t != 0u) { c -= 4; n = t; }
    t = n << 2; if (t != 0u) { c -= 2; n = t; }
    t = n << 1; if (t != 0u) { c -= 1; }

    return c;
}

偶數模為特殊形式2^k。因為x (mod 2^k) = x & (2^k - 1),可以大量簡化模運算:

/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod 2 ^ k). (<paramref name="modulo"/> = 2 ^k)
/// </summary>
public static uint ModPowPowerOf2(uint value, uint exponent, uint modulo)
{
    modulo--;
    value &= modulo;

    var result = BigInteger.One;
    var bits = Convert.ToString(exponent, 2).ToCharArray();
    for (var i = 0; i < bits.Length; ++i)
    {
        result = (result * result) & modulo;
        if (bits[i] == '1')
        {
            result = (result * value) & modulo;
        }
    }

    return (uint)result;
}

已知:

a ^ e (mod oddmod) = MonModPowOdd(a, e, oddm);

a ^ e (mod evenmod) = ModPowPowerOf2(a, e, evenm)。

因為gcd(oddm, 2) = 1。而evenm = 2^k只有兩個因子:1和2。所以gcd(oddm, evenm) = 1。

這樣就可以用中國餘數定理求得最終結果:a ^ e (mod m)

/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod <paramref name="modulo"/>) <paramref name="modulo"/> is even
/// </summary>
private static uint MonModPowEven(uint value, uint exponent, uint modulo)
{
    //modulo = oddmodulo * evenmodulo
    var k = modulo.NumberOfTrailingZeros();
    var evenmodulo = 1u << k;
    var oddmodulo = (modulo >> k);

    var even = ModPowPowerOf2(value, exponent, evenmodulo);
    if (oddmodulo == 1) return even;
    //china remainder theorem
    var odd = MonModPowOdd(value, exponent, oddmodulo);
    var oddinverse = ModInverse(oddmodulo, evenmodulo);
    var eveninverse = ModInverse(evenmodulo, oddmodulo);
    var result = new BigInteger(evenmodulo) * odd * eveninverse;
    result += new BigInteger(oddmodulo) * even * oddinverse;
    result %= modulo;

    return (uint)result;
}


/// <summary>
/// <paramref name="value"/> ^ -1 (mod <paramref name="modulo"/>)
/// </summary>
public static uint ModInverse(uint value, uint modulo)
{
    if (modulo == 1) return 0;

    var mod = Mod(value, modulo);
    if (mod == 1) return 1;

    var result = ExtendedEuclid(value, modulo);
    if (result[2] != 1)
        throw new ArgumentOutOfRangeException("mod inverse does not exist");

    return (uint)result[0];
}


/// <summary>
/// 返回{x, y, gcd}, 使得<paramref name="a"/>x + <paramref name="b"/>y = gcd(<paramref name="a"/>, <paramref name="b"/>)
/// </summary>
public static long[] ExtendedEuclid(long a, long b)
{
    if (a < 0 || b < 0)
        throw new ArgumentOutOfRangeException("a and b should all be nonegative");

    if (a == 0) return new long[] { 0, 1, b };
    if (b == 0) return new long[] { 1, 0, a };

    var aZeros = a.NumberOfTrailingZeros();
    var bZeros = b.NumberOfTrailingZeros();
    var k = Math.Min(aZeros, bZeros);
    a >>= k; b >>= k;

    var u = new BigInteger[] { 1, 0, a };
    var v = new BigInteger[] { 0, 1, b };
    //from here on we maintain
    //u[0]a + u[1]b = u[2]
    //v[0]a + v[1]b = v[2]
    while (u[2].IsEven)
    {
        u[2] >>= 1;
        if (u[0].IsEven && u[1].IsEven)
        {
            u[0] >>= 1; u[1] >>= 1;
        }
        else
        {
            u[0] = (u[0] + b) >> 1; u[1] = (u[1] - a) >> 1;
        }
    }

    while (u[2] != v[2])
    {
        if (v[2].IsEven)
        {
            v[2] >>= 1;
            //Commentary: note that here, since v[2] is even
            //(i)  if v[0], v[1] are both odd then so are a, b
            //(ii) if v[0] is odd and v[1] even then a must be even, so b is odd
            //(iii)if v[1] is odd and v[0] even then b must be even, so a is odd
            //so for each of (i), (ii) and (iii) v[0] + b and v[1] - a are even
            if (v[0].IsEven && v[1].IsEven)
            {
                v[0] >>= 1; v[1] >>= 1;
            }
            else
            {
                v[0] = (v[0] + b) >> 1; v[1] = (v[1] - a) >> 1;
            }
        }
        else if (v[2] < u[2])
        {
            Swap(ref u[0], ref v[0]);
            Swap(ref u[1], ref v[1]);
            Swap(ref u[2], ref v[2]);
        }
        else
        {
            v[0] -= u[0];
            v[1] -= u[1];
            v[2] -= u[2];
        }
    }

    return new long[] { (long)v[0], (long)v[1], (long)v[2] << k };
}

最後就可以得到MonModPow程式碼:

/// <summary>
/// <paramref name="value"/> ^ <paramref name="exponent"/> (mod <paramref name="modulo"/>)
/// </summary>
public static uint MonModPow(uint value, uint exponent, uint modulo)
{
    if (modulo == 1)
        return 0U;

    if (exponent == 0)
        return 1U;

    if (value == 1)
        return 1U;
    if (value == 0)
        return 0U;

    return IsEven(modulo)
           ? MonModPowEven(value, exponent, modulo)
           : MonModPowOdd(value, exponent, modulo);
}