再观BZOJ2219数论之神...
题解:
当时A掉的时候其实并不会...
于是现在膜别人的题解终于会啦!
求方程$x^a\equiv{b}\pmod{n}$在$[0,n-1]$内的解的个数.
我们考虑通过对于两边取指标,转化为线性同余方程问题.
但是发现$n$不一定存在原根,我们将n进行质因数分解.
原方程的解的个数等于每一个子方程$x^a\equiv{b}\pmod{p_i^{q_i}}$解的个数的乘积.
我们考虑中国剩余定理,就能发现这样的乘法原理显然成立.
于是我们只考虑每一个子方程的解的个数.
我们考虑以下几种情况:
(1)$b$存在指标,即$(b,p_i)=1$
我们对两边取指标得到$a\times{ind(x)}\equiv{ind(b)}\pmod{\phi(p_i^{q_i})}$
由线性同余方程的知识可知,如果$(a,\phi(p_i^{q_i}))\not|ind(b)$,则无解;
否则解的个数为$(a,\phi(p_i^{q_i}))$.
(2)$b$不存在指标,且$b=0$
我们可以得到$p_i|x$,不妨令$x=cp_i^{x'}((c,p_i)=1)$,可以得到$ax'\geq{q_i},即x'\geq{\lceil\frac{q_i}{a}\rceil}$.
所以我们有只需满足$p_i^{\lceil\frac{q_i}{a}\rceil}|x$,于是解的个数就为$p_i^{q_i-\lceil\frac{q_i}{a}\rceil}$.
(3)$b$不存在指标,且$b\not=0$
不妨令$b=cp_i^{b'}((c,p_i)=1)$.
若$b'\%a\not=0$,显然无解.
否则我们将等式两边同除$p_i^{b'}$得到:
$(\frac{x}{p_i^{\frac{b'}{a}}})^a\equiv{c}\pmod{p_i^{q_i-b'}}$
令$x'=\frac{x}{p_i^{\frac{b'}{a}}}$,由于满足(2)的条件,我们可以求出$x'$的解的个数.
注意$x'$的范围是$[0,p_i^{q_i-b'})$,但原式的$\frac{x}{p_i^{\frac{b'}{a}}}$范围是$[0,p_i^{q_i-\frac{b'}{a}})$.
这相当于我们求出的解是在每个长度为$p_i^{q_i-b'}$的区间中的解的个数.
所以我们需要将解的个数乘以$p^{b'-\frac{b'}{a}}$.
于是我们就彻底解决了这个问题!
代码:
#include<cmath> #include<cstdio> #include<cctype> #include<climits> #include<cstring> #include<cstdlib> #include<iostream> #include<algorithm> using namespace std; typedef long long ll; ll pow(ll x,ll y,ll p){ ll re=1,t=x; for(;y;y>>=1,t=t*t%p) if(y&1) re=re*t%p; return re; } ll pow(ll x,ll y){ ll re=1,t=x; for(;y;y>>=1,t*=t) if(y&1) re*=t; return re; } ll gcd(ll a,ll b){ return(!b)?a:gcd(b,a%b); } inline void exgcd(ll a,ll b,ll&d,ll&x,ll&y){ if(!b){ d=a; x=1; y=0; } else{ exgcd(b,a%b,d,y,x); y-=x*(a/b); } } inline ll inv(ll a,ll n){ ll d,x,y; exgcd(a,n,d,x,y); return(x%n+n)%n; } inline ll get_g(ll p){ static int pr[110]; int top=0; ll t=p-1; for(int i=2;i*i<=p-1;++i){ if(t%i==0){ pr[++top]=i; while(t%i==0) t/=i; } } if(t!=1) pr[++top]=t; for(int i=1;;++i){ bool ok=1; for(int j=1;j<=top;++j) if(pow(i,(p-1)/pr[j],p)==1){ ok=0; break; } if(ok) return i; } } struct Hashset{ static const int mod=13131; int head[mod],next[30010],f[30010],v[30010],ind; inline void reset(){ ind=0; memset(head,-1,sizeof head); } inline void insert(int x,int _v){ int ins=x%mod; for(int j=head[ins];j!=-1;j=next[j]) if(f[j]==x){ v[j]=min(v[j],_v); return; } f[ind]=x; v[ind]=_v; next[ind]=head[ins]; head[ins]=ind++; } inline int operator[](int x){ int ins=x%mod; for(int j=head[ins];j!=-1;j=next[j]) if(f[j]==x) return v[j]; return -1; } }M; ll BSGS(ll c,ll a,ll b,ll p){ ll m=ceil(sqrt(p)),t=1,ct=c; M.reset(); for(int i=0;i<m;++i,(t*=a)%=p,(ct*=a)%=p) M.insert(ct,i); t=inv(t,p); for(int i=0;i*m<p;++i,(b*=t)%=p) if(M[b]!=-1) return i*m+M[b]; return -1; } ll EXBSGS(ll a,ll b,ll p){ ll c=1,tim=0,_gcd; while((_gcd=gcd(a,p))!=1){ if(b%_gcd) return -1; p/=_gcd; ++tim; (c*=(a/_gcd))%=p; (b/=_gcd)%=p; } ll ans=BSGS(c,a,b,p); if(ans!=-1) ans+=tim; return ans; } ll solve(ll a,ll b,ll p,ll t){ b%=pow(p,t); if(b==0) return pow(p,t-((t-1)/a+1)); int _b=0; while(b%p==0) b/=p,++_b; ll phi=pow(p,t)*(p-1)/p; ll g=get_g(p); if(_b==0){ ll ind_b=EXBSGS(g,b,pow(p,t)); if(ind_b%gcd(a,phi)) return 0; return gcd(a,phi); } else{ if(_b%a) return 0; else return solve(a,b,p,t-_b)*pow(p,_b-_b/a); } } int main(){ #ifndef ONLINE_JUDGE freopen("tt.in","r",stdin); #endif int T; cin>>T; ll A,B,K,n,_n; int cnt; while(T--){ scanf("%lld%lld%lld",&A,&B,&K); n=_n=K<<1^1; ll ans=1; for(int i=2;i*i<=n;++i){ if(_n%i==0){ cnt=0; while(_n%i==0) ++cnt,_n/=i; ans*=solve(A,B,i,cnt); } } if(_n!=1) ans*=solve(A,B,_n,1); printf("%lld\n",ans); } return 0; }
Sep 02, 2015 01:10:02 PM
看个博客还有注册个号差评。。
Sep 02, 2015 06:31:24 PM
@ww140142: 那是为了要评论没办法啊。。。
Sep 02, 2015 06:32:54 PM
并不是这个意思,我是说那个网址进去要在OJ注册账号= =
Sep 02, 2015 08:05:59 PM
@ww140142: 好吧0.0