原根&离散对数相关
command_block
2020-03-16 08:12:21
- ## 离散对数问题
求解模方程$a^x=c\pmod p$
著名的数论难题,目前还没有$O(poly(\log))$的做法。
如果$p$是质数,这就是经典`BSGS`.
考虑费马小定理 : $a^x=a^{x\mod (p-1)}$,所以幂是有循环节的。
我们考虑分块,即$a^{kT+b}$,
预处理$a^{0...(T-1)}$作为$b$,此时要求$a^{kT}=ca^{-b}\pmod p$.
将$ca^{-b}$放入`Hash`表中。
然后预处理$a^{kT}$并在`Hash`表中查询即可。
取$T=\sqrt{p}$即可得到复杂度为$O(\sqrt{p})$,`Hash`表可能带`\log`.
- 模板题 : [P3846 [TJOI2007]可爱的质数](https://www.luogu.com.cn/problem/P3846)
由于$m$是质数,这就是个经典离散对数问题。质数真的非常可爱啊qwq。
```cpp
#include<algorithm>
#include<cstdio>
#include<map>
#define ll long long
using namespace std;
ll mod,T;
ll powM(ll a,int t=mod-2)
{
ll ret=1;
while(t){
if (t&1)ret=ret*a%mod;
a=a*a%mod;t>>=1;
}return ret;
}
ll x,y;
map<int,int> sav;
int main()
{
scanf("%lld%lld%lld",&mod,&x,&y);
while(T*T<=mod)T++;
ll buf=powM(x),pro=1;
for (int i=0;i<T;i++){
sav[pro*y%mod]=i;
pro=pro*buf%mod;
}buf=powM(pro);pro=1;
for (int i=0;i<T;i++){
if (sav.count(pro)){
printf("%d\n",i*T+sav[pro]);
return 0;
}pro=pro*buf%mod;
}puts("no solution");
return 0;
}
```
附送 :
- [P4884 多少个1？](https://www.luogu.com.cn/problem/P4884)
- [P2485 [SDOI2011]计算器](https://www.luogu.com.cn/problem/P2485) (一些特判)
- 在递推式上的扩展 : [stO Itst!](https://www.luogu.com.cn/blog/pengSiJin/solution-p4000)
- 当模数不是质数 : [P4195 【模板】exBSGS/Spoj3105 Mod](https://www.luogu.com.cn/problem/P4195)
由鸽笼原理容易得知,每个数的幂的循环节不大于$O(p)$,当$a\perp p$的时候,直接使用普通`BSGS`。
问题在于,当$a^x=c\pmod p$中$a\not\perp p$的时候,我们就不能求逆元了。
考虑求$d=(a,p)$,让方程两边除以$d$,这样能产生互质。
问题在于可能$d\!\!\not |\ c$,由于$a$的幂膜上$p$之后一定残余$d$的倍数,可得$c≠1$时无解。判掉。
现在我们成功除了一下,变成$a^{x-1}\frac{a}{d}=\frac{c}{d}\pmod{\frac{p}{d}}$
此时,有$\frac{a}{d}\perp \frac{p}{d}$,我们就可以求逆了。这里需要`EXgcd`算法。
$a^{x-1}=\frac{c}{d}(\frac{a}{d})^{-1}\pmod{\frac{p}{d}}$
注意,我们习惯$\frac{a}{(a,b)}\perp\frac{b}{(a,b)}$,但是这里$a$不一定和$\frac{p}{d}$互质,于是我们就进入了子问题。
如此重复直到$d=1$为止,此时我们就可以正常求逆了。
注意可能在转化过程中就已有$a=c$,这时要及时停下来。
每次$p$至少减半,所以这部分复杂度为$O(\log p)$可以不计。
注意还原子问题以及其他细节,具体可以看代码。
```cpp
#include<algorithm>
#include<cstdio>
#include<map>
#define ll long long
using namespace std;
ll gcd(ll a,ll b)
{return !b ? a : gcd(b,a%b);}
void exgcd(ll a,ll b,ll &x,ll &y){
if (a==1&&b==0){x=1;y=0;return ;}
exgcd(b,a%b,y,x);y-=x*(a/b);
}
ll inv(ll a,int mod){
ll x,y;exgcd(a,mod,x,y);
return (x+mod)%mod;
}
map<int,int> sav;
ll BSGS(ll a,ll c,ll mod)
{
sav.clear();
int T=1;while(T*T<mod)T++;
ll buf=1,xp=inv(a,mod);
for (int i=0;i<T;i++){
if (!sav.count(buf*c%mod))
sav[buf*c%mod]=i;
buf=buf*xp%mod;
}xp=inv(buf,mod);buf=1;
for (int i=0;i<=T;i++){
if (sav.count(buf))
return i*T+sav[buf];
buf=buf*xp%mod;
}return -1000;
}
ll solve(ll a,ll c,int mod)
{
if (a==c)return 1;
ll d=gcd(a,mod);
if (d==1)return BSGS(a,c,mod);
if (c%d)return -1000;
mod/=d;
return solve(a%mod,(c/d)*inv(a/d,mod)%mod,mod)+1;
}
ll mlog(ll a,ll c,int mod)
{
if (c==1)return 0;
if (!a)return !c ? 1 : -1;
return solve(a,c,mod);
}
int main()
{
ll a,c;int mod;
while(1){
scanf("%lld%d%lld",&a,&mod,&c);
if (!mod)break;
a%=mod;c%=mod;
ll ret=mlog(a,c,mod);
printf(ret<0
? "No Solution\n"
: "%lld\n",ret
);
}return 0;
}
```
- ## 原根和阶
- **阶** : $a$ 在模 $p$ 意义下的阶是 : 最小的整数 $k$ 使得 $a^k=1\pmod p$.
对于$a\not\perp p$的情况,阶为$∞$,或认为不存在。
说白了就是幂的最小循环节,根据欧拉定理这个上界是$\varphi(p)$.
记$δ_p(a)$为$a$在模$p$下的阶。
- **原根**
满足$δ_p(g)=\varphi(p)$的$g$(达到上界),我们称之为$p$的原根。
模板题 : [P6091 【模板】原根](https://www.luogu.com.cn/problem/P6091)
形如$2,4,2p^c,p^c(p\text{是奇质数})$的数才有原根。
一旦我们求出了其中一个原根$g$,我们可以这样构造出其余的原根 :
所有的$g^k$当$\big(k\perp\varphi(p)\big)$,这样恰好构造出$\varphi(\varphi(p))$个,而且能够证明没有遗漏。
**理解** : 当$k\perp\varphi(p)$时,在$\varphi(p)$的环上每次走$k$步能够遍历整个环。
这也告诉我们,所有$t\perp p$的$t$都能表示为原根的幂次,这一点很重要。
还能推出,原根是相当稠密的,另一个已被证明的结论是,最小原根的大小是$O(p^{0.25})$的。
现在问题变成了求出最小原根,我们可以从小到大枚举。
如果按照定义求阶,每次的复杂度是$O(p)$的,显然不优。
- **引理** : 模$p$的阶如果存在,一定是$\varphi(p)$的约数。
**证明** : 设这个数有阶,则被表示成$g^k$的形式。
能够在$\varphi(p)$大小的环上得到一个$\frac{\varphi(p)}{\gcd(k,\varphi(p))}$大小的环。
这就是它的阶,而且显然为$\varphi(p)$的约数。
所以,我们每次计算阶的时候,只需要验证$\varphi(p)$的约数即可。
当然,我们不需要确切地计算出阶,我们只需要查看阶是否等于$\varphi(p)$而已。
注意到当$a^k=1\pmod p$,那么$a^{ck}=1\pmod p$。
我们可以取出$\varphi(p)$的质因数$p_1,p_2...p_m$,分别判定$\dfrac{\varphi(p)}{p_i}$即可。
因为这能够不遗漏地覆盖所有其它约数的倍数,而又达不到$\varphi(p)$本身。可以视为高维空间的可重复差分。
复杂度是$O(p^{0.25}\omega(p)\log p+\sqrt{p})$
求出原根之后,构造并排序,总复杂度$O(Tp\log p)$.
讲道理,直接求最小原根不就好了,哪有题目需要同时使用很多个原根啊?
注意特判`2`.
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 1005000
#define ll long long
#define pf printf
using namespace std;
inline int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
int mod;
ll powM(ll a,int t=mod-2)
{
ll ret=1;
while(t){
if (t&1)ret=ret*a%mod;
a=a*a%mod;t>>=1;
}return ret;
}
int gcd(int a,int b)
{return !b ? a : gcd(b,a%b);}
int getphi(int n)
{
int ret=n,p=2;
while(p*p<=n){
if (n%p==0){
while(n%p==0)n/=p;
ret=ret/p*(p-1);
}p++;
}if (n>1)ret=ret/n*(n-1);
return ret;
}
int d[105],tn;
void getft(int n)
{
int sav=n,p=2;
while(p*p<=n){
if (n%p==0){
while(n%p==0)n/=p;
d[++tn]=sav/p;
}p++;
}if (n>1)d[++tn]=sav/n;
}
bool check(int n)
{
for (int i=1;i<=tn;i++)
if (powM(n,d[i])==1)
return 0;
return 1;
}
int ans[MaxN];
void solve()
{
mod=read();
int dr=read(),phi=getphi(mod),g=0;
if (mod==2){
puts("1");
if (dr==1)printf("1");
puts("");
return ;
}
tn=0;getft(phi);
for (int i=1;i<=200;i++)
if (gcd(i,mod)==1&&check(i))
{g=i;break;}
tn=0;
if (g){
for (int i=1;i<phi;i++)
if (gcd(i,phi)==1)
ans[++tn]=powM(g,i);
sort(ans+1,ans+tn+1);
}
printf("%d\n",tn);
for (int i=dr;i<=tn;i+=dr)
printf("%d ",ans[i]);
puts("");
}
int main()
{
int T=read();
while(T--)solve();
return 0;
}
```
- ## 例题
- [P5605 小A与两位神仙](https://www.luogu.com.cn/problem/P5605)
这道题让我发现我的数论在这个分支有多么渣,所以我才来写了这篇文章。
**题意** : 给出模数$m$,保证其是奇质数的幂。
有$T$组$x,y$,保证都与$m$互质,询问是否存在自然数$a$使得$x^a=y\pmod m$
$m\leq 10^{18},T\leq 2\times 10^4$,时限$\texttt{3s}$,允许`__int128`.
设$m$的原根为$g$,由于$x,y$与$m$互质,可以分别设为$g^b,g^c$.
方程变为 $g^{ab}=g^c\pmod {m}$,也即$ab=c\pmod {\varphi(m)}$
暴力`BSGS`求离散对数就能拿到部分分了。
根据斐蜀定理,这个方程有解的充要条件是$(b,\varphi(m))|(c,\varphi(m))$
观察$g^b$的幂组成的集合,相当于在$\varphi(m)$的环上,每次跳$b$步,则有$\frac{\varphi(m)}{\gcd(b,\varphi(m))}$个元素。
这个集合的大小正是$g^b$的阶$δ_m(x)$,这就能够得到$(b,\varphi(m))=\dfrac{\varphi(m)}{δ_m(x)}$
简单变化一下,就能得到原来的条件等价于$δ_m(y)|δ_m(x)$,现在问题就是快速求阶。
前面介绍过,阶必定是$\varphi(m)$的约数。由此我们可以$O(d(\varphi(m))\log m)$求解。
在$10^{18}$以内$d_{\max}≈10^5$,就算有素数的限制,也能够造出$2\times10^4$,看起来会超时的样子。
注意到当$a^k=1\pmod m
$,那么$a^{ck}=1\pmod m$,这能够查看阶是否是某个数的约数。
初始的阶为$\varphi(m)$,对于其每个**质因数**,逐次降低次数来判定,就能摸到阶在这个质数上的最高次数。
需要`Pollard-Rho`来求出$\varphi(m)$以及其所有质因数。
每次求阶的复杂度为$O(\omega(m)\log m)$,总的复杂度为$O(\sqrt[4]{m}+T\omega(m)\log m)$
```cpp
#include<algorithm>
#include<cstdio>
#include<cmath>
#define ll long long
#define sf scanf
using namespace std;
inline ll mul(ll a,ll b,ll m){
ll r=((long double)a/m*b+0.5);
r=a*b-r*m;
return r<0?r+m:r;
}
ll powM(ll a,ll t,ll m)
{
ll ans=1;
while(t){
if (t&1)ans=mul(ans,a,m);
a=mul(a,a,m);t>>=1;
}return ans;
}
bool mr(ll n,ll a)
{
ll t=n-1;
while(!(t&1))t>>=1;
ll buf=powM(a,t,n);
if (buf==1||buf==n-1)return 0;
while((t<<=1)<n-1){
buf=mul(buf,buf,n);
if (buf==n-1)return 0;
}return 1;
}
const int testp[8]={2,3,5,7,13,19,61,233};
bool ptest(ll n)
{
if (n<2)return 0;
for (int i=2;i*i<=min(n,10000ll);i++)
if (n%i==0)return 0;
if (n<=10000)return 1;
for (int i=0;i<8;i++)
if (mr(n,testp[i]))return 0;
return 1;
}
inline ll gcd(ll a,ll b){
if(a==0)return b;
if(a<0)a=-a;
ll t;
while(b){t=a%b;a=b;b=t;}
return a;
}
ll sav[150];int tot;
const int lim=128;
ll prho(ll n,ll c)
{
ll x1=(c+1)%n,x2=(mul(x1,x1,n)+c)%n,buf;
tot=0;
while(1){
buf=mul(x1-x2,buf,n);
sav[++tot]=x1-x2;
if (tot==lim){
if (gcd(buf,n)>1)break;
tot=0;
}
x1=mul(x1,x1,n)+c;
x2=mul(x2,x2,n)+c;
x2=mul(x2,x2,n)+c;
}
for (int i=1;i<=tot;i++){
buf=gcd(sav[i],n);
if (buf>1)return buf;
}return n;
}
ll p[105];int tn;
void solve(ll n)
{
if (ptest(n)){p[++tn]=n;return ;}
ll sav=prho(n,rand()%(n-1)+1);
while(sav==n)sav=prho(n,rand()%(n-1)+1);
solve(sav);solve(n/sav);
}
ll getphi(ll n)
{
tn=1;
for (int i=60;i;i--){
ll v=pow(n,1.0/i);
if (v-1&&!powM(v-1,i,n)){p[1]=v-1;break;}
if (v &&!powM(v ,i,n)){p[1]=v ;break;}
if (v+1&&!powM(v+1,i,n)){p[1]=v+1;break;}
}if (!p[1])p[1]=n;
ll sav=p[1]-1;
for (int i=2;i<=100;i++){
if (sav%i==0){
p[++tn]=i;
while(sav%i==0)sav/=i;
}
}
if (sav>1)solve(sav);
sav=p[1]-1;
for (int i=2;i<=tn;i++){
if (sav%p[i]==0){
while(sav%p[i]==0)sav/=p[i];
}else p[i]=0;
}return n/p[1]*(p[1]-1);
}
ll mod,phi;
ll ord(ll u)
{
ll ans=1;
for (int i=1;i<=tn;i++)if (p[i]){
ll sav=phi;int c=0;
while(sav%p[i]==0){sav/=p[i];c++;}
for (ll x=1;c;x*=p[i],c--)
if (powM(u,sav*x,mod)>1)
ans*=p[i];
}return ans;
}
int T;
int main()
{
sf("%lld%d",&mod,&T);
phi=getphi(mod);
for (int i=1;i<=T;i++){
ll x,y;sf("%lld%lld",&x,&y);
puts(ord(x)%ord(y)==0 ? "Yes" : "No");
}return 0;
}
```