【数学相关】模板

更新日志

2018/9/17 绝大多数模板
2018/9/18 康托展开与康托逆展开(组合)
2018/9/25 矩阵树定理(组合)
2018/10/9 01字典树(进制)
2018/10/10 对矩阵树定理模板的修改
2018/10/11 数论分块(数论)
2018/10/12 对快速傅里叶变换模板的修改
2018/10/16 欧拉函数(数论)
2018/10/16 可持久化01字典树(进制)
2018/10/17 对斐波那契取模模板的修改
2018/10/24 快速莫比乌斯变换/子集反演(反演)
2018/11/03 发现上面那个东西出锅了，打不开，先删了，反正NOIp也不考

杂项

long long mul(long long a,long long b,long long mod)
{
    long long res=0;
    for(;b;a=(a+a)%mod,b>>=1)
        if(b&1) res=(res+a)%mod;
    return res;
}

龟速乘法

1 int FP(int a,int k,int mod)
2 {
3     int ans=1;
4     for(;k;a=1LL*a*a%mod,k>>=1)
5         if(k&1) res=1LL*res*a%mod;
6     return res;
7 }

快速幂

float InvSqrt(float x)
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f375a86- (i>>1);
    x = *(float*)&i;
    x = x*(1.5f-xhalf*x*x);
    return x;
}

^(-0.5)

数论

ll phi(ll x)
{
    ll k=x;
    ll ans=x;
    for(ll i=2;i*i<=k;i++)
    {
        if (k%i==0)
        {
             ans=ans/i*(i-1);
             while (k%i==0) k/=i;
        } 
    }
    if(k>1) ans=ans/k*(k-1);
    return ans%mod;
}

欧拉函数

 1 int p[N];
 2 bool jud[N];
 3 bool vis[N];
 4 int n,m,cnt;
 5 
 6 void prime()
 7 {
 8     int i,j;
 9     for(int i=2;i<N;i++)
10     {
11         if(!vis[i])
12         {
13             p[++cnt]=i;
14             jud[i]=true;
15         }
16         for(int j=1;j<=cnt;j++)
17         {
18             int v=i*p[j];
19             if(v>N) break;
20             vis[v]=true;
21             if(i%p[j]==0) continue;
22         }
23     }
24     return;
25 }

欧拉筛素数

int T,tmp,phi[N],p[666666],cnt=0;
bool vis[N];

void getphi()
{
    phi[1]=1;
    vis[1]=true;
    for(int i=2;i<N;i++)
    {
        if(!vis[i])
        {
            p[++cnt]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt;j++)
        {
            int v=i*p[j];
            if(v>=N) break;
            vis[v]=true;
            if(i%p[j]==0)
            {
                phi[v]=phi[i]*p[j];
                break;
            }
            else phi[v]=phi[i]*(p[j]-1);
        }
    }
}

欧拉函数线性筛

int mbs[N],p[N],n;
bool vis[N];

void mobius()
{
    int num=0;
    vis[1]=true;
    mbs[1]=sum[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!vis[i])
        {
            p[num++]=i;
            mbs[i]=-1;
        }
        for(int j=0;j<=num&&i*p[j]<=N;j++)
        {
            vis[i*p[j]]=true;
            if(i%p[j]==0)
            {
                mbs[i*p[j]]=0;
                break;
            }
            mbs[i*p[j]]=-mbs[i];
        }
    }
}

莫比乌斯函数线性筛

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

long long inv(long long a,long long n)
{
    long long d,x,y;
    exgcd(a,n,d,x,y);
    return d==1?(x+n)%n:-1;
}

单点逆元

int inv[N];
void inv(int n, int p)
{
    inv[1]=1;
    for(int i=2;i<=n;i++)
    inv[i]=1LL*(p-p/i)*inv[p%i]%p;
}

线性逆元

int gcd(int x,int y)
{
    if(y<0) y=-y;
    if(x<0) x=-x;
    return y?gcd(y,x%y):x;
}

最大公约数

    for(int i=1;i<=n;i++)
    {
        scanf("%d",&tmp);
        ans=gcd(ans,tmp);
    }

裴蜀定理

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

扩展欧几里德算法

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

int get(int a,int b,int d,int x,int y);
{
    exgcd(b,mod,d,x,y);
    if(d!=1) return -1;
    else return 1LL*a*((x%mod+mod)%mod)%mod);
}

有理数取余

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

int crt(int *a,int *w,int len)
{
    int i,d,x,y,m,n,ret;
    ret=0;
    n=1;
    for(int i=0;i<len;i++) n*=w[i];
    for(int i=0;i<len;i++)
    {
        m=n/w[i];
        d=exgcd(w[i],m,x,y);
        ret=(ret+y*m*a[i])%n;
    }
    return (n+ret%n)%n;
}

中国剩余定理

long long mul(long long a,long long b,long long mod)
{
    long long res=0;
    while(b)
    {
        if(b&1) res=(res+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return res;
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b) 
    {
        x=1;y=0;
        return a;
    }
    long long gcd=exgcd(b,a%b,x,y);
    long long tmp=x;
    x=y;
    y=tmp-a/b*y;
    return gcd;
}

long long excrt()
{
    long long x,y,k;
    long long m=b[1],ans=a[1];
    for(int i=2;i<=n;i++)
    {
        long long aa=m,bb=b[i],c=(a[i]-ans%bb+bb)%bb;
        long long gcd=exgcd(aa,bb,x,y),bg=bb/gcd;
        if(c%gcd) return -1;
        x=mul(x,c/gcd,bg);
        ans+=x*m;
        m*=bg;
        ans=(ans%m+m)%m;
    }
    return (ans%m+m)%m;
}

扩展中国剩余定理

int bsgs(int a,int b,int c,int d)
{
    int len=ceil(sqrt(c+0.5));
    int tmp=b,sum=d;
    T.clear();
    for(int i=1;i<=len;i++)
    {
        T[tmp=1LL*tmp*a%c]=i;
    }
    tmp=FP(a,len,c);
    for(int i=1;i<=len;i++)
    {
        sum=1LL*sum*tmp%c;
        if(T.count(sum))
        {
            return i*len-T[sum];
        }
    }
    return -1; 
}

BSGS

int gcd(int a,int b)
{
    if(!b) return a;
    return gcd(b,a%b);
}

int FP(int a,int k,int mod)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*a*ans%mod;
    return ans;
}

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

int bsgs(int a,int b,int c,int d)
{
    int len=ceil(sqrt(c+0.5));
    int tmp=b,sum=d;
    T.clear();
    for(int i=1;i<=len;i++)
    {
        T[tmp=1LL*tmp*a%c]=i;
    }
    tmp=FP(a,len,c);
    for(int i=1;i<=len;i++)
    {
        sum=1LL*sum*tmp%c;
        if(T.count(sum))
        {
            return i*len-T[sum];
        }
    }
    return -1; 
}

int exbsgs(int a,int b,int c)
{
    if(c==1) return 0;
    a%=c;
    b%=c;
    if(b==1) return 0;
    if(a==0)
    {
        if(b==0) return 1;
        return -1;
    }
    int d=1,cnt=0;
    for(int g=gcd(a,c);g!=1;g=gcd(a,c))
    {
        if(b%g) return -1;
        c/=g;
        b/=g;
        d=1LL*d*(a/g)%c;
        cnt++;
        if(b==d) return cnt;
    }
    a%=c;
    int ans=bsgs(a,b,c,d);
    if(ans==-1) return -1;
    else return ans+cnt;
}

扩展BSGS

bool MR(long long n)
{
    if(n==2) return true;
    if(n<2||!(n&1)) return false;
    for(int i=0;i<Times;i++)
    {
        long long a=((long long)rand())%(n-2)+2;
        long long m=n-1,k=0;
        for(;(m&1)==0;k++,m>>=1);
        long long x=FP(a,m,n);
        for(int j=1;j<=k;j++)
        {
            long long y=multi(x,x,n);
            if(y==1&&x!=1&&x!=n-1) return false;
            x=y;
        }
        if(x!=1) return false;
    }
    return true;
}

Miller-Rabin

void read(long long &x)
{
    long long f=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')f=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=x*10+s-'0';s=getchar();}
    x*=f;
}

long long getmax(long long a,long long b)
{
    if(a>b) return a;
    return b;
}

long long getmin(long long a,long long b)
{
    if(a<b) return a;
    return b;
}

long long gcd(long long a,long long b)
{
    return b?gcd(b,a%b):a;
}

long long multi(long long a,long long b,long long mod)
{
    long long ans=0;
    a%=mod;
    for(;b;b>>=1,a=(a+a)%mod)
        if(b&1) ans=(ans+a)%mod;
    return ans;
}

long long FP(long long a,long long b,long long mod)
{
    long long ans=1;
    a%=mod;
    for(;b;b>>=1,a=multi(a,a,mod))
        if(b&1) ans=multi(ans,a,mod);
    return ans;
}

bool MR(long long n)
{
    if(n==2) return true;
    if(n<2||!(n&1)) return false;
    for(int i=0;i<Times;i++)
    {
        long long a=((long long)rand())%(n-2)+2;
        long long m=n-1,k=0;
        for(;(m&1)==0;k++,m>>=1);
        long long x=FP(a,m,n);
        for(int j=1;j<=k;j++)
        {
            long long y=multi(x,x,n);
            if(y==1&&x!=1&&x!=n-1) return false;
            x=y;
        }
        if(x!=1) return false;
    }
    return true;
}

long long PR(long long n,long long c)
{
    long long i=1,k=2;
    long long x=rand()%(n-1)+1;
    long long y=x;
    while(114514)
    {
        i++; 
        x=(multi(x,x,n)+c)%n;
        long long d=gcd(y-x+n,n);
        if(1!=d&&d!=n) return d;
        if(y==x) return n;
        if(i==k) y=x,k<<=1;
    }
}

void find(long long n,int c)
{
    if(n<=1) return;
    if(MR(n))
    {
        minn=getmin(n,minn);
        return;
    }
    long long p=PR(n,c--);    
    find(n/p,c);
    find(p,c);
}

long long get(long long n)
{
    if(MR(n)) return n;
    long long k=-1;
    minn=n;
    find(n,16381);
    k=getmax(get(minn),k);
    k=getmax(get(n/minn),k);
    return k;
}

Pollard-Rho

for(int l=1,r;l<=n;l=r+1)
{
    r=n/(n/l);
    ans+=(r-l+1)*(n/l);
}

数论分块

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;

const int mod=1000000007;
long long n;
long long q=1,p=1;
long long a1=1,a2=1;

struct Matrix
{
    long long T[4][4];    
}E,P;

Matrix mul(Matrix X,Matrix Y)
{
    Matrix Z;
    memset(Z.T,0,sizeof(Z.T));
    for(register int i=1;i<=2;i++)
    for(register int j=1;j<=2;j++)
    for(register int k=1;k<=2;k++)
    {
        Z.T[i][j]=(1LL*Z.T[i][j]%mod+1LL*X.T[i][k]*Y.T[k][j]%mod)%mod;
    }
    return Z;
}

void pow(long long k)
{
    for(;k;k>>=1,E=mul(E,E))
        if(k&1) P=mul(P,E);
    return;
}

int main()
{
    scanf("%lld",&n);
    if(n<=1) return puts("1"),0;
    P.T[1][1]=a2;
    P.T[1][2]=a1;
    E.T[1][1]=p;
    E.T[2][1]=q;
    E.T[1][2]=1;
    pow(n-2);
    cout << P.T[1][1]%mod << endl;
    return 0;
}

斐波那契数取模

组合数学

int N,P,cnt;
long long ans=1;
int prime[MAXN],maxp[MAXN],c[MAXN];
bool sign[MAXN];

void add(int x,int f)
{
    while(x!=1){    
        c[maxp[x]]+=f;
        x/=prime[maxp[x]];
    }
}

int main()
{
    scanf("%d%d",&N,&P);
    for(int i=2;i<=N*2;i++)
    {
        if(!sign[i]) prime[++cnt]=i,maxp[i]=cnt;
        for(int j=1;j<=cnt;j++)
        {
            if(i*prime[j]>2*N) break;
            sign[prime[j]*i]=true;
            maxp[i*prime[j]]=j;
            if(!(i%prime[j])) break;
        }
    }
    for (int i=N+1;i<=2*N;i++) add(i,1);//?
    for (int i=1;i<=N+1;i++) add(i,-1);//?
    for(int i=1;i<=cnt;i++) while(c[i]--) ans=(ans*prime[i])%P;
    return printf("%lld\n",ans),0;
}

卡特兰数取模

int lucas(int x,int y)
{
    if(x<y) return 0;
    if(x<mod) return (1LL*b[x]*a[y]*a[x-y])%mod;
    return (1LL*lucas(x/mod,y/mod)*lucas(x%mod,y%mod))%mod;
}

卢卡斯定理

namespace Mathmatics
{
    long long FastPower(long long a,long long k,long long mod)
    {
        long long ans=1;
        for(;k;k>>=1,a=1LL*a*a%mod)
            if(k&1) ans=1LL*ans*a%mod;
        return ans;
    }
    
    long long ExGCD(long long a,long long b,long long &x,long long &y)
    {
        if(!b)
        {
            x=1;
            y=0;
            return a;
        }
        long long ans=ExGCD(b,a%b,x,y);
        long long t=x;
        x=y;
        y=t-a/b*y;
        return ans;
    }
    
    long long Inv(long long n,long long mod)
    {
        long long x,y;
        ExGCD(n,mod,x,y);
        x+=mod;
        return x>mod?x-mod:x;
    }
    
    long long CRT(long long a,long long mod)
    {
        return a*Inv(p/mod,mod)%p*(p/mod)%p; 
    }
    
    long long fac(long long n,long long pi,long long pk)
    {
        if(!n) return 1;
        long long ans=1;
        for(long long i=2;i<=pk;i++)
        {
            if(i%pi)
            {
                ans=ans*i%pk;
            }
        }
        ans=FastPower(ans,n/pk,pk);
        for(long long i=2;i<=n%pk;i++)
        {
            if(i%pi)
            {
                ans=ans*i%pk;
            }
        }
        return ans*fac(n/pi,pi,pk)%pk;
    }
        
    long long C(long long n,long long m,long long pi,long long pk)
    {
        long long a=fac(n,pi,pk);
        long long d1=fac(m,pi,pk);
        long long d2=fac(n-m,pi,pk);
        long long k=0;
        for(long long i=n;i;i/=pi) k+=i/pi;
        for(long long i=m;i;i/=pi) k-=i/pi;
        for(long long i=n-m;i;i/=pi) k-=i/pi;
        return a*Inv(d1,pk)%pk*Inv(d2,pk)%pk*FastPower(pi,k,pk)%pk;
    }

    long long ExLucas(long long n,long long m)
    {
        long long ans=0;
        long long tmp=p,pk=0;
        static int lim=sqrt(p)+5;
        for(int i=2;i<=lim;i++)
        {
            if(tmp%i==0)
            {
                pk=1;
                while(tmp%i==0)
                {
                    pk*=i;
                    tmp/=i;
                }
                ans=(ans+CRT(C(n,m,i,pk),pk))%p;
            }
        }
        if(tmp>1) ans=(ans+CRT(C(n,m,tmp,tmp),tmp))%p;
        return ans;
    }
}

扩展卢卡斯定理

void getfac()
{
    fac[0]=fac[1]=1;
    for(int i=2;i<=20;i++)
    fac[i]=fac[i-1]*i;
}

long long cantor(long long *A,int N)
{
    long long res=0;
    for(int i=0;i<N;i++)
    {
        int cnt=0;
        for(int j=i+1;j<N;j++)
        {
            if(A[j]<A[i]) cnt++;
        }
        res+=1LL*fac[N-i-1]*cnt;
    }
    return res;
}

void decantor(long long x,int N)
{
    memset(vis,0,sizeof(vis));
    for(int i=1,j=1;i<=N;i++)
    {
        long long cur=x/fac[N-i];
        for(j=1;j<=N;j++)
        {
            if(!vis[j])
            {
                if(!cur) break;
                else cur--;
            }
        }
        printf("%d ",j);
        vis[j]=true;
        x=x%fac[N-i];
    }
    puts("");
}

int N,K;

int main()
{
    getfac();
    scanf("%d%d",&N,&K);
    while(K--)
    {
        scanf("%s",opt);
        if(opt[0]=='P')
        {
            scanf("%lld",&a[0]);
            decantor(a[0]-1,N);
        }
        else if(opt[0]=='Q')
        {
            for(int i=0;i<N;i++)
            scanf("%lld",&a[i]);
            long long ans=0;
            ans=cantor(a,N)+1;
            printf("%lld",ans);
            puts("");
        }
    }
    return 0;
}

康托展开与逆展开

int Gauss()
{
    int ans=1;
    n--;
    for(int i=1;i<=n;i++)
    for(int j=i+1;j<=n;j++)
    {
        while(T[j][i])
        {
            int q=T[i][i]/T[j][i];
            for(int k=i;k<=n;k++)
            {
                T[i][k]=(1LL*T[i][k]+mod-1LL*T[j][k]*q%mod)%mod;
                swap(T[i][k],T[j][k]);
            }
            ans*=-1;
        }
        if(T[i][i]==0) return printf("0\n"),0;
    }  
    for(int i=1;i<=n;i++)
    ans=1LL*ans*T[i][i]%mod;
    
    printf("%d\n",(ans+mod)%mod);
    return 0;
}

无向图矩阵树定理

多项式算法

#include<bits/stdc++.h>
using namespace std;

const int N=20000010;
const double PI=acos(-1.0);

struct Complex
{
    double real;
    double imag;
};

Complex operator + (Complex a,Complex b)
{
    double X=a.real+b.real;
    double Y=a.imag+b.imag;
    return (Complex) {X,Y};
}

Complex operator - (Complex a,Complex b)
{
    double X=a.real-b.real;
    double Y=a.imag-b.imag;
    return (Complex) {X,Y};
}

Complex operator * (Complex a,Complex b)
{
    double X=a.real*b.real-a.imag*b.imag;
    double Y=a.real*b.imag+a.imag*b.real;
    return (Complex) {X,Y};
}

Complex F[N];
Complex G[N];
int n,m,rev[N];
int len=1,tt;

void FFT(Complex *T,int deg,int opt)
{
    for(int i=0;i<deg;i++)
    if(i<rev[i]) swap(T[i],T[rev[i]]);
    for(int step=1;step<deg;step<<=1)
    {    
        Complex Wn;
        double k=PI/step;
        Wn.real=cos(k);
        Wn.imag=sin(k)*opt;
        for(int j=0;j<deg;j+=(step<<1))    
        {    
            Complex w;
            w.real=1.0;
            w.imag=0.0;
            for(int k=0;k<step;k++,w=w*Wn)
            {
                Complex X=T[j+k];
                Complex Y=w*T[j+k+step];
                T[j+k]=X+Y;
                T[j+k+step]=X-Y;
            }
        }
    }
}

void getlen()
{
    int deg=n+m;
    while(len<=deg) len<<=1,tt++;
}

void getrev()
{
    for(int i=0;i<len;i++)
    rev[i]=((rev[i>>1]>>1)|((i&1)<<(tt-1)));
}

void calc()
{    
    FFT(F,len,1);
    FFT(G,len,1);
    for(int i=0;i<=len;i++) F[i]=F[i]*G[i];
    FFT(F,len,-1);
}

int main()
{
    scanf("%d%d",&n,&m);
    getlen();
    getrev();
    for(int i=0;i<=n;i++) scanf("%lf",&F[i].real);
    for(int i=0;i<=m;i++) scanf("%lf",&G[i].real);
    calc();
    for(int i=0;i<=m+n;i++)
    printf("%d ",(int)(F[i].real/len+0.5));
    return 0;
}

快速傅里叶变换

void NTT(int *a,int n,int opt)
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int step=0;(1<<step)<n;++step)
    {
        int wn=(~opt?g[step]:ng[step]);
        int tt=1<<step;
        for(int j=0;j<n;j+=tt<<1)
        for(int w=1,l=0;l<tt;l++,w=1LL*wn*w%mod)
        {
            int &x=a[j+l];
            int &y=a[j+l+tt];
            int t=1LL*w*y%mod;
            y=(x-t<0)?x-t+mod:x-t;
            x=(x+t>=mod)?x+t-mod:x+t;
        }
    }
    if(!(~opt)) for(int i=0;i<n;i++) a[i]=1LL*a[i]*inv[n]%mod;
}

快速数论变换

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
const int N=400010;
const int M=30000;
const double PI=acos(-1.0);
int read()
{
    int X=0,ww=0; char ch=0;
    while(!isdigit(ch)) {ww|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return ww?-X:X;
}
void write(int xx)
{
    if(xx<0) putchar('-'),xx=-xx;
    if(xx>9) write(xx/10);
    putchar(xx%10+'0');
}

int n,m,mod,l,kk;
int rev[N],ans[N];

complex<double> a1[N],a2[N],a[N],b1[N],b2[N],w[N];

void FFT(complex<double> *c,int jud)
{
    for(int i=0;i<l;i++) if(i>rev[i]) swap(c[i],c[rev[i]]);
    for(int step=1;step<l;step<<=1)
    {
        for(int j=0;j<l;j+=(step<<1))
        {
            for(int k=j;k<j+step;k++)
            {
                complex<double>W=w[l/step*(k-j)]; W.imag()*=jud;
                complex<double>x=c[k];
                complex<double>y=W*c[k+step];
                c[k]=x+y;
                c[k+step]=x-y;
            }
        }
    }
}

void work(complex<double>*f,complex<double>*g,int s)
{
    for(int i=0;i<l;i++) a[i]=f[i]*g[i];
    FFT(a,-1);
    for(int i=0;i<=m+n;i++) (ans[i]+=((long long)(a[i].real()/l+0.5)%mod*s+mod)%mod)%=mod;
}

int tmp;
int main()
{
    
    n=read(),m=read(),mod=read();
    for(int i=0;i<=n;i++)
    {
        tmp=read();
        a1[i]=tmp/M;
        b1[i]=tmp%M;
    }
    for(int i=0;i<=m;i++)
    {
        tmp=read();
        a2[i]=tmp/M;
        b2[i]=tmp%M;
    }
    for(l=1;l<=m+n;l<<=1,kk++);
    for(int i=0;i<l;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(kk-1));
    for(int i=0;i<l;i++) w[i].real()=cos(PI*i/l),w[i].imag()=sin(PI*i/l);
    FFT(a1,1);FFT(b1,1);
    FFT(a2,1);FFT(b2,1);
    work(a1,a2,M*M%mod);
    work(a1,b2,M%mod);
    work(a2,b1,M%mod);
    work(b1,b2,1);
    for(int i=0;i<=m+n;i++) printf("%d ",ans[i]);
    return puts(""),0;
}

任意模NTT

void cdq(int l,int r)
{
    if(l==r) return;
    int mid=(l+r)/2;
    cdq(l,mid);
    
    int deg=r-l+1;
    for(int i=0;i<=deg*2;i++) A[i]=B[i]=0;
    for(int i=0;i<=mid-l;i++) A[i]=F[i+l];
    for(int i=1;i<deg;i++) B[i]=G[i];
    deg=NTT(A,deg,1);
    deg=NTT(B,deg,1);
    for(int i=0;i<deg;i++) A[i]=1LL*A[i]*B[i]%mod;
    deg=NTT(A,deg,-1);
    for(int i=mid+1;i<=r;i++)
        F[i]=(F[i]+A[i-l]>=mod)?F[i]+A[i-l]-mod:F[i]+A[i-l];
        
    cdq(mid+1,r);
}

分治FFT

#include<bits/stdc++.h>
#define N 400010

using namespace std;

const double PI=acos(-1.0);

const int INF=0x3f3f3f3f;
const int mod=1000000007;
const int M=sqrt(mod);

struct Complex
{
    double real;
    double imag;
};

Complex operator + (Complex a,Complex b)
{
    return (Complex){a.real+b.real,a.imag+b.imag};
}
Complex operator - (Complex a,Complex b)
{
    return (Complex){a.real-b.real,a.imag-b.imag};
}
Complex operator * (Complex a,Complex b)
{
    return (Complex){a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real};
}

Complex A1[N],A2[N];
Complex B1[N],B2[N];
Complex W[N];

int a[N],b[N],c[N],d[N],rev[N],n;

int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) {w|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return w?-X:X;
}
void write(int x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}

int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

#define DFT 1
#define IDFT -1

void FFT(Complex *f,int deg,int opt)
{
    for(int i=0;i<deg;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int step=1;step<deg;step<<=1)
      for(int j=0;j<deg;j+=(step<<1))
        for(int k=0;k<step;k++)
        {
            Complex w=W[deg/step*k]; w.imag*=opt;
            Complex X=f[j+k];
            Complex Y=f[step+j+k]*w;
            f[j+k]=(X+Y);
            f[j+k+step]=(X-Y);
        }
    if(opt==IDFT) for(int i=0;i<deg;i++) f[i].real/=1.0*deg;
}

void work(int *r,int i)
{
    r[i]=0;
    r[i]=(r[i]+1LL*(long long)(A1[i].real+0.5)%mod*M%mod*M%mod)%mod;
    r[i]=(r[i]+1LL*(long long)(B1[i].real+0.5)%mod*M%mod)%mod;
    r[i]=(r[i]+1LL*(long long)(A2[i].real+0.5)%mod)%mod;
    r[i]=(r[i]+mod)%mod;
}

void mul(int *f,int *g,int deg,int *ret)
{
    for(int i=0;i<(deg<<1);i++) A1[i]=A2[i]=B1[i]=B2[i]=Complex {0,0};
    for(int i=0;i<deg;i++)
    {
        f[i]%=mod;
        g[i]%=mod;
        A1[i].real=f[i]/M*1.0;
        A2[i].real=f[i]%M*1.0;
        B1[i].real=g[i]/M*1.0;
        B2[i].real=g[i]%M*1.0;
    }
    int len=1,tt=0; for(len=1;len<=deg;len<<=1) ++tt;
    for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(tt-1));
    for(int i=1;i<len;i<<=1)
      for(int j=0;j<i;j++)
        W[len/i*j].real=cos(PI*j/i),W[len/i*j].imag=sin(PI*j/i);
    FFT(A1,len,DFT);
    FFT(A2,len,DFT);
    FFT(B1,len,DFT);
    FFT(B2,len,DFT);
    for(int i=0;i<len;i++)
    {
        Complex temp=A1[i]*B1[i];
        B1[i]=B1[i]*A2[i];
        A2[i]=A2[i]*B2[i];
        B2[i]=B2[i]*A1[i];
        A1[i]=temp;
        B1[i]=B1[i]+B2[i];
    }
    FFT(A1,len,IDFT);
    FFT(A2,len,IDFT);
    FFT(B1,len,IDFT);
    for(int i=0;i<deg;i++) work(ret,i);
}

void inv(int deg,int *a,int *b)
{
    if(deg==1)
    {
        b[0]=FP(a[0],mod-2);
        return;
    }
    inv((deg+1)>>1,a,b);
    mul(a,b,deg,c);
    mul(c,b,deg,d);
    for(int i=0;i<deg;i++) b[i]=(1LL*(2LL*b[i]-d[i])%mod+mod)%mod;
}

int main()
{
//    freopen("in.txt","r",stdin);
    n=read();
    for(int i=0;i<n;i++) a[i]=read();
    int K=1; for(;K<n;K<<=1);
    inv(K,a,b);
    for(int i=0;i<n;i++) write(b[i]),putchar(' ');
    return puts(""),0;
}

多项式求逆

#include<bits/stdc++.h>
using namespace std;
const int N=400010;
const int INF=0x3f3f3f3f;
const int mod=998244353;
int G=3;
int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) {w|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return w?-X:X;
}
void write(int x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}

int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

int rec[N],a[N],b[N],q[N],rev[N],ff[N],gg[N];

int NTT(int *a,int tot,int opt)
{
    int L=0,len=1;
    for(;len<tot;len<<=1,L++);
    for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1)
    {
        int ng=FP(G,(mod-1)/(i<<1));
        if(!(~opt)) ng=FP(ng,mod-2);
        for(int j=0;j<len;j+=i<<1)
        {
            int w=1;
            for(int k=0;k<i;k++,w=1LL*w*ng%mod)
            {
                int x=a[j+k],y=1LL*w*a[j+k+i]%mod;
                a[j+k]=(x+y)%mod;
                a[j+k+i]=(x-y+mod)%mod;
            }
        }
    }
    if(!(~opt)){
        int inv=FP(len,mod-2);
        for(int i=0;i<len;i++) a[i]=1LL*a[i]*inv%mod;
    }
    return len;
}
void inv(int deg,int *a,int *b)
{
    if(deg==1)
    {
        b[0]=FP(a[0],mod-2);
        return;
    }
    inv((deg+1)>>1,a,b);
    for(int i=0;i<deg;i++) ff[i]=a[i],gg[i]=b[i];
    int tmp=NTT(ff,deg<<1,1);
    NTT(gg,deg<<1,1);
    for(int i=0;i<tmp;i++) gg[i]=1LL*((2LL-1LL*ff[i]*gg[i]%mod)+mod)%mod*gg[i]%mod;
    NTT(gg,deg<<1,-1);
    for(int i=0;i<deg;i++) b[i]=gg[i];
    for(int i=0;i<tmp;i++) ff[i]=gg[i]=0;
}
int main()
{
//    freopen("in.txt","r",stdin);
    int n=read(),m=read();
    for(int i=0;i<=n;i++) a[i]=read();
    for(int i=0;i<=n;i++) rec[i]=a[i];
    reverse(a,a+n+1);
    for(int i=n-m+1;i<=n;i++) a[i]=0;
    for(int i=0;i<=m;i++) b[i]=read();
    reverse(b,b+m+1);
    inv(n-m+1,b,q);
    int leng=NTT(a,2*n-2*m+1,1);
    NTT(q,leng,1);
    for(int i=0;i<leng;i++) q[i]=1LL*a[i]*q[i]%mod;
    NTT(q,leng,-1);
    reverse(q,q+n-m+1);
    for(int i=n-m+1;i<leng;i++) q[i]=0;
    for(int i=0;i<n-m+1;i++) write(q[i]),putchar(' ');
    puts("");
    for(int i=0;i<=n;i++) a[i]=rec[i];
    reverse(b,b+m+1);
    leng=NTT(q,n+1,1);
    NTT(b,leng,1);
    for(int i=0;i<leng;i++) q[i]=1LL*q[i]*b[i]%mod;
    NTT(q,leng,-1);
    for(int i=0;i<m;i++) write((1LL*(a[i]-q[i]+mod)%mod)),putchar(' ');
    return puts(""),0;
}

多项式除法

// luogu-judger-enable-o2
#include <bits/stdc++.h>
using namespace std;

const int N=4000010;
const int mod=998244353;
const int G=3;

int len,tt,rev[N],g[22];

int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

void getl(int ll)
{
    int tt=(int)floor(log(ll)/log(2)+0.3)-1;
    for(int i=1;i<ll;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<tt);
}

int pls(int a,int b)
{
    return (a+=b)>=mod?a-mod:a;
}

void NTT(int *a,int opt)
{
    for(int i=1;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int step=1,rr=1;(step<<1)<=len;step<<=1,++rr)
    {
        int wn=g[rr];
        for(int j=0;j<len;j+=(step<<1))
        for(int k=0,w=1;k<step;++k,w=1LL*w*wn%mod)
        {
            int t=1LL*w*a[k+j+step]%mod;
            a[j+k+step]=pls(a[j+k],mod-t);
            a[j+k]=pls(a[j+k],t);
        }
    }
    if(!(~opt))
    {
        reverse(a+1,a+len);
        for(int i=0;i<len;i++) a[i]=1LL*a[i]*FP(len,mod-2)%mod;
    }
}

int raq[N],rbq[N];

void mul(int *a,int *b) //V
{
    memset(raq,0,sizeof(raq));
    memset(rbq,0,sizeof(rbq));
    for(int i=0;i<(len>>1);i++) raq[i]=a[i],rbq[i]=b[i];
    NTT(raq,1);
    NTT(rbq,1);
    for(int i=0;i<len;i++) raq[i]=1LL*raq[i]*rbq[i]%mod;
    NTT(raq,-1);
    for(int i=0;i<len;i++) a[i]=raq[i];
}

int rai[N],rbi[N];

void inv(int *a,int *b,int deg)
{
    if(deg==1)
    {
        b[0]=FP(a[0],mod-2);
        return;
    }
    inv(a,b,deg>>1);
    len=deg<<1;
    getl(len);
    for(int i=0;i<(len>>1);i++) rai[i]=a[i],rbi[i]=b[i];
    NTT(rai,1);
    NTT(rbi,1);
    for(int i=0;i<len;i++) rbi[i]=1LL*rbi[i]*rbi[i]%mod*rai[i]%mod;
    NTT(rbi,-1);
    for(int i=0;i<(len>>1);i++) b[i]=((2LL*b[i]%mod-rbi[i])%mod+mod)%mod;
    for(int i=0;i<=len;i++) rai[i]=rbi[i]=0;
}

int A[N],B[N];

void drvt(int *a,int deg)
{
    for(int i=1;i<=deg;i++)
        a[i-1]=1LL*a[i]*i%mod;
    a[deg]=0;
}

void itgt(int *a,int deg)
{
    for(int i=deg;i;i--)
        a[i]=1LL*a[i-1]*FP(i,mod-2)%mod;
    a[0]=0;
}

int rln[N];

void getln(int *a,int deg)
{
    len=2;
    memset(rln,0,sizeof(rln));
    while(len<=deg) len<<=1;
    inv(a,rln,len);
    drvt(a,deg);
    while(len<=(deg<<2)) len<<=1;
    getl(len);
    mul(a,rln);
    itgt(a,deg);
    for(int i=deg+1;i<=len;i++) a[i]=0;
}

int T[N],K[N];

void getexp(int *a,int deg)
{
    int t=0,z;
    memset(K,0,sizeof(K));
    K[0]=1;
    len=z=2;
    while(z<=(deg<<1))
    {
        z<<=1;
        t^=1;
        for(int i=0;i<(z>>1);i++) T[i]=K[i];
        getln(T,(z>>1)-1);
        for(int i=0;i<(z>>1);i++) T[i]=pls(a[i]+(i==0),mod-T[i]);
        getl(len=z);
        mul(K,T);
        for(int i=z;i<=(z<<1);i++) K[i]=T[i]=0;
    }
    for(int i=0;i<=deg;i++) a[i]=K[i];
}

int F[N],n;

int main()
{
    for(int i=1,j=2;i<=20;++i,j<<=1)
        g[i]=FP(G,(mod-1)/j);
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d",&F[i]);
    getexp(F,n);
    for(int i=0;i<n;i++) printf("%d ",F[i]);
    return puts(""),0;
}

多项式指数

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;

const int N=400010;
const int mod=998244353;

int n,m,len,tt,rev[N];
int A[N],B[N],ra[N],rb[N],F[N],G[N];

void getl(int ll)
{
    for(len=1,tt=0;len<=ll;len<<=1) tt++;tt--;
    for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<tt);
}

int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

void NTT(int *a,int opt)
{
    for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int step=1;step<len;step<<=1)
    {
        int g=FP(3,(mod-1)/(step<<1));
        if(!(~opt)) g=FP(g,mod-2);
        for(int j=0,p=step<<1;j<len;j+=p)
            for(int k=0,w=1;k<step;k++,w=1LL*w*g%mod)
            {
                int X=a[j+k];
                int Y=1LL*a[j+k+step]*w%mod;
                a[j+k]=(X+Y)%mod;
                a[j+k+step]=((X-Y)%mod+mod)%mod;
            }
    }
    if(!(~opt)) for(int i=0;i<len;i++) a[i]=1LL*a[i]*FP(len,mod-2)%mod;
}

void inv(int *a,int *b,int deg)
{
    if(deg==1)
    {
        b[0]=FP(a[0],mod-2);
        return;
    }
    inv(a,b,deg>>1);
    getl(deg);
    for(int i=0;i<deg;i++) ra[i]=a[i],rb[i]=b[i];
    for(int i=deg;i<len;i++) ra[i]=rb[i]=0;
    NTT(ra,1);
    NTT(rb,1);
    for(int i=0;i<len;i++) ra[i]=1LL*ra[i]*rb[i]%mod*rb[i]%mod;
    NTT(ra,-1);
    for(int i=0;i<deg;i++) b[i]=((2LL*b[i]%mod-ra[i])%mod+mod)%mod;
}

void drvt(int *a,int *b,int ll)
{
    for(int i=1;i<ll;i++)
        b[i-1]=1LL*i*a[i]%mod;
    b[ll-1]=0;
}

void itgt(int *a,int *b,int ll)
{
    for(int i=1;i<ll;i++)
        b[i]=1LL*a[i-1]*FP(i,mod-2)%mod;
    b[0]=0;
}

void getln(int *f,int *g,int s)
{
    drvt(f,A,s);
    inv(f,B,s);
    for(len=1,tt=0;len<=s;len<<=1) tt++;
    for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(tt-1));
    NTT(A,1);
    NTT(B,1);
    for(int i=0;i<len;i++) A[i]=1LL*A[i]*B[i]%mod;
    NTT(A,-1);
    itgt(A,g,s);
}

int main()
{
//    freopen("i.txt","r",stdin);
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d",&F[i]);
    for(m=1;m<=n;m<<=1); getln(F,G,m);
    for(int i=0;i<n;i++) printf("%d ",G[i]);
    return puts(""),0;
}

多项式对数

#include<bits/stdc++.h>
using namespace std;

const int N=400010;
const int mod=998244353;

int n,m,rev[N],G=3;
int A[N],iv[N],sq[N],F[N];

int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

int NTT(int *a,int tot,int opt)
{
    int len=1,tt=0;
    for(;len<tot;len<<=1,++tt);tt--;
    for(int i=0;i<len;i++)
    {
        rev[i]=(rev[i>>1]>>1|((i&1)<<tt));
        if(i<rev[i]) swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<len;step<<=1)
    {
        int g=FP(G,(mod-1)/(step<<1));
        if(!(~opt)) g=FP(g,mod-2);
        for(int j=0;j<len;j+=step<<1)
        {
            int w=1;
            for(int k=0;k<step;k++,w=1LL*w*g%mod)
            {
                int X=a[j+k];
                int Y=1LL*w*a[j+k+step]%mod;
                a[j+k]=(X+Y)%mod;
                a[j+k+step]=((X-Y)%mod+mod)%mod;
            }
        }
    }
    if(!(~opt)) for(int i=0;i<len;i++) a[i]=1LL*a[i]*FP(len,mod-2)%mod;
    return len;
}

void inv(int *a,int *b,int deg)
{
    if(deg==1)
    {
        memset(b,0,sizeof(b[0])*N);
        b[0]=FP(a[0],mod-2);
        return;
    }
    inv(a,b,(deg+1)>>1);
    for(int i=0;i<deg;i++) A[i]=a[i];
    int len=NTT(A,deg<<1,1);
    NTT(b,deg<<1,1);
    for(int i=0;i<len;i++) b[i]=1LL*b[i]*(2LL-1LL*A[i]*b[i]%mod+mod)%mod;
    NTT(b,deg<<1,-1);
    for(int i=deg;i<len;i++) b[i]=0;
    for(int i=0;i<len;i++) A[i]=0;
}

void getsq(int *f,int *g,int deg)
{
    if(deg==1)
    {
        memset(g,0,sizeof(g));
        g[0]=1;
        return;
    }
    getsq(f,g,(deg+1)>>1);
    inv(g,iv,deg);
    for(int i=0;i<deg;i++) A[i]=f[i];
    int len=NTT(A,deg<<1,1);
    NTT(g,deg<<1,1);
    NTT(iv,deg<<1,1);
    for(int i=0;i<len;i++) g[i]=1LL*(1LL*g[i]*g[i]%mod+A[i])*iv[i]%mod*FP(2,mod-2)%mod;
    NTT(g,deg<<1,-1);
    for(int i=deg;i<len;i++) g[i]=0;
    for(int i=0;i<len;i++) A[i]=0;
}

int tmp;

int main()
{
    scanf("%d %d",&n,&m);
    for(int i=1;i<=n;i++)
    {
        scanf("%d",&tmp);
        if(tmp<=m) F[tmp]=1;
    }
    m++;
    for(int i=0;i<m;i++) F[i]=(mod-4LL*F[i]%mod)%mod;
    F[0]=(F[0]+1)%mod;
    getsq(F,sq,m);
    sq[0]=(sq[0]+1)%mod;
    
    inv(sq,iv,m);
    for(int i=1;i<m;i++)
    {
        tmp=2LL*iv[i]%mod;
        printf("%d\n",tmp);
    }
    return 0;
}

多项式开根

进制/反演

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define N 1000010
using namespace std;

const int mod=998244353;

int read()
{
    int X=0,ww=0; char ch=0;
    while(!isdigit(ch)) {ww|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return ww?-X:X;
}

int FP(int a,int k)
{
    int ans=1;
    for(;k;a=1LL*a*a%mod,k>>=1)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

int add(int a,int b)
{
    if((a+=b)>=mod) a-=mod;
    return a;
}
int sub(int a,int b)
{
    if((a-=b)<0) a+=mod;
    return a;
}

void FMT_or(int *f,int deg,int jud)
{
    for(int step=0;step<deg;step++)
      for(int j=0;j<(1<<deg);j++)
        if(j>>step&1)
            f[j]=add(f[j],jud^1?mod-f[j^1<<step]:f[j^1<<step]);
    return;
}

void FMT_and(int *f,int deg,int jud)
{
    for(int step=0;step<deg;step++)
      for(int j=(1<<deg)-1;~j;j--)
        if(~j>>step&1)
            f[j]=add(f[j],jud^1?mod-f[j|1<<step]:f[j|1<<step]);
    return;
}

void FWT(int *f,int deg,int jud)
{
    for(int step=1;(step<<1)<=(1<<deg);step<<=1)
      for(int j=0;j<(1<<deg);j+=(step<<1))
        for(int k=0;k<step;k++)
        {
            int X=f[j+k];
            int Y=f[j+k+step];
            f[j+k]=1LL*(jud^1?FP(2,mod-2):1)*add(X,Y)%mod;
            f[j+k+step]=1LL*(jud^1?FP(2,mod-2):1)*sub(X,Y)%mod;
        }
    return;
}

int n,len;
int F1[N],G1[N],F2[N],G2[N],F3[N],G3[N],H1[N],H2[N],H3[N];

void Or()
{
    FMT_or(F1,n,1);
    FMT_or(G1,n,1);
    for(int i=0;i<len;i++) H1[i]=1LL*F1[i]*G1[i]%mod;
    FMT_or(H1,n,-1);
    for(int i=0;i<len;i++) printf("%d ",H1[i]);
    puts("");
}

void And()
{
    FMT_and(F2,n,1);
    FMT_and(G2,n,1);
    for(int i=0;i<len;i++) H2[i]=1LL*F2[i]*G2[i]%mod;
    FMT_and(H2,n,-1);
    for(int i=0;i<len;i++) printf("%d ",H2[i]);
    puts("");
}

void Xor()
{
    FWT(F3,n,1);
    FWT(G3,n,1);
    for(int i=0;i<len;i++) H3[i]=1LL*F3[i]*G3[i]%mod;
    FWT(H3,n,-1);
    for(int i=0;i<len;i++) printf("%d ",H3[i]);
    puts("");
}

int main()
{
    scanf("%d",&n);
    len=1<<n;
    for(int i=0;i<len;i++) F1[i]=F2[i]=F3[i]=read();
    for(int i=0;i<len;i++) G1[i]=G2[i]=G3[i]=read();
    Or();
    And();
    Xor();
    return 0;
}

快速沃尔什变换

#include<bits/stdc++.h>
#define N 70
using namespace std;
long long a[N],ans,tmp;
int n;

int main()
{
    scanf("%d",&n);
    while(n--)
    {
        scanf("%lld",&tmp);
        for(int x=60;~x;x--)
        {
            if (tmp&(1LL<<x))
            {
                if(!a[x])
                {
                    a[x]=tmp;
                    break;
                }
                else tmp^=a[x];
            }
        }
    }
    for(int x=60;~x;x--)
    {
        if((ans^a[x])>ans) ans^=a[x];
    }
    return printf("%lld\n",ans),0;
}

线性基

int Gray(int x)
{
    return x^(x>>1);
}

格雷码

int Trie[32*MAXN][2];
long long val[MAXN*32];
void insert(long long d)
{
    int root=0;
    for(int i=32;~i;i--)
    {
        int id=(d>>i)&1;
        if(!Trie[root][id])
        Trie[root][id]=++tot;
        root=Trie[root][id];
    }
    val[root]=d;
} 
long long query(long long d)
{
    int root=0;
    for(int i=32;~i;i--)
    {
        int id=(d>>i)&1;
        if(Trie[root][id^1])
        root=Trie[root][id^1];
        else root=Trie[root][id];
    }
    return val[root];
}

01Trie

namespace PT
{
    struct Trie_node
    {
        int root;
        int nxt[2];
        int cnt;
    }Trie[MAXN*24];
    
    void insert(int x)
    {
        int tmp=Trie[cnt].root;
        Trie[++cnt].root=tot+1;
        for(int i=1<<24;i;i>>=1)
        {
            Trie[++tot].cnt=Trie[tmp].cnt+1;
            bool c=x&i;
            Trie[tot].nxt[c]=tot+1;
            Trie[tot].nxt[!c]=Trie[tmp].nxt[!c];
            tmp=Trie[tmp].nxt[c];
        }
        Trie[++tot].cnt=Trie[tmp].cnt+1;
    }
    
    int query(int l,int r,int x)
    {
        if(l>r) return 0;
        l=Trie[l-1].root;
        r=Trie[r].root;
        int ans=0;
        for(int i=1<<24;i;i>>=1)
        {
            bool c=x&i;
            int t1=Trie[Trie[r].nxt[!c]].cnt;
            int t2=Trie[Trie[l].nxt[!c]].cnt;
            if(t1-t2)
            {
                ans+=i;
                l=Trie[l].nxt[!c];
                r=Trie[r].nxt[!c];
            }
            else
            {
                l=Trie[l].nxt[c];
                r=Trie[r].nxt[c];
            }
        }
        return ans;
    }
}
using namespace PT;

可持久化01Trie

矩阵

inline Matrix mul(Matrix X,Matrix Y)
{
    Matrix Z;
    memset(Z.T,0,sizeof(Z.T));
    for(register int i=1;i<=n;i++)
    for(register int j=1;j<=n;j++)
    for(register int k=1;k<=n;k++)
    {
        Z.T[i][j]=(1LL*Z.T[i][j]%mod+1LL*X.T[i][k]*Y.T[k][j]%mod)%mod;
    }
    return Z;
}

矩阵乘法

inline Matrix mul(Matrix X,Matrix Y)
{
    Matrix Z;
    memset(Z.T,0,sizeof(Z.T));
    for(register int i=1;i<=n;i++)
    for(register int j=1;j<=n;j++)
    for(register int k=1;k<=n;k++)
    {
        Z.T[i][j]=(1LL*Z.T[i][j]%mod+1LL*X.T[i][k]*Y.T[k][j]%mod)%mod;
    }
    return Z;
}

Matrix pow(Matrix M,long long k)
{
    Matrix ans=E;
    for(;k;k>>=1,M=mul(M,M))
        if(k&1) ans=mul(ans,M);
    return ans;
}

矩阵快速幂

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;

const int N=1010;
const int mod=1000000007;

int n;
long long T[N][N];

inline int FP(int a,int k)
{
    int ans=1;
    for(;k;k>>=1,a=1LL*a*a%mod)
        if(k&1) ans=1LL*ans*a%mod;
    return ans;
}

int main()
{
    scanf("%d",&n);
    for(register int i=1;i<=n;i++)
    {
        for(register int j=1;j<=n;j++)
        {
            scanf("%lld",&T[i][j]);
        }
        T[i][i+n]=1;
    }
    for(register int i=1;i<=n;i++)
    {        
        for(register int j=i;j<=n;j++)
        {
            if(T[j][i])
            {
                for(register int k=1;k<=n+n;k++)
                swap(T[i][k],T[j][k]);
                break;
            }
        }
        if(!T[i][i])
        {
            puts("No Solution");
            return 0;
        }
        int inv=FP(T[i][i],mod-2);
        for(register int j=i;j<=n+n;j++)
            T[i][j]=1LL*T[i][j]*inv%mod;
        for(register int j=1;j<=n;j++)
        {
            if(j!=i)
            {
                inv=T[j][i];
                for(register int k=i;k<=n+n;k++)
                T[j][k]=(1LL*T[j][k]-1LL*inv*T[i][k]%mod+mod)%mod;
            }
        }
    }
    for(register int i=1;i<=n;i++,puts(""))
    for(register int j=n+1;j<=n+n;j++) printf("%lld ",T[i][j]);
    return 0;
}

矩阵求逆

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;

const int mod=1000000007;
long long n,T;

struct Matrix
{
    long long T[4][4];    
}E,P;

inline Matrix mul(Matrix X,Matrix Y)
{
    Matrix Z;
    memset(Z.T,0,sizeof(Z.T));
    for(register int i=1;i<=3;i++)
    for(register int j=1;j<=3;j++)
    for(register int k=1;k<=3;k++)
    {
        Z.T[i][j]=(1LL*Z.T[i][j]%mod+1LL*X.T[i][k]*Y.T[k][j]%mod)%mod;
    }
    return Z;
}

Matrix pow(Matrix M,long long k)
{
    Matrix ans=E;
    for(;k;k>>=1,M=mul(M,M))
        if(k&1) ans=mul(ans,M);
    return ans;
}

int main()
{
    E.T[1][1]=1;
    E.T[3][1]=1;
    E.T[1][2]=1;
    E.T[2][3]=1;
    scanf("%lld",&T);
    while(T--)
    {
        scanf("%lld",&n);
        if(n==1)
        {
            printf("1\n");
            continue;
        }
        memset(P.T,0,sizeof(P.T));
        P=pow(E,n-2);
        printf("%lld\n",P.T[1][1]);
    }
    return 0;
}

矩阵加速

#include<bits/stdc++.h>
using namespace std;
const int N=110;
const double eps=1e-6;
#define p jud
#define a mx

double mx[N][N];
bool jud;
int n;

void gauss()
{
     for(int i=1;i<=n;i++)
    {
        int r=i;
  for(int j=i+1;j<=n;j++)
        if(fabs(mx[r][i])<fabs(mx[j][i])) r=j;
  if(r!=i) for(int j=1;j<=n+1;j++) swap(mx[i][j],mx[r][j]);
  if(fabs(mx[i][i])<eps) {jud=false; return;}
        for(int j=1;j<=n;j++)
        if(j!=i)
        {
            double q=mx[j][i]/mx[i][i];
            for(int k=i+1;k<=n+1;k++)
                mx[j][k]-=mx[i][k]*q;
        }
    }
    for(int i=1;i<=n;i++)
    mx[i][n+1]/=mx[i][i];
}

int main()
{
    jud=true;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n+1;j++)
        scanf("%lf",&mx[i][j]);
    gauss();
    if(!jud) return printf("No Solution\n"),0;
    for(int i=1;i<=n;i++) printf("%.2lf\n",mx[i][n+1]);
    return 0;
}

高斯消元法

posted @ 2018-09-17 14:45 _ConveX 阅读(271) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

_ConveX

吉林一中第一蒟蒻

【数学相关】模板

数论

组合数学

多项式算法

进制/反演

矩阵

公告