[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