主要有三个点:
1、展开模乘的式子得到运算逻辑
2、处理模加中的溢出问题 (如果 a + b = c 溢出,则 a + b mod m = c mod m + 2^64 mod m)?
3、如何在模数和被模数相差十分大时,快速算出 a mod m
#include<stdio.h>
//2^64-1
static unsigned long long maxNum = -1;
//2^63
static unsigned long long smaxNum = (unsigned long long) 1 << 63;
//得到二进制下数a的第i位的值
unsigned long long bitof(unsigned long long a, unsigned long long i)
{
return (a >> i) & 1;
}
//模运算
//a和m相差很大时,仍能较快算出 a mod m
unsigned long long mod(unsigned long long a, unsigned long long m)
{
unsigned long long m1 = m;
while (a >= m)
{
//m1 >= smaxNum 左移后会溢出
while (m1 < smaxNum && m1 < a)
{
m1 = m1 << 1;
}
while (m1 > a && m1 > m)
{
m1 = m1 >> 1;
}
a -= m1;
}
return a;
}
//模加运算
//( a + b ) mod m = ( a mod m + b mod m ) mod m
//注意 a + b 的溢出
unsigned long long plusmod(unsigned long long a, unsigned long long b, unsigned long long m)
{
a = mod(a,m);
b = mod(b,m);
unsigned long long sum = a + b;
//a + b 何时发生溢出,很容易错!!
while (a != 0 && b != 0 && b - 1 >= maxNum - a)
{
a = mod(sum,m);
b = mod(maxNum,m) + mod(1,m);
sum = a + b;
}
return mod(sum,m);
}
//模乘运算
// a * b mod m = a * (bnbn-1...b1) mod m
// = ( ∑ a * bi * 2^(i-1) mod m ) mod m
unsigned long long multimod(unsigned long long a, unsigned long long b , unsigned long long m)
{
unsigned long long res = 0;
int i = 1;
for(; i < 64; ++i)
{
if(bitof(a,i) == 0) continue;
int j = 0;
unsigned long long b1 = b;
for (; j < i; ++j)
{
b1 = plusmod(b1,b1,m);
}
res = plusmod(res,b1,m);
}
if(bitof(a,0) == 1) res = plusmod(res,b,m);
return res;
}
|