代码实现一些数学操作

· · 个人记录

又是一发学习笔记。

这次和以往有点不同，以前的学习笔记是纯数学干货，这次是代码。

备注：exgcd在Chapter III “快速”类别里。不要找不到就说这里没有了然后骂街

Chapter I 筛法

埃筛 

#define maxn 1000000
bool isp[maxn+1];
void eratos(int n){
    memset(isp,true,sizeof(isp));
    isp[0]=isp[1]=false;
    for(int i=2;i<=n;i++){
        if(isp[i]){
            for(int j=i*i;j<=n;j+=i) isp[j]=false;
        }
    }
}

欧拉筛 

bool isp[10001];int p[10001];
void Eulersieve(int n){
    memset(isp,true,sizeof(isp));
    p[0]=0;
    for(int i=2;i<=n;i++){
        if(isp[i]) p[++p[0]]=i;
        for(int j=1;j<=p[0]&&i*p[j]<=n;j++){
            isp[i*p[j]]=false;
            if(i%p[j]==0) break;
        }
    }
}

Chapter III 欧拉函数

欧拉函数 

int Eulerphi(int n){
    int res=n;
    for(int i=2;i*i<=n;i++){
        if(n%i==0){
            res=res/i*(i-1);
            for(;n%i==0;n/=i);
        }
    }
    if(n!=1) res=res/n*(n-1);
    return res;
}

范围求欧拉函数 

int euler[2<<30];
void EulerphiPLUS(int n){
    for(int i=0;i<n;i++) euler[i]=i;
    for(int i=2;i<n;i++){
        if(euler[i]==i){
            for(int j=i;j<n;j+=i) euler[j]=euler[j]/i*(i-1);
        }
    }
}

Chapter III 快速

拓展gcd

由于快速求出一元二次方程，所以也归类在“快速”里

int exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }else{
        int d=exgcd(b,a%b,y,x);
        y-=a/b*x;
        return d;
    }
}

快速幂


void ksm(int a,int b,int p){
    int ans=1;
    while(b!=0){
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans%p;
}

Chapter IV 逆元

线性逆元

注意：需要p为质数

int inv[1000000];
void getinv(int x,int p){
    inv[1]=1;
    for(int i=2;i<=x;i++)
        inv[i]=(p-p/i)*(inv[p%i])%p;
}

利用exgcd求逆元，不需要m为质数

前置：exgcd
int get_inverse(int a,int m){
    int x,y,d;
    d=exgcd(a,m,x,y);
    if(d!=1) return 0;
    else return (x+m)%m;
}

阶乘逆元

前置：快速幂

注意：需要p为质数

#define N 10005
int fac[N];
int finv[N];//阶乘的逆元
void init(int n,int p){
    fac[0]=1;
    for(int i=1;i<=n;i++){
        fac[i]=fac[i-1]*i%p;
    }
    finv[n]=ksm(fac[n],p-2,p);
    for(int i=n;i>=1;i--){
        finv[i-1]=finv[i]*i%p;
    }
}

黑科技阶乘逆元

同样需要让p为质数

finv[1]=1;
for(int i=2;i<=n;i++){
    finv[i]=(mod-mod/i)*finv[mod%i]%mod;
}
fac[0]=finv[0]=1;
for(int i=1;i<=n;i++){
    fac[i]=fac[i-1]*i%mod;
    finv[i]=finv[i]*finv[i-1]%mod;
}

本质上是线性的逆元，直接乘起来。

Chapter V 组合

p为质数时求组合取模

前置：快速幂
#define lll long long
lll C(lll a,lll b,lll p){
    if(a==b||b==0)
        return 1;
    if(b>a-b){
        b=a-b;
    }
    lll x=1;
    lll y=1;
    for(lll i=(a-b+1);i<=a;i++){
        x=x*i%p;
    }
    for(lll i=1;i<=b;i++){
        y=y*i%p;
    }
    return (x*ksm(y,p-2,p))%p;
}

无预处理组合数

前置：快速幂

注意：需要p为质数

#define lll long long
lll C(lll a,lll b,lll p){
    if(a==b||b==0)
        return 1;
    if(b>a-b){
        b=a-b;
    }
    lll x=1;
    lll y=1;
    for(lll i=(a-b+1);i<=a;i++){
        x=x*i%p;
    }
    for(lll i=1;i<=b;i++){
        y=y*i%p;
    }
    return (x*ksm(y,p-2,p))%p;
}

预处理组合数

注意：x,y 要很小

#define mod 998244353
int q;int n,m,k; 
int c[N][N];
void init(int x,int y){
    for(int i=0;i<=a;i++){
        for(int j=0;j<=b;j++){
            if(!j)
                c[i][j]=1;
            else
                c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
        }
    }
}

lucas定理（原）

前置：组合数

注意：需要p为质数

LL Lucas(LL n,LL m,LL p){
    if(n<m)
        return 0;
    if(!n)
        return 1;
    return Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}

中国剩余定理

前置：快速幂

注意：需要p为质数

int b[N];//模数 
void crt(int n){
    int val;
    fo(i,1,n)
        val=(val+a[i]*(mod/b[i])%mod*ksm(mod/b[i],b[i]-2,b[i]))%mod;
    return val;
}