素数测试+大数的因素分解 的概率算法

版权声明：转载时请以超链接形式标明文章原始出处和作者信息及本声明 http://www.blogbus.com/jingtianzi-logs/43281804.html

素数在密码届被广泛使用。

如何判断一个素是素数，迄今为止最有效的方法是筛数法,做个表（打表法）即素数的倍数一定是合数。筛选法的效率很高，但是遇到大素数就无能为力了。
米勒罗宾算法是一个相当著名的判断是否是素数的算法，能够很大概率的判断一个数是否是素数（接近100%，在一定范围内能达到100%）。在工程方面被广泛应用。

 

米勒罗宾算法的核心是课上提到的费马小定理：
即如果一个数n 它的 a^(n-1)%n !=1 (0<a<n) 则它一定是合数。

如果 a^(n-1)%n ==1 (0<a<n) 则它可能是合数可能是素数。

概率算法的概率就在这个 a上体现。

/*

1 随机取一个 a

2 如果 它不满足 a^(n-1)%n ==1
3 则它一定是 合数

4 退出

5 如果它满足 a^(n-1)%n ==1

6 则它是一个素数的概率是1/2

7 回到 1

*/

如果对这个过程重复50次，每次都没说它是合数，那这个数是素数的概率是多少那？ 它只有(1/2)^50可能不是素数！
很显然就像是投硬币那样如果你投了1000次，每次都投到正面 ，那你就会怀疑这个硬币2面都是正面！

米勒罗宾算法的算法很简单，但是由于需要计算a^(n-1)%n这类的问题，如果不经行优化，对于大数计算时间复杂度仍旧很高。

下面Witness(long a,long n)这个函数就是用java实现的a^(n-1)%n的比较快速的大数计算。
Witness计算的核心是"模运算"

？？？

/*

把a 分解成 t个2 和 奇数u的乘积 ，比如 144 可以分解成4个2 和一个9的乘积
于是根据"模运算"，我们可以把计算a^(n-1)%n,分解成计算 4次(a^2)迭代计算和a^9%n的计算
对于4次(a^2)迭代计算，只需一个for,执行 t次 对于a^9%n的计算可参考Modular_Exponentiation(long a, long b,long n)实现

具体方法为

/*

定义一个数d为结果，初始值为1
将9化为2进制1001(2^0方位是第一位)共四位 则做四次操作 d=d*d%n 如果当前位为1则再做 d=d*a%n
*/

*/
以上算法其实就是第一次课堂练习最后道题的算法！ 当然整个算法还有许多细节：
1比如2个64位整数的乘积超过64位，会产生溢出。

解决这个问题的方法有3个

1）使用汇编写大数乘法和模

2）使用数组编写大数
3）使用JAVA，Java有大数类，方便实惠！

所以这次开始用c++写，后来改成了java 2随机数的产生，尽量产生足够长的奇数，效果好！

3实现可以做张素数表，加速计算，我就实现打了张1900个素数的大表！大家不要鄙视打表，在计算机中很常用的。
当然米勒罗宾算法只能判断它是否是素数，但如果知道了它不是素数后，则么样才能得到它的因数那？
正常的算法是O(logn)的，当然有更快的算法"POLLARD启发式搜索" 这个算法我看了一个好长时间才看懂的。核心是中国剩余定理！相当巧妙！
这个算法是启发式的，不一定能在有限时间内找到，需要写的很好才不会遇到偶尔死循环的情况。
要优化可以考虑优化gcd的计算。这个算法我写的方法还有很多漏洞，相当不完善，自己功力也不够解释清楚。 大家有兴趣的话还是自己看书吧，祝能看懂！
当然还有留下点以为，对于分解因数有更巧妙的办法吗？

素数判定随机化算法 根据费马小定理+随机化.注意pro函数之所以可以直接写成-=而不是%=,是因为x,y<n而且过程中的和不可能>=2*n,所以采取-

/* 
*jundge whether a number is primer or not. 
*/ 
#include <iostream> 
using namespace std; 
typedef unsigned __int64 llong; 
llong mod_pro(llong x,llong y,llong n) 
{ 
    llong ret=0,tmp=x%n; 
    while(y) 
        { 
            if(y&0x1)if((ret+=tmp)>n)ret-=n; 
            if((tmp<<=1)>n)tmp-=n; 
            y>>=1; 
        } 
    return ret; 
} 
llong mod(llong a,llong b,llong c) 
{ 
    llong ret=1; 
    while(b) 
        { 
            if(b&0x1)ret=mod_pro(ret,a,c); 
            a=mod_pro(a,a,c); 
            b>>=1; 
        } 
    return ret; 
} 
llong ran() 
{ 
    llong ret=rand(); 
    return ret*rand(); 
} 
bool is_prime(llong n,int t) 
{ 
    if(n<2)return false; 
    if(n==2)return true; 
    if(!(n&0x1))return false; 
    llong k=0,m,a,i; 
    for(m=n-1;!(m&1);m>>=1,k++); 
    while(t--) 
        { 
            a=mod(ran()%(n-2)+2,m,n); 
            if(a!=1) 
                { 
                    for(i=0;i<k&&a!=n-1;i++) 
                        a=mod_pro(a,a,n); 
                     
                    if(i>=k)return false; 
                } 
        } 
    return true; 
} 
int main() 
{ 
    llong n; 
    while(scanf("%I64u",&n)!=EOF) 
        if(is_prime(n,3)) 
            cout<<"It is a prime number.\n"; 
        else 
            cout<<"It is not a prime number.\n"; 
    return 0; 
}