Poj 3696 miller-rabin pollard-rho解法
挺好的题目
（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;
}
