Poj 3696 millerrabin pollardrho解法

挺好的题目

(1)888…8可以表示为8*10^0+8*10^1+…+8*10^(n-1)

可看出,抽出每一项组成等比数列,根据等比数列和公式,得,8/9(10^n-1)

所以8/9(10^n-1)=kL,k为正整数,即L|(8/9(10^n-1))

(2)令gcd(8,L)=s,故(10^n-1)*8/s=9kL/s

8,9互素,故8/s与9L/s互素,

故10^n-1=0(mod 9L/s)

    10^n=1(mod 9L/s)

若10和9L/s不互素,则10^n(mod 9L/s)必不为1,

所以10和9L/s互素(欧拉定理)

(3)10和9L/s互素,即10在模9L/s的乘法群内,10^n构成循环子群,设其规模为s(即,要求的最小的n),

10^s=1(mod 9L/s)

    故s|phi(9L/s)(拉格朗日定理),

(4)所以要求s,只需从小到大枚举求出phi(9L/s)的每个约数,有满足10^s=1(mod 9L/s),输出s

(5)我写的代码纯粹是为了练 miller-rabin 和pollard,看网上的代码打印素数表也行

 

代码:

#include <stdio.h>

#include <iostream>

#include <stdlib.h>

#include <time.h>

#include <cmath>

#include <algorithm>

#include <string.h>

using namespace std;

typedef long long ll;

 

ll prfactor[1000],prcnt[1000],factor[1000],ans,primenum,factornum;

 

bool cmp(ll a,ll b){return a<b;}

 

 

ll gcd(ll a,ll b)

{

if(b==0) return a;

else return gcd(b,a%b);

}

 

ll mulmod(ll a,ll b,ll n) //a*b%n

{

ll ret=0;

a%=n;

while(b>=1)

{

if(b&1)

{

ret+=a;

if(ret>=n) ret-=n;

}

a<<=1;

if(a>=n) a-=n;

b>>=1;

}

return ret;

}

 

ll exmod(ll a,ll b,ll n)

{

ll ret=1;

a%=n;

while(b>=1)

{

if(b&1)

ret=mulmod(ret,a,n);

a=mulmod(a,a,n);

b>>=1;

}

return ret;

}

 

bool witness(ll a,ll n)

{

int i,t=0;//n-1=m=2^t*u,t为2的指数

ll m=n-1,x,y;

while(m%2==0) {m>>=1;t++;}

x=exmod(a,m,n);//x=a^u

for(i=1;i<=t;i++)

{

y=exmod(x,2,n);

if(y==1 && x!=1 && x!=n-1)

return 1;

x=y;

}

if(y!=1) return 1;

return 0;

}

 

bool millerrabin(ll n,int times=10)

{

ll a;int i;

if(n==2) return 1;

if(n==1 || n%2==0) return 0;

srand(time(NULL));

for(i=1;i<=times;i++)

{

a=rand()%(n-1)+1;

if(witness(a,n)) return 0;

}

 

return 1;

 

}

 

ll rho(ll n,int c)

{

ll i,k,x,y,d;

srand(time(NULL));

i=1;k=2;

y=x=rand()%n;

while(1)

{

i++;

x=(mulmod(x,x,n)+c)%n;

d=gcd(y-x,n);

if(d>1 && d<n) return d;

if(y==x) break;

if(i==k){y=x;k*=2;}

}

return n;

}

 

void pollard(ll n,int c)

{

if(n<=1) return;

if(millerrabin(n)) {prfactor[primenum++]=n;return;}

ll m=n;

while(m>=n)

m=rho(m,c--);

pollard(m,c);

pollard(n/m,c);

}

/*

void dfs(ll i,ll x,ll k)

{

if(i>primenum) return ;

if(x>ans && x<=k) ans=x;

dfs(i+1,x,k);

x*=prfactor[i];

if(x>ans && x<=k) ans=x;

dfs(i+1,x,k);

}*/

 

ll eu(ll n)

{

for(int i=0;i<=primenum-1;i++)

n=n/prfactor[i]*(prfactor[i]-1);

return n;

}

 

void dfs(ll val,ll d)//dfs求一个数的所有因子

{

if(d==primenum)

{

factor[factornum++]=val;

return;

}

ll sum=1;

for(int i=1;i<=prcnt[d];i++)

{

sum*=prfactor[d];

dfs(val*sum,d+1);

}

dfs(val,d+1);

}

 

int main()

{

ll a,tmp,s,q;

int i,j,tt=0;

freopen("in.txt","r",stdin);

while(scanf("%lld",&a)!=EOF)

{

if(a==0)

break;

a=a/gcd(a,8)*9;

if(gcd(a,10)!=1)

{

printf("Case %d: 0\n",++tt);

continue;

}

 

primenum=0;

factornum=0;

pollard(a,107);

sort(prfactor,prfactor+primenum,cmp);

int cnt1=0,cnt2=0;

memset(prcnt,0,sizeof(prcnt));

memset(factor,0,sizeof(factor));

for(i=0;i<=primenum-1;i++)

{

while(i<primenum-1&&prfactor[i]==prfactor[i+1])

{

prcnt[cnt2]++;

i++;

}

prcnt[cnt2]++;

prfactor[cnt1++]=prfactor[i];

cnt2++;

}

primenum=cnt1;

 

ll phi=eu(a);

 

primenum=0;

factornum=0;

pollard(phi,107);

sort(prfactor,prfactor+primenum,cmp);

cnt1=0,cnt2=0;

memset(prcnt,0,sizeof(prcnt));

memset(factor,0,sizeof(factor));

for(i=0;i<=primenum-1;i++)

{

while(i<primenum-1&&prfactor[i]==prfactor[i+1])

{

prcnt[cnt2]++;

i++;

}

prcnt[cnt2]++;

prfactor[cnt1++]=prfactor[i];

cnt2++;

}

primenum=cnt1;

 

dfs(1,0);

sort(factor,factor+factornum,cmp);

for(i=0;i<factornum;i++)

if(exmod(10,factor[i],a)==1)

{printf("Case %d: %lld\n",++tt,factor[i]);break;}

 

}

return 1;

}

 

原文地址:https://www.cnblogs.com/inpeace7/p/2397601.html