数论相关
1.逆元
引入
我们知道在同余系里面加，减，乘都是可以直接完成。
但是不能除。
我们要想在模 \(p\) 意义下除以一个数，就可以乘这个数的逆元。
这个数记作 \(x^{-1}\).
满足 \(x \cdot x^{-1} \equiv 1 \pmod p\)
线性求 \(1\sim n\) 逆元：
$i^{-1} = - (p \bmod i)^{-1}\cdot \left \lfloor \frac{p}{i} \right \rfloor $
证明：设 \(t=\left \lfloor \frac{p}{i} \right \rfloor\),\(r=p \bmod i\).
则
\(t*i+r\equiv 0 \pmod p\).
同时乘 \(i^{-1}\)
\(t+r*i^{-1} \equiv 0 \pmod p\).
故
\(i^{-1}\equiv -t * r^{-1} \pmod p\)
2.裴蜀定理/扩展欧几里得
引入
裴蜀定理
\(d=\gcd(x,y)\)，则必有 \(ax+by=d\).
推论： \(a\perp b\) 的充要条件是 \(ax+by=d\).
扩展欧几里得算法
用于求解 \(ax+by=d\) 问题。
code
void exgcd(LL a,LL b,LL &x,LL &y) {
if(b==0) {x=1; y=0; return ;}
exgcd(b,a%b,x,y);
int z=x; x=y; y=z-(a/b)*y;
}
这求出的是一组特解，设它们是 \(x_0,y_0\).
通解为：
\(x=x_0+k\cdot \frac{y}{d}\).
\(y=y_0+k\cdot \frac{x}{d}\).
$k \in \mathbb{Z} $.
对于任意形式的不定方程 \(ax+by=c\)，只要把 \(x,y\) 同乘 \(\frac{c}{d}\) 即可。
若 \(d \not \mid c\)，则无解。
扩欧还可以求解逆元。
只要求 \(x^{-1}*x+k*p=1\) 即可。
其中 \(x,p\) 为已知项，需求 \(x^{-1},k\).
应用
P3986 斐波那契数列
设 \(f\) 为原本斐波那契数列， \(g\) 为题目的数列。
考虑矩阵乘法：
我们得到一个重要性质：
\(g_i = a\cdot f_{i-2} + b\cdot f_{i-1}\)
于是我们只要解 \(k=a\cdot f_{i-2} + b\cdot f_{i-1}\) 中 \((a,b)\) 的个数即可。
枚举 \(f\) ，由于斐波那契呈指数增长，所以项数是 \(\log\) 级别。
3.欧拉定理
引入
欧拉定理
\(a^{\varphi(p)}\equiv 1 \pmod p\)
需要满足 \(a\perp p\).
即 \(a^b\equiv a^{b \mod \varphi(p)}\pmod p\)
说明 \(a^b \mod p\) 是以 \(\varphi(p)\) 为循环的。（不一定是最小循环节）
证明：link
欧拉定理可以求逆元。 \(x^{\varphi(p)-1}\equiv x^{-1} \pmod p\) .
扩展欧拉定理
若 \(b> \varphi(p)\) 则，\(a^b\equiv a^{b\mod \varphi(p)+\varphi(p)}\pmod p\)
否则 \(a^b\equiv a^b\pmod p\)
对任意 \(a,p\) 都满足。
即 \(a^b \mod p\) 是经过 \(\varphi(p)\) 才进入循环的。
应用
P4139 上帝与集合的正确用法
求 \(2^{2^{2...}} \mod p\)
设上面这个式子为 \(x\).
则 \(x=2^{(2^{2...} \mod \varphi(p)+\varphi(p))}\)
在设 \(y=2^{2...} \mod \varphi(p)\)
我们再求 \(y\).
以此类推，我们最终发现，模数会为 1.
那后面就不用算了。
可以证明 \(p=\varphi(p)\)，最多只会更新 \(\log p\) 次。
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e7+10;
int t,p,phi[N],prime[N],tot;
bool v[N];
void init() {
phi[1]=1;
for(int i=2; i<N; i++) {
if(!v[i]) prime[++tot]=i,phi[i]=i-1;
for(int j=1; j<=tot&&i*prime[j]<N; j++) {
v[prime[j]*i]=1;
if(i%prime[j]==0) {phi[i*prime[j]]=phi[i]*prime[j]; break;}
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
}
int power(int a,int b,int p) {
int res=1;
for(; b; b>>=1) {
if(b&1) res=1ll*res*a%p;
a=1ll*a*a%p;
}
return res;
}
int solve(int p) {
if(p==1) return 0;
return power(2,solve(phi[p])+phi[p],p);
}
int main() {
init();
scanf("%d",&t);
for(; t; t--) scanf("%d",&p),printf("%d\n",solve(p));
return 0;
}
P3934 [Ynoi2016] 炸脖龙 I
求区间 \([l,r]\)，中， \(a_l^{a_{l+1}^{...a_r}}\)
带修改。
考虑扩展欧拉定理。
跟上题类似。
迭代到 \(p=\varphi(p)\)， \(p=1\) 时即可。
注意我们不能直接套用 \(a^b\equiv a^{b\mod \varphi(p)+\varphi(p)}\pmod p\).
因为 \(b\) 可能不大于 \(p\).
所以加个特判。
由于我们只用单点查询，所以修改用一个差分树状数组即可。
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
const LL N=5e5+5;
const LL M=2e7+5;
LL n,m;
LL a[N],c[N];
LL phi[M],prime[N*10],tot;
LL st[N];
bool isp[M],flag;
void add(LL p,LL x){
for(; p<=n; p+=p&-p) c[p]+=x;
}
LL ask(LL p) {
LL res=0;
for(; p; p-=p&-p) res+=c[p];
return res;
}
void init() {
memset(isp,true,sizeof(isp));
phi[1]=1;
for(LL i=2; i<M; i++){
if(isp[i]) prime[++tot]=i,phi[i]=i-1;
for(LL j=1; j<=tot&&i*prime[j]<M; j++){
isp[i*prime[j]]=false;
if(i%prime[j]==0) {
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
}
LL Qpow(LL a,LL b,LL mod) {
LL res=1;
flag=false;
if(a>=mod) flag=true,a%=mod;
for(; b; b>>=1) {
if(b&1) {
res=res*a;
if(res>=mod) flag=true,res%=mod;
}
a=a*a;
if(a>=mod) flag=true,a%=mod;
}
return res;
}
signed main() {
init();
scanf("%lld%lld",&n,&m);
for(LL i=1; i<=n; i++) scanf("%lld",&a[i]);
for(LL i=1; i<=n; i++) add(i,a[i]-a[i-1]);
for(; m; m--) {
LL opt,x,y,p;
scanf("%lld%lld%lld%lld",&opt,&x,&y,&p);
if(opt==1) {
add(x,p);
add(y+1,-p);
} else {
LL now=x;
st[x]=p;
p=phi[p];
while(p>1&&now<y) {
st[++now]=p;
p=phi[p];
}
LL res=1;
for(LL i=now; i>=x; i--) {
res=Qpow(ask(i),res,st[i]);
if(flag) res+=st[i];
}
printf("%lld\n",res%st[x]);
}
}
return 0;
}
P3747 [六省联考 2017] 相逢是问候
套一个小清新线段树。
需要用到光速乘。
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
const int N=5e4+5,M=2e4;
int n,m,P,c,mod[60],pw[60][M+5],pwr[60][M+5],cnt;
int phi(int n) {
int ret=n;
for(int i=2; i*i<=n; i++){
if(n%i) continue;
ret=1ll*ret*(i-1)/i;
while(n%i==0) n/=i;
}
if(n>1) ret=1ll*ret*(n-1)/n;
return ret;
}
int Mod(LL x,int p) {
return x>mod[p]?(x%mod[p]+mod[p]):x;
}
int Fpow(int b,int i) {
return Mod(1ll*pw[i][b/M]*pwr[i][b%M],i);
}
int calc(int x,int y,int p) {
if(!y) return Mod(x,p);
if(p==cnt) return 1;
return Mod(Fpow(calc(x,y-1,p+1),p),p);
}
void init() {
mod[0]=P;
while(mod[cnt]>1) {
cnt++; mod[cnt]=phi(mod[cnt-1]);
}
mod[++cnt]=1;
for(int i=0; i<=cnt; i++) {
pw[i][0]=pwr[i][0]=1;
for(int j=1; j<=M; j++) pwr[i][j]=Mod(1ll*pwr[i][j-1]*c,i);
for(int j=1; j<=M; j++) pw[i][j]=Mod(1ll*pw[i][j-1]*pwr[i][M],i);
}
}
int a[N][60],dat[N<<2],t[N<<2];
void pushup(int p) {
t[p]=min(t[p<<1],t[p<<1|1]);
dat[p]=(dat[p<<1]+dat[p<<1|1])%P;
}
void build(int p,int l,int r) {
if(l==r) return dat[p]=a[l][0],void();
int mid=(l+r)>>1;
build(p<<1,l,mid);
build(p<<1|1,mid+1,r);
pushup(p);
}
void modify(int p,int l,int r,int L,int R) {
if(t[p]>=cnt) return;
if(l==r) return dat[p]=a[l][++t[p]],void();
int mid=(l+r)>>1;
if(L<=mid) modify(p<<1,l,mid,L,R);
if(R>mid) modify(p<<1|1,mid+1,r,L,R);
pushup(p);
}
int query(int p,int l,int r,int L,int R) {
if(L<=l&&r<=R) return dat[p];
int mid=(l+r)>>1,ret=0;
if(L<=mid) ret=(ret+query(p<<1,l,mid,L,R))%P;
if(R>mid) ret=(ret+query(p<<1|1,mid+1,r,L,R))%P;
return ret;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
cin>>n>>m>>P>>c;
init();
for(int i=1; i<=n; i++) {
cin>>a[i][0];
for(int j=1; j<=cnt; j++) a[i][j]=calc(a[i][0],j,0)%P;
}
build(1,1,n);
for(int i=1,op,l,r; i<=m; i++){
cin>>op>>l>>r;
if(op) printf("%d\n",query(1,1,n,l,r));
else modify(1,1,n,l,r);
}
return 0;
}
4.Lucas 定理
引入
我们求解组合数时，如果模数较小，我们就比较难求。
因为我们要用阶乘，而超过模数的阶乘会变成 0.
Lucas 定理： \(C_{n}^{m}=C_{\left \lfloor n/p \right \rfloor }^{\left \lfloor m/p \right \rfloor} \cdot C_{n \bmod p}^{m \bmod p}\)
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
const int N=1e5+10;
LL n,m,p,t,inv[N],mul[N],mul2[N];
LL C(LL a,LL b) {
if(a<b) return 0;
return mul[a]*mul2[a-b]%p*mul2[b]%p;
}
LL lucas(LL a,LL b) {
if(a<p&&b<p) return C(a,b);
return C(a%p,b%p)*lucas(a/p,b/p)%p;
}
signed main() {
scanf("%lld",&t);
for(; t; t--) {
scanf("%lld%lld%lld",&n,&m,&p);
inv[0]=mul[0]=mul2[0]=1;
inv[1]=mul[1]=mul2[1]=1;
for(int i=2; i<p; i++) inv[i]=((-(p/i)*inv[p%i])%p+p)%p;
for(int i=2; i<p; i++) mul[i]=mul[i-1]*i%p;
for(int i=2; i<p; i++) mul2[i]=mul2[i-1]*inv[i]%p;
printf("%lld\n",lucas(n+m,n));
}
return 0;
}
注意特判：递归过程中可能会出现 \(n<m\) 返回 0，那么整个式子就是 0.
证明：link
应用
P4345 [SHOI2015]超能粒子炮·改
设 \(f(n,k)\) 为 \(\sum ^k_{i=0} C_n^i\)。
求 \(f(n,k)\bmod p,p=2333\).
我们尝试用 Lucas 定理形式拆解。
\(C_{n}^{m}=C_{n/p}^{m/p} \cdot C_{n \bmod p}^{m \bmod p}\)
\(f(n,k) = C_{n/p}^{0/p}*C_{n \bmod p}^{0\bmod p}+C_{n/p}^{1/p}*C_{n \bmod p}^{1\bmod p}+....\)
我们发现可以按 \(C_{n/p}^{i/p}\) 分组.
则 \(f(n,k)= \sum_j^{k/p-1} (C_{n/p}^j * f(n\bmod p,p-1)) + C_{n/p}^{k/p}*f(n\bmod p,k\bmod p)\)
\(=f(n/p,k/p-1)*f(n\bmod p,p-1) + C_{n/p}^{k/p}*f(n\bmod p,k\bmod p)\)
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
const LL p=2333;
LL t,n,k;
LL f[p+10][p+10],c[p+10][p+10];
LL C(LL n,LL m) {
if(n<p&&m<p) return c[n][m];
return C(n%p,m%p)*C(n/p,m/p)%p;
}
LL solve(LL n,LL k) {
if(!n) return 1;
if(!k) return 1;
if(k<0) return 0;
if(n<p&&k<p) return f[n][k];
return (solve(n%p,p-1)*solve(n/p,k/p-1)+C(n/p,k/p)*solve(n%p,k%p))%p;
}
signed main() {
scanf("%lld",&t);
for(int i=0; i<p; i++) f[i][0]=c[i][0]=1;
for(int i=1; i<p; i++)
for(int j=1; j<p; j++)
c[i][j]=(c[i-1][j]+c[i-1][j-1])%p,f[i][j]=(f[i][j-1]+c[i][j])%p;
for(; t; t--) {
scanf("%lld%lld",&n,&k);
printf("%lld\n",solve(n,k));
}
return 0;
}
5.中国剩余定理
引入
中国剩余定理
可以求解同余方程组：
\( \begin{cases} x \equiv a_1 \pmod {p_1} \\ x \equiv a_2 \pmod {p_2} \\ x \equiv a_3 \pmod {p_3} \\ \end{cases} \)
要求 \(p\) 两两互质。
设 \(M=\prod p,M_i=M/p_i\)
求出 \(t_i=M_i^{-1} \pmod {p_i}\)
\(ans= \sum a_i t_i M_i\)
扩展中国剩余定理：
\(p\) 不互质也可以。
考虑一一合并所有问题。
已知
\( \begin{cases} x \equiv a_1 \pmod {p_1} \\ x \equiv a_2 \pmod {p_2} \\ \end{cases} \)
则 \(x=k_1p_1 + a_1=k_2p_2 + a_2\).
移项 \(k_1p_1-k_2p_2 = -a_1+a_2\).
那么我们发现这是经典不定方程形式，扩欧求解出 \(k_1,k_2\)，或判断无解。
若有解，利用 \(k_1,k_2\) 求出 \(x^*\) 为一个特解。
则 \(x \equiv x^* \pmod {lcm (p_1,p_2)}\)
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
LL n,p1,c1,p2,c2;
LL mul(LL a,LL b,LL mod) {
LL res=0;
for(; b; b>>=1) {
if(b&1) res=(res+a)%mod;
a=(a+a)%mod;
}
return res;
}
LL gcd(LL a,LL b) {
if(b==0) return a;
return gcd(b,a%b);
}
void exgcd(LL a,LL b,LL &x,LL &y) {
if(b==0) {x=1; y=0; return ;}
exgcd(b,a%b,x,y);
int z=x; x=y; y=z-(a/b)*y;
}
signed main() {
ios::sync_with_stdio(false);
cin>>n;
c1=0; p1=1;
for(int i=1; i<=n; i++) {
cin>>p2>>c2;
LL d=gcd(p1,p2),lc=p1/d*p2,k,y;
c2=(c2-c1%p2+p2)%p2;
exgcd(p1,p2,k,y);
k=(k%(p2/d)+p2/d)%(p2/d);
k=mul(k,c2/d,p2);
c1=(c1+mul(p1,k,lc))%lc;
p1=lc;
}
cout<<c1<<endl;
return 0;
}
应用
P2480 [SDOI2010]古代猪文
求 \(g^{\sum d|n C_n^d} \bmod 999911659\).
考虑欧拉定理：由于 999911659 是质数，所以求 \(\sum d|n C_n^d \bmod 999911658\).
我们直接 Lucas 会爆炸，所以观察发现 \(999911658 = 2*3*4697*35617\)
所以我们对这四个质数一一做一遍 Lucas ，再用中国剩余定理合并即可。
code
#include<bits/stdc++.h>
using namespace std;
const int N=36666;
using LL=long long;
LL p[4]={2,3,4679,35617},c[4];
LL inv[N],mul1[N],mul2[N];
LL n,g;
LL M[4],M0=999911658,t[4],ans;
LL mod=999911659;
vector<LL> s;
LL power(LL a,LL b,LL p) {
LL res=1;
for(; b; b>>=1) {
if(b&1) res=res*a%p;
a=a*a%p;
}
return res;
}
LL C(LL n,LL m,LL p) {
if(n<m) return 0;
return mul1[n]*mul2[n-m]%p*mul2[m]%p;
}
LL Lucas(LL n,LL m,LL p) {
if(n<p&&m<p) return C(n,m,p);
return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}
signed main() {
scanf("%lld%lld",&n,&g);
if(g%mod==0) return puts("0"),0;
for(int d=1; d*d<=n; d++) {
if(n%d==0) {
s.push_back(d);
if(n/d!=d) s.push_back(n/d);
}
}
for(int k=0; k<4; k++) {
inv[0]=mul1[0]=mul2[0]=1;
inv[1]=mul1[1]=mul2[1]=1;
for(int i=2; i<p[k]; i++) inv[i]=-(p[k]/i)*inv[p[k]%i]%p[k];
for(int i=2; i<p[k]; i++) mul1[i]=mul1[i-1]*i%p[k];
for(int i=2; i<p[k]; i++) mul2[i]=mul2[i-1]*inv[i]%p[k];
for(int i=0; i<(int)s.size(); i++) {
c[k]=(c[k]+Lucas(n,s[i],p[k]))%p[k];
}
c[k]=(c[k]%p[k]+p[k])%p[k];
M[k]=M0/p[k];
t[k]=power(M[k],p[k]-2,p[k]);
ans=ans+M[k]*t[k]*c[k];
}
ans=(ans%M0+M0)%M0;
printf("%lld\n",power(g,ans,mod));
return 0;
}
P4774 [NOI2018] 屠龙勇士
题意复杂，难以概括。
先用平衡树维护出每次攻击的武器。
然后跑一个扩展中国剩余定理即可。
记得用龟速乘或“骆可强”快速乘防止溢出。
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
const int N=1e5+10;
multiset<LL> st;
int T,n,m;
LL s[N],p[N],b[N],z[N],a[N];
LL Mul(LL a,LL b,LL p) {
LL res=0;
for(; b; b>>=1) {
if(b&1) res=(res+a)%p;
a=(a+a)%p;
}
return res;
}
LL gcd(LL a,LL b) {
if(b==0) return a;
return gcd(b,a%b);
}
void exgcd(LL a,LL b,LL &x,LL &y) {
if(b==0) {x=1; y=0; return ;}
exgcd(b,a%b,x,y);
LL z=x; x=y; y=z-(a/b)*y;
}
void solve() {
st.clear();
cin>>n>>m;
for(int i=1; i<=n; i++) cin>>s[i];
for(int i=1; i<=n; i++) cin>>p[i];
for(int i=1; i<=n; i++) cin>>b[i];
for(int i=1; i<=m; i++) cin>>z[i],st.insert(z[i]);
for(int i=1; i<=n; i++) {
auto it=st.upper_bound(s[i]);
if(it==st.begin()) a[i]=*st.begin(),st.erase(st.begin());
else a[i]=*(--it),st.erase(it);
st.insert(b[i]);
}
LL p0=1,c0=0; // x = p0*k + c0
for(int i=1; i<=n; i++) { // a*x = s (mod p)
LL fa=a[i],fp=p[i],fs=s[i];
LL d=gcd(fa,fp),x,y;
if(fs%d!=0) return printf("-1\n"),void();
fs/=d;
exgcd(fa/=d,fp/=d,x,y);
x=(x%fp+fp)%fp;
fs=Mul(fs,x,fp); // x = fs (mod fp)
// x = p0 * k + c0= fp * y + fs -> p0*k-fp*y=fs-c0
d=gcd(p0,fp);
if((fs-c0)%d!=0) return printf("-1\n"),void();
exgcd(p0,fp,x,y);
x=x*((fs-c0)/d); y=y*((fs-c0)/d);
LL lc=p0/d*fp;
c0=(Mul(x,p0,lc)+c0)%lc;
p0=lc;
}
LL idx=1;
for(int i=2; i<=n; i++) {
if(ceil((double)s[i]/a[i])>ceil((double)s[idx]/a[idx]))
idx=i;
}
while(c0<ceil((double)s[idx]/a[idx])) c0+=p0;
printf("%lld\n",c0);
}
int main() {
ios::sync_with_stdio(false);
cin>>T;
for(; T; T--) solve();
return 0;
}
6.BSGS
引入
BSGS:
求解高次同余方程 \(a^x \equiv b \pmod p (a \perp p)\).
注意到 \(x<p\).
设 $t=\left \lceil \sqrt p \right \rceil $
我们知道任意 \(x<p\)，都可以表示为 \(x=it-j (i,j\le t)\)。
则 \(a^x\equiv a^{it-j}\equiv b \pmod p\).
那么 \(a^{it} \equiv b * a^j \pmod p\)
预处理出所有 \(b * a^j\)，将它们存在哈希表里，或 STL map 里。
再处理出所有 \(a^{it}\)，判断是否存在有 \(b*a^j = a^{it}\).
答案为 \(it-j\).
code
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
LL p,a,b;
unordered_map<LL,LL> hsh;
LL Qpow(LL a,LL b,LL p) {
LL res=1;
for(; b; b>>=1) {
if(b&1) res=res*a%p;
a=a*a%p;
}
return res;
}
LL bsgs(LL p,LL a,LL b) {
hsh.clear();
b%=p;
LL t=sqrt(p)+1;
for(int i=0; i<t; i++) {
LL s=b*Qpow(a,i,p)%p;
hsh[s]=i;
}
a=Qpow(a,t,p);
for(int i=1; i<=t; i++) {
LL val=Qpow(a,i,p);
int j=hsh.find(val)==hsh.end()?-1:hsh[val];
if(j>=0&&i*t-j>=0) return i*t-j;
}
return -1;
}
signed main() {
scanf("%lld%lld%lld",&p,&a,&b);
LL ans=bsgs(p,a,b);
if(ans==-1) puts("no solution");
else printf("%lld\n",ans);
return 0;
}
笔者发现， BSGS 似乎与 meet in middle 些许相似。
扩展 BSGS:
处理 \(a\not \perp p\) 情况.
\(p,b\) 不断除以 \(gcd(a,p)\)，直到 \(a\perp p\).
若期间 \(b\) 不能整除，则无解。
设总共除了 \(d\)，除了 \(cnt\) 次.
则转换为 \(a^x * a^{cnt} = \frac{b}{d} \pmod {\frac{p}{d}}\)
答案为 \(x+cnt\).
code
#include<bits/stdc++.h>
using namespace std;
int p,a,b;
map<int,int> hsh;
int Qpow(int a,int b,int p) {
int res=1;
for(; b; b>>=1) {
if(b&1) res=1ll*res*a%p;
a=1ll*a*a%p;
}
return res;
}
int gcd(int a,int b) {
if(b==0) return a;
return gcd(b,a%b);
}
int bsgs(int a,int b,int p,int ad) {
hsh.clear();
int t=sqrt(p)+1,s=1;
for(int i=0; i<t; i++,s=1ll*s*a%p) hsh[1ll*s*b%p]=i;
int tmp=s; s=ad;
for(int i=0; i<=t; i++,s=1ll*s*tmp%p)
if(hsh.find(s)!=hsh.end())
if(i*t-hsh[s]>=0) return i*t-hsh[s];
return -1;
}
int exbsgs(int a,int b,int p) {
a%=p; b%=p;
if(b==1||p==1) return 0;
int cnt=0,d,ad=1;
for(; (d=gcd(a,p))!=1; ) {
if(b%d) return -1;
cnt++; b/=d; p/=d;
ad=(1ll*ad*a/d)%p;
if(ad==b) return cnt;
}
int ans=bsgs(a,b,p,ad);
if(ans==-1) return -1;
return ans+cnt;
}
signed main() {
for(; ; ) {
scanf("%d%d%d",&a,&p,&b);
if(!a&&!b&&!p) break;
int ans=exbsgs(a,b,p);
if(ans==-1) puts("No Solution");
else printf("%d\n",ans);
}
return 0;
}
应用
P3306 [SDOI2013] 随机数生成器
矩阵乘法的 BSGS.
code
#include<bits/stdc++.h>
using namespace std;
int p,a,b,x1,u;
struct Mtr {
int n,m;
int a[3][3];
bool operator < (const Mtr Q) const {
return a[1][1]<Q.a[1][1];
}
} X,F,S,I;
map<Mtr,int> hsh;
Mtr Mul(Mtr A,Mtr B) {
Mtr C; C.n=A.n; C.m=B.m;
for(int i=1; i<=C.n; i++)
for(int j=1; j<=C.m; j++)
C.a[i][j]=0;
for(int i=1; i<=A.n; i++)
for(int j=1; j<=B.m; j++)
for(int k=1; k<=A.m; k++)
C.a[i][j]=(C.a[i][j]+1ll*A.a[i][k]*B.a[k][j])%p;
return C;
}
void solve() {
hsh.clear();
scanf("%d%d%d%d%d",&p,&a,&b,&x1,&u);
if(x1==u) {
printf("1\n");
return ;
}
if(a==0) {
if(x1==u) printf("1\n");
else if(b==u) printf("2\n");
else printf("-1\n");
return ;
}
if(a==1&&b==0) {
if(x1==u) printf("1\n");
else printf("-1\n");
return ;
}
I.n=I.m=2; I.a[1][1]=I.a[2][2]=1;
X.n=X.m=2; X.a[1][1]=a; X.a[1][2]=b; X.a[2][1]=0; X.a[2][2]=1;
F.n=2; F.m=1; F.a[1][1]=x1; F.a[2][1]=1;
S.n=2; S.m=1; S.a[1][1]=u; S.a[2][1]=1;
int t=sqrt(p)+1;
Mtr G=I;
for(int i=0; i<t; i++) {
Mtr J=Mul(G,S);
hsh[J]=i;
G=Mul(G,X);
}
Mtr T=I;
for(int i=0; i<=t; i++) {
Mtr J=Mul(T,F);
if(hsh.find(J)!=hsh.end()) {
if(-hsh[J]+i*t>=0) {
printf("%d\n",-hsh[J]+i*t+1);
return ;
}
}
T=Mul(T,G);
}
printf("-1\n");
return ;
}
int main() {
int T; scanf("%d",&T);
for(; T; T--) solve();
return 0;
}
