zoukankan      html  css  js  c++  java
  • [模板] Miller_Rabin和Pollard_Rho

    Miller_Rabin

    快速($O(slogn)$,s为尝试次数)地测试一个数是否是质数

    首先有费马小定理$a^{p-1}=1 (mod p)$当p为质数时成立,所以可以随机选择a来以这个式子作为一定的判断依据,但并不是所有合数都不满足这个式子,甚至存在合数对所有的a都不满足这个式子

    然后有二次探测定理:若p是质数,$a^2=1 (mod p)$成立当且仅当$a=1 (mod p)$或$a=p-1 (mod p)$

    证明:移项可得$(a-1)(a+1)=0 (mod p)$,于是有$p|(a-1)(a+1)$,若p为质数,则只能有$p|(a-1)$或者$p(a+1)$

    所以,如果有$a^2=1 (mod p)$但a不等于1或者p-1的话,证明p一定是合数

    然而当我们随便选一个a,很难保证它的平方就等于1...但我们还有费马小定理!

    所以我们可以考虑把p-1表示成$s*2^t$的形式,然后从$a^s$一路平方上去,并进行判断

    这个a可以考虑选择前10个质数,在1e18的范围内应该都是对的

    Pollard_Rho

    快速(期望$O(n^{1/4})$)地找到n的某个非平凡因子

    首先,用Miller_Rabin来判断n是否是质数(这也是P_R和M_R唯一的关联)

    然后我们考虑随机一个数x,然后计算$gcd(x,n)$,如果结果不是1或者n,就找到了n的一个非平凡因子

    然而纯随机肯定是不行的..

    考虑一种生成随机数的方式:$x_i=x_{i-1}^2+c (mod n)$,$c$和$x_0$为随机的常数,这样生成出的数一定会成一个"ρ"形,而且根据生日悖论(期望随机$sqrt{n}$次就能得到在n的范围内相等的两个数),环长和前面一段的长度都是期望$sqrt{n}$的

    更进一步的是,如果考虑a是n的一个最小质因数,那么$x_i$在模a意义下,会形成一个期望$sqrt{a}$的ρ

    这样的话,我随机在这个模n的ρ形上取两个数做gcd,期望只要$sqrt{a}$也就是$n^{1/4}$次,就能找到一个n的因数

    于是可以使用a走一步,b走两步的方法判圈,同时计算$gcd(a-b,n)$

    但注意到一次gcd需要log的时间,所以我们可以采用每127(我也不知道这个玄学的数是哪来的)个计算一次gcd的方法来优化,就是说,每次连续算127个$a-b$并乘起来,然后算一下gcd,如果发现gcd大于1了,再重新算一下这127个看看是哪个大于1了

    例题

    luogu4718 (求n的最大质因子)
     1 #include<bits/stdc++.h>
     2 #include<tr1/unordered_map>
     3 #define CLR(a,x) memset(a,x,sizeof(a))
     4 #define MP make_pair
     5 #define fi first
     6 #define se second
     7 using namespace std;
     8 typedef long long ll;
     9 typedef unsigned long long ull;
    10 typedef long double ld;
    11 typedef pair<int,int> pa;
    12 const int maxn=233;
    13 
    14 inline ll rd(){
    15     ll x=0;char c=getchar();int neg=1;
    16     while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();}
    17     while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
    18     return x*neg;
    19 }
    20 
    21 int pri[10]={2,3,5,7,11,13,17,19,23,29};
    22 
    23 inline ll mul(ll a,ll b,ll p){
    24     return ((a*b-(ll)((ld)a*b/p)*p)%p+p)%p;
    25 }
    26 
    27 inline ll fpow(ll a,ll b,ll p){
    28     ll r=1;
    29     while(b){
    30         if(b&1) r=mul(r,a,p);
    31         a=mul(a,a,p),b>>=1;
    32     }return r;
    33 }
    34 
    35 inline bool judge(ll p){
    36     if(p==2) return 1;
    37     else if(p==1||p%2==0) return 0;
    38     ll s=0,t=p-1;
    39     while(!(t&1)) s++,t>>=1;
    40     for(int i=0;i<10&&pri[i]<p;i++){
    41         ll x=fpow(pri[i],t,p);
    42         for(int j=1;j<=s;j++){
    43             ll r=mul(x,x,p);
    44             if(r==1&&x!=1&&x!=p-1) return 0;
    45             x=r;
    46         }
    47         if(x!=1) return 0;
    48     }return 1;
    49 }
    50 
    51 inline ll gcd(ll a,ll b){return a?gcd(b%a,a):b;}
    52 inline ll rand(ll r){
    53     return (((ll)rand()<<30)|rand())%(r-1)+1;
    54 }
    55 
    56 inline ll calc(ll x){
    57     if(judge(x)) return x;
    58     while(1){
    59         ll a=rand(x-1),c=rand(x-1),b=(mul(a,a,x)+c)%x;
    60         // printf("~%lld %lld %lld
    ",a,b,c);
    61         bool bl=0;
    62         while(a!=b){
    63             ll s=1;ll aa=a,bb=b;
    64             if(!bl){
    65                 for(int i=1;i<=127&&aa!=bb;i++){
    66                     s=mul(s,aa-bb,x);
    67                     aa=(mul(aa,aa,x)+c)%x;
    68                     bb=(mul(bb,bb,x)+c)%x,bb=(mul(bb,bb,x)+c)%x;
    69                 }
    70             }
    71             ll g=gcd(s+x,x);
    72             if(g>1) bl=1;
    73             else a=aa,b=bb;
    74             if(bl){
    75                 g=gcd(a-b+x,x);
    76                 if(g>1) return max(calc(g),calc(x/g));
    77                 a=(mul(a,a,x)+c)%x;
    78                 b=(mul(b,b,x)+c)%x,b=(mul(b,b,x)+c)%x;
    79             }
    80             
    81             
    82         }
    83     }
    84 }
    85 
    86 int main(){
    87     srand(1919810);
    88     for(int T=rd();T;T--){
    89         ll x=rd();
    90         ll re=calc(x);
    91         if(re==x) puts("Prime");
    92         else printf("%lld
    ",re);
    93     }
    94     return 0;
    95 }
  • 相关阅读:
    互联网 DBA 需要做那些事(转)
    mysql_connect和mysql_pconnect区别(转)
    Redis应用案例,查找某个值的范围(转)
    PHP 正则表达式常用函数使用小结
    PHP转换UTF-8和GB2312的URL编码(转)
    PHP 打印调用函数入口地址(堆栈)
    php CI框架nginx 配置
    apache部署多个项目
    Apache+php在windows下的安装和配置
    appium测试之获取appPackage和appActivity
  • 原文地址:https://www.cnblogs.com/Ressed/p/11076543.html
Copyright © 2011-2022 走看看