对正整数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的二进制长度为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 = 1UL;
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 oddm) = MonModPowOdd(a, e, oddm);
a ^ e (mod 2^k) = ModPowPowerOf2(a, e, 2^k)。
因为gcd(oddm, 2) = 1。而2^k只有一个素因子:2。所以gcd(oddm, 2^k) = 1。
这样就可以用中国余数定理求得最终结果:a ^ e (mod 2^k * oddm) = 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 = evenmodulo * odd * eveninverse;
result += oddmodulo * even * oddinverse;
result %= modulo;
return (uint)result;
}
/// <summary>
/// <paramref name="value"/> ^ -1 (mod <paramref name="modulo"/>) if gcd(<paramref name="value"/>, <paramref name="modulo"/>) = 1。others return -1
/// </summary>
public static long ModInverse(long value, long modulo)
{
if (modulo <= 0)
throw new ArgumentOutOfRangeException("modulo should be positive");
if (modulo == 1) return 0;
var mod = Mod(value, modulo);
if (mod == 1) return 1;
var result = ExtendedEuclid(mod, modulo);
return result[2] == 1 ? Mod(result[0], modulo) : -1;
}
/// <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) return new long[] { 0, 1, b };
if (b == 0) return new long[] { 1, 0, a };
var anegtive = a < 0;
var bnegtive = b < 0;
if (anegtive) a = -a;
if (bnegtive) b = -b;
var aZeros = a.NumberOfTrailingZeros();
var bZeros = b.NumberOfTrailingZeros();
var k = Math.Min(aZeros, bZeros);
a >>= k; b >>= k;
var u = new long[] { 1, 0, a };
var v = new long[] { 0, 1, b };
var t = new long[3];
while (v[2] != 0)
{
var q = u[2] / v[2];
t[0] = u[0] - v[0] * q;
t[1] = u[1] - v[1] * q;
t[2] = u[2] - v[2] * q;
u[0] = v[0]; u[1] = v[1]; u[2] = v[2];
v[0] = t[0]; v[1] = t[1]; v[2] = t[2];
}
if (anegtive) u[0] = -u[0];
if (bnegtive) u[1] = -u[1];
u[2] <<= k;
return u;
}
public static void Swap<T>(ref T left, ref T right)
{
var temp = left;
left = right;
right = temp;
}
最后就可以得到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);
}