【算法总结】多项式相关

【快速傅里叶变换】

〖相关资料〗

《虚数的图解》        《虚数的意义》        《FFT学习笔记》

《从多项式乘法到快速傅里叶变换》        《Fast Fourier Transform》

《【快速傅里叶变换】【FFT】【WikiOI】【P3132】【高精度练习之超大整数乘法】》

〖模板代码〗

  [FFT]

const double pi=acos(-1);
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};};
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
{
    int n=1;while(n<n1+n2)n<<=1;
    fft(c1,n,1);fft(c2,n,1);
    for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
    fft(c1,n,-1);
    for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
}

  [MTT]

#include<cstdio>
#include<algorithm> 
#include<cstring>
#include<cstdlib>
#include<cmath>
#define double long double
#define LL long long
using namespace std;
const int N=524288+5;
const double pi=acos(-1);
int n,n1,n2,mod,qmod=32768;
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};};
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
cpx a1[N],a2[N],b1[N],b2[N];
cpx c1[N],c2[N],c3[N];
cpx w,omg[N],omgi[N];
struct mtt{LL s[N];}a,b,ans,zero;
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        for(cpx *p=a;p!=a+n;p+=l)
            for(int i=0;i<m;i++)
            {
                if(f==1)w=omg[n/l*i];
                else w=omgi[n/l*i];
                cpx t=p[m+i]*w;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
    }
    if(f==-1)for(int i=0;i<n;i++)a[i].r/=n;
}
mtt operator * (mtt x,mtt y)
{
    for(int i=0;i<n;i++)
    {
        a1[i].r=x.s[i]/qmod;a1[i].i=0;
        a2[i].r=x.s[i]%qmod;a2[i].i=0;
        b1[i].r=y.s[i]/qmod;b1[i].i=0;
        b2[i].r=y.s[i]%qmod;b2[i].i=0;
    }
    fft(a1,n,1);fft(a2,n,1);
    fft(b1,n,1);fft(b2,n,1);
    for(int i=0;i<n;i++)
    {
        c1[i]=a1[i]*b1[i];
        c2[i]=a2[i]*b1[i]+a1[i]*b2[i];
        c3[i]=a2[i]*b2[i];
    }
    fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1);
    mtt ret=zero;
    for(int i=0;i<n;i++)
        ret.s[i]=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod
                +(LL)(c2[i].r+0.5)%mod*qmod%mod
                +(LL)(c3[i].r+0.5)%mod)%mod;
    return ret;
}
int main()
{
    n1=read()+1;n2=read()+1;mod=read();
    n=1;while(n<n1+n2)n<<=1;
    for(int i=0;i<n;i++)
    {
        omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i));
        omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i));
    }
    for(int i=0;i<n1;i++)a.s[i]=read();
    for(int i=0;i<n2;i++)b.s[i]=read();
    ans=a*b;
    for(int i=0;i<n1+n2-1;i++)printf("%lld ",ans.s[i]);
    return 0;
}

〖相关题目〗

1.【bzoj2179】FFT快速傅立叶

题意：给出两个n位10进制整数x和y，输出x*y。

分析：FFT裸题

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;
const double pi=acos(-1);
const int N=131072+5;
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
int n,ans[N];
char ch[N];
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
{
    int n=1;while(n<n1+n2)n<<=1;
    fft(c1,n,1);fft(c2,n,1);
    for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
    fft(c1,n,-1);
    for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
}
int main()
{
    scanf("%d",&n);
    scanf("%s",ch);
    for(int i=0;i<n;i++)a[i].r=ch[n-i-1]-'0';
    scanf("%s",ch);
    for(int i=0;i<n;i++)b[i].r=ch[n-i-1]-'0';
    mul(a,n,b,n,ans);int m=n*2-2;
    for(int i=0;i<=m;i++)
        if(ans[i]>=10)
        {
            ans[i+1]+=ans[i]/10;ans[i]%=10;
            if(i==m)m++;
        }
    for(int i=m;i>=0;i--)printf("%d",ans[i]);
    return 0;
}

2.【bzoj2194】快速傅立叶之二

题意：计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n。

分析：Heret1cの博客

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;
const double pi=acos(-1);
const int N=262144+5;
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
int n,ans[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
{
    int n=1;while(n<n1+n2)n<<=1;
    fft(c1,n,1);fft(c2,n,1);
    for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
    fft(c1,n,-1);
    for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
}
int main()
{
    n=read();
    for(int i=0;i<n;i++)a[i]=read(),b[n-i-1]=read();
    mul(a,n,b,n,ans);
    for(int i=n-1;i<n*2-1;i++)printf("%d\n",ans[i]);
    return 0;
}

3.【bzoj3527】[Zjoi2014]力

题意：见原题

分析：Mychaelの博客

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;
const double pi=acos(-1);
const int N=262144+5;
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}f1[N],f2[N],g[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
int n;double x,ans1[N],ans2[N];
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    {
        scanf("%lf",&x);
        f1[i].r=f2[n-i-1].r=x;
    }
    for(int i=1;i<n;i++)g[i].r=1.0/i/i;
    int mx=1;while(mx<n*2)mx<<=1;
    fft(f1,mx,1);fft(f2,mx,1);fft(g,mx,1);
    for(int i=0;i<mx;i++)f1[i]=f1[i]*g[i];
    for(int i=0;i<mx;i++)f2[i]=f2[i]*g[i];
    fft(f1,mx,-1);fft(f2,mx,-1);
    for(int i=0;i<n*2-1;i++)ans1[i]=f1[i].r/mx;
    for(int i=0;i<n*2-1;i++)ans2[i]=f2[i].r/mx;
    for(int i=0;i<n;i++)
        printf("%.3lf\n",ans1[i]-ans2[n-i-1]);
    return 0;
}

4.【bzoj3160】万径人踪灭

题意：给定一个由'a'和'b'构成的字符串，求不连续回文子序列的个数。

分析：PoPoQQQの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long
using namespace std;
const int N=262144;
const int M=1e6+5;
const int mod=1e9+7;
const double pi=acos(-1);
int len,ans;
int bin[M],tmp[N],p[N];
char ch[M],s[N];
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
int manacher()
{
    int n=0;s[++n]='#';
    for(int i=1;i<=len;i++)s[++n]=ch[i],s[++n]='#';
    int mxr=0,id=0,sum=0;
    for(int i=1;i<=n;i++)
    {
        if(i<mxr)p[i]=min(p[2*id-i],mxr-i);else p[i]=1;
        for(;i+p[i]<=n&&s[i-p[i]]==s[i+p[i]];p[i]++);
        if(i+p[i]>mxr)mxr=i+p[i],id=i;
        sum=(sum+p[i]/2)%mod;
    }
    return sum;
}
int main()
{
    scanf("%s",ch+1);len=strlen(ch+1);
    bin[0]=1;for(int i=1;i<=M-5;i++)bin[i]=bin[i-1]*2%mod;
    int n=1;while(n<len*2)n<<=1;
    for(int i=1;i<=len;i++)a[i].r=(ch[i]=='a');
    fft(a,n,1);for(int i=0;i<n;i++)b[i]=a[i]*a[i];
    for(int i=0;i<n;i++)a[i].r=a[i].i=0;
    for(int i=1;i<=len;i++)a[i].r=(ch[i]=='b');
    fft(a,n,1);for(int i=0;i<n;i++)b[i]=b[i]+a[i]*a[i];
    fft(b,n,-1);for(int i=0;i<n;i++)tmp[i]=(LL)(b[i].r/n+0.5);
    for(int i=0;i<n;i++)ans=(ans+bin[(tmp[i]+1)>>1]-1)%mod;
    printf("%d",(ans-manacher()+mod)%mod);
    return 0;
}

5.【bzoj3992】[SDOI2015]序列统计

题意：见原题

分析：L_0_Forever_LFの博客

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#define LL long long
using namespace std;
const double pi=acos(-1);
const int N=2e4+5;
const int mod=1004535809;
int n,ci,m,x,num,tmp,qmod;
int lg[N],b[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}
w,a1[N],b1[N],a2[N],b2[N],c1[N],c2[N],c3[N],omg[N],omgi[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        for(cpx *p=a;p!=a+n;p+=l)
        {
            for(int i=0;i<m;i++)
            {
                if(f==1)w=omg[n/l*i];
                else w=omgi[n/l*i];
                cpx t=p[m+i]*w;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
    if(f==-1)for(int i=0;i<n;i++)a[i].r/=n;
}
struct mtt{LL s[N];}zero,st,ans;
int power(int a,int b,int mod)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int find(int n)
{
    int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0;
    for(int i=2;i<=mx;i++)
        if(num%i==0)
        {
            b[++tot]=i;
            while(num%i==0)num/=i;
        }
    if(num!=1)b[++tot]=num;
    for(int i=2;i<=n;i++)
    {
        bool ok=true;
        for(int j=1;j<=tot;j++)
            if(power(i,(n-1)/b[j],n)==1){ok=false;break;}
        if(ok)return i;
    }
    return 0;
}
void pre()
{
    int g=find(m),sum=1;
    for(int i=0;i<m-1;i++)
    {
        lg[sum]=i;
        sum=1ll*sum*g%m;
    }
}
int tim=0; 
mtt operator * (mtt x,mtt y)
{
    for(int i=0;i<n;i++)
    {
        a1[i].r=x.s[i]/qmod;a1[i].i=0;
        a2[i].r=x.s[i]%qmod;a2[i].i=0;
        b1[i].r=y.s[i]/qmod;b1[i].i=0;
        b2[i].r=y.s[i]%qmod;b2[i].i=0;
    }
    fft(a1,n,1);fft(a2,n,1);
    fft(b1,n,1);fft(b2,n,1);
//  printf("-----------%d-------\n",n);
    tim++;
    for(int i=0;i<n;i++)
    {
        c1[i]=a1[i]*b1[i];
        c2[i]=a2[i]*b1[i]+a1[i]*b2[i];
        c3[i]=a2[i]*b2[i];
//      if(tim<=10)printf("%.3lf %.3lf %.3lf\n",c1[i].r,c2[i].r,c3[i].r);
    }
    fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1);
    mtt ret=zero;
    for(int i=0;i<n;i++)
    {
        LL temp=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod
                +(LL)(c2[i].r+0.5)%mod*qmod%mod
                +(LL)(c3[i].r+0.5)%mod)%mod;
        ret.s[i%(m-1)]=(ret.s[i%(m-1)]+temp)%mod;
    }
    return ret;
}
void work(mtt a,int b)
{
    bool flag=false;
    while(b)
    {
        if(b&1)
        {
            if(!flag)ans=a,flag=true;
            else ans=ans*a;
        }
        a=a*a;b>>=1;
//      puts("");
//      for(int i=0;i<n;i++)printf("%lld ",a.s[i]);
//      puts("");
    }
}
int main()
{
    ci=read();m=read();x=read();num=read();
    pre();n=1;while(n<m+m)n<<=1;
    for(int i=0;i<n;i++)
    {
        omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i));
        omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i));
    }
    for(int i=0;i<n;i++)zero.s[i]=0;
    st=zero;
    for(int i=1;i<=num;i++)
    {
        tmp=read();
        if(tmp)st.s[lg[tmp]]=1;
    }
//  for(int i=0;i<n;i++)printf("%d %d\n",i,lg[i]);
//  for(int i=0;i<n;i++)printf("%lld ",st.s[i]);
//  puts("");
    qmod=sqrt(1.0*mod);
    work(st,ci);
    printf("%lld",(ans.s[lg[x]]%mod+mod)%mod);
    return 0;
}

6.【bzoj4827】[Hnoi2017]礼物

题意：见原题

分析：cdcqの博客

#include<cstdio>
#include<algorithm> 
#include<cstring>
#include<cmath>
#define LL long long
using namespace std;
const int N=262144+5;
const double pi=acos(-1);
int n,m,anss,sum,x[N],y[N],ans[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
}
void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
{
    int n=1;while(n<n1+n2)n<<=1;
    fft(c1,n,1);fft(c2,n,1);
    for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
    fft(c1,n,-1);
    for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
}
int main()
{
    n=read();m=read();
    for(int i=0;i<n;i++)x[i]=read();
    for(int i=0;i<n;i++)a[i].r=x[n-i-1];
    for(int i=0;i<n;i++)y[i]=read(),b[i].r=b[i+n].r=y[i],sum+=x[i]-y[i];
    mul(a,n,b,n*2-1,ans);
    for(int i=n-1;i<=2*n-2;i++)anss=max(anss,ans[i]);
    anss=-anss*2;int mx=0x3f3f3f3f;
    for(int c=0;c<=m;c++)mx=min(mx,sum*c*2+n*c*c);
    anss+=mx;
    for(int i=0;i<n;i++)anss+=x[i]*x[i]+y[i]*y[i];
    printf("%d",anss);
    return 0;
}

7.【bzoj4259】残缺的字符串

题意：见原题

分析：Clarisの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=524288+5;
const double pi=acos(-1);
int n,n1,n2,cnt,out[N],a[N],b[N];
char s1[N],s2[N];
struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};};
cpx A[N],B[N],ans[N];
cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cpx *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
        for(cpx *p=a;p!=a+n;p+=l)
        {
            cpx w=cpx(1,0);
            for(int i=0;i<m;i++)
            {
                cpx t=p[m+i]*w;w=w*nw;
                p[m+i]=p[i]-t;p[i]=p[i]+t;
            }
        }
    }
    if(f==-1)for(int i=0;i<n;i++)a[i].r/=n;
}
int main()
{
    scanf("%d%d%s%s",&n1,&n2,s1,s2);
    n=1;while(n<n1||n<n2)n<<=1;
    for(int i=0;i<n1;i++)if(s1[i]!='*')a[n1-i-1]=s1[i]-'a'+1;
    for(int i=0;i<n2;i++)if(s2[i]!='*')b[i]=s2[i]-'a'+1;
    for(int i=0;i<n;i++)A[i]=cpx(a[i]*a[i]*a[i],0);
    for(int i=0;i<n;i++)B[i]=cpx(b[i],0);
    fft(A,n,1);fft(B,n,1);
    for(int i=0;i<n;i++)ans[i]=A[i]*B[i];
    for(int i=0;i<n;i++)A[i]=cpx(a[i]*a[i],0);
    for(int i=0;i<n;i++)B[i]=cpx(b[i]*b[i],0);
    fft(A,n,1);fft(B,n,1);
    for(int i=0;i<n;i++)ans[i]=ans[i]-A[i]*B[i]*cpx(2,0);
    for(int i=0;i<n;i++)A[i]=cpx(a[i],0);
    for(int i=0;i<n;i++)B[i]=cpx(b[i]*b[i]*b[i],0);
    fft(A,n,1);fft(B,n,1);
    for(int i=0;i<n;i++)ans[i]=ans[i]+A[i]*B[i];
    fft(ans,n,-1);
    for(int i=n1-1;i<n2;i++)
        if(ans[i].r<0.5)out[++cnt]=i-n1+2;
    printf("%d\n",cnt);
    for(int i=1;i<=cnt;i++)printf("%d ",out[i]);
    return 0;
}

【生成函数】

〖相关题目〗

1.【bzoj3028】食物

题意：给出n种重量为1的食物以及它们的限制，求出食物个数总和恰好为n的方案数。

分析：maijingの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int mod=1e4+7;
int n,ans;
char ch[505];
int main()
{
    scanf("%s",ch+1);
    n=strlen(ch+1);
    for(int i=1;i<=n;i++)
        ans=(ans*10+ch[i]-'0')%mod;
    printf("%d\n",ans*(ans+1)%mod*(ans+2)%mod*1668%mod);
    return 0;
}

【快速数论变换】

〖模板代码〗

  [NTT]

const int N=131072+5;
const int mod=998244353;
int qmod(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
     int k=0;while((1<<k)<n)k++;
     for(int i=0;i<n;i++)
     {
         int t=0;
         for(int j=0;j<k;j++)
             if(i&(1<<j))t|=(1<<(k-j-1));
         if(i<t)swap(a[i],a[t]);
     }
     for(int l=2;l<=n;l<<=1)
     {
         int m=l>>1,nw=qmod(3,(mod-1)/l);
         if(f==-1)nw=qmod(nw,mod-2);
         for(int *p=a;p!=a+n;p+=l)
         {
             int w=1;
             for(int i=0;i<m;i++)
             {
                 int t=1ll*p[m+i]*w%mod;
                 p[m+i]=(p[i]-t+mod)%mod;
                 p[i]=(p[i]+t)%mod;
                 w=1ll*w*nw%mod;
             }
         }
     }
     if(f==-1)
     {
         int inv=qmod(n,mod-2);
         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
     }
}
int main()
{
    n=1;while(n<len1+len2)n<<=1;
    ntt(a,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod;
    ntt(c,n,-1);
    return 0;
}

  [多项式求逆]

int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l);
        if(f==-1)nw=power(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void pinv(int *a,int *b,int n)
{
    if(n==1){b[0]=power(a[0],mod-2);return;}
    pinv(a,b,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
    ntt(c,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}

  [多项式开根]

int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l);
        if(f==-1)nw=power(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void pinv(int *a,int *b,int n)
{
    if(n==1){b[0]=power(a[0],mod-2);return;}
    pinv(a,b,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
    ntt(c,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)
        b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}
void psqrt(int *a,int *b,int *invb,int n)
{
    if(n==1){b[0]=1;return;}
    psqrt(a,b,invb,n>>1);
    memset(invb,0,4*n);pinv(b,invb,n);
    for(int i=0;i<n;i++)c[i]=a[i];
    for(int i=n;i<n+n;i++)c[i]=0;
    n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1);
    for(int i=0;i<n;i++)
        b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}

〖相关题目〗

1.【bzoj2179】FFT快速傅立叶

题意：给出两个n位10进制整数x和y，输出x*y。

分析：NTT裸题

#include<cstdio>
#include<algorithm>
#include<cstring>
#define LL long long
using namespace std;
const int N=131072+5;
const int mod=998244353;
int n,len,a[N],b[N],c[N];
char s[N];
int qmod(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
     int k=0;while((1<<k)<n)k++;
     for(int i=0;i<n;i++)
     {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
     }
     for(int l=2;l<=n;l<<=1)
     {
        int m=l>>1,nw=qmod(3,(mod-1)/l);
        if(f==-1)nw=qmod(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
     }
     if(f==-1)
     {
        int inv=qmod(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
     }
}
int main()
{
    scanf("%d",&len);
    n=1;while(n<len*2)n<<=1;
    scanf("%s",s+1);
    for(int i=1;i<=len;i++)a[len-i]=s[i]-'0';
    scanf("%s",s+1);
    for(int i=1;i<=len;i++)b[len-i]=s[i]-'0';
    ntt(a,n,1);ntt(b,n,1);
//  for(int i=0;i<n;i++)printf("%d %d\n",a[i],b[i]);
    for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod;
    ntt(c,n,-1);
    for(int i=0;i<n;i++)
        if(c[i]>=10)
        {
            c[i+1]+=c[i]/10;c[i]%=10;
            if(i==n-1)n++;
        }
    n--;while(!c[n]&&n)n--;
    for(int i=n;i>=0;i--)printf("%d",c[i]);
    return 0;
}

2.【51nod1135】原根

题意：给出1个质数P，找出P最小的原根。

分析：无

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
int mod,num,mx,tot,b[1005];
int qmod(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int main()
{
    scanf("%d",&mod);
    mx=(int)(sqrt(mod)+0.5);num=mod-1;
    for(int i=2;i<=mx;i++)
        if(num%i==0)
        {
            b[++tot]=i;
            while(num%i==0)num/=i;
        }
    if(num!=1)b[++tot]=num;
    for(int i=2;i<=mod;i++)
    {
        bool ok=true;
        for(int j=1;j<=tot;j++)
            if(qmod(i,(mod-1)/b[j])==1){ok=false;break;}
        if(ok){printf("%d",i);return 0;}
    }
    return 0;
}

3.【bzoj3992】[SDOI2015]序列统计

题意：见原题

分析：WerKeyTom_FTDの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cstdlib>
using namespace std;
const int N=2e4+5;
const int mod=1004535809;
int n,m,x,num,tmp;
int lg[N],b[N],c[N],ans[N],f[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int power(int a,int b,int mod)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int find(int n)
{
    int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0;
    for(int i=2;i<=mx;i++)
        if(num%i==0)
        {
            b[++tot]=i;
            while(num%i==0)num/=i;
        }
    if(num!=1)b[++tot]=num;
    for(int i=2;i<=n;i++)
    {
        bool ok=true;
        for(int j=1;j<=tot;j++)
            if(power(i,(n-1)/b[j],n)==1){ok=false;break;}
        if(ok)return i;
    }
    return 0;
}
void pre()
{
    int g=find(m),sum=1;
    for(int i=0;i<m-1;i++)
    {
        lg[sum]=i;
        sum=1ll*sum*g%m;
    }
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l,mod);
        if(f==-1)nw=power(nw,mod-2,mod);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2,mod);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void mul(int *a,int *b,int n)
{
    for(int i=0;i<n;i++)c[i]=b[i];
    ntt(a,n,1);ntt(c,n,1);
    for(int i=0;i<n;i++)a[i]=1ll*a[i]*c[i]%mod;
    ntt(a,n,-1);
    for(int i=m-1;i<n;i++)
        a[i%(m-1)]=(a[i%(m-1)]+a[i])%mod,a[i]=0;
}
void work(int *a,int b)
{
    int n=1;while(n<m+m)n<<=1;
    ans[0]=1;
    while(b)
    {
        if(b&1)mul(ans,a,n);
        mul(a,a,n);b>>=1;
    }
}
int main()
{
    n=read();m=read();x=read();num=read();
    pre();
    for(int i=1;i<=num;i++)
    {
        tmp=read();
        if(tmp)f[lg[tmp]]=1;
    }
    work(f,n);
    printf("%d",ans[lg[x]]);
    return 0;
}

4.【bzoj4555】[Tjoi2016&Heoi2016]求和

题意：见原题

分析：zltttttの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cstdlib>
using namespace std;
const int N=262144+5;
const int mod=998244353;
int n,m,ans;
int fac[N],inv[N],f[N],g[N],c[N];
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l);
        if(f==-1)nw=power(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void pinv(int *a,int *b,int n)
{
    if(n==0){b[0]=power(a[0],mod-2);return;}
    pinv(a,b,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
    ntt(c,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}
int main()
{
    scanf("%d",&m);
    n=1;while(n<=m)n<<=1;
    fac[0]=1;
    for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
    inv[n]=power(fac[n],mod-2);
    for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod;
    for(int i=1;i<=n;i++)f[i]=1ll*(mod-inv[i])*2%mod;
    f[0]=1;pinv(f,g,n);
    for(int i=0;i<=m;i++)ans=(ans+1ll*g[i]*fac[i]%mod)%mod;
    printf("%d",ans);
    return 0;
}

5.【bzoj3456】城市规划

题意：求n个点的无向简单连通图个数。

分析：PoPoQQQの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=262144+5;
const int mod=1004535809;
int n,mx;
int fac[N],inv[N],c[N];
int A[N],B[N],C[N];
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l);
        if(f==-1)nw=power(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void pinv(int *a,int *b,int n)
{
    if(n==1){b[0]=power(a[0],mod-2);return;}
    pinv(a,b,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
    ntt(c,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)
        b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}
int main()
{
    scanf("%d",&n);
    mx=1;while(mx<=n)mx<<=1;
    fac[0]=1;
    for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
    inv[n]=power(fac[n],mod-2);
    for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod;
    for(int i=0;i<=n;i++)
        C[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i]%mod;
    pinv(C,A,mx);
//  for(int i=1;i<=n;i++)printf("%d\n",A[i]);
    for(int i=1;i<=n;i++)
        B[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i-1]%mod;
    mx<<=1;ntt(A,mx,1);ntt(B,mx,1);
    memset(C,0,sizeof(C));
    for(int i=0;i<mx;i++)C[i]=1ll*A[i]*B[i]%mod;
    ntt(C,mx,-1);
    printf("%lld",1ll*C[n]*fac[n-1]%mod);
    return 0;
}

6.【bzoj3625】[Codeforces Round #250]小朋友和二叉树

题意：见原题

分析：wzq_QwQの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=262144+5;
const int mod=998244353;
int n,num,m,x,inv;
int f[N],g[N],invg[N],c[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void ntt(int *a,int n,int f)
{
    int k=0;while((1<<k)<n)k++;
    for(int i=0;i<n;i++)
    {
        int t=0;
        for(int j=0;j<k;j++)
            if(i&(1<<j))t|=(1<<(k-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1,nw=power(3,(mod-1)/l);
        if(f==-1)nw=power(nw,mod-2);
        for(int *p=a;p!=a+n;p+=l)
        {
            int w=1;
            for(int i=0;i<m;i++)
            {
                int t=1ll*p[m+i]*w%mod;
                p[m+i]=(p[i]-t+mod)%mod;
                p[i]=(p[i]+t)%mod;
                w=1ll*w*nw%mod;
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void pinv(int *a,int *b,int n)
{
    if(n==1){b[0]=power(a[0],mod-2);return;}
    pinv(a,b,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
    ntt(c,n,1);ntt(b,n,1);
    for(int i=0;i<n;i++)
        b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}
void psqrt(int *a,int *b,int *invb,int n)
{
    if(n==1){b[0]=1;return;}
    psqrt(a,b,invb,n>>1);
    memset(invb,0,4*n);pinv(b,invb,n);
    for(int i=0;i<n;i++)c[i]=a[i];
    for(int i=n;i<n+n;i++)c[i]=0;
    n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1);
    for(int i=0;i<n;i++)
        b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod;
    ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
}
int main()
{
    num=read();m=read();inv=power(2,mod-2);
    for(int i=1;i<=num;i++)x=read(),f[x]=1;
    for(int i=1;i<=m;i++)f[i]=(-f[i]*4%mod+mod)%mod;
    n=1;while(n<=m)n<<=1;
    f[0]=1;psqrt(f,g,invg,n);
    g[0]=(g[0]+1)%mod;
    memset(f,0,sizeof(f));
    pinv(g,f,n);
    for(int i=1;i<=m;i++)printf("%d\n",f[i]*2%mod);
    return 0;
}

【拉格朗日插值】

〖相关知识〗

给定 $n+1$ 个点 $(x_0,y_0),\cdots (x_n,y_n)$ ，求穿过它们的 $n$ 次函数。

令 $\ell _j(x)=\prod _{i=0,i\neq j}^{n}\frac{x-x_i}{x_j-x_i}$ ，拉格朗日插值多项式 $L(x)=\sum_{j=0}^{n}y_j \ell _j(x)$。

〖相关题目〗

1.【bzoj4559】成绩比较

题意：见原题

分析：cdsszjjの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int mod=1e9+7;
const int N=105;
int n,m,k;
int u[N],r[N],g[N],w[N];
int c[N][N],f[N][N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int lag(int u,int r)
{
    memset(g,0,sizeof(g));
    for(int x=1;x<=n+1;x++)
    {
        for(int i=1;i<=x;i++)
            g[x]=(g[x]+1ll*power(x-i,r-1)*power(i,n-r)%mod)%mod;
        if(u==x)return g[x];
    }
    for(int i=1;i<=n+1;i++)
    {
        w[i]=1;
        for(int j=1;j<=n+1;j++)
            if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
        w[i]=power(w[i],mod-2);
    }
    int sum=0;
    for(int i=1;i<=n+1;i++)
    {
        int tmp=1ll*w[i]*g[i]%mod;
        for(int j=1;j<=n+1;j++)
            if(i!=j)tmp=1ll*tmp*(u-j)%mod;
        sum=(sum+tmp)%mod;
    }
    return sum;
}
int C(int n,int m)
{
    if(n<0||m<0||n<m)return 0;
    return c[n][m];
}
int main()
{
    n=read();m=read();k=read();
    for(int i=1;i<=m;i++)u[i]=read();
    for(int i=1;i<=m;i++)r[i]=read();
    for(int i=0;i<=n;i++)c[i][0]=1;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=i;j++)
            c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
    f[0][0]=1;
    for(int i=1;i<=m;i++)
    {
        int tmp=lag(u[i],r[i]);
        for(int j=0;j<=n-1;j++)
        {
            for(int w=0;w<=j;w++)
                f[i][j]=(f[i][j]+1ll*C(n-1-w,j-w)*C(w,r[i]-1-j+w)%mod*f[i-1][w]%mod)%mod;
            f[i][j]=1ll*f[i][j]*tmp%mod;
        }
    }
    printf("%d",f[m][n-1-k]);
    return 0;
}

2.【bzoj2655】calc

题意：见原题

分析：ez_ywwの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const int N=505;
int m,n,mod,ans;
int w[N],f[N*2][N];
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int main()
{
    scanf("%d%d%d",&m,&n,&mod);
    f[0][0]=1;
    for(int i=1;i<=2*n;i++)
    {
        f[i][0]=f[i-1][0];
        for(int j=1;j<=n;j++)
            f[i][j]=(1ll*f[i-1][j-1]*i%mod*j%mod+f[i-1][j])%mod;
    }
    if(m<=2*n){printf("%d",f[m][n]);return 0;}
    for(int i=0;i<=2*n;i++)
    {
        w[i]=1;
        for(int j=0;j<=2*n;j++)
            if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
        w[i]=power(w[i],mod-2);
    }
    for(int i=0;i<=2*n;i++)
    {
        int sum=1ll*f[i][n]*w[i]%mod;
        for(int j=0;j<=2*n;j++)
            if(i!=j)sum=1ll*sum*(m-j)%mod;
        ans=(ans+sum)%mod;
    }
    printf("%d\n",ans);
    return 0;
}

3.【bzoj3453】tyvj 1858 XLkxc

题意：见原题

分析：DZYOの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iostream>
using namespace std;
const int N=150;
const int mod=1234567891;
int T,k,a,n,d,m,sum,ans;
int g[N],w[N],f[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
void work()
{
    k=read();a=read();n=read();d=read();
    for(int i=1;i<=k+3;i++)g[i]=power(i,k);
    for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod;
    for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod;
    for(int i=1;i<=k+3;i++)
    {
        w[i]=1;
        for(int j=1;j<=k+3;j++)
            if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
        w[i]=power(w[i],mod-2);
    }
    for(int x=0;x<=k+4;x++)
    {
        f[x]=0;m=(a+1ll*x*d%mod)%mod;
        for(int i=1;i<=k+3;i++)
        {
            sum=1ll*g[i]*w[i]%mod;
            for(int j=1;j<=k+3;j++)
                if(i!=j)sum=(1ll*sum*(m-j)%mod+mod)%mod;
            f[x]=(1ll*f[x]+sum)%mod;
        }
    }
    for(int i=1;i<=k+4;i++)f[i]=(1ll*f[i-1]+f[i])%mod;
    for(int i=0;i<=k+4;i++)
    {
        w[i]=1;
        for(int j=0;j<=k+4;j++)
            if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
        w[i]=power(w[i],mod-2);
    }
    ans=0;
    for(int i=0;i<=k+4;i++)
    {
        sum=1ll*f[i]*w[i]%mod;
        for(int j=0;j<=k+4;j++)
            if(i!=j)sum=(1ll*sum*(n-j)%mod+mod)%mod;
        ans=(1ll*ans+sum)%mod;
    }
    printf("%d\n",ans);
}
int main()
{
    T=read();
    while(T--)work();
    return 0;
}