简单数论
2021/4/18 10:56:07
本文主要是介绍简单数论,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!
exgcd
求解形如ax+by=c的方程的一组解
则:
\[x=x_0+\frac{b}{(a,b)}t,t\in Z \]void exgcd(ll a,ll b) { if(!b) { x=1,y=0; return; } ll k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; }
同余方程
合并方程(见小本本)
模板
ybtoj H. 中国剩余定理
#include<bits/stdc++.h> #define ll long long using namespace std; ll gcd(ll a,ll b) { return !b?a:gcd(b,a%b); } ll x,y; void exgcd(ll a,ll b) { if(!b) { x=1,y=0; return; } ll k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } inline ll ksj(ll a,ll b,ll p) { ll ret=0; while(b) { if(b&1) ret=(ret+a)%p; a=(a+a)%p,b>>=1; } return ret; } int n; int main() { while(scanf("%d",&n)!=EOF) { ll qp,qa; scanf("%lld%lld",&qp,&qa); bool fl=0; for(int i=2;i<=n;i++) { ll p,a; scanf("%lld%lld",&p,&a); if(!fl) { if(qa<a) { swap(qp,p),swap(qa,a); } ll g=gcd(qp,p); if((qa-a)%g) { fl=1; continue; } p/=g; exgcd(qp/g,p); x=(ksj(x,((qa-a)/g),p)+p)%p; ll t=qp*p; qa=(qa-ksj(qp,x,t)+t)%t; qp=t; } if(fl) puts("-1"); else printf("%lld\n",qa); } return 0; }
BSGS
求解最小的x形如
\[A^x=B(\% p) \]若(A,p)=1,则只需要判断\(x\in [0,\phi(p))\)
可以用分块的思想,将其分成\(a\times \sqrt n+b,b\in [0,\sqrt n)\)
然后用hash存储b,枚举a,用exgcd判断
若(A,p)!=1,则需要使(A,p)=1
可将方程写成:
\[A^x+py=B \]两边一直同除(A,p)即可,若(A,p)不能整除B,则无解
注意:如果存在\(A=0\),需要特判\(A=0\)的情况(因为exgcd时会使解x为0,可以找到对应的\(b\),而实际上答案只有两种\(0^0=1,0\),无解)
模板
ybtoj C. 【例题3】扩展BSGS
#include<bits/stdc++.h> #define ll long long using namespace std; int gcd(int a,int b) { return !b?a:gcd(b,a%b); } ll x,y; void exgcd(int a,int b) { if(!b) { x=1,y=0; return; } int k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } unordered_map<int,int>mp; int a,p,b; int bsgs() { int g=gcd(a,p),c=1,num=0; for(int g=gcd(a,p);g>1;g=gcd(a,p)) { if(b%g) return -1; b/=g,p/=g; c=(ll)c*(a/g)%p; num++; } int m=sqrt(p); mp.clear(); int t=1; for(int i=0;i<m;i++,t=(ll)t*a%p) { if(!mp.count(t)) { mp[t]=i; } } for(int i=0;(i+1)*m<=p;i++,c=(ll)c*t%p) { exgcd(c,p); x=((ll)x*b%p+p)%p; if(mp.count(x)) { return i*m+num+mp[x]; } } return -1; } int main() { while(scanf("%d%d%d",&a,&p,&b),a!=0||b!=0||p!=0) { int ans=bsgs(); if(ans==-1) puts("No Solution"); else printf("%d\n",ans); } return 0; }
例题1
Luogu P3306 [SDOI2013] 随机数生成器
求出通项公式(见本本)后BSGS(注意分类讨论)
WA了2次:
- 忘记特判a=1的情况(除0)
- 忘记特判a=0的情况(使不可能变为可能)
#include<bits/stdc++.h> #define ll long long using namespace std; int a,p,b,x1,Y1,t; inline int ksm(ll a,int b) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } ll x,y; void exgcd(int a,int b) { if(!b) { x=1,y=0; return; } int k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } unordered_map<int,int>mp; int bsgs() { int m=sqrt(p); mp.clear(); int t=1; for(int i=1;i<=m;i++) { t=(ll)t*a%p; if(!mp.count(t)) mp[t]=i; } int c=1; for(int i=0;i*m<=p;i++,c=(ll)c*t%p) { exgcd(c,p); x=((ll)x*b%p+p)%p; if(mp.count(x)) { return i*m+mp[x]; } } return -1; } int main() { int T; scanf("%d",&T); for(int e=1;e<=T;e++){ scanf("%d%d%d%d%d",&p,&a,&b,&x1,&t); if(a==0) { if(x1==t) { puts("1"); } else if(b==t) { puts("2"); } else { puts("-1"); } continue; } if(a==1) { t=(t+b-x1)%p; if(!b) { if(!t) puts("1"); else puts("-1"); continue; } exgcd(b,p); x=(x%p*t%p+p-1)%p+1; printf("%lld\n",x); continue; } int tmp=(ll)b*ksm(a-1,p-2)%p; Y1=(tmp+x1)%p; if(!Y1) { if((tmp+t)%p) { puts("-1"); } else puts("1"); continue; } b=((ll)(t+tmp)*a%p*ksm(Y1,p-2)%p+p)%p; int ans=bsgs(); printf("%d\n",ans); } return 0; }
Lucas
-
若p是质数,
\[C_n^m=C_{n\%p}^{m\%p}\times C_{n/p}^{m/p} \] -
反之,需要exLucas(质因子分解后扩展合并方程,将p提出即可,具体见小本本)
模板
Luogu P4720 【模板】扩展卢卡斯
Wa了一次:\(\%p^a\) 手贱写成了\(\%p\)
#include<bits/stdc++.h> #define ll long long using namespace std; ll n,m; inline int ksm(ll a,ll b,int p) { ll ret=1; while(b ){ if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } const int N=2e6+5; int jc[N]; inline int ask(ll x,int p,int pa) { if(x==0) return 1; return (ll)ksm(jc[pa],x/pa,pa)*ask(x/p,p,pa)%pa*jc[x-(x/pa)*pa]%pa; } inline int work(int p,int pa) { ll k1=0,k2=0,k3=0; for(ll i=p;i<=n;i*=p) { k1+=n/i; } for(ll i=p;i<=m;i*=p) { k2+=m/i; } for(ll i=p;i<=n-m;i*=p) { k3+=(n-m)/i; } jc[0]=1; for(int i=1;i<=pa;i++) { jc[i]=jc[i-1]; if(i%p) jc[i]=(ll)jc[i]*i%pa; } int phi=pa/p*(p-1); return (ll)ask(n,p,pa)*ksm(ask(n-m,p,pa),phi-1,pa)%pa*ksm(ask(m,p,pa),phi-1,pa)%pa*ksm(p,k1-k2-k3,pa)%pa; } ll x,y; void exgcd(int a,int b) { if(!b) { x=1,y=0; return; } int k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } int main() { int t; scanf("%lld%lld%d",&n,&m,&t); int qa=0,qp=0; for(int p=2;p<=sqrt(t);p++) { if(t%p==0) { int pa=1; for(;t%p==0;t/=p) pa=pa*p; if(!qp) { qp=pa,qa=work(p,pa); } else { int a=work(p,pa); exgcd(pa,qp); x=(x*(a-qa)%qp+qp)%qp; qp=qp*pa; qa=((-x*pa+a)%qp+qp)%qp; } } } if(t>1) { int p=t; if(!qp) { qp=p,qa=work(p,p); } else { int a=work(p,p); exgcd(p,qp); x=(x*(a-qa)%qp+qp)%qp; qp=qp*p; qa=((-x*p+a)%qp+qp)%qp; } } printf("%d\n",qa); return 0; }
例题1
Luogu P2480 [SDOI2010]古代猪文
求解:求\(G^{\sum{d|n}\ C_n^d}\ mod\ 999911659\)
需要求\(\sum_{d|n}C_{n}^{d} mod 999911658\)
套用exLucas也行,但发现:999911658=2*3*4679*35617
所以只需要对这4个用Lucas求解,然后合并方程
#include<bits/stdc++.h> #define ll long long const int p=999911659; using namespace std; inline int ksm(ll a,int b,int p) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } const int N=40000; struct A{ int p,inv[N],jc[N]; inline int C(int x,int y) { if(x>y) return 0; if(x<p&&y<p) { // printf("%lld %lld %lld\n",jc[y],inv[x],inv[y-x]); return (ll)jc[y]*inv[x]%p*inv[y-x]%p; } return (ll)C(x/p,y/p)*C(x%p,y%p)%p; } inline void init(int x) { p=x; inv[0]=inv[1]=1; for(int i=2;i<x;i++) { inv[i]=(ll)(p-p/i)*inv[p%i]%p; } jc[0]=1; for(int i=1;i<x;i++) { jc[i]=(ll)i*jc[i-1]%p; inv[i]=(ll)inv[i]*inv[i-1]%p; } } }c[5]; ll x,y; void exgcd(int a,int b) { if(!b) { x=1,y=0; return; } int k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } inline int ask(int m,int n) { int qp=c[1].p,qa=c[1].C(m,n); for(int i=2;i<=4;i++) { int p=c[i].p,a=c[i].C(m,n); if(qa<a) { swap(qp,p),swap(qa,a); } exgcd(qp,p); x=(x*(qa-a)%p+p)%p; int t=qp*p; qa=(qa-qp*x%t+t)%t; qp=t; } return qa; } int n,g; int main() { c[1].init(2),c[2].init(3),c[3].init(4679),c[4].init(35617); scanf("%d%d",&n,&g); int m=sqrt(n); ll ans=1; for(int i=1;i<=m;i++) { if(n%i==0) { ans=ans*ksm(g,ask(i,n),p)%p; if(i*i!=n) { ans=ans*ksm(g,ask(n/i,n),p)%p; } } } printf("%lld\n",ans); return 0; }
原根和阶
-
阶:最小的整数\(r\)满足\(a^r=0(\% p)\),记作\(ord_p(a)\)
当且仅当(a,p)互质时才存在阶
性质: 若\((a,p)=1,r|\phi(p)\)(之前证过了)
-
原根:满足阶=\(\phi(p)\)的数,
- 可以用\(<1,g^1,g^2,\cdots,g^{\phi(p)-1}>\)表示\(<1,2,\cdots,p-1>\)(所有与\(p\)互质的数)
- 当\(p\)为2,4,奇质数的幂次,及2*奇质数的幂次时原根一定存在,且最小的原根\(\leq p^{\frac{1}{4}}\)
- 如果存在原根,原根个数为\(\phi(\phi(p))\)
所以求原根只需要判断对于所有\(x\)是\(\phi(p)\)的约数(不包括\(\phi(p)\)),满足\(a^x!=1\)即可,更进一步,\(\phi(p)\)除以其质因子都满足上诉条件即可
模板
ybtoj E. 【例题5】最小原根
#include<bits/stdc++.h> #define ll long long using namespace std; int p; inline int ksm(ll a,int b) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } int main() { scanf("%d",&p); for(int x=2;x<=p-1;x++) { bool fl=0; int m=sqrt(p-1); for(int i=1;i<m;i++) { if((p-1)%i==0) { if(ksm(x,i)==1) { fl=1; break; } } } if(!fl) { for(int i=m;i>1;i--) { if((p-1)%i==0) { if(ksm(x,(p-1)/i)==1) { fl=1; break; } } } if(!fl) { printf("%d\n",x); return 0; } } } return 0; }
N次剩余
求解:\(x^k=n\)在\(\% p\)意义下的所有解,只考虑\(p\)为素数的情况
需要应用之前素数原根的性质
-
首先n=0没法用原根表示,所以需要特判掉
-
然后不妨设\(x=g^a,n=g^t\),则\(g^{ak}=g^t(\% p)\)
则:
\[ak=t(\%\phi(p)) \]
用exgcd和BSGS求出即可
模板
ybtoj F. 【例题6】N次剩余
忘记特判0了。。
#include<bits/stdc++.h> #define ll long long using namespace std; ll K,p,b; vector<ll>V; inline ll mo(ll x) { return x>=p?x-p:x; } ll gcd(ll a,ll b) { return !b?a:gcd(b,a%b); } inline ll ksj(ll a,ll b,ll p) { ll ret=0; while(b) { if(b&1) ret=(ret+a)%p; a=(a+a)%p,b>>=1; } return ret; } inline ll ksm(ll a,ll b) { ll ret=1; while(b) { if(b&1) ret=ksj(ret,a,p); a=ksj(a,a,p),b>>=1; } return ret; } inline ll get() { int m=sqrt(p-1); for(ll x=2;x<=p-1;x++) { bool fl=0; for(int i=1;i<m;i++) { if((p-1)%i==0) { if(ksm(x,i)==1) { fl=1; break; } } } if(!fl) { for(int i=m;i>1;i--) { if((p-1)%i==0) { if(ksm(x,(p-1)/i)==1) { fl=1; break; } } } if(!fl) return x; } } } ll x,y; void exgcd(ll a,ll b) { if(!b) { x=1,y=0; return; } ll k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } unordered_map<ll,int>mp; inline ll ask(ll a) { int m=sqrt(p-1); ll t=1; mp.clear(); for(int i=0;i<m;i++,t=ksj(t,a,p)) { if(!mp.count(t)) mp[t]=i; } ll c=1; for(int i=0;(ll)m*(i+1)<=p;i++,c=ksj(c,t,p)) { exgcd(c,p); x=mo(ksj(x,b,p)+p); if(mp.count(x)) return (ll)m*i+mp[x]; } } int main() { int id=0; while(scanf("%lld%lld%lld",&K,&p,&b)!=EOF) { id++,printf("case%d:\n",id); ll G=get(),g=gcd(K,p-1),t=ask(G); if(t%g) { puts("-1"); continue; } t/=g; ll tt=(p-1)/g; exgcd(K/g,tt); x=(ksj(x,t,tt)+tt)%tt; V.clear(); int S=0; for(;x<p;x+=tt) { V.push_back(ksm(G,x)),S++; } sort(V.begin(),V.end()); printf("%lld\n",V[0]); for(int i=1;i<S;i++) { if(V[i]!=V[i-1])printf("%lld\n",V[i]); } } return 0; }
例题 1
ybtoj N. 数论之神
考虑拆解\(2k+1=p_1^{a1}p_2^{a2}\cdots p_c^{a_c}\)
则解的个数为\(x^A=B(\% p_i^{a_i})\) 的解的个数的乘积(feel一下)
接下来分类讨论\(x^A=B(\% p_i^{a_i})\)的解的个数
Wa了一次:跳出搞错了
#include<bits/stdc++.h> #define ll long long using namespace std; int gcd(int a,int b) { return !b?a:gcd(b,a%b); } inline int ksm(int a,int b) { int ret=1; while(b) { if(b&1) ret=ret*a; a=a*a,b>>=1; } return ret; } inline int ksmp(ll a,int b,int p) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } int A,B,K,phi,p,b; inline int get(int x,int num) { int m=sqrt(phi); for(int i=2;i<x;i++) { bool fl=0; for(int j=2;j<=m;j++) { if(phi%j==0) { if(ksmp(i,j,x)==1||ksmp(i,phi/j,x)==1) { fl=1; break; } } } if(!fl) return i; } return 0; } ll x,y; void exgcd(int a,int b) { if(!b) { x=1,y=0; return; } int k=a/b; exgcd(b,a-k*b); int t=x; x=y,y=t-k*y; } map<int,int>mp; int bsgs(int a,int b,int p) { int m=sqrt(p-1); mp.clear(); int t=1; for(int i=0;i<m;i++,t=(ll)t*a%p) { if(!mp.count(t)) { mp[t]=i; } } int c=1; for(int i=0;i*m<p;i++,c=(ll)c*t%p) { exgcd(c,p); x=((ll)x*b%p+p)%p; if(mp.count(x)) { return i*m+mp[x]; } } } int main() { int T; scanf("%d",&T); while(T--) { scanf("%d%d%d",&A,&B,&K); int ans=1,t=K*2+1; for(p=2;p<=sqrt(t);p++) { if(ans&&t%p==0) { int pa=1; int a=0; for(;t%p==0;t/=p) pa=pa*p,a++; int b=B%pa; if(b%pa==0) { a-=(a-1)/A+1; ans=ans*ksm(p,a); continue; } if(b%p==0) { int t1=0; for(;b%p==0;b/=p)t1++; if(t1%A) { ans=0; continue; } else { ans=ans*ksm(p,t1-t1/A); pa=ksm(p,a-t1); } } phi=pa/p*(p-1); int g=get(pa,a); if(!g) { ans=0; continue; } b=bsgs(g,b,pa); if(b%gcd(A,phi)) { ans=0;continue; } ans*=gcd(A,phi); } } if(ans&&t>1) { p=t; int pa=1; int a=0; for(;t%p==0;t/=p) pa=pa*p,a++; int b=B%pa; if(b==0) { a-=(a-1)/A+1; ans=ans*ksm(p,a); printf("%d\n",ans); continue; } if(b%p==0) { int t1=0; for(;b%p==0;b/=p)t1++; if(t1%A) { ans=0; printf("%d\n",ans); continue; } else { ans=ans*ksm(p,t1-t1/A); pa=ksm(p,a-t1); } } phi=pa/p*(p-1); int g=get(pa,a); b=bsgs(g,b,pa); if(b%gcd(A,phi)) { ans=0; printf("%d\n",ans); continue; } ans*=gcd(A,phi); } printf("%d\n",ans); } return 0; }
二次剩余
求解:\(x^2=n\)在\(\% p\) 意义下,p为奇素数
法1:当作N次剩余
法2:直接背诵:
-
若存在\(n\)的二次剩余,则\(n^{\frac{p-1}{2}}=1\),否则\(n^{\frac{p-1}{2}}=-1\)
-
求解:先找出\(a\),使得\(w=a^2-n\)无二次剩余,答案即为\((a+\sqrt w)^{\frac{p+1}{2}}\)
其中可以设虚根\(i=\sqrt w\) ,然后用类似虚数的运算即可,最后得到的是实数
模板
Luogu P5491 【模板】二次剩余
WA了一次:虚根是\(\sqrt w\) 不是\(i\),所以实部中的一个乘积需要乘上\(w\)
#include<bits/stdc++.h> #define ll long long using namespace std; int p; inline int rnd() { return rand()*rand()%p; } inline int ksm(ll a,int b) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } struct cp{ll x,y; }; int w; inline cp mul(cp i,cp j) { return (cp){((i.x*j.x%p+i.y*j.y%p*w%p)%p+p)%p,(i.x*j.y%p+i.y*j.x%p)%p}; } inline int cpksm(cp a,int b) { cp ret=(cp){1,0}; while(b) { if(b&1) ret=mul(ret,a); a=mul(a,a),b>>=1; } return ret.x; } int ask(int x) { x%=p; if(ksm(x,p-1>>1)==p-1) return -1; int a=rnd(); while(w=(((ll)a*a-x)%p+p)%p,ksm(w,p-1>>1)!=p-1) { a=rnd(); } return cpksm((cp){a,1},p+1>>1); } int main() { int T; scanf("%d",&T); while(T--) { int x; scanf("%d%d",&x,&p); if(!x) { puts("0"); continue; } int t1=ask(x); if(t1==-1) puts("Hola!"); else { int t2=p-t1; if(t1>t2) swap(t1,t2); if(t2!=t1) printf("%d %d\n",t1,t2); else printf("%d\n",t1); } } }
例题1
ybtoj Fib数列
见小本本
#include<bits/stdc++.h> #define ll long long const int p=1e9+9; using namespace std; inline int rnd() { return rand()*rand()%p; } inline int ksm(ll a,int b) { ll ret=1; while(b) { if(b&1) ret=ret*a%p; a=a*a%p,b>>=1; } return ret; } struct cp{ll x,y; }; int w; inline cp mul(cp i,cp j) { return (cp){((i.x*j.x%p+i.y*j.y%p*w%p)%p+p)%p,(i.x*j.y%p+i.y*j.x%p)%p}; } inline int cpksm(cp a,int b) { cp ret=(cp){1,0}; while(b) { if(b&1) ret=mul(ret,a); a=mul(a,a),b>>=1; } return ret.x; } int ask(int x) { x%=p; if(!x) return 0; if(ksm(x,p-1>>1)==p-1) return -1; int a=rnd(); while(w=(((ll)a*a-x)%p+p)%p,ksm(w,p-1>>1)!=p-1) { a=rnd(); } return cpksm((cp){a,1},p+1>>1); } ll x,y; void exgcd(ll a,ll b) { if(!b) { x=1,y=0; return; } ll k=a/b; exgcd(b,a-k*b); ll t=x; x=y,y=t-k*y; } unordered_map<ll,int>mp; int a,b,n; inline int bsgs() { int m=sqrt(p-1); ll t=1; mp.clear(); for(int i=0;i<m;i++,t=t*a%p) { if(!mp.count(t)) mp[t]=i; } ll c=1; for(int i=0;(ll)m*(i+1)<=p;i++,c=c*t%p) { exgcd(c,p); x=(x*b%p+p)%p; if(mp.count(x)) return (ll)m*i+mp[x]; } return p+1; } int main() { int m=ask(5); a=(m+1)/2; scanf("%d",&n); int tmp=ask(5LL*n%p*n%p+4),ans=p+1; if(tmp!=-1) { b=((ll)m*n%p+tmp)*ksm(2,p-2)%p; int t=bsgs(); if(!(t&1)) ans=min(ans,t); b=(((ll)m*n%p-tmp)*ksm(2,p-2)%p+p)%p; t=bsgs(); if(!(t&1)) ans=min(ans,t); } tmp=ask((5LL*n%p*n%p-4+p)%p); if(tmp!=-1) { b=((ll)m*n%p+tmp)*ksm(2,p-2)%p; int t=bsgs(); if(t&1) ans=min(ans,t); b=(((ll)m*n%p-tmp)*ksm(2,p-2)%p+p)%p; t=bsgs(); if(t&1) ans=min(ans,t); } if(ans>p) puts("-1"); else printf("%d\n",ans); return 0; }
这篇关于简单数论的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!
- 2024-06-19《2023版Java工程师》课程升级公告
- 2024-06-15matplotlib作图不显示3D图,怎么办?
- 2024-06-1503-Loki 日志监控
- 2024-06-1504-让LLM理解知识 -Prompt
- 2024-06-05做软件测试需要懂代码吗?
- 2024-06-0514-ShardingSphere的分布式主键实现
- 2024-06-03为什么以及如何要进行架构设计权衡?
- 2024-05-31全网首发第二弹!软考2024年5月《软件设计师》真题+解析+答案!(11-20题)
- 2024-05-31全网首发!软考2024年5月《软件设计师》真题+解析+答案!(21-30题)
- 2024-05-30【Java】百万数据excel导出功能如何实现