大步小步
BSGS算法
普通版
在以前我们学了用拓欧求解线性同余方程 $ax\equiv 1\pmod b$ ，现在我们来思考高次同余方程，考虑如下问题：
求解关于 $x$ 的方程 $a^x\equiv b\pmod p$ ，其中 $gcd(a,p)=1$
不难想到欧拉定理 $a^{\varphi(p)}\equiv 1\pmod p$ ，设 $a^t\equiv a^{x\mod \varphi(p)}\pmod p$ ，易得 $t\in[0,\varphi(p)-1]$ ，我们只需找到 $t$ 使 $x=t$ 时方程成立即可，为了代码简洁（懒得求欧拉函数）不 妨扩大 $t$ 的范围，使 $t\in[0,p]$
直接枚举 $t$ 计算显得过于暴力，BSGS算法使用一个“以空间换时间”的思想，用分块加hash表处理：
令 $k=\lceil\sqrt{p}\rceil$ （代码中为简洁一般是 $k=\sqrt{n}+1$ ），将 $[0,p]$ 每 $k$ 个分为一段，则可以把 $t$ 表示为 $t=kx-y,x\in[1,k],y\in[0,k-1]$ (这里的 $x$ 已经不是上面方程的 $x$ 了，而是新设的未知数，不用 $+y$ 而用 $-y$ 是为了在后面移项后不出现逆元) ，这样， $t$ 可以取遍 $[1,p]$ 中所有值，剩下的 $0$ 特判即可
那么问题转化为求出一组 $x,y$ 使 $t=kx-y,a^t\equiv b\pmod p$ ，很暴力的思路是，枚举每一个 $x$ 判断是否存在 $y$ ，这里就体现出用 $-y$ 的好处，我们来变个形：
$$
a^t\equiv b\pmod p\Rightarrow a^{kx-y}\equiv b\pmod p\Rightarrow a^{kx}\equiv b\times a^y\pmod p
$$
可以提前将 $b\times a^y$ 存到hash表里，枚举 $x$ 判定即可，如果有多个 $y$ 的 $b\times a^y$ 相同，我们保留较大的（这样减了以后 $t$ 最小）
时间复杂度为 $O(\sqrt{p})$
BSGS板子
代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
| #include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50000+5,INF=0x3f3f3f3f;
const int H=20051221;
struct Hash{
struct Edge{
int key,val,ne;
}e[N<<2];
int h[H],idx=0;
void clear(){
for(int i=0;i<H;++i) h[i]=-1;
}
int find(int key){
int _key=(key%H+H)%H;
for(int i=h[_key];i!=-1;i=e[i].ne)
if(e[i].key==key) return e[i].val;
return -1;
}
int change(int key,int val){
int _key=(key%H+H)%H;
for(int i=h[_key];i!=-1;i=e[i].ne)
if(e[i].key==key) return e[i].val=val;
return -1;
}
void insert(int key,int val){
if(change(key,val)!=-1) return ;
int _key=(key%H+H)%H;
e[idx].key=key,e[idx].val=val,e[idx].ne=h[_key],h[_key]=idx++;
}
}ha;
int bsgs(int a,int b,int p){
if(1%p==b%p) return 0;
int k=sqrt(p)+1;
ha.clear();
for(int i=0,j=b%p;i<k;++i){
ha.insert(j,i);
j=(LL)j*a%p;
}
int ak=1;
for(int i=0;i<k;++i) ak=(LL)ak*a%p;
for(int i=1,j=ak,_h;i<=k;++i){
_h=ha.find(j);
if(_h!=-1) return (LL)i*k-_h;
j=(LL)j*ak%p;
}
return -INF;
}
int main(){
int a,p,b;
scanf("%d%d%d",&p,&a,&b);
int ans=bsgs(a,b,p);
ans<0?printf("no solution\n"):printf("%d\n",ans);
return 0;
}
|
拓展版
考虑更一般的情况：
$a^x\equiv b\pmod p$ ，其中 $a,p$ 不一定互质，求解 $x$ ，理所当然想办法转化为普通版，过程如下：
若 $a^0\equiv b\pmod p$ ， $x=0$ ，结束
若 $x\ge1$ 且 $gcd(a,p)==1$ 用普通版BSGS
若 $gcd(a,p)=d,d>1$ ，则 $a^x\equiv b\pmod p\Leftrightarrow a^x+kp=b$
明显，若 $d\nmid b$ 原方程无解
否则有 $d|b$ ，可以所有数同时除以 $d$ ，化为 $\frac{a}{d}\times a^{x-1}\equiv\frac{b}{d}\pmod {\frac{p}{d}}$ 此时 $gcd(\frac{a}{d},\frac{p}{d})=1$ ，可以求 $\frac{a}{b}$ 的逆元，于是化为 $a^{x-1}\equiv\frac{b}{d}\times(\frac{a}{d})^{-1}\pmod {\frac{p}{d}}$ ，让 $a’=a,x’=x-1,b’=\frac{b}{d}\times(\frac{a}{d})^{-1},p’=\frac{p}{d}$ ，再次回到第一步求解 $a’^{x’}\equiv b’\pmod {p’}$
拓展BSGS
本题有多组数据，每次清空hash表还不如每次重开一个map快（unordered_map也行）
代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
| #include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50000+5,INF=0x3f3f3f3f;
int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1,y=0;
return a;
}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
int bsgs(int a,int b,int p){
if(1%p==b%p) return 0;
int k=sqrt(p)+1;
map<int,int>ha;
for(int i=0,j=b%p;i<k;++i){
ha[j]=i;
j=(LL)j*a%p;
}
int ak=1;
for(int i=0;i<k;++i) ak=(LL)ak*a%p;
for(int i=1,j=ak;i<=k;++i){
if(ha.find(j)!=ha.end()) return i*k-ha[j];
j=(LL)j*ak%p;
}
return -INF;
}
int exbsgs(int a,int b,int p){
b=(b%p+p)%p;
if(1%p==b%p) return 0;
int x,y;
int d=exgcd(a,p,x,y);
if(d>1){
if(b%d) return -INF;
exgcd(a/d,p/d,x,y);
return exbsgs(a,(LL)b/d*x%(p/d),p/d)+1;
}
return bsgs(a,b,p);
}
int main(){
int a,p,b;
while(scanf("%d%d%d",&a,&p,&b),a||p||b){
int ans=exbsgs(a,b,p);
ans<0?printf("No Solution\n"):printf("%d\n",ans);
}
return 0;
}
|
ps：据说把那两个 return 0 去掉就可以求次大解
例题
随机数生成器
本题的特判较多，我们慢慢分析
首先是 $a,b$ 可以等于 $0$ ，先特判掉，保证 $a,b>0$
再将问题转化：给定 $p,a,b,x_1,t$ 和递推公式 $x_{i+1}=a\times x_i+b\pmod p$ 求一个最小的 $n$ 使得 $x_n=t$
当然，我们的第一步应该是找到 $x_n$ 的通项公式，我要开始变形了（以下均默认 $\mod p$）：
$$
\begin{align}
&\text{设}c\in\mathbb{R}\text{且使得}x_i+c=a(x_{i-1}+c)\\
\Rightarrow&x_i=ax_{i-1}+c(a-1)\\
&\text{又}\because x_i=ax_{i-1}+b\\
\Rightarrow&c=\frac{b}{a-1}(a\ne1)\\
\Rightarrow&x_i+\frac{b}{a-1}=a(x_{i-1}+\frac{b}{a-1})(a\ne1)\\
&\text{等比数列通项公式得}x_i+\frac{b}{a-1}=a^{i-1}(x_1+\frac{b}{a-1})(a\ne1)\\
\Rightarrow&x_n+\frac{b}{a-1}=a^{n-1}(x_1+\frac{b}{a-1})(a\ne1)\\
&\text{又}\because x_n=t\\
\Rightarrow&t+\frac{b}{a-1}=a^{n-1}(x_1+\frac{b}{a-1})(a\ne1)\\
&\text{特判掉}x_1+\frac{b}{a-1}=0\\
\Rightarrow&a^{n-1}=\frac{t+\frac{b}{a-1}}{x_1+\frac{b}{a-1}}(a\ne1)
\end{align}
$$
特判 $a=1$ 的情况，然后，由于 $p$ 问题就转化为 $a^x\equiv b\pmod p$ 的形式，解出 $x$ 即是 $n-1$
以上，用了三个特判，我们把问题转换为BSGS
代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
| #include<bits/stdc++.h>
#define LL long long
using namespace std;
const int INF=0x3f3f3f3f;
int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1,y=0;
return a;
}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
int inv(int a,int p){
int x,y;
exgcd(a,p,x,y);
return (x+p)%p;
}
int bsgs(int a,int b,int p){
if(1%p==b%p) return 0;
unordered_map<int,int>ha;
int k=sqrt(p)+1;
for(int i=0,j=b%p;i<k;++i){
ha[j]=i;
j=(LL)j*a%p;
}
int ak=1;
for(int i=0;i<k;++i) ak=(LL)ak*a%p;
for(int i=1,j=ak;i<=k;++i){
if(ha.find(j)!=ha.end()) return (LL)i*k-ha[j];
j=(LL)j*ak%p;
}
return -INF;
}
int main(){
int T;
scanf("%d",&T);
int p,a,b,x1,t;
while(T--){
scanf("%d%d%d%d%d",&p,&a,&b,&x1,&t);
if(a==0){
if(x1==t) puts("1");
else if(b==t) puts("2");
else puts("-1");
}
else if(a==1){
if(b==0) t==x1?puts("1"):puts("-1");
else{
int x;
x=((LL)inv(b,p)*(t-x1)%p+p)%p;
printf("%d\n",x+1);
}
}
else{
int C=(LL)b*inv(a-1,p)%p;
int A=(x1+C)%p;
if(A==0){
int u=(-C+p)%p;
t==u?puts("1"):puts("-1");
}
else{
int B=(t+C)%p;
int ans=bsgs(a,(LL)B*inv(A,p)%p,p);
ans<0?puts("-1"):printf("%d\n",ans+1);
}
}
}
return 0;
}
|