原根&离散对数相关

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]可爱的质数

由于m是质数,这就是个经典离散对数问题。质数真的非常可爱啊qwq。

#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？
P2485 [SDOI2011]计算器 (一些特判)
在递推式上的扩展 : stO Itst!
当模数不是质数 : P4195 【模板】exBSGS/Spoj3105 Mod

由鸽笼原理容易得知,每个数的幂的循环节不大于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)可以不计。

注意还原子问题以及其他细节,具体可以看代码。

#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 【模板】原根

形如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.

#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与两位神仙

这道题让我发现我的数论在这个分支有多么渣,所以我才来写了这篇文章。

题意 : 给出模数m,保证其是奇质数的幂。

有T组x,y,保证都与m互质,询问是否存在自然数a使得x^a=y\pmod m

设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

初始的阶为$\varphi(m)$,对于其每个**质因数**,逐次降低次数来判定,就能摸到阶在这个质数上的最高次数。需要`Pollard-Rho`来求出$\varphi(m)$以及其所有质因数。每次求阶的复杂度为$O(\omega(m)\log m)$,总的复杂度为$O(\sqrt[4]{m}+T\omega(m)\log m)

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

原根&离散对数相关

离散对数问题

原根和阶

例题