费马小定理，积模分解公式，高效幂，快速幂模 及 米勒-拉宾检验
有N为任意正整数，P为素数，且N不能被P整除（显然N和P互质），则有：
N^P%P=N(即：N的P次方除以P的余数是N) 或 (N^(P-1))%P=1
互相变形：
原式可化为：
(N^P-N)%P=0
(N*(N^(P-1)-1))%P=0
所以，N*(N^(P-1)-1)是N和P的公倍数
又因为 N与P互质，而互质数的最小公倍数为它们的乘积
所以一定存在正整数M使得等式成立：N*(N^(P-1)-1)=M*N*P
所以N^(P-1)-1=M*P
因为M是整数
所以(N^(P-1)-1)%P=0
即：
N^(P-1)%P=1
2.积模分解公式
引理， 若：X%Z=0，则(X+Y)%Z=Y%Z
设有X、Y和Z三个正整数，则必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z
证明：
1. 当X和Y都比Z大时，必有整数A和B使下面的等式成立：
X=Z*I+A（1）
Y=Z*J+B（2）除模运算的性质
将（1）和（2）代入(X*Y)modZ得：
((Z*I+A)(Z*J+B))%Z
所以 (Z*(Z*I*J+I*A+I*B)+A*B)%Z（3）
所以 Z*(Z*I*J+I*A+I*B)能被Z整除
概据引理，（3）式可化简为：(A*B)%Z
又因为：A=X%Z，B=Y%Z，代入上式，成为原式右边。
2. 当X比Z大而Y比Z小时：
X=Z*I+A
代入(X*Y)%Z得：
(Z*I*Y+A*Y)%Z
根据引理，转化得：(A*Y)%Z
因为A=X%Z，又因为Y=Y%Z，代入上式，即得到原式右边。
同理，当X比Z小而Y比Z大时，原式也成立。
3. 当X比Z小，且Y也比Z小时，X=X%Z，Y=Y%Z，所以原式成立。
3.快速乘方算法
可参见http://www.cnblogs.com/bl4nk/archive/2011/04/20/2022998.html 中的power函数
unsigned Power(unsigned n, unsigned p) //n^p
{
unsigned ans = 1;
while (p > 1)
{
if (( p & 1 )!=0) ans *= n;
n *= n;
p /= 2;
}
return n * ans;
}
4."蒙格马利”快速幂模算法
__int64 Montgomery(__int64 n, __int64 p, __int64 m)
{ // 快速计算 (n ^ p) % m 的值，与power算法极类似
__int64 r = n % m; // 这里的r可不能省
__int64 k = 1;
while (p > 1)
{
if ((p & 1)!=0)
{
k = (k * r) % m; // 直接取模
}
r = (r * r) % m; // 同上
p /= 2;
}
return (r * k) % m; // 还是同上
}
蒙格马利极速版：
unsigned Montgomery(unsigned n,unsigned p,unsigned m)
{ //快速计算(n^p)%m的值
unsignedk=1;
n%=m;
while(p!=1)
{
if(0!=(p&1))k=(k*n)%m;
n=(n*n)%m;
p>>=1;
}
return(n*k)%m;
}
5.素数判断
根据费马小定理，对于两个互质的素数N和P，必有：N^(P-1)%P=1
对于N，从素数表中取出任意的素数对其进行费马测试，如果取了很多个素数，N仍未测试失败，那么则认为N是素数。当然，测试次数越多越准确，但一般来讲
50次就足够了。费马测试并不完全可靠。
现在我们发现了重要的一点，费马定理是素数的必要条件而非充分条件。这种不是素数，但又能通过费马测试的数字还有不少，数学上把它们称为Carmichael数，现在
数学家们已经找到所有10 ^
16以内的Carmichael数，最大的一个是9585921133193329。我们必须寻找更为有效的测试方法。数学家们通过对费马小定理的研究，并加以扩展，
总结出了多种快速有效的素数测试方法，目前最快的算法是米勒-拉宾检验算法。
米勒-拉宾检验 wikipedia : Miller–Rabin primality test
米勒-拉宾检验是一个不确定的算法，只能从概率意义上判定一个数可能是素数，测试用的素数表越多，准确率越高，但并不能确保。算法流程如下:
1.选择T个随机数A，并且有A<N成立。
2.找到R和M，使得N=2*R*M+1成立。
快速得到R和M的方式：N用二进制数B来表示，令C=B-1。因为N为奇数（素数都是奇数），所以C的最低位为0，从C的最低位的0开始向高位统计，一直到遇到第一个1。这时0的个数即为R，M为B右移R位的值。
3.如果A^M%N=1，则通过A对于N的测试，然后进行下一个A的测试
4.如果A^M%N!=1，那么令i由0迭代至R，进行下面的测试
5.如果A^((2^i)*M)%N=N-1则通过A对于N的测试，否则进行下一个i的测试
6.如果i=r，且尚未通过测试，则此A对于N的测试失败，说明N为合数。
7.进行下一个A对N的测试，直到测试完指定个数的A
通过验证得知，当T为素数，并且A是平均分布的随机数，那么测试有效率为1 / ( 4 ^ T )。如果T > 8那么测试失误的机率就会小于10^(-5)，这对于一般的应用是足够了。如果需要求的素数极大，或着要求更高的保障度，可以适当调高T的值。
下面是代码：
bool RabbinMillerTest( unsigned n )
{
if (n<2)
{ // 小于2的数即不是合数也不是素数
throw 0;
}
const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数
for(int i=0;i<nPrimeListSize;++i)
{// 按照素数表中的数对当前素数进行判断
if (n/2+1<=g_aPrimeList[i])
{// 如果已经小于当前素数表的数，则一定是素数
return true;
}
if (0==n%g_aPrimeList[i])
{// 余数为0则说明一定不是素数
return false;
}
}
// 找到r和m，使得n = 2^r * m + 1;
int r = 0, m = n - 1; // ( n - 1 ) 一定是合数
while ( 0 == ( m & 1 ) )
{
m >>= 1; // 右移一位
r++; // 统计右移的次数
}
const unsigned nTestCnt = 8; // 表示进行测试的次数
for ( unsigned i = 0; i < nTestCnt; ++i )
{ // 利用随机数进行测试，
int a = g_aPrimeList[ rand() % nPrimeListSize ];
if ( 1 != Montgomery( a, m, n ) )
{
int j = 0;
int e = 1;
for ( ; j < r; ++j )
{
if ( n - 1 == Montgomery( a, m * e, n ) )
{
break;
}
e <<= 1;
}
if (j == r)
{
return false;
}
}
}
return true;
}
