小数论专题basic(待更新)

扩展GCD

给定 a,b a , b ,求使得 ax+by=gcd(a,b) a x + b y = g c d ( a , b ) 的一组解(x,y)。
我们知道, gcd(a,b)=gcd(b,a%b) g c d ( a , b ) = g c d ( b , a % b )
所以通过辗转相除法,可以得到这些方程:
a=bm+n,m=a div b,n=a%b a = b ∗ m + n , m = a   d i v   b , n = a % b
m=nm+n m = n ∗ m ′ + n ′
通过这样的轮换,最后得到 gcd(a,b)x+0y=gcd(a,b) g c d ( a , b ) ∗ x + 0 ∗ y = g c d ( a , b ) ,返回 x=1,y=0 x = 1 , y = 0
ax+by=gcd(a,b) a x + b y = g c d ( a , b )
bx+a%by=gcd(a,b) b x ′ + a % b ∗ y ′ = g c d ( a , b )
合并同类项,最后得出 x=y,y=xyab x = y ′ , y = x ′ − y ′ ∗ ⌊ a b ⌋
自己推一下式子。

int extgcd(int a,int b,int &x,int &y){
    if (!b){
        x=1;y=0;
        return a;
    }
    int c=extgcd(b,a%b,x,y),t;
    t=x;
    x=y;
    y=t-y*(a/b);
    return c;
}

矩阵快速幂

差点忘了矩乘怎么打…….

struct matrix{
    int mat[N][N];
    void clear(){
        memset(mat,0,sizeof(mat));
    }
    void unit(){
        clear();
        for(int i=1;i<=n;i++)mat[i][i]=1;
    }
};matrix a,b,c;
matrix mul(matrix a,matrix b){
    matrix c;c.clear();
    for(int k=1;k<=n;k++)
        for(int i=1;i<=n;i++)if(a.mat[i][k])//这个优化很强大
            for(int j=1;j<=n;i++)if(b.mat[k][j])
                 c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
    return c;
}
matrix mksm(matrix a,int p){
    matrix b;b.unit();
    for(;p;p/=2,a=mul(a,a))
        if(p&1) b=mul(b,a);//注意顺序,矩乘没有交换律
    return b;
}

线筛求φ(N)和μ(N)

通过线筛质数来求。

void getprime(){
    int i,j;
    for(i=2;i<=N;i++){
        if(!bz[i]){
            prime[++prime[0]]=i;
            phi[i]=i-1;
            mu[i]=-1;
        }
        fo(j,1,prime[0]){
            if(i*prime[j]>N)break;
            bz[i*prime[j]]=0;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            } else{
                mu[i*prime[j]]=-mu[i];
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }
}

中国剩余定理

就是给定这么几个方程:
Aa1(mod p1) A ≡ a 1 ( m o d   p 1 )
Aa2(mod p2) A ≡ a 2 ( m o d   p 2 )
……
Aan(mod pn) A ≡ a n ( m o d   p n )
求最小的A。(一般任意两个p之间互质,如果不互质自己对号入座)
考虑 A=p1k1+a1=p2k2+a2 A = p 1 ∗ k 1 + a 1 = p 2 ∗ k 2 + a 2
移项,得 p1k1p2k2=a2a1 p 1 ∗ k 1 − p 2 ∗ k 2 = a 2 − a 1 ,即 p1k1a2a1+p2k2a2a1=1 p 1 ∗ k 1 a 2 − a 1 + p 2 ∗ − k 2 a 2 − a 1 = 1 ,那么用 extgcd e x t g c d 求出 k1a2a1 k 1 a 2 − a 1 k2a2a1 − k 2 a 2 − a 1 就大功告成,合并了两个同余方程。以此类推,将n个同余方程合并到一个的时候,所得到的 an a n ′ 就是最小的A。

void get(){
    int i,ans,last,l1,p1,p2,a1,a2;
    last=a[1];l1=p[1];
    fo(i,2,p[0]){
        p1=l1,p2=p[i];
        a1=last,a2=a[i];
        extgcd(p1,p2,k1,k2);
        k1=k1*(a2-a1);
        l1=l1*p2;
        last=k1*p1+a1;
    }
    ans=(last%P+P)%P;//P为lcm(p1,p2,...,pn)
    return ans;
}

Lucas定理

C(n,m)%p C ( n , m ) % p
Lucas(n,m,p)=Lucas(np,mp,p)C(n%p,m%p,p) L u c a s ( n , m , p ) = L u c a s ( ⌊ n p ⌋ , ⌊ m p ⌋ , p ) ∗ C ( n % p , m % p , p )
证明:设 n=ap+b,m=cp+d n = a p + b , m = c p + d ,现证明 Ccp+dap+b=CcaCdb%p) C a p + b c p + d = C a c ∗ C b d % p )
xZ ∀ x ∈ Z (1+x)ap+b=Σap+bi=0xiCin ( 1 + x ) a p + b = Σ i = 0 a p + b x i ∗ C n i
(1+x)p1+xp (mod p) ( 1 + x ) p ≡ 1 + x p   ( m o d   p )
(1+x)ap+b[(1+x)p]a(1+x)b(1+xp)a(1+x)b[Σai=0Ciaxip][Σbi=0Cibxi] ( 1 + x ) a p + b ≡ [ ( 1 + x ) p ] a ∗ ( 1 + x ) b ≡ ( 1 + x p ) a ∗ ( 1 + x ) b ≡ [ Σ i = 0 a C a i ∗ x i p ] ∗ [ Σ i = 0 b C b i ∗ x i ]
那么二项式中 xcp+dCcaCdbCcp+dap+b(modp) x c p + d 项 的 系 数 ≡ C a c ∗ C b d ≡ C a p + b c p + d ( m o d p )
∴得证。

(LL)即为long long
LL lucas(LL n,LL m,LL p){
    if(!m)return 1;
    return lucas(n/p,m/p,p)*C(n%p,m%p);
}

高斯消元

运用加减消元法进行消元,一次消一个,总复杂度 O(n3) O ( n 3 )

void gauss(){
    int i,j,k;
    double temp;
    fo(i,1,n){
        j=i;
        while(!mat[j][i]&&j<=n)j++;
        if(j==n+1)continue;
        if(j!=i)fo(k,1,n+1)swap(mat[i][k],mat[j][k]);
        fo(j,1,n)
            if(j!=i && mat[j][i]){
                temp=mat[i][i]/mat[j][i];
                fo(k,1,n+1)mat[j][k]=mat[j][k]*temp-mat[i][k];
            }
    }
    fd(i,n,1)
        if(mat[i][i]!=0){
            fo(j,i+1,n)mat[i][n]-=mat[i][j]*x[j];
            x[i]=mat[i][n]/mat[i][i];
        }
}

FFT

这里就不详细讲了。
虚数的定义。

struct Cpx{
    DB r,i;
    Cpx(){}
    Cpx(DB R,DB I){r=R,i=I;}
    Cpx operator +(const Cpx&x)const{return Cpx(r+x.r,i+x.i);}
    Cpx operator -(const Cpx &x)const{return Cpx(r-x.r,i-x.i);}
    Cpx operator *(const Cpx &x)const{return Cpx(r*x.r-i*x.i,r*x.i+i*x.r);}
};
void FFT(Cpx *y,int opz){
    int i,h,k;
    fo(i,0,m-1)if(i<Rev[i])swap(y[i],y[Rev[i]]);
    for(h=2;h<=m;h<<=1){
        Cpx wn(cos(2*Pi/h),opz*sin(2*Pi/h));
        for(i=0;i<m;i+=h){
            Cpx w(1,0);
            fo(k,i,i+(h>>1)-1){
                Cpx u=y[k],t=w*y[k+(h>>1)];
                y[k]=u+t;
                y[k+(h>>1)]=u-t;
                w=w*wn;
            }
        }
    }
    if(opz==-1){
        fo(i,0,m-1)y[i].r/=m;
    }
}

fo(i,0,m-1)Rev[i]=Rev[i>>1]>>1|(i&1)<<len-1;//求i的反过来的二进制。
//意思是,假设目前要加一个位,已经知道不加位的反过来的。如果最后一位是0,反过来后第一位是0.
//最后一位是1同理。

NTT(有原根版)

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#define N 200010
#define M 525000 
#define LL long long
#define mo 924844033
#define P(a) putchar(a)
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
LL i,j,k,l,n,m,ans;
LL u,v,L,len; 
LL jc[N],ny[N],fa[N];
LL siz[N],f[N],h[N];
LL f1[M],f2[M];
LL w[M],_w[M],Rev[M];
LL Ans[N];
LL read(){
    LL fh=1,rs=0;char ch;
    while((ch<'0'||ch>'9')&&(ch^'-'))ch=getchar();
    if(ch=='-')fh=-1,ch=getchar();
    while(ch>='0'&&ch<='9')rs=(rs<<3)+(rs<<1)+(ch^'0'),ch=getchar();
    return fh*rs;
}
LL ksm(LL x,LL y){
    LL rs=1;
    for(;y;y>>=1,x=(x*x)%mo)
    if(y&1)rs=(rs*x)%mo;
    return rs;
}
void pre_(){
    int i;
    jc[0]=jc[1]=ny[0]=ny[1]=1;
    fo(i,2,N-10)jc[i]=(jc[i-1]*i)%mo;
    ny[N-10]=ksm(jc[N-10],mo-2);
    fd(i,N-11,2)ny[i]=(ny[i+1]*(i+1))%mo; 
}
void NTT(LL *y,LL *w){
    LL i,h,k;
    fo(i,0,L-1)if(i<Rev[i])swap(y[i],y[Rev[i]]);
    for(h=2;h<=L;h<<=1){
        for(i=0;i<L;i+=h){
            fo(k,i,i+(h>>1)-1){
                LL u=y[k],t=(w[L/h*(k-i)]*y[k+(h>>1)])%mo;
                y[k]=(u+t)%mo;
                y[k+(h>>1)]=((u-t)%mo+mo)%mo;
            }
        }
    }
}
int main(){
    n=read();
    //...
    len=0;
    for(L=1;L<(n+1)*2;L<<=1)len++;
    w[0]=w[L]=1;
    w[1]=ksm(5,(mo-1)/L);
    fo(i,1,L-1) w[i]=w[i-1]*1ll*w[1]%mo;
    fo(i,0,L) _w[i]=w[L-i];
    fo(i,0,L-1)Rev[i]=Rev[i>>1]>>1|(i&1)<<len-1;
    fo(i,0,n)f1[i]=f[n-i];
    fo(i,n+1,L-1)f1[i]=0;
    fo(i,0,n)f2[i]=h[i];
    fo(i,n+1,L-1)f2[i]=0;
    NTT(f1,w);
    NTT(f2,w);
    fo(i,0,L-1)f1[i]=f1[i]*f2[i]%mo;
    NTT(f1,_w);
    LL inv=ksm(L,mo-2);
    fo(i,0,L-1)f1[i]=(f1[i]*inv)%mo;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值