2019 BUAA Spring Training 8 题解

Solved	A	B	C	D	E	F
6/6	Ø	Ø	Ø	Ø	Ø	Ø

O for passing during the contest
Ø for passing after the contest
! for attempted but failed
· for having not attempted yet

比赛链接

A HDU 4609 - 3-idiots

题目描述

给定一堆长度，任选三个问构成三角形的概率。

解题思路

建桶，求自身卷积，得到的便是每一个长度由原数组组合出来的方案数。去除自身选择自身、两个相互选择后得到不重复的方案数。

对初始数组排序，从小到大取，每次取该边为最长边，即需要满足剩下两个边都在$i$之前出现且其和大于$a_i$。于是先求和再去重（两个都在$i$后面，一个在前一个在后，一个是$i$另一个任取）即可。

AC代码

点击

#include<bits/stdc++.h>
typedef long long ll;
struct complex{
    double x,y;
    complex(double xx=0,double yy=0){x=xx;y=yy;}
    complex operator+(const complex a)const{return {x+a.x,y+a.y};}
    complex operator-(const complex a)const{return {x-a.x,y-a.y};}
    complex operator*(const complex a)const{return {x*a.x-y*a.y,x*a.y+y*a.x};}
};
#define N 463010
complex a[N],wn1,wnk,t;
int m,n,mx,limit;
int r[N];
void fft(complex *F,int sign){
    int i,j,len;
    for(i=0;i<limit;i++)
        if(i<r[i])std::swap(F[r[i]],F[i]);
    for(len=1;len<limit;len<<=1){
        wn1=complex(cos(acos(-1)/len),sign*sin(acos(-1)/len));
        for(j=0;j<limit;j+=(len<<1)){
            wnk=complex(1,0);
            for(i=j;i<j+len;i++){
                t=F[i+len]*wnk;
                F[i+len]=F[i]-t;
                F[i]=F[i]+t;
                wnk=wnk*wn1;
            }
        }
    }
    if(sign==-1)for(i=0;i<limit;i++)F[i].x/=limit;
}
ll b[N],s[N],p[N];
int main(){
    int i,n,T,x;
    scanf("%d",&T);
    while(T--){
        int mx=0;
        scanf("%d",&n);
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        for(i=0;i<n;i++){
            scanf("%d",&x),b[x]++;
            p[i]=x;if(x>mx)mx=x;
        }
        mx<<=1;limit=1;
        while(limit<=mx)limit<<=1;
        for(i=0;i<=limit;i++)a[i].x=b[i];
        for(i=0;i<=limit;i++)r[i]=(r[i>>1]>>1)|((i&1)?limit>>1:0);
        fft(a,1);
        for(i=0;i<=limit;i++)a[i]=a[i]*a[i];
        fft(a,-1);
        for(i=0;i<=limit;i++)s[i]=round(a[i].x);
        s[0]=0;
        for(i=0;i<n;i++)s[p[i]<<1]--;
        for(i=0;i<=limit;i++)s[i]>>=1;
        for(i=1;i<=limit;i++)s[i]+=s[i-1];
        ll ans=0;
        std::sort(p,p+n);
        for(i=0;i<n;i++){
            ans+=s[limit]-s[p[i]];
            ans-=1LL*i*(n-i-1);
            ans-=n-1;
            ans-=1LL*(n-i-1)*(n-i-2)/2;
        }
        printf("%.7f\n",ans*1.0/(1LL*n*(n-1)*(n-2)/6));
    }
    return 0;
}

B HDU 4914 - Linear recursive sequence

题目描述

$f[i]=0,i\leq 0$

$f[i]=a\times f[i-p]+b\times f[i-q],i>0$

给定$n,a,b,p,q(1≤n≤10^9,0≤a,b≤10^9,1≤p<q≤10^4)$，求$f[n]$。多组数据。

解题思路

一开始以为是线性递推，学了一下午打了个板子上去$T$掉了。

原理就是$k$阶矩阵$M$的任意幂次都能用$M$小于$k$次的幂次来线性表示，故可以通过特征多项式构造这个系数数列，利用快速幂+$FFT$解决。

TLE代码

点击

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
#define M 65536
struct comp{
    double x,y;
    comp():x(0),y(0){}
    comp(const double &_x,const double &_y):x(_x),y(_y){}
};
comp operator+(const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
comp conj(const comp &a){return comp(a.x,-a.y);}
double PI=acos(-1);
comp w[M+5];
int rev[M+5];
int lim,mx,mod;
void fft(comp *a,int n){
    int i,j,k,lyc;
    for(i=0;i<n;i++)if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(i=2,lyc=n>>1;i<=n;i<<=1,lyc>>=1)
        for(j=0;j<n;j+=i){
            comp *l=a+j,*r=a+j+(i>>1),*p=w;
            for(k=0;k<(i>>1);k++){
                comp tmp=*r**p;
                *r=*l-tmp,*l=*l+tmp;
                ++l,++r,p+=lyc;
            }
        }
}
void fft_prepare(){
    int i;
    for(i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)<<(mx-1));
    for(i=0;i<lim;i++)w[i]=comp(cos(2*PI*i/lim),sin(2*PI*i/lim));
}
void conv(int *x,int *y,int *z,int n,int m,int mod){
    int i,j;
    static comp ta[M+5],tb[M+5],a[M+5],b[M+5];
    static int X[M+5],Y[M+5];
    for(i=0;i<n;i++)X[i]=(x[i]+mod)%mod;
    for(;i<lim;i++)X[i]=0;
    for(i=0;i<m;i++)Y[i]=(y[i]+mod)%mod;
    for(;i<lim;i++)Y[i]=0;
    for(i=0;i<lim;i++){
        ta[i]=comp(X[i]&32767,X[i]>>15);
        tb[i]=comp(Y[i]&32767,Y[i]>>15);
    }
    fft(ta,lim);fft(tb,lim);
    for(i=0;i<lim;i++){
        j=(lim-i)%lim;//0的特判
        comp da,db,dc,dd;
          da=(ta[i]+conj(ta[j]))*comp(0.5,0);
        db=(ta[i]-conj(ta[j]))*comp(0,-0.5);
        dc=(tb[i]+conj(tb[j]))*comp(0.5,0);
        dd=(tb[i]-conj(tb[j]))*comp(0,-0.5);
        a[j]=da*dc+da*dd*comp(0,1);
        b[j]=db*dc+db*dd*comp(0,1);
    }
    fft(a,lim);fft(b,lim);
    for(i=0;i<lim;i++){
        ll da,db,dc,dd;
        da=(ll)(a[i].x/lim+0.5)%mod;
        db=(ll)(a[i].y/lim+0.5)%mod;
        dc=(ll)(b[i].x/lim+0.5)%mod;
        dd=(ll)(b[i].y/lim+0.5)%mod;
        z[i]=(((da+((db+dc)<<15)+(dd<<30))%mod)+mod)%mod;
    }
}
int qpow(ll x,int p){
    int ans=1;
    for(;p;p>>=1,x=x*x%mod)if(p&1)ans=ans*x%mod;
    return ans;
}
void polyinv(int *A,int *ans,int n){
    static int B[2][M+5];
    int bas,cur,i;
    memset(B,0,sizeof(B));
    B[0][0]=qpow(A[0],mod-2);
    for(mx=2,bas=1,lim=4,cur=1;bas<n*2;mx++,lim<<=1,cur^=1,bas<<=1){
        fft_prepare();
        memset(B[cur],0,sizeof(B[cur]));
        for(i=0;i<bas;i++)B[cur][i]=(ll)B[cur^1][i]*2%mod;
        conv(B[cur^1],B[cur^1],B[cur^1],bas,bas,mod);
        conv(B[cur^1],A,B[cur^1],bas,bas,mod);
        for(i=0;i<bas;i++)B[cur][i]=((ll)B[cur][i]-B[cur^1][i]+mod)%mod;
    }
    for(i=0;i<n;i++)ans[i]=(B[cur^1][i]+mod)%mod;
}
int flag;
void polydiv(int *a,int *b,int *H,int *R,int n,int m){
    int i;
    static int rF[M+5],rG[M+5],t1[M+5];
    memset(t1,0,sizeof(int)*n);
    for(i=0;i<=n;i++)rF[n-i]=a[i];
    if(!flag){
        flag=1;
        for(i=0;i<=m;i++)rG[m-i]=b[i];
        polyinv(rG,rG,n-m+1);
    }
    mx=0,lim=1;while(lim<n<<1)lim<<=1,mx++;
    fft_prepare();
    conv(rF,rG,H,n,n,mod);
    reverse(H,H+n-m+1);
    for(i=n-m+1;i<=n;i++)H[i]=0;
    mx=0,lim=1;while(lim<n+m)lim<<=1,mx++;
    fft_prepare();
    conv(b,H,t1,m,n,mod);
    for(i=0;i<m;i++)R[i]=((ll)a[i]-t1[i]+mod)%mod;
    for(;i<=n;i++)R[i]=0;
}
int a[M+5],g[M+5];
int D[M+5];
int res[M+5],bas[M+5],n,k;
void mtmul(int *a,int *b){
    mx=0,lim=1;
    while(lim<k*2)lim<<=1,mx++;
    fft_prepare();
    conv(a,b,a,k,k,mod);
    polydiv(a,g,D,a,k*2-2,k);
}
void qpow(int x){
    res[0]=bas[1]=1;
    for(;x;x>>=1,mtmul(bas,bas))
        if(x&1)mtmul(res,bas);
}
int main(){
    int i,aa,bb,p,q;
    mod=119;
    while(~scanf("%d%d%d%d%d",&n,&aa,&bb,&p,&q)){
        memset(res,0,sizeof(res));
        memset(bas,0,sizeof(bas));
        memset(g,0,sizeof(g));
        memset(a,0,sizeof(a));
        flag=0;
        k=q+1;n+=q;
        g[k]=118;
        g[k-p]=aa%mod;g[k-q]=bb%mod;
        for(i=0;i<k;i++)a[i]=1;
        qpow(n);
        ll ans=0;
        for(i=0;i<k;i++)(ans+=(ll)a[i]*res[i])%=mod;
        printf("%lld\n",ans);
    }
    return 0;
}

~~过后，不久，~~就发现没理解透这个板子，其实这个矩阵第一行就两个数……

~~然后就可以暴力了~~

AC代码

点击

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
#define M 65536
struct comp{
    double x,y;
    comp():x(0),y(0){}
    comp(const double &_x,const double &_y):x(_x),y(_y){}
};
comp operator+(const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int rev[M+5];
int lim,mx,mod;
void fft(comp *F,int sign){
    int i,j,len;
    for(i=0;i<lim;i++)if(i<rev[i])swap(F[rev[i]],F[i]);
    for(len=1;len<lim;len<<=1){
        comp wn1=comp(cos(acos(-1)/len),sign*sin(acos(-1)/len));
        for(j=0;j<lim;j+=(len<<1)){
            comp wnk=comp(1,0);
            for(i=j;i<j+len;i++){
                comp t=F[i+len]*wnk;
                F[i+len]=F[i]-t;
                F[i]=F[i]+t;
                wnk=wnk*wn1;
            }
        }
    }
    if(sign==-1)for(i=0;i<lim;i++)F[i].x/=lim;
}
void fft_prepare(){
    int i;
    for(i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)<<(mx-1));
}
int a[M+5],g[M+5];
int res[M+5],bas[M+5],n,k,p,q,aa,bb;
void mtmul(int *a,int *b){
    int i;
    static comp a1[M+5],a2[M+5];
    memset(a1,0,sizeof(a1));
    memset(a2,0,sizeof(a2));
    for(i=0;i<q;i++)a1[i].x=a[i],a2[i].x=b[i];
    fft(a1,1);fft(a2,1);
    for(i=0;i<lim;i++)a1[i]=a1[i]*a2[i];
    fft(a1,-1);
    for(i=0;i<lim;i++)a[i]=(int)((ll)(a1[i].x+0.5)%mod);
    for(i=2*q-2;i>=q;i--)if(a[i]){
        (a[i-q]+=a[i]*bb)%=mod;
        (a[i-p]+=a[i]*aa)%=mod;
    }
}
void qpow(int x){
    res[0]=bas[1]=1;
    for(;x;x>>=1,mtmul(bas,bas))
        if(x&1)mtmul(res,bas);
}
int main(){
    int i;
    mod=119;
    while(~scanf("%d%d%d%d%d",&n,&aa,&bb,&p,&q)){aa%=mod,bb%=mod;
        memset(res,0,sizeof(res));
        memset(bas,0,sizeof(bas));
        a[0]=1;
        for(i=1;i<q;i++)a[i]=(i-p>=0?(a[i-p]*aa+bb):aa+bb)%mod;
        if(n<q){
            printf("%d\n",a[n]);
            continue;
        }
        lim=1;mx=0;while(lim<=(q-1)<<1)lim<<=1,mx++;
        fft_prepare();
        qpow(n);
        int ans=0;
        for(i=0;i<q;i++)ans=(ans+res[i]*a[i])%mod;
        printf("%d\n",ans);
    }
    return 0;
}

C HDU 5322 - Hope

题目描述

定义一个排列$a_n$的价值如下：$a_i$向右边第一个大于它的数连边，所有连通分量的大小相乘后平方。定义 $f(n)$为$n$的全排列的价值之和，多次询问求 $f(n)$。

解题思路

先根据当前最大值的位置推出来$DP$式子：$dp[n]=\sum_{i=1}^ndp[n-i]i^2\frac{(n-1)!}{(n-i)!}=(n-1)!\sum_{i=0}^{n-1}i^2\frac{dp[i]}{i!}$。

然后就可以快乐分治$FFT$求解了。对于每一个区间$[l,r]$，先计算出左半边，后计算出左半边对于右半边$dp$的贡献（可用卷积表示），再分治解决右半边即可。

AC代码

点击

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
#define M 524290
struct comp{
    double x,y;
    comp():x(0),y(0){}
    comp(const double &_x,const double &_y):x(_x),y(_y){}
};
comp operator+(const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
comp conj(const comp &a){return comp(a.x,-a.y);}
double PI=acos(-1);
comp w[M+5];
int rev[M+5];
int lim,mx,mod;
void fft(comp *a,int n){
    int i,j,k,lyc;
    for(i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(i=2,lyc=n>>1;i<=n;i<<=1,lyc>>=1)
        for(j=0;j<n;j+=i){
            comp *l=a+j,*r=a+j+(i>>1),*p=w;
            for(k=0;k<(i>>1);k++){
                comp tmp=*r**p;
                *r=*l-tmp,*l=*l+tmp;
                ++l,++r,p+=lyc;
            }
        }
}
void fft_prepare(){
    int i;
    for(i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)<<(mx-1));
    for(i=0;i<lim;i++)w[i]=comp(cos(2*PI*i/lim),sin(2*PI*i/lim));
}
void conv(int *x,int *y,int *z,int n,int m,int mod){
    int i,j;
    mx=0;lim=1;while(lim<n+m+1)lim<<=1,mx++;
    fft_prepare();
    static comp ta[M+5],tb[M+5],a[M+5],b[M+5];
    static int X[M+5],Y[M+5];
    for(i=0;i<n;i++)X[i]=(x[i]+mod)%mod;
    for(;i<lim;i++)X[i]=0;
    for(i=0;i<m;i++)Y[i]=(y[i]+mod)%mod;
    for(;i<lim;i++)Y[i]=0;
    for(i=0;i<lim;i++){
        ta[i]=comp(X[i]&32767,X[i]>>15);
        tb[i]=comp(Y[i]&32767,Y[i]>>15);
    }
    fft(ta,lim);fft(tb,lim);
    for(i=0;i<lim;i++){
        j=(lim-i)%lim;
        comp da,db,dc,dd;
  		da=(ta[i]+conj(ta[j]))*comp(0.5,0);
        db=(ta[i]-conj(ta[j]))*comp(0,-0.5);
        dc=(tb[i]+conj(tb[j]))*comp(0.5,0);
        dd=(tb[i]-conj(tb[j]))*comp(0,-0.5);
        a[j]=da*dc+da*dd*comp(0,1);
        b[j]=db*dc+db*dd*comp(0,1);
    }
    fft(a,lim);fft(b,lim);
    for(i=0;i<lim;i++){
        ll da,db,dc,dd;
        da=(ll)(a[i].x/lim+0.5)%mod;
        db=(ll)(a[i].y/lim+0.5)%mod;
        dc=(ll)(b[i].x/lim+0.5)%mod;
        dd=(ll)(b[i].y/lim+0.5)%mod;
        z[i]=(((da+((db+dc)<<15)+(dd<<30))%mod)+mod)%mod;
    }
}
int dp[M+5]={1},fac[M+5]={1},inv[M+5]={1};
ll qp(ll a,int p){
    ll ans=1;
    for(;p;p>>=1,a=a*a%mod)if(p&1)ans=ans*a%mod;
    return ans;
}
int a[M+5],b[M+5];
void solve(int l,int r){
    if(l==r)return;
    int i,mid=(l+r)>>1;
    solve(l,mid);
    for(i=0;i<=mid-l;i++)a[i]=(ll)dp[i+l]*inv[i+l]%mod;
    conv(a,b,a,mid-l+1,r-l+1,mod);
    for(i=mid+1;i<=r;i++)dp[i]=(dp[i]+(ll)a[i-l]*fac[i-1])%mod;
    solve(mid+1,r);
}
int main(){
    int i,n;
    mod=998244353;
    for(i=1;i<M;i++)
        fac[i]=(ll)fac[i-1]*i%mod,inv[i]=qp(fac[i],mod-2),
        b[i]=(ll)i*i%mod;
    solve(0,100000);
    while(~scanf("%d",&n))printf("%d\n",dp[n]);
    return 0;
}

D HDU 4656 - Evaluation

题目描述

$x_k=b\times c^{2k}+d $
$F(x)=a_0x_0+a_1x_1+a_2 x_2+…+a_{n-1 }x_{n-1} $

给定$n,a_i,b,c,d$，求序列$F(x_i)$。

解题思路

神仙题$orzzzzzzz$

~~$FFT$——从入门到放弃.jpg~~

整整连学带调写了一天半，终于过了

我们先来化一下式子：

$F(x_k)=\sum_{i=0}^{n-1}a_ix_k^i$

$=\sum_{i=0}^{n-1}a_i\sum_{j=0}^{i}b^jc^{2kj}d^{i-j}C_j^i$

$=\sum_{j=0}^{n-1}\frac{b^jc^{2kj}}{j!}\sum_{i=j}^{n-1}a_ii!\frac{d^{i-j}}{(i-j)!}$

设$A(i)=a_ii!,B(i)=\frac{d^i}{i!},p(j)=\sum_{i=j}^{n-1}A(i)B(i-j)$

显然可将$A$翻转：$A’(i)=A(n-i)$

则$p(n-j)=\sum_{i=1}^{n-j}A’(i)B(n-i-j)$，设$p’(i)=p(n-i)$，则$p’(j)=\sum_{i=1}^{j}B(j-i)$，可以用卷积求出。

$F(x_k)=\sum_{j=0}^{n-1}\frac{b^jc^{2kj}}{j!}p(j)$

$=\sum_{j=0}^{n-1}\frac{b^jc^{k^2+j^2-(k-j)^2}}{j!}p(j)$

$=c^{k^2}\sum_{j=0}^{n-1}\frac{q(j)}{c^{(k-j)^2}}$，其中$q(j)=p(j)\frac{b^jc^{j^2}}{j!}$。

于是设$g(i)=\frac1{i^2}$，$\frac{F(x_k)}{c^{k^2}}=\sum_{j=0}^{n-1}q(j)g(k-j)$，为了避免负数情况，将$g(i)$移动$n$（参考$ZJOI2014$力），即$g’(i)=g(i-n)$，最后$[n,2n-1]$区间即为所求卷积。

其中需要对给定质数求模，故用拆系数$FFT$进行运算，再用共轭优化$DFT$次数…

详情参考毛啸2016集训队论文《再探快速傅里叶变换》。

来说几个坑点吧：拆系数的时候注意正负号，注意零的特判，$c$的$i^2$次幂求的时候用$qp(c,i^2%(mod-1))$，以及这题卡空间……

MLE代码

点击

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define N 524290
#define M N
struct comp{
    double x,y;
    comp():x(0),y(0){}
    comp(const double &_x,const double &_y):x(_x),y(_y){}
};
comp operator+(const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
comp conj(const comp &a){return comp(a.x,-a.y);}
double PI=acos(-1);
comp w[M];
int rev[M];
int A[M],B[M],C[M],lim,mx,mod=1000003;
void fft(comp *a,int n){
    int i,j,k,lyc;
    for(i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(i=2,lyc=n>>1;i<=n;i<<=1,lyc>>=1)
        for(j=0;j<n;j+=i){
            comp *l=a+j,*r=a+j+(i>>1),*p=w;
            for(k=0;k<(i>>1);k++){
                comp tmp=*r**p;
                *r=*l-tmp,*l=*l+tmp;
                ++l,++r,p+=lyc;
            }
        }
}
void fft_prepare(){
    int i;
    for(i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)<<(mx-1));
    for(i=0;i<lim;i++)w[i]=comp(cos(2*PI*i/lim),sin(2*PI*i/lim));
}
comp a[M+5],b[M+5],dfta[M+5],dftb[M+5],dftc[M+5],dftd[M+5];
void conv(int *x,int *y,int *z){
    int i,j;
    for(i=0;i<lim;i++)(x[i]+=mod)%=mod,(y[i]+=mod)%=mod;
    for(i=0;i<lim;i++){
        a[i]=comp(x[i]&32767,x[i]>>15);
        b[i]=comp(y[i]&32767,y[i]>>15);
    }
    fft(a,lim);fft(b,lim);
    for(i=0;i<lim;i++){
        j=(lim-i)%lim;//0的特判
        comp da,db,dc,dd;
        da=(a[i]+conj(a[j]))*comp(0.5,0);
        db=(a[i]-conj(a[j]))*comp(0,-0.5);
        dc=(b[i]+conj(b[j]))*comp(0.5,0);
        dd=(b[i]-conj(b[j]))*comp(0,-0.5);
        dfta[j]=da*dc;
        dftb[j]=da*dd;
        dftc[j]=db*dc;
        dftd[j]=db*dd;
    }
    for(i=0;i<lim;i++){
        a[i]=dfta[i]+dftb[i]*comp(0,1);
        b[i]=dftc[i]+dftd[i]*comp(0,1);
    }
    fft(a,lim);fft(b,lim);
    for(i=0;i<lim;i++){
        ll da,db,dc,dd;
        da=(ll)(a[i].x/lim+0.5)%mod;
        db=(ll)(a[i].y/lim+0.5)%mod;
        dc=(ll)(b[i].x/lim+0.5)%mod;
        dd=(ll)(b[i].y/lim+0.5)%mod;
        z[i]=(da+((db+dc)<<15)+(dd<<30))%mod;
    }
}
int arr[N],fac[N]={1},inv[N]={1},bpow[N]={1},dpow[N]={1};
int qp(int x,int p){
    int ans=1;
    for(;p;p>>=1,x=(ll)x*x%mod)if(p&1)ans=(ll)ans*x%mod;
    return ans;
}
int Ai[N],Bi[N],P[N],p[N],ans[N],m[N],cpow[N]={1};
int main(){
    int i,n,b,c,d;
    scanf("%d%d%d%d",&n,&b,&c,&d);
    while((1<<mx)<4*n)mx++;
    lim=1<<mx;
    fft_prepare();
    for(i=1;i<n;i++)fac[i]=(ll)fac[i-1]*i%mod,inv[i]=qp(fac[i],mod-2);
    for(i=1;i<n;i++)dpow[i]=(ll)dpow[i-1]*d%mod,bpow[i]=(ll)bpow[i-1]*b%mod,cpow[i]=qp(c,(ll)i*i%(mod-1));
    for(i=0;i<n;i++)scanf("%d",&arr[i]);
    for(i=0;i<=n;i++)Ai[i]=(ll)arr[n-i]*fac[n-i]%mod,Bi[i]=(ll)dpow[i]*inv[i]%mod;
    conv(Ai,Bi,P);
    for(i=0;i<n;i++)p[i]=(ll)P[n-i]*bpow[i]%mod*inv[i]%mod*cpow[i]%mod;
    for(i=0;i<n+n;i++)m[i]=qp(i-n>=0?cpow[i-n]:cpow[n-i],mod-2);
    conv(m,p,ans);
    for(i=0;i<n;i++){
        ans[i+n]=(ll)ans[i+n]*cpow[i]%mod;
        printf("%d\n",(ans[i+n]+mod)%mod);
    }
    return 0;
}

疯狂删减数组，勉强卡过$65536kb$内存限制。

AC代码

点击

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define N 524290
#define M N
struct comp{
    double x,y;
    comp():x(0),y(0){}
    comp(const double &_x,const double &_y):x(_x),y(_y){}
};
comp operator+(const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
comp conj(const comp &a){return comp(a.x,-a.y);}
double PI=acos(-1);
comp w[M];
int rev[M];
int lim,mx,mod=1000003;
void fft(comp *a,int n){
    int i,j,k,lyc;
    for(i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(i=2,lyc=n>>1;i<=n;i<<=1,lyc>>=1)
        for(j=0;j<n;j+=i){
            comp *l=a+j,*r=a+j+(i>>1),*p=w;
            for(k=0;k<(i>>1);k++){
                comp tmp=*r**p;
                *r=*l-tmp,*l=*l+tmp;
                ++l,++r,p+=lyc;
            }
        }
}
void fft_prepare(){
    int i;
    for(i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)<<(mx-1));
    for(i=0;i<lim;i++)w[i]=comp(cos(2*PI*i/lim),sin(2*PI*i/lim));
}
comp a[M],b[M],ta[M],tb[M];
void conv(int *x,int *y,int *z){
    int i,j;
    for(i=0;i<lim;i++)(x[i]+=mod)%=mod,(y[i]+=mod)%=mod;
    for(i=0;i<lim;i++){
        ta[i]=comp(x[i]&32767,x[i]>>15);
        tb[i]=comp(y[i]&32767,y[i]>>15);
    }
    fft(ta,lim);fft(tb,lim);
    for(i=0;i<lim;i++){
        j=(lim-i)%lim;//0的特判
        comp da,db,dc,dd;
        da=(ta[i]+conj(ta[j]))*comp(0.5,0);
        db=(ta[i]-conj(ta[j]))*comp(0,-0.5);
        dc=(tb[i]+conj(tb[j]))*comp(0.5,0);
        dd=(tb[i]-conj(tb[j]))*comp(0,-0.5);
        a[j]=da*dc+da*dd*comp(0,1);
        b[j]=db*dc+db*dd*comp(0,1);
    }
    fft(a,lim);fft(b,lim);
    for(i=0;i<lim;i++){
        ll da,db,dc,dd;
        da=(ll)(a[i].x/lim+0.5)%mod;
        db=(ll)(a[i].y/lim+0.5)%mod;
        dc=(ll)(b[i].x/lim+0.5)%mod;
        dd=(ll)(b[i].y/lim+0.5)%mod;
        z[i]=(da+((db+dc)<<15)+(dd<<30))%mod;
    }
}
int arr[N],fac[N]={1},inv[N]={1};
int qp(int x,int p){
    int ans=1;
    for(;p;p>>=1,x=(ll)x*x%mod)if(p&1)ans=(ll)ans*x%mod;
    return ans;
}
int Ai[N],Bi[N],P[N],cpow[N]={1};
int main(){
    int i,n,b,c,d;
    scanf("%d%d%d%d",&n,&b,&c,&d);
    while((1<<mx)<4*n)mx++;
    lim=1<<mx;
    fft_prepare();
    for(i=1;i<n;i++)fac[i]=(ll)fac[i-1]*i%mod,inv[i]=qp(fac[i],mod-2);
    for(i=1;i<n;i++)cpow[i]=qp(c,(ll)i*i%(mod-1));
    for(i=0;i<n;i++)scanf("%d",&arr[i]);
    for(i=0;i<=n;i++)Ai[i]=(ll)arr[n-i]*fac[n-i]%mod,Bi[i]=(ll)qp(d,i)*inv[i]%mod;
    conv(Ai,Bi,P);
    Ai[n]=0;
    for(i=0;i<n;i++)Ai[i]=(ll)P[n-i]*qp(b,i)%mod*inv[i]%mod*cpow[i]%mod;
    for(i=0;i<n+n;i++)Bi[i]=qp(i-n>=0?cpow[i-n]:cpow[n-i],mod-2);
    conv(Ai,Bi,P);
    for(i=0;i<n;i++){
        P[i+n]=(ll)P[i+n]*cpow[i]%mod;
        printf("%d\n",(P[i+n]+mod)%mod);
    }
    return 0;
}

E CF662C - Binary Table

题目描述

给定一个$n\times m(1\leq m\leq 20,1\leq n\leq 100000)$的$01$矩阵，每次可以反转任意一列一行，求最终$1$的最少个数。

解题思路

看到$m$范围极小，考虑枚举每一行的翻转状态$S$。

由于每一行反转后，后面的所有行答案便已经确定，考虑求出$ans[S]$。

那么可以设$a[i]$为$i$这个状态出现的次数，$b[i]$为$i$这个状态中对应的$1$的个数和$0$的个数中的较小值，则$ans[s]=\sum_{j\oplus k=s}a[j]\times b[k]$，用$FWT$板子即可求出。

AC代码

点击

#include<bits/stdc++.h>
#define N (1<<20)+10
typedef long long ll;
ll a[N],b[N],c[N];
void fwtxor(int f,int l,ll a[]){
    int i,j,k;ll t;
    if(f==-1)f=2;
    for(i=1;i*2<=l;i<<=1)
        for(j=0;j<l;j+=i*2)
            for(k=0;k<i;k++)
                t=a[j+k],
                a[j+k]=(t+a[j+k+i])/f,
                a[j+k+i]=(t-a[j+k+i])/f;
}
char ch[22][100010];
int main(){
    int i,j,n,m;
    scanf("%d%d",&n,&m);
    for(i=0;i<n;i++)scanf("%s",ch[i]);
    for(j=0;j<m;j++){
        int now=0;
        for(i=0;i<n;i++)now=(now<<1)+ch[i][j]-'0';
        a[now]++;
    }
    for(i=1;i<1<<n;i++)b[i]=b[i/2]+(i&1);
    for(i=0;i<1<<n;i++)b[i]=std::min(b[i],n-b[i]);
    ll ans=n*m;
    n=1<<n;
    fwtxor(1,n,a);fwtxor(1,n,b);
    for(i=0;i<n;i++)c[i]=a[i]*b[i];
    fwtxor(-1,n,c);
    for(i=0;i<n;i++)ans=std::min(ans,c[i]);
    printf("%I64d",ans);
    return 0;
}

F CF914G - Sum the Fibonacci

题目描述
给定一个集合$S$，求出使得$(s_a|s_b)&(s_c)&(s_d\oplus s_e)==2^p,s_a&s_b==0$的$fib(s_a|s_b)\times fib(s_c)\times fib(s_d\oplus s_e)$的和，其中$s_a,s_b,s_c,s_d,s_e$均为$S$的子集。$s_i<2^{17}$。

解题思路
非常有意思的$FMT+FWT$练习题。先求出每个数出现的次数$s_i$，求出$s_a,s_b$的子集卷积，$s_d,s_e$的对称差（异或）卷积，分别乘上$fib$再进行（或）变换，乘起来再做逆运算即可。

既然是写博，就要自己在复习一遍，故下面是对$FMT$和子集卷积的理解。

首先是$FMT$。

定义$H_S=\sum_{L\subseteq S}\sum_{R\subseteq S}f_L\times g_R(L\bigcup R\subseteq S)$（①），我们要求出这个卷积$H$。
定义莫比乌斯变换为$f_S^{‘}=\sum_{T\subseteq S}f_T$。
对于①式，左右两端取莫比乌斯变换：
$H_S^{‘}=\sum_{L\subseteq S}\sum_{R\subseteq S}f_L\times g_R$
故可以进行变换，卷积运算。变换过程：
设 $f_S^{‘}(i)=\sum_{T\subseteq S}f_T((S-T)\subseteq(1,2,…,i))$
那么，对于任意不包含${i}$的集合$S$，有$f_S^{‘}(i)=f_S^{‘}(i-1)$，$f_{S\bigcup {i}}^{‘}(i)=f_{S\bigcup {i}}^{‘}(i-1)+f_S^{‘}(i-1)$，这样就可以$O(n\times 2^n)$解决了。反演同理。

对于子集卷积：

定义$H_S=\sum_{L\subseteq S}f_L\times g_{S-L}$，我们现在要求出这个卷积$H$。因为需要满足交集为空，所以增设一维表示集合大小。设$f_{i,S}$为集合中有$i$个元素、集合表示为$S$的值，那么求卷积的时候$h_{i,S}=\sum_{j=0}^i f_{j,S}\times g_{i-j,S}$（从$f$中选出$j$个元素，从$g$中选出$i-j$个元素，只有他们互不相交的情况下才能够保证$h_{i,S}$合理）。（非常玄学，怎么想出来的）
那么最后的结果就是$h_{cnt_i,i}$这个数组了。

AC代码
点击
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
#define N 18
#define M (1<<17)
#define P M+10
ll a[P],b[P],w=1000000007,subset[P],fib[P]={0,1};
ll Xor[P];
ll div(int x,ll p){
    if(x==1)return p;
    if(p&1)return (p+w)/2%w;
    else return p/2;
}
void fwtand(ll a[],int l,int f){
    int i,j,k;
    for(i=1;i*2<=l;i<<=1)
        for(j=0;j<l;j+=i*2)
            for(k=0;k<i;k++)
                (a[j+k]+=a[j+k+i]*f)%=w;
}
void fwtxor(ll a[],int l,int f){
    int i,j,k;ll t;
    if(f==-1)f=2;
    for(i=1;i*2<=l;i<<=1)
        for(j=0;j<l;j+=i*2)
            for(k=0;k<i;k++)
                t=a[j+k],
                (a[j+k]=div(f,t+a[j+k+i]))%=w,
                (a[j+k+i]=div(f,t-a[j+k+i]))%=w;
}
void fmt(ll a[],int n,int f){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<1<<n;j++)if(j&(1<<i))a[j]=(a[j]+f*a[j^(1<<i)])%w;
}
ll f[20][P],cnt[P];
void subsetconvolution(int n,ll a[],ll c[]){
    int i,j,k;
    for(i=0;i<1<<n;i++)cnt[i]=cnt[i>>1]+(i&1);
    for(i=0;i<1<<n;i++)f[cnt[i]][i]=a[i];
    for(i=0;i<=n;i++)fmt(f[i],n,1);
    for(i=n;i>=0;i--){
        for(j=0;j<1<<n;j++){
            ll x=0;
            for(k=0;k<=i;k++)(x+=f[k][j]*f[i-k][j])%=w;
            f[i][j]=x;
        }
    }
    for(i=0;i<=n;i++)fmt(f[i],n,-1);
    for(i=0;i<1<<n;i++)c[i]=f[cnt[i]][i];
}
int main(){
    int i,n,x;
    for(i=2;i<=M;i++)fib[i]=(fib[i-1]+fib[i-2])%w;
    scanf("%d",&n);
    for(i=0;i<n;i++)
        scanf("%d",&x),b[x]++,a[x]++;
    subsetconvolution(17,a,subset);
    fwtxor(a,M,1);
    for(i=0;i<M;i++)Xor[i]=a[i]*a[i]%w;
    fwtxor(Xor,M,-1);
    for(i=0;i<M;i++){
        subset[i]=subset[i]*fib[i]%w;
        b[i]=b[i]*fib[i]%w;
        Xor[i]=Xor[i]*fib[i]%w;
    }
    fwtand(subset,M,1);
    fwtand(b,M,1);
    fwtand(Xor,M,1);
    for(i=0;i<M;i++)b[i]=b[i]*subset[i]%w*Xor[i]%w;
    fwtand(b,M,-1);
    ll ans=0;
    for(i=1;i<M;i<<=1)
        ans=(ans+b[i])%w;
    printf("%I64d",(ans+w)%w);
    return 0;
}