好文章,整理收藏。

1.费马小定理:

有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.快速乘方算法

可参见 中的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;
}

在2^31-1以下,有一些伪素数可以通过{2,3,5,7,11}的测试,

已知的有 29341 = 13 * 37 * 61 , 3215031751 = 151 * 751 * 28351