【模板整合计划】数论数学—数论
【模板整合计划】数论数学—数论
一：【质数】
1.【暴力判】
【模板】素数、コンテスト、素数 \(\text{[AT807]}\)
#include<cstdio>
#include<cmath>
int x;
inline bool judge(int n){
if(n<4)return 1;
if(n%2==0)return 0;
int half=sqrt(n);
for(int i=3;i<=half;i+=2)
if(n%i==0)return 0;
return 1;
}
int main(){
scanf("%d",&x);
puts(judge(x)?"YES":"NO");
}
2.【Miller Rabin】
#include<cstdio>
#define LL long long
#define Re register LL
LL n,prime[10]={2,3,5,7,11,13,17,19,23,29};
inline LL cf(Re x,Re k,Re P){
Re s=0;
while(k){
if(k&1)(s+=x)%=P;
(x<<=1)%=P,k>>=1;
}
return s%P;
}
inline LL mi(Re x,Re k,Re P){
Re s=1;
while(k){
if(k&1)s=cf(s,x,P)%P;
x=cf(x,x,P)%P,k>>=1;
}
return s%P;
}
inline bool Miller_Rabin(Re x){
Re s=0,t=x-1;
if(x==2)return 1;
if(x<2||!(x&1))return 0;
while(!(t&1))s++,t>>=1;
for(Re i=0;i<10&&prime[i]<x;++i){
Re a=prime[i],b=mi(a,t,x);
for(Re j=1,k;j<=s;++j){
k=cf(b,b,x);
if(k==1&&b!=1&&b!=x-1)return 0;
b=k;
}
if(b!=1)return 0;
}
return 1;
}
int main(){
while(scanf("%lld",&n)!=EOF)puts(Miller_Rabin(n)?"Y":"N");
}
【模板】\(\text{PON - Prime or Not [SP288]}\)
3.【埃氏筛】
#include<algorithm>
#include<cstdio>
#define Re register int
using namespace std;
const int N=1e7+3;
int n,x,T,cnt,pri[N/3];bool pan[N];
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
inline void get_pri(Re N){
pan[0]=pan[1]=1;
for(Re i=2;i<=N;i++)
if(!pan[i]){
pri[++cnt]=i;
for(Re j=i;j<=N/i;j++)pan[i*j]=1;
}
}
int main(){
in(n),in(T),get_pri(n);
while(T--)in(x),puts(pan[x]?"No":"Yes");
}
4.【欧拉筛】
#include<cstdio>
#include<ctime>
const int N=1e6;
int cnt,i,j,pri[N/3];bool pan[N+1];
inline void get_pri(){
pan[0]=pan[1]=1;
for(i=2;i<=N;++i){
if(!pan[i])pri[++cnt]=i;
for(j=1;j<=cnt&&i*pri[j]<=N;++j){
pan[i*pri[j]]=1;
if(i%pri[j]==0)break;
}
}
}
int main(){
get_pri();
for(j=0,i=1;i<=cnt;++i)printf("%d ",pri[i]);
}
二：【约数】
1.【欧几里得 (gcd)】
#include<algorithm>
#include<cstdio>
int a,b;
inline int gcd(Re a,Re b){return (!b)?a:(a%b==0?b:gcd(b,a%b));}
//来自某大佬无比魔性的诡异写法（至今未看懂）：
inline int GCD(int a,int b){while(b^=a^=b^=a%=b);return a;}
int main(){
while(scanf("%d%d",&a,&b))
printf("%d\n",gcd(a,b));
}
2.【欧拉函数（埃氏筛）】
【模板】\(\text{ETF - Euler Totient Function}\) \(\text{[SP4141]}\)
#include<cstdio>
const int N=1e6+10;
int cnt,n=N-10,phi[N],pan[N],pri[N/3];
inline void get_phi(){
pan[0]=pan[1]=1,phi[1]=1;
for(int i=2;i<=n;i++)if(!phi[i]){
for(int j=i;j<=n;j+=i){
if(!phi[j])phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}
}
int main(){
get_phi();
for(int i=1;i<=n;++i)printf("%d: %d\n",i,phi[i]);
}
3.【欧拉函数（欧拉筛）】
【模板】\(\text{ETF - Euler Totient Function}\) \(\text{[SP4141]}\)
#include<cstdio>
const int N=1e6;
int cnt,i,j,pan[N+1],pri[N/3],phi[N+1];
inline void get_phi(){
pan[0]=pan[1]=1,phi[1]=1;
for(i=2;i<=N;++i){
if(!pan[i])pri[++cnt]=i,phi[i]=i-1;
for(j=1;j<=cnt&&i*pri[j]<=N;++j){
pan[i*pri[j]]=1;
if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
else{phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
}
int main(){
get_phi();
for(i=1;i<=N;++i)printf("%d: %d\n",i,phi[i]);
}
三：【同余】
1.【逆元】
(1).【递推】
#include<cstdio>
#define LL long long
#define Re register int
const int N=3e6+3;
int n,P,inv[N];
int main(){
// freopen("123.txt","r",stdin);
scanf("%d%d",&n,&P);
inv[1]=1;
for(Re i=2;i<=n;++i)inv[i]=(LL)(P-P/i)*inv[P%i]%P;
for(Re i=1;i<=n;++i)printf("%d\n",inv[i]);
}
(2).【快速幂】
这道题卡快速幂，会 \(\text{TLE}\) 两个点。。。。
#include<cstdio>
#define LL long long
#define Re register int
int n,P;
inline int inv(Re x,Re P){
Re s=1,k=P-2;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s%P;
}
int main(){
scanf("%d%d",&n,&P);
for(Re i=1;i<=n;++i)printf("%d\n",inv(i,P));
}
3.【扩展欧拉定理】
【模板】扩展欧拉定理（\(a^b \mod m\)） \([U55950]\)
#include<cstdio>
#include<cctype>
#define LL long long
#define Re register LL
LL i,a,b,m,x,s,phi,flag;char c;
inline LL cf(Re x,Re k){
Re s=0;
while(k){
if(k&1)(s+=x)%=m;
(x<<=1)%=m,k>>=1;
}
return s;
}
int main(){
scanf("%lld%lld",&a,&m);phi=x=m;
for(i=2;i*i<=m;++i)
if(!(x%i)){phi-=phi/i;while(!(x%i))x/=i;}
if(x>1)phi-=phi/x;
while(!isdigit(c=getchar()));
while(isdigit(c))flag|=((b=(b<<1)+(b<<3)+(c^48))>=phi),b%=phi,c=getchar();
b+=flag?phi:0;x=1;
while(b){if(b&1)x=cf(x,a)%m;a=cf(a,a)%m,b>>=1;}
printf("%lld",x);
}
4.【扩展欧几里得 (Exgcd)】
#include<cstdio>
#define Re register int
int a,b,x,y;
inline void exgcd(Re a,Re b,Re &x,Re &y){
if(!b){x=1,y=0;return;}
exgcd(b,a%b,x,y);Re x0=x,y0=y;
x=y0,y=x0-a/b*y0;return;
}
int main(){
scanf("%d%d",&a,&b),exgcd(a,b,x,y);
printf("%d\n",(x%b+b)%b);
}
5.【中国剩余定理 (CRT)】
【模板】猜数字 \(\text{[TJOI2009] [P3868]}\)
#include<cstdio>
#define LL long long
#define Re register LL
const int N=13;
LL n,ans,lcm=1,a[N],m[N];
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
inline void exgcd(Re a,Re b,Re &x,Re &y){
if(!b){x=1,y=0;return;}
exgcd(b,a%b,x,y);Re x0=x,y0=y;
x=y0,y=x0-a/b*y0;return;
}
inline LL cf(Re x,Re k,Re P){
Re s=0;
while(k){
if(k&1)(s+=x)%=P;
(x+=x)%=P,k>>=1;
}
return s;
}
int main(){
// freopen("123.txt","r",stdin);
in(n);
for(Re i=1;i<=n;++i)in(m[i]),in(a[i]),lcm*=m[i];
for(Re i=1;i<=n;++i){
Re x,y,M=lcm/m[i];exgcd(M,m[i],x,y);//M[i]*inv+m[i]*y=1
x=(x%m[i]+m[i])%m[i],(ans+=cf(cf(a[i],M,lcm),x,lcm))%=lcm;
}
printf("%lld\n",ans);
}
6.【扩展中国剩余定理 (EXCRT)】
【模板】扩展中国剩余定理（\(\text{EXCRT}\)）\(\text{[P4777]}\)
#include<cstdio>
#define LL long long
#define Re register LL
const int N=1e5+3;
LL n,ans,lcm=1,a[N],m[N];
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
inline LL exgcd(Re a,Re b,Re &x,Re &y){
if(!b){x=1,y=0;return a;}
Re gcd=exgcd(b,a%b,x,y);Re x0=x,y0=y;
x=y0,y=x0-a/b*y0;return gcd;
}
inline LL cf(Re x,Re k,Re P){
Re s=0;
while(k){
if(k&1)(s+=x)%=P;
(x+=x)%=P,k>>=1;
}
return s;
}
int main(){
// freopen("123.txt","r",stdin);
in(n);
for(Re i=1;i<=n;++i)in(m[i]),in(a[i]),lcm*=m[i];
ans=a[1]%m[1],lcm=m[1];
for(Re i=2;i<=n;++i){//ans+lcm*k=a[i](mod m[i]) -> lcm*k=a[i]-ans(mod m[i]) -> lcm*k+m[i]*y=a[i]-ans
Re c=((a[i]-ans)%m[i]+m[i])%m[i],x,y;
Re gcd=exgcd(lcm,m[i],x,y);
Re k=cf(x,c/gcd,m[i]/gcd);
ans+=lcm*k,lcm*=m[i]/gcd,ans=(ans%lcm+lcm)%lcm;
}
printf("%lld\n",ans);
}
7.【拔山盖世 / 大步小步 / BSGS (Baby Steps Giant Steps)】
【模板】\(\text{BSGS}\) / 可爱的质数 \(\text{[TJOI2007] [P3846]}\)
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<map>
#define LL long long
#define Re register int
using namespace std;
const int N=1e5+3;
int x,y,m,P;map<int,int>B;
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
//x^k = y (mod P)
//k = am-b (m=sqrt(P)+1, a\in[1,m], b\in[0,m-1])
//(x^m)^a = y x^b (mod P)
inline int mi(Re x,Re k){
Re s=1;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s;
}
int main(){
// freopen("123.txt","r",stdin);
in(P),in(x),in(y),m=sqrt(P)+1;
if(y==1)return !puts("0");
for(Re b=0,s=y;b<m;++b)B[s]=b,s=(LL)s*x%P;
Re tmp=mi(x,m);
for(Re a=1,s=1;a<=m;++a){
s=(LL)s*tmp%P;
if(B.find(s)!=B.end()){printf("%d\n",a*m-B[s]);return 0;}
}
puts("no solution");
}
8.【扩展 BSGS (EXBSGS)】
【模板】扩展 \(\text{BSGS [P4195]}\)
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<map>
#define LL long long
#define Re register int
using namespace std;
int x,y,z;map<int,int>B;
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
inline int mi(Re x,Re k,Re P){
Re s=1;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s;
}
inline int gcd(Re a,Re b){return !b?a:gcd(b,a%b);}
inline void exBSGS(Re x,Re y,Re P){
if(y==1){puts("0");return;}
Re d=gcd(x,P),tmp=1,t=0;
while(d!=1){
if(y%d){puts("No Solution");return;}
++t,y/=d,P/=d,tmp=(LL)tmp*(x/d)%P;
if(y==tmp){printf("%d\n",t);return;}
d=gcd(x,P);
}
Re m=sqrt(P)+1;B.clear();
for(Re b=0,s=y;b<m;++b)B[s]=b,s=(LL)s*x%P;
Re s=tmp;tmp=mi(x,m,P);
for(Re a=1;a<=m;++a){
s=(LL)s*tmp%P;
if(B.find(s)!=B.end()){printf("%d\n",a*m-B[s]+t);return;}
}
puts("No Solution");
}
int main(){
// freopen("123.txt","r",stdin);
while(1){
in(x),in(y),in(z);
if(!x&&!y&&!z)break;
exBSGS(x,z,y);
}
}
9.【二次剩余】
#include<algorithm>
#include<cstring>
#include<cstdio>
#define LL long long
#define Re register int
using namespace std;
int n,T,P,W,det;
inline void in(Re &x){
int f=0;x=0;char c=getchar();
while(c<'0'||c>'9')f|=c=='-',c=getchar();
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
x=f?-x:x;
}
inline int mi(Re x,Re k){
Re s=1;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s;
}
struct CP{
int x,y;CP(Re X=0,Re Y=0){x=X,y=Y;}
inline CP operator*(const CP &O)const{return CP(((LL)x*O.x%P+(LL)y*O.y%P*W%P)%P,((LL)x*O.y%P+(LL)y*O.x%P)%P);}
inline CP operator*=(const CP &O){return *this=*this*O;}
};
inline CP mi_(CP x,Re k){
CP s=CP(1,0);
while(k){
if(k&1)s*=x;
x*=x,k>>=1;
}
return s;
}
inline LL Sqrt(Re n){
if(!n)return 0;Re x;
while(1){
x=rand()%P,W=((LL)x*x%P-n+P)%P;
if(mi(W,(P-1)/2)==P-1)break;
}
return mi_(CP(x,1),(P+1)/2).x;
}
int main(){
// freopen("123.txt","r",stdin);
in(T);
while(T--){
in(n),in(P);
if(mi(n,(P-1)/2)==P-1)puts("Hola!");//无解
else{
det=Sqrt(n);
if(det)printf("%d %d\n",min(det,P-det),max(det,P-det));//两解
else printf("%d\n",det);//单解
}
}
}
10.【N 次剩余】
【模板】\(\text{N}\) 次剩余 \(\text{[P5668]}\)
还不会，先咕着。
四：【类欧几里得算法】
【模板】类欧几里得算法 \(\text{[P5170]}\)
还不会，先咕着。
五：【各类反演及数论筛法】
1.【莫比乌斯函数(Mobius)】
(1).【埃氏筛】
#include<cstdio>
const int N=1e7+3;
int i,j,n,cnt,pan[N],miu[N];
inline void get_miu(){
for(i=1;i<=n;++i)miu[i]=1;
for(i=2;i<=n;++i)
if(!pan[i]){
miu[i]=-1;
for(j=i<<1;j<=n;j+=i)pan[j]=1,miu[j]*=(j/i%i==0)?0:-1;
}
}
int main(){
n=N;
get_miu();
for(i=1;i<=n;++i)printf("%d: %d\n",i,miu[i]);
}
(2).【线性筛】
#include<cstdio>
const int N=1e7+3;
int i,j,n,cnt,pri[N>>1],pan[N],miu[N];
inline void get_miu(){
pan[1]=1,miu[1]=1;
for(i=2;i<=n;++i){
if(!pan[i])pri[++cnt]=i,miu[i]=-1;
for(j=1;j<=cnt&&i*pri[j]<=n;++j){
pan[i*pri[j]]=1;
if(i%pri[j])miu[i*pri[j]]=-miu[i];
else{miu[i*pri[j]]=0;break;}
}
}
}
int main(){
n=N;
get_miu();
for(i=1;i<=n;++i)printf("%d: %d\n",i,miu[i]);
}
2.【杜教筛】
(1).【map 记忆化】
【模板】 杜教筛（\(\text{Sum}\)）\(\text{[P4213]}\)
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<map>
#define Re register int
#define LL long long
using namespace std;
const int N23=1664511+3,N=2147483647;
int x,T,cnt,pan[N23],pri[N23],phi[N23],miu[N23];
int TO,Smiu[N23];LL Sphi[N23];
map<int,int>Smiu_;map<int,LL>Sphi_;
inline void in(Re &x){
Re fu=0;x=0;char ch=getchar();
while(ch<'0'||ch>'9')fu|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=fu?-x:x;
}
inline void get_(Re N){
miu[1]=phi[1]=pan[1]=1;
for(Re i=2;i<=N;++i){
if(!pan[i])pri[++cnt]=i,miu[i]=-1,phi[i]=i-1;
for(Re j=1;j<=cnt&&pri[j]<=N/i;++j){
pan[i*pri[j]]=1;
if(i%pri[j])miu[i*pri[j]]=-miu[i],phi[i*pri[j]]=phi[i]*phi[pri[j]];
else{miu[i*pri[j]]=0,phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
for(Re i=1;i<=N;++i)Smiu[i]=Smiu[i-1]+miu[i],Sphi[i]=Sphi[i-1]+phi[i];
}
inline LL get_Sphi(Re n){
if(n<=TO)return Sphi[n];
if(Sphi_[n])return Sphi_[n];
LL ans=(LL)n*(n+1)>>1;
for(Re l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(LL)(r-l+1)*get_Sphi(n/l);
}
return Sphi_[n]=ans;
}
inline int get_Smiu(Re n){
if(n<=TO)return Smiu[n];
if(Smiu_[n])return Smiu_[n];
Re ans=1;
for(Re l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(r-l+1)*get_Smiu(n/l);
}
return Smiu_[n]=ans;
}
int main(){
// freopen("123.txt","r",stdin);
in(T),get_(TO=N23-3);
while(T--)in(x),printf("%lld %d\n",get_Sphi(x),get_Smiu(x));
}
(2).【数组记忆化】
【模板】 杜教筛（\(\text{Sum}\)）\(\text{[P4213]}\)
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define Re register int
#define LL long long
using namespace std;
const int N23=1664511+3,N13=1303,N=2147483647;
int x,T,cnt,pan[N23],pri[N23],phi[N23],miu[N23];
int TO,vis1[N13],vis2[N13],Smiu[N23],Smiu_[N13];
LL Sphi[N23],Sphi_[N13];
inline void in(Re &x){
Re fu=0;x=0;char ch=getchar();
while(ch<'0'||ch>'9')fu|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=fu?-x:x;
}
inline void get_(Re N){
miu[1]=phi[1]=pan[1]=1;
for(Re i=2;i<=N;++i){
if(!pan[i])pri[++cnt]=i,miu[i]=-1,phi[i]=i-1;
for(Re j=1;j<=cnt&&pri[j]<=N/i;++j){
pan[i*pri[j]]=1;
if(i%pri[j])miu[i*pri[j]]=-miu[i],phi[i*pri[j]]=phi[i]*phi[pri[j]];
else{miu[i*pri[j]]=0,phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
for(Re i=1;i<=N;++i)Smiu[i]=Smiu[i-1]+miu[i],Sphi[i]=Sphi[i-1]+phi[i];
}
inline LL get_Sphi(Re n){
if(n<=TO)return Sphi[n];
Re nn=N/n;
if(vis1[nn])return Sphi_[nn];
LL ans=(LL)n*(n+1)>>1;
for(Re l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(LL)(r-l+1)*get_Sphi(n/l);
}
vis1[nn]=1;
return Sphi_[nn]=ans;
}
inline int get_Smiu(Re n){
if(n<=TO)return Smiu[n];
Re nn=N/n;
if(vis2[nn])return Smiu_[nn];
Re ans=1;
for(Re l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(r-l+1)*get_Smiu(n/l);
}
vis2[nn]=1;
return Smiu_[nn]=ans;
}
int main(){
// freopen("123.txt","r",stdin);
in(T),get_(TO=N23-3);
while(T--){
in(x),printf("%lld %d\n",get_Sphi(x),get_Smiu(x));
memset(vis1,0,sizeof(vis1));//注意对于每次询问的n，n/i都不同
memset(vis2,0,sizeof(vis2));//所以要清空
}
}
3.【最值反演 (Min-Max容斥)】
【模板】 \(\text{Card Collector [Hdu4336]}\)
#include<algorithm>
#include<cstring>
#include<cstdio>
#define LD double
#define LL long long
#define Re register int
using namespace std;
const int N=23,M=1048576+3;
int n,V,cnt[M];LD ans,p[N],Min[M];
inline void in(Re &x){
int f=0;x=0;char ch=getchar();
while(ch<'0'||ch>'9')f|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=f?-x:x;
}
int main(){
// freopen("123.txt","r",stdin);
while(~scanf("%d",&n)){
V=(1<<n)-1,ans=0;
for(Re i=1;i<=n;++i)scanf("%lf",&p[i]);
for(Re s=1;s<=V;++s){
Min[s]=0,cnt[s]=cnt[s>>1]+(s&1);
for(Re i=1;i<=n;++i)if(s&(1<<i-1))Min[s]+=p[i];
Min[s]=1.0/Min[s];
}
for(Re t=1;t<=V;++t)ans+=(cnt[t]&1)?Min[t]:-Min[t];
printf("%lf\n",ans);
}
}
