关键算法:
Miller-Rabbin判断素数。http://www.dxmtb.com/blog/miller-rabbin/
Pollard-Rho找出整数的任意一个因子。http://www.dxmtb.com/blog/pollard-rho/
先跟着他讲的写了一遍。太神了。
找因子的时候因为找出来的是任意一个。所以还要继续二分直到找到最小的。
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int S=30,C=240;
typedef long long ll;
ll An;
ll mul(ll a,ll b,ll c){
ll r; for(r=0;b;a=(a<<1)%c,b>>=1) if(b&1) r=(r+a)%c; return r;
}ll pow(ll a,ll b,ll c){
ll r; for(r=1;b;a=mul(a,a,c),b>>=1) if(b&1) r=mul(r,a,c); return r;
}int miller_rabbin(ll n){
if(n==2) return 1; if(n<2||!(n&1)) return 0;
int t=0; ll a,x,y,u=n-1;
while(!(u&1)) t++, u>>=1;
for(int i=0;i<S;i++){
a=rand()%(n-1)+1;
x=pow(a,u,n);
for(int j=0;j<t;j++){
y=mul(x,x,n);
if(y==1&&x!=1&&x!=n-1) return 0;
x=y;
}if(x!=1) return 0;
}return 1;
}ll gcd(ll a,ll b){
return (!b)?a:gcd(b,a%b);
}ll pollard_rho(ll n,ll c){
ll x,y,d,i=1,k=2;
x=rand()%(n-1)+1;y=x;
while(++i){
x=(mul(x,x,n)+c)%n;
d=gcd(y-x,n);
if(1<d&&d<n) return d;
if(x==y) return n;
if(i==k) y=x,k<<=1;
}
}void find(ll n,ll c){
ll m;
if(n==1) return;
if(miller_rabbin(n)){ An=min(An,n); return;}
m=n; while(m==n) m=pollard_rho(n,c--);
find(m,c); find(n/m,c);
}int main(){
int tc; scanf("%d",&tc); ll n;
while(tc--){
scanf("%lld",&n);
if(miller_rabbin(n)) printf("Prime\n");
else{ An=1ll<<60ll; find(n,1LL*C); printf("%lld\n",An);}
}return 0;
}
给定N和G。求Ans=G^(sigma(C(n,i)|i是n的约数))。Ans对999911659取模。
很高兴地发现它是质数。
其实求组合数的话。费马小定理乱搞除法取模。然后用lucas定理。很快就可以求出sigma(C(n,i)|i是n的约数)。
关键是求Ans的时候,这个指数太大了怎么破!
欢迎这个式子登场:a^b mod c=a^(b mod phi(c) + phi(c))。
注意到phi(质数)=质数-1。那么我们要模的就是999911658了。
可是lucas只对质数适用啊。
很高兴的分解质因数,发现999911658=2*3*4679*35617。
那么只需要模这4个数。然后用中国剩余定理合并。
这道题就被水了。一开始95分。发现指数里没有+phi(c),导致指数是0,挂了一个点。
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=999911659ll,mod2=999911658ll;
ll ft[]={2,3,4679,35617},ft2[]={499955829,333303886,213702,28074};
ll p[4],q[4],n,g,m=0;
ll power(ll a,ll b,ll c){
ll r=1; while(b){ if(b&1) r=r*a%c; b>>=1; a=a*a%c;}
return r;
}ll calc(ll a,ll b,ll c){
if(a<b) return 0; if(b>a-b) b=a-b; ll r=1,ca=1,cb=1;
for(ll i=0;i<b;i++) ca=ca*(a-i)%c,cb=cb*(b-i)%c;
return ca*power(cb,c-2,c)%c;
}ll lucas(ll a,ll b,ll c){
ll r=1; while(a&&b&&r) r=r*calc(a%c,b%c,c)%c,a/=c,b/=c;
return r;
}ll ext_gcd(ll a,ll b,ll&x,ll&y){
if(!b){ x=1; y=0; return a;}
ll xx,yy,d=ext_gcd(b,a%b,xx,yy);
x=yy; y=xx-a/b*yy; return d;
}ll equ(ll a,ll b,ll c){
ll x,y,d=ext_gcd(a,c,x,y); x*=b/d;
while(x>=c/d) x-=c/d; while(x<0) x+=c/d;
return x;
}ll solve(ll a,ll b){
for(int i=0;i<4;i++) p[i]=lucas(a,b,ft[i]);
for(int i=0;i<4;i++) q[i]=equ(ft2[i],1,ft[i]);
ll r=0; for(int i=0;i<4;i++) r=(r+(p[i]*q[i]%mod2)*ft2[i]%mod2)%mod2;
return r;
}int main(){
cin>>n>>g;
for(ll i=1;i*i<=n;i++) if(n%i==0){
m=(m+solve(n,i))%mod2;
if(i*i!=n) m=(m+solve(n,n/i))%mod2;
}cout<<power(g,m+mod2,mod)<<endl;
return 0;
}
观察到n的范围极小,m的范围也才1000000。那就可以枚举答案了。
枚举答案M。满足对于任何(i,j),都满足C[i]+x*P[i]=C[j]+x*P[j](mod M)的最小正整数解x大于min(L[i],L[j]),或者不存在x。
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;ll An=0,c[20],p[20],l[20]; int n;
ll ext_gcd(ll a,ll b,ll&x,ll&y){
if(!b){ x=1; y=0; return a;}
ll xx,yy,d=ext_gcd(b,a%b,xx,yy);
x=yy; y=xx-a/b*yy; return d;
}ll equ(ll a,ll b,ll c){
ll x,y,d=ext_gcd(a,c,x,y); if(b%d) return -1;
x*=b/d; while(x>c/d)x-=c/d;
while(x<0) x+=c/d; return x;
}int check(ll m){
for(int i=1;i<n;i++) for(int j=i+1;j<=n;j++){
ll d; if(p[i]>=p[j])d=equ(p[i]-p[j],c[j]-c[i],m);
else d=equ(p[j]-p[i],c[i]-c[j],m);
if(d!=-1&&d<=min(l[i],l[j]))return 1;
}return 0;
}int main(){
scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%lld%lld%lld",&c[i],&p[i],&l[i]),An=max(An,c[i]);
for(;check(An);An++); printf("%lld\n",An); return 0;
}