1467: Pku3243 clever Y - #include <Seter> - using namespace Orz;

1467: Pku3243 clever Y

Seter posted @ 2012年2月19日 16:55 in BZOJ with tags POJ NumberTheory , 2452 阅读

http://www.zybbs.org/JudgeOnline/problem.php?id=1467

RunID User Problem Result Memory Time Language Code Length Submit Time
206331 Seter 1467 Accepted 1432 kb 52 ms C/Edit 1599 B 2012-02-17 20:55:23

终于弄明白扩展小步大步算法了囧,我的数学真是太差了。。

AC大神的代码很好懂,但是解释就……我反正看不懂。

题意:解方程A^x=B(mod C)。


首先假设C是素数。那么特判掉A是C的倍数的情况后,(A,C)=1。

由于C是素数,根据XX小定理,可以证明A^C=A(mod C),于是如果存在x,最大可能的值就是C-1。

弄一个M=sqrt(C),令x=iM+j(0<=j<M),式子转化为(A^M)^i*A^j=B(mod C)

再令p=(A^M)^i mod C,q=A^j mod C,则pq=B(mod C)。另外,(A,C)=1可以推出(p,C)=1。

我们从0开始枚举i直到i*M超过C-1,则相当于知道了p,要求q。

此时根据XX定理,[0,C)中有且仅有一个q,可以用exgcd算,至于exgcd...自己查!

这样我们就知道了q,问题转化为求是否存在A^j=q(mod C)。咦?这个和原题一样啊,难道递归?

注意到j<M=sqrt(C),已经比较小了,暴力即可!

当然,这个暴力的意思是,将A^0...A^(M-1) mod C用hash等保存下来,然后每次去里面找。

伪代码大概像这样:

 

IF A mod C = 0 特判
M = sqrt(C)
Insert (A^0...A^(M-1) mod C) into HASHMAP
FOR i = 0 to C/M do
p = (A^M)^i mod C
q = Inval(p,B,C) //求p对C的逆元然后乘上B,求逆元用exgcd
IF q is in HASHMAP return i * M + FIRST(q) //FIRST表示对于q,最早插入时的指数
done
return FAIL

当C不是素数的时候能否直接套用呢?当然不可以……最直接的问题就是,不一定存在逆元!

考虑一个G'同时是ABC的因数,令B'=B/G',C'=C/G',则当x不等于0时,(A/G')*A^(x-1)=B'(mod  C')

这样多弄几次,C的因数就越来越少,直到(A,C)=1。

那么如何选取G'呢?AC大神告诉你:不断取G'=gcd(A,C'),直到G'=1。如果任意一个G'不是B'的因数则一定无解。

假设取了r次G',然后所有G'的积是G,则问题变为(A^r/G)*A^(x-r)=B'(mod C'),令D=A^r/G(一定是整数),则由于此时(A,C')=1,所以(D,C')=1。

那么上面算法中只要变成求(D*p)对于C'的逆元就可以了,返回的答案还要加上r。

还有一点小问题,就是这样得出的答案是大于等于r的。但是即使每次G'=2,r最大也只有log2(C),那么这些再暴力求解就可以了。

 

#include <stdio.h>
#include <math.h>
#include <string.h>
#define LL long long
#define maxn 40000
#define hashp 50021
typedef struct _Hash{struct _Hash*next;int v,t}Hash,*PHash;
Hash st[maxn],*h[hashp],*stp,*p;
inline exgcd(x,y,a,b)int*a,*b,x,y;{
	int aa=*b=0,bb=*a=1,ta,tb;
	if(!x||!y)return -1;
	while(y){
		ta=aa;tb=bb;
		aa=*a-x/y*aa;bb=*b-x/y*bb;
		*a=ta;*b=tb;
		ta=x;x=y;y=ta%y;
	}return x;
}inline getint(){
	int re=0,ch=getchar();
	while(ch<'0'||ch>'9')if((ch=getchar())<=0)return -1;
	for(;ch>='0'&&ch<='9';ch=getchar())re=re*10+ch-48;
	return re;
}inline PHash find(x,t,c)int x,t;_Bool c;{
	int i=(x%hashp+hashp)%hashp;PHash r=h[i];
	for(;r;r=r->next)if(r->v==x)return r;
	if(c)return 0;stp->next=h[i];stp->v=x;stp->t=t;return h[i]=stp++;
}inline gcd(x,y){int t;while(y)t=x,x=y,y=t%y;return x;}
inline solve(a,b,c){
	int i=-1,tmp,r=0,M,x;LL buf=1,D=1,K=1;
	for(stp=st;++i!=50;buf=buf*a%c)if(buf==b)return i;
	while((tmp=gcd(a,c))!=1){
		if(b%tmp)return -1;
		r++;c/=tmp;b/=tmp;D=D*a/tmp%c;
	}M=sqrt(c);memset(h,0,hashp<<2);
	for(buf=1,i=-1;++i!=M;K=K*a%c,buf=buf*a%c)find((int)buf,i,0);
	for(i=-1;(LL)++i*M<=c;D=D*K%c){
		exgcd((int)D,c,&tmp,&x);
		tmp=(LL)tmp*b%c;if(tmp<0)tmp+=c;
		if(p=find(tmp,0,1))return i*M+p->t+r;
	}return -1;
}main(){
	int a,b,c;
	while(a=getint(),c=getint(),b=getint())b>=c||(a=solve(a,b,c))<0?puts("No Solution"):printf("%d\n",a);
	return 0;
}
By Seter

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter
Host by is-Programmer.com | Power by Chito 1.3.3 beta | © 2007 LinuxGem | Design by Matthew "Agent Spork" McGee