离散对数的应用
Ew_Cors
·
·
算法·理论
前言
这里主要是讲离散对数在题目中的应用，一些基础的部分可以去前面的日报或者我的另一篇专栏上看。
请确保在阅读前已经了解了一些离散对数的相关知识。
应用方法
首先忽略简单板子题和不会的板子题。
离散对数和普通对数一样，都是一种把乘法转化为加法的运算。
那么我们可以根据这一点，将题目要求的乘法或者幂次运算转化为更易求的加法或乘法运算，从而使用诸如矩阵快速幂，卷积以及其他数据结构或者算法进行维护和查询。
如果题目中所有幂次操作的底数都是固定的，我们也可以对这些底数求以原根为底的离散对数，通过预处理原根的幂而做到 O(1) 光速幂，从而优化时间复杂度。
下面的三道练手题都是基于上面的两点而使用离散对数的。
写的代码马蜂极其丑陋，除非实在搞不懂之外，大家不要看！！！
应用
Problem：Codeforces-1106F
给你一个长为 k 的数组 b 和一个递推式（i>k）：
f_i=\left(\prod_{j = 1}^{k}f_{i-j}^{b_j}\right)\bmod p
已知 f_1=f_2=\dots=f_{k-1}=1，求出任意一个满足 f_n=m 的 f_k 的值，或输出无解。
---
Solution：
对于 $f$ 的每一项，只和前面的 $k$ 项有关，这让我们想起矩阵加速递推。
但是题目的递推式是连乘形式，怎样把它转化成 $\sum ab$ 的形式呢？
离散对数！
对 $f$ 序列取离散对数，那么 $f_i=\sum_{j = 1}^{k}b_jf_{i-j}$。
直接使用矩阵快速幂：
令矩阵 $B=\begin{bmatrix}b_1&b_2&\cdots&b_{k-1}&b_k\\1&0&\cdots&0&0\\0&1&\cdots&0&0\\\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&\cdots&1&0\end{bmatrix}$。
则 $\log_gm\equiv f_n\equiv\sum_{i=1}^kB^{n-k}_{1,i}f_{k-i+1}\equiv B^{n-k}_{1,1}f_k\pmod{\varphi(p)}$。
扩欧求解即可。
时间复杂度 $O(k^3\log n)$。
---
Code：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int Mod=998244353,Phi=Mod-1,sqrtP=31597,G=3;
int qp(int a,int b,int mod){int ans=1;while(b){if(b&1)ans=1LL*ans*a%mod;a=1LL*a*a%mod;b>>=1;}return ans;}
namespace BSGS{
const int P=2e6+3;
int head[P+5],key[P+5],val[P+5],nxt[P+5],cnt;
void insert(int x,int id){
int hs=x%P;
key[++cnt]=x;val[cnt]=id;nxt[cnt]=head[hs];head[hs]=cnt;
}
int query(int x){
for(int i=head[x%P];i;i=nxt[i])if(key[i]==x)return val[i];
return -1;
}
void init(){
long long now=1,bnow=1;int lim=1e3;
for(int i=0;i<lim;i++)now=now*G%Mod;
for(int i=1;i<=Mod/lim+1;i++){
bnow=bnow*now%Mod;insert(bnow,i*lim);
}
}
int dlog(long long x){
for(int i=0;;x=x*G%Mod,i++){
int ret=query(x);if(~ret)return ret-i;
}
return -1;
}
}
int k,n,m,b[105];
struct Matrix{int c[105][105];};
Matrix operator *(Matrix a,Matrix b){
Matrix c;memset(c.c,0,sizeof c.c);
for(int i=1;i<=k;i++)for(int j=1;j<=k;j++)for(int l=1;l<=k;l++)
c.c[i][j]=(c.c[i][j]+1LL*a.c[i][l]*b.c[l][j]%Phi)%Phi;
return c;
}
Matrix ans,bs;
void exgcd(long long a,long long b,long long &x,long long &y){
if(b==0)return x=1,y=0,void();
exgcd(b,a%b,x,y);long long t=x;x=y;y=t-y*(a/b);
}
map<int,int>mm;
int main(){
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
BSGS::init();
cin>>k;for(int i=1;i<=k;i++)cin>>b[i];
cin>>n>>m;m=(BSGS::dlog(m)%Phi+Phi)%Phi;
for(int i=1;i<=k;i++)bs.c[1][i]=b[i],ans.c[i][i]=1;
for(int i=2;i<=k;i++)bs.c[i][i-1]=1;
int x=n-k;
while(x){if(x&1)ans=ans*bs;bs=bs*bs;x>>=1;}
int sum=ans.c[1][1];
mm[0]=0;
int now=1LL*sum*sqrtP%Phi,dnow=1;
for(int i=1;i<=Phi/sqrtP+1;i++){dnow=(dnow+now)%Phi;mm[(dnow-m+Phi)%Phi]=i*sqrtP;}
now=1;
for(int i=0;i<sqrtP;i++,now=(now+sum)%Phi)if(mm.count(now)){cout<<qp(G,(mm[now]-i+Phi)%Phi,Mod);return 0;}
cout<<"-1";
return 0;
}
```
---
Problem：[HDU-7352](https://acm.hdu.edu.cn/showproblem.php?pid=7352)
时限 11.5s。
给你一个长为 $n$ 的数组 $a$。$q$ 组询问，每次询问一段区间 $l,r$，令 $b_l,b_{l+1},\dots,b_r$ 为 $a_l,a_{l+1},\dots,a_r$ 排序后的值，你需要求出：
$$\prod_{i=l}^{r-1}b_i^{b_{i+1}}\mod p$$
$1\leq n\leq 10^6,1\leq a_i\leq 10^9,1\leq q\leq 5\times 10^3,p=10^9+7$。
---
Solution：
每次询问都要排序，因此我会的普通的数据结构没办法高效做到这一点。
观察到当 $l,r$ 移动时询问较好维护，于是考虑莫队。
使用普通莫队的话，复杂度需要带上查询前驱后继的 $\log$，于是使用回滚莫队和链表，这样就可以做到 $O(\log p)$ 移动左右指针了。这里的 $\log$ 是快速幂。
现在的时间复杂度是 $O(n\sqrt m\log p)$，算一下 $2\times 10^9$。
有点危险，~~但感觉能过，自信~~测极限数据，60s。
这不是常数能解决的问题，我们只能对复杂度进行优化。
考虑对 $a$ 数组做一遍离散对数，这样更新的时候复杂度就是 $O(1)$ 的了。
但是有两个问题：
> 1. 预处理是 $O(n\sqrt p)$ 的，如何优化？
对于 $x<\left\lceil\sqrt p\right\rceil$ 的离散对数，我们可以通过离散对数是完全积性函数这个性质进行线性筛，复杂度 $O(\pi(\sqrt p)\times \sqrt p)$。
其中 $\pi(p)$ 为 $1\sim p$ 内质数个数，$\pi(p)\approx\dfrac p{\ln p}$，于是化简后时间复杂度 $O(\dfrac p{\sqrt{\ln p}})$。
而对于 $x\geq \left\lceil\sqrt p\right\rceil$ 的离散对数，我们对 $p$ 做 $x$ 的带余除法：$p=qx+r$。
$$qx+r\equiv 0\pmod p$$
$$qx\equiv -r\pmod p$$
$$\log_g qx\equiv \log_g -r\pmod{\varphi(p)}$$
$$\log_gq+\log_gx\equiv \log_g-1+\log_gr\pmod{\varphi(p)}$$
注意到：$\log_g-1\equiv\dfrac{\log_g1}2\equiv\dfrac{\varphi(p)}2\pmod{\varphi(p)}$，
$$\log_g x\equiv\log_gr+\dfrac{\varphi(p)}2-\log_gq\pmod{\varphi(p)}$$
由 $q,r$ 的定义可以知道：$q<\left\lceil\sqrt p\right\rceil,r<x$。因此我们可以往下递归。但是我们没有办法保证 $r$ 规模的缩减速度，这是假的。
观察到：$p=qx+r=(q+1)x+r-x$，因此也有：
$$\log_gx\equiv\log_g(x-r)-\log_g(q+1)\pmod{\varphi(p)}$$
而 $\min(r,x-r)\leq \dfrac x2$，因此可以直接递归了，单次复杂度 $O(\log x)$。
预处理部分从而被我们优化至 $O(\dfrac p{\sqrt{\ln p}}+n\log p)$。
> 2. 你求完离散对数之后不还是要快速幂吗，这 $\log$ 还是不能省啊。
有两种解决办法：
第一种：我们在莫队时，只求离散对数，在得到答案时再进行快速幂，莫队复杂度 $O(n\sqrt m+m\log p)$。
第二种：对于 $0<x\leq \sqrt p$，我们预处理 $g^i$ 和 $g^{i\sqrt p}$ 次方，求幂时直接查表即可，预处理复杂度多一个根号。
其中第二种理论复杂度更优。~~虽然没多大区别就是了~~
使用以上算法，可以在 $O(\dfrac p{\sqrt{\ln p}}+n(\sqrt m+\log p))$ 的时间复杂度内通过此题。
---
Code：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int N=1e6,Mod=1e9+7,Phi=Mod-1,G=5,sqrtP=31624,V=1e9;
namespace BSGS{
const int P=2e6+3;
int head[P+5],key[N+5],val[N+5],nxt[N+5],cnt;
void insert(int x,int id){
int hs=x%P;
key[++cnt]=x;val[cnt]=id;nxt[cnt]=head[hs];head[hs]=cnt;
}
int query(int x){
for(int i=head[x%P];i;i=nxt[i])if(key[i]==x)return val[i];
return -1;
}
void init(){
long long now=1,bnow=1;int lim=1e3;
for(int i=0;i<lim;i++)now=now*G%Mod;
for(int i=1;i<=V/lim+1;i++){
bnow=bnow*now%Mod;insert(bnow,i*lim);
}
}
int dlog(long long x){
for(int i=0;;x=x*G%Mod,i++){
int ret=query(x);if(~ret)return ret-i;
}
return -1;
}
}
namespace Gexp{
int bG[1<<15|5],gG[1<<15|5];
void init(){
bG[0]=gG[0]=1;
for(int i=1;i<=(1<<15);i++)bG[i]=1LL*bG[i-1]*G%Mod;
for(int i=1;i<=(1<<15);i++)gG[i]=1LL*gG[i-1]*bG[1<<15]%Mod;
}
long long getexp(int x){return 1LL*bG[x&((1<<15)-1)]*gG[x>>15]%Mod;}
}
const int dlg_1=500000003;
int dlg[sqrtP+5],pri[sqrtP+5],pcnt;
bool isprime[sqrtP+5];
void dlog_prepare(){
BSGS::init();
dlg[1]=0;
for(int i=2;i<=sqrtP;i++){
if(!isprime[i])dlg[i]=BSGS::dlog(i),pri[++pcnt]=i;
for(int j=1;j<=pcnt&&pri[j]*i<=sqrtP;j++){
isprime[i*pri[j]]=1,dlg[i*pri[j]]=(dlg[i]+dlg[pri[j]])%Phi;
if(i%pri[j]==0)break;
}
}
}
long long get_dlog(int x){
if(x<=sqrtP)return dlg[x];
int q=Mod/x,r=Mod-q*x;
return r<x-r?get_dlog(r)+dlg_1-dlg[q]:get_dlog(x-r)-dlg[q+1];
}
int n,m,a[N+5];
int loga[N+5];
void get_loga(){
dlog_prepare();
for(int i=1;i<=n;i++)
loga[i]=(get_dlog(a[i])%Phi+Phi)%Phi;
}
struct Moalgo_NoAdd{int l,r,id;}q[5005];int B;
bool operator <(Moalgo_NoAdd a,Moalgo_NoAdd b){return a.l/B==b.l/B?a.r>b.r:a.l>b.l;}
vector<pair<int,int>>lne;
int pre[N+5],nxt[N+5];
int now=1;
void build(int l,int r){
now=1;
for(int i=l;i<=r;i++)lne.push_back(make_pair(a[i],i));
auto it=lne.end()-(r-l+1);sort(it,lne.end());
inplace_merge(lne.begin(),it,lne.end());
int p=0;
for(auto i:lne){pre[i.second]=p,nxt[p]=i.second;p=i.second;}
pre[0]=lne.back().second,nxt[lne.back().second]=0;
for(int i=0;i<(int)lne.size()-1;i++)
now=1LL*now*Gexp::getexp(1LL*loga[lne[i].second]*lne[i+1].first%Phi)%Mod;
}
void del(int x){
if(pre[x])now=1LL*now*Gexp::getexp(Phi-1LL*a[x]*loga[pre[x]]%Phi)%Mod;
if(nxt[x])now=1LL*now*Gexp::getexp(Phi-1LL*a[nxt[x]]*loga[x]%Phi)%Mod;
if(pre[x]&&nxt[x])now=1LL*now*Gexp::getexp(1LL*a[nxt[x]]*loga[pre[x]]%Phi)%Mod;
pre[nxt[x]]=pre[x];nxt[pre[x]]=nxt[x];
}
void add(int x){pre[nxt[x]]=x;nxt[pre[x]]=x;}
int ans[5005];
int main(){
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
cin>>n>>m;for(int i=1;i<=n;i++)cin>>a[i];
get_loga();Gexp::init();
for(int i=1;i<=m;i++)cin>>q[i].l>>q[i].r,q[i].id=i;
B=max(1.0,n/sqrt(m));sort(q+1,q+m+1);
for(int bl=n/B*B,br=n,p=1;br;br=bl-1,bl-=B){
bl=max(bl,1);build(bl,br);int l=bl,r=n;
while(p<=m&&q[p].l>=bl){
while(r>q[p].r)del(r--);
int tmp=now;
while(l<q[p].l)del(l++);
ans[q[p].id]=now;
while(l>bl)add(--l);
now=tmp;p++;
}
}
for(int i=1;i<=m;i++)cout<<ans[i]<<'\n';
return 0;
}
```
---
Problem：[HDU-7361](https://acm.hdu.edu.cn/showproblem.php?pid=7361)
时限 13s。
给你一个长为 $n$ 的数组 $a$（从 $1$ 开始）和长为 $m$ 的 $0/1$ 数组 $b$（从 $0$ 开始，并且 $0/1$ 数组这个条件不重要），你需要求出：
$$\sum_{i=1}^{n-m+1}\prod_{j=i}^{i+m-1}a_j^{b_{j-i}}\mod p$$
多测，$1\leq T\leq 10^3,1\leq m\leq n\leq 5\times 10^5,0\leq a_i<p=1062125569,\sum n\leq 10^6$。
---
Solution：
这道题没有上一道题莫队的性质：在低复杂度内移动询问左右端点。
但是 $a,b$ 数组固定，这让我们想到离散对数。
对数组 $a$ 求离散对数并枚举 $i$，询问转化为：
$$\sum_{j=i}^{i+m-1}b_{j-i}a_j\mod\varphi(p)$$
咦？ $\sum_{j=i}^{i+m-1}b_{j-i}a_j$ 看起来有点眼熟。
我们把 $b$ 数组扩充为长为 $n$ 的数组并将新增部分赋值为 $0$，反转，从而转化为：
$$\sum_{j=i}^{i+n-1}b_{n+i-j}a_j\mod\varphi(p)$$
将 $b$ 数组其他地方都赋值为 $0$，可以发现这个式子就是将 $a,b$ 数组卷积后的第 $n+i$ 项。
$\varphi(p)$ 不是质数，所以需要用到任意模数 NTT，这道题就做完了。
一些实现细节跟上面那道题差不多。
时间复杂度 $O(\dfrac p{\sqrt{\ln p}}+n(\log n+\log p))$。
---
Code：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int N=5e5,Mod=1062125569,Phi=Mod-1,G=11,sqrtP=32592,V=Mod;
namespace BSGS{
const int P=2e6+3;
int head[P+5],key[P+5],val[P+5],nxt[P+5],cnt;
void insert(int x,int id){
int hs=x%P;
key[++cnt]=x;val[cnt]=id;nxt[cnt]=head[hs];head[hs]=cnt;
}
int query(int x){
for(int i=head[x%P];i;i=nxt[i])if(key[i]==x)return val[i];
return -1;
}
void init(){
long long now=1,bnow=1;int lim=1e3;
for(int i=0;i<lim;i++)now=now*G%Mod;
for(int i=1;i<=V/lim+1;i++){
bnow=bnow*now%Mod;insert(bnow,i*lim);
}
}
int dlog(long long x){
for(int i=0;;x=x*G%Mod,i++){
int ret=query(x);if(~ret)return ret-i;
}
return -1;
}
}
namespace Gexp{
int bG[1<<15|5],gG[1<<15|5];
void init(){
bG[0]=gG[0]=1;
for(int i=1;i<=(1<<15);i++)bG[i]=1LL*bG[i-1]*G%Mod;
for(int i=1;i<=(1<<15);i++)gG[i]=1LL*gG[i-1]*bG[1<<15]%Mod;
}
long long getexp(int x){return 1LL*bG[x&((1<<15)-1)]*gG[x>>15]%Mod;}
}
const int dlg_1=Phi/2;
int dlg[sqrtP+5],pri[sqrtP+5],pcnt;
bool isprime[sqrtP+5];
void dlog_prepare(){
BSGS::init();
dlg[1]=0;
for(int i=2;i<=sqrtP;i++){
if(!isprime[i])dlg[i]=BSGS::dlog(i),pri[++pcnt]=i;
for(int j=1;j<=pcnt&&pri[j]*i<=sqrtP;j++){
isprime[i*pri[j]]=1,dlg[i*pri[j]]=(dlg[i]+dlg[pri[j]])%Phi;
if(i%pri[j]==0)break;
}
}
}
long long get_dlog(int x){
if(x<=sqrtP)return dlg[x];
int q=Mod/x,r=Mod-q*x;
return r<x-r?get_dlog(r)+dlg_1-dlg[q]:get_dlog(x-r)-dlg[q+1];
}
int qp(int a,int b,int mod){int ans=1;while(b){if(b&1)ans=1LL*ans*a%mod;a=1LL*a*a%mod;b>>=1;}return ans;}
namespace MTT{
const int Mod1=998244353,Mod2=1004535809,Mod3=469762049,G=3;
struct Int{
int v1,v2,v3;
Int():v1(0),v2(0),v3(0){}
Int(int _v):v1(_v),v2(_v),v3(_v){}
Int(int _v1,int _v2,int _v3):v1(_v1),v2(_v2),v3(_v3){}
};
Int operator +(Int a,Int b){return Int((a.v1+b.v1)%Mod1,(a.v2+b.v2)%Mod2,(a.v3+b.v3)%Mod3);}
Int operator -(Int a,Int b){return Int((a.v1+Mod1-b.v1)%Mod1,(a.v2+Mod2-b.v2)%Mod2,(a.v3+Mod3-b.v3)%Mod3);}
Int operator *(Int a,Int b){return Int(1LL*a.v1*b.v1%Mod1,1LL*a.v2*b.v2%Mod2,1LL*a.v3*b.v3%Mod3);}
int RealMod;
__int128 MN=(__int128)Mod1*Mod2*Mod3,M1,M2,M3,t1,t2,t3;
void InitMod(int Mod){
RealMod=Mod;
M1=471892799929712641LL;t1=850356001LL;
M2=468937312667959297LL;t2=395249030LL;
M3=1002772198720536577LL;t3=354521948LL;
}
int Real(Int a){
return (((1LL*a.v1*M1%MN*t1%MN)+
(1LL*a.v2*M2%MN*t2%MN)+
(1LL*a.v3*M3%MN*t3%MN))%MN)%RealMod;
}
int rev[1<<20|5];
void revid(int N){
for(int w=1,d=N>>1;w<N;w<<=1,d>>=1)
for(int p=0;p<w;p++)rev[p|w]=rev[p]|d;
}
Int pG[1<<20|5];
void initpG(int N){
const Int g(qp(G,(Mod1-1)/N,Mod1),qp(G,(Mod2-1)/N,Mod2),qp(G,(Mod3-1)/N,Mod3));
pG[0]=Int(1);for(int i=1;i<N;i++)pG[i]=pG[i-1]*g;
}
void NTT(Int *F,int N){
for(int i=0;i<N;i++)if(rev[i]>i)swap(F[rev[i]],F[i]);
for(int w=1;w<N;w<<=1){
for(int t=0;t<N;t+=2*w){
Int nG(1);
for(int i=0;i<w;i++,nG=nG*pG[N/2/w]){
const Int x=F[t|i],y=nG*F[t|i|w];
F[t|i]=x+y,F[t|i|w]=x-y;
}
}
}
}
Int F[1<<20|5],H[1<<20|5];
void Mul(int *f,int *h,int lenf,int lenh){
int N=1;while(N<=lenf+lenh)N<<=1;revid(N);initpG(N);
for(int i=0;i<N;i++)F[i]=Int(f[i]),H[i]=Int(h[i]);
NTT(F,N);NTT(H,N);for(int i=0;i<N;i++)F[i]=F[i]*H[i];
NTT(F,N);reverse(F+1,F+N);Int inv(qp(N,Mod1-2,Mod1),qp(N,Mod2-2,Mod2),qp(N,Mod3-2,Mod3));
for(int i=0;i<N;i++)f[i]=Real(F[i]*inv);
}
}
int F[1<<20|5],H[1<<20|5],n,m;
int b[N+5];
void init(){Gexp::init();dlog_prepare();MTT::InitMod(Phi);}
void work(){
cin>>n>>m;n--;
for(int i=0,x;i<=n;i++)cin>>x,F[i]=(get_dlog(x)%Phi+Phi)%Phi,F[i+n]=H[i]=H[i+n]=0;
for(int i=1;i<=m;i++)cin>>b[i];
for(int i=1;i<=m;i++)H[b[m]-b[i]]=1;
MTT::Mul(F,H,n,b[m]);
long long sum=0;
for(int i=b[m];i<=n;i++)sum=(sum+Gexp::getexp(F[i]))%Mod;
cout<<sum<<'\n';
}
int main(){
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
init();
int T;cin>>T;while(T--)work();
return 0;
}
```
---
## 结束
尽管离散对数的题目大多都是板题，但学了总比没学好。
并且离散对数肯定也不局限于我提到的方法，大家也可以集思广益多多讨论 ~~，虽然我觉得是没有~~。
总之谢谢大家！