# [template integration plan] NB number theory

## [prime number]

### [violent sentence]

Prime, コ ン テ ス, prime \ ([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");
}```

### [Miller Rabin]

Prime, コ ン テ ス, prime \ ([AT807] \)

```#include<cstdio>
#define LL long long
int prime[10]={2,3,5,7,11,13,17,19,23,29};
inline int slow_cf(int a,int b,int mod){
int x=0;
while(b){
if(b&1)x=(x+a)%mod;
a=a*2%mod;b>>=1;
}
return x%mod;
}
inline int quick_pow(int x,int k,int mod){
int s=1;
while(k){
if(k&1)s=slow_cf(s,x,mod)%mod;
x=slow_cf(x,x,mod)%mod;k>>=1;
}
return s%mod;
}
inline bool Miller_Rabin(int x){
int i,j,k,s=0,t=x-1;
if(x==2)return 1;
if(x<2||!(x&1))return 0;
while(!(t&1))s++,t>>=1;
for(i=0;i<10&&prime[i]<x;++i){
int a=prime[i];
int b=quick_pow(a,t,x);
for(j=1;j<=s;++j){
k=slow_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(){
int n;
scanf("%d",&n);
puts(Miller_Rabin(n)?"YES":"NO");
}```

### [sieve]

```#include<cstdio>
#define R register int
int pri[10005],pan[10005],gs;
inline void creat_prime(){
for(R i=2;i<=10000;i++)
if(pan[i]==0){
pri[++gs]=i;
for(R j=2;i*j<=10000;j++)pan[i*j]=1;
}
}
inline int sakura(){
creat_prime();
for(R i=1;i<=gs;i++)printf("%d\n",pri[i]);
}
int QAQWQ=sakura();
int main(){}```

### [screen sieve]

```#include<cstdio>
#include<ctime>
const int N=1e6;
int cnt,i,j,pri[N/3];bool pan[N+1];
inline void get_pri(){
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;
}
}
int main(){
get_pri();
for(j=0,i=1;i<=cnt;++i)printf("%d ",pri[i]);
}```

## [number]

### [GCD Euclid]

```#include<algorithm>
#include<cstdio>
int a,b;
inline int gcd(int n,int m){return n%m?gcd(m,n%m):m;}
//The weird writing method from some big guy's incomparable magic (not understood yet):
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));
}```

### Euler function (wire screen)

```#include<cstdio>
const int N=1e6;
int cnt,i,j,pan[N+1],pri[N/3],phi[N+1];
inline void get_phi(){
for(i=2;i<=N;++i){
if(!pan[i])pri[++cnt]=i,phi[i]=i-1;
for(j=1;j<=cnt;++j){
if(i*pri[j]>N)break;
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]);
}```

### Euler function

```#include<cstdio>
const int N=1e6;
int cnt,n=N,phi[N/3],pan[N],pri[N];
inline void get_phi(){
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]);
}```

## [congruence]

### Inverse element (recurrence)

Multiplicative inverse \ ([P3811] \)

```#include<cstdio>
#define LL long long
int i,n,P,inv[3000003];
int main(){
scanf("%d%d",&n,&P);
inv[1]=1;
for(i=2;i<=n;++i)inv[i]=(LL)(P-P/i)*inv[P%i]%P;
for(i=1;i<=n;++i)printf("%d\n",inv[i]);
}```

### Inverse element (fast power)

```#include<cstdio>
#define LL long long
int a,P;
int quick_pow(int x,int k){
LL s=1,a=x;
while(k){
if(k&1)s=(s*a)%P;
a=(a*a)%P,k>>=1;
}
return s%P;
}
int main(){
scanf("%d%d",&a,&P);
printf("%d",quick_pow(a,P-2));
}```

### [extended Euler theorem]

Extended Euler theorem (\ (a ^ B \)% \ (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);
}```

### [EXGCD extended Euclid]

(solving linear congruence equation)

Congruence equation \ ([P1082] \)

```#include<cstdio>
#define LL long long
LL a,b,x,y,gcd;
inline LL exgcd(LL a,LL b,LL &x,LL &y){
if(!b){x=1;y=0;return a;}
gcd=exgcd(b,a%b,x,y);
int t=x;x=y;y=t-a/b*y;
return gcd;
}
int main(){
scanf("%lld%lld",&a,&b);
exgcd(a,b,x,y);
printf("%lld",(x+b)%b);
}```

### CRT China remainder theorem

(solution of Congruence Equations) (P pairwise coprime)

Guess the number \ ([TJOI2009]\) \([P3868] \)

```#include<cstdio>
#define LL long long
#define Re register LL
#define F() for(i=1;i<=k;++i)
LL i,k,x,y,fu,ans,lcm=1,a[13],m[13];char ch;
inline void in(Re &x){
fu=x=0;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 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;
}
inline void exgcd(Re &x,Re &y,Re a,Re b){
if(!b){x=1,y=0;return;}
exgcd(x,y,b,a%b);Re X=x;
x=y,y=X-a/b*y;
}
int main(){
in(k);
F()in(a[i]);
F()in(m[i]),lcm*=m[i],(a[i]+=m[i])%=m[i];
F(){
Re M=lcm/m[i];exgcd(x,y,M,m[i]);
(x+=m[i])%=m[i];(ans+=cf(cf(M,x,lcm),a[i],lcm))%=lcm;
}
printf("%lld",(ans+lcm)%lcm);
}```

### [EXCRT extends China's residual theorem]

(solution of Congruence Equations) (\ (P \) is not necessarily two coprime)

Extended Chinese Remainder Theorem (\ (EXCRT \) \ ([P4777] \)

```#include<cstdio>
#define LL long long
#define Re register LL
const int N=1e5+3;
LL a,b,c,x,y,i,n,M,fu,gcd,ans,tmp,m[N],ai[N];char ch;
inline void in(Re &x){
fu=x=0;ch=getchar();
while(ch<'0'||c>'9')fu|=ch=='-',ch=getchar();
while(ch>='0'&&c<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=fu?-x:x;
}
inline void cf(Re &s,Re k,Re P){
Re a=s;s=0;
while(k){
if(k&1)(s+=a)%=P;
(a<<=1)%=P,k>>=1;
}
}
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 t=x;x=y;y=t-a/b*y;
return gcd;
}
int main(){
in(n);
for(i=1;i<=n;++i)in(m[i]),in(ai[i]);
ans=ai[1],M=m[1];
for(i=2;i<=n;++i){
a=M,b=m[i],c=(ai[i]-ans%b+b)%b;
gcd=exgcd(a,b,x,y),tmp=b/gcd;
cf(x,c/gcd,tmp);
ans+=x*M,M*=tmp,((ans%=M)+=M)%=M;
}
printf("%lld",ans);
}```

## Matrix multiplication

### [ A × B ]

```#include<cstdio>
const int mod=1e9+7;
int i,j,k,n,p,m,a[505][505],b[505][505],c[505][505];
int main(){
scanf("%d%d%d",&n,&p,&m);
for(i=1;i<=n;i++)
for(j=1;j<=p;a[i][j]=(a[i][j]+mod)%mod,j++)
scanf("%d",&a[i][j]);
for(i=1;i<=p;i++)
for(j=1;j<=m;b[i][j]=(b[i][j]+mod)%mod,j++)
scanf("%d",&b[i][j]);
for(i=1;i<=n;i++)
for(j=1;j<=m;j++)
for(k=1;k<=p;k++)
c[i][j]=(c[i][j]+(long long)a[i][k]*b[k][j]%mod)%mod;
for(i=1;i<=n;printf("\n"),i++)
for(j=1;j<=m;j++)
printf("%d ",c[i][j]);
}```

### [sequence recursive acceleration]

Matrix acceleration (sequence) \ ([P1939] \)

(structure packing)

```//f[1]=f[2]=f[3]=1,f[n]=f[n-1]+f[n-3]
//[0 0 1]
//[1 0 0]
//[0 1 1]
#include<cstring>
#include<cstdio>
#define LL long long
const int P=1e9+7;
struct QAQ{
int a[4][4];
void CL(){memset(a,0,sizeof(a));a[1][3]=a[2][1]=a[3][2]=a[3][3]=1;}
QAQ operator * (QAQ &x){
int i,j,k;QAQ ans;memset(ans.a,0,sizeof(ans.a));
for(i=1;i<4;i++)
for(j=1;j<4;j++)
for(k=1;k<4;k++)
(ans.a[i][j]+=(LL)a[i][k]*x.a[k][j]%P)%=P;
return ans;
}
};
int ans,n,T;
int sovle(int k){
if(k<1)return 1;
QAQ s,x;s.CL();x.CL();
while(k){
if(k&1)s=(s*x);
x=(x*x);k>>=1;
}
ans=0;
for(k=1;k<4;k++)(ans+=s.a[k][2])%=P;
return ans;
}
int main(){
scanf("%d",&T);
while(T--){
scanf("%d",&n);
printf("%d\n",sovle(n-3));
}
}```

## [combination count]

### [combination number (recursive)]

```#include<cstdio>
#define R register int
#define LL long long
LL c[1005][1005],a,b;
int main(){
for(R i=0;i<=1000;i++)c[0][i]=c[i][i]=1;
for(R i=1;i<=1000;i++)
for(R j=1;j<i;j++)
c[j][i]=c[j-1][i-1]+c[j][i-1];
scanf("%lld%lld",&a,&b);
printf("%lld\n",c[a][b]);
}```

### Number of combinations

```#include<cstdio>
#define LL long long
LL C[1005][1005],a,b;
inline LL cc(LL m,LL n){
if(m>n)return 0;
if(!m||m==n)return C[m][n]=1;
if(C[m][n])return C[m][n];
return C[m][n]=cc(m,n-1)+cc(m-1,n-1);
}
int main(){
scanf("%lld%lld",&a,&b);
printf("%lld\n",cc(a,b));
}```

### Lucas theorem

```#include<cstdio>
#define LL long long
LL a,b,mod,jc[1000005];
LL niv(LL x,LL p){
LL ans=1,k=p-2;
while(k){
if(k&1)ans=(ans*x)%p;
x=(x*x)%p;k>>=1;
}
return ans%p;
}
inline LL c(LL m,LL n,LL p){
if(n<m)return 0;
if(n<p)return jc[n]*niv(jc[m],p)%p*niv(jc[n-m],p)%p;
return c(m/p,n/p,p)*c(m%p,n%p,p)%p;
}
int main(){
scanf("%lld%lld%lld",&a,&b,&mod);
jc[0]=1;for(LL i=1;i<=1000000;i++)jc[i]=jc[i-1]*i%mod;
printf("%lld\n",c(a,b,mod));
}```

## [inclusion exclusion principle and Mobius function]

### Mobius function

```#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]);
}```

### [Mobius function (wire screen))

```#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]);
}```

Tags: PHP

Posted on Sun, 03 Nov 2019 15:45:14 -0800 by CoB-Himself