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

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

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;
}

• inline
• block