【日記】Montgomery Moduluar Exponentiation
阿新 • • 發佈：2018-11-06
對正整數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);
}