【WC2016模拟】几何

Description

n<=60000,T<=5
时限0.8S

Solution

忽略掉题面最开始三个字
显然这题分为两部分
第一部分是求出Dp[i]表示i-多面体的选择方案数。
第二部分是把Dp[1]~Dp[n]组合起来。

第二部分显然可以用分治FFT来搞（求n个一次多项式的乘积），我们来看第一部分
考虑Dp[n]，枚举棱上选了多少条边，$Dp[n]=\sum_{k=0}^{4}C_4^k3^k\sum_{i=n}^{6n-12+k}C_{6n-12}^{j-k}$
考虑把后边那个部分写成减法的形式，$Dp[n]=\sum_{k=0}^{4}C_4^k3^k(2^{6n-12}-\sum_{i=0}^{n-1}C_{6n-12}^{j-k})$
观察到k很小，我们不妨直接求出$\sum_{i=0}^{n-1}C_{6n-12}^{j}$，剩下的加一些组合数即可。
设$Bad(n)=\sum_{i=0}^{n+1}C_{6n}^{i}$，显然我们要求的东西就是$Bad(n-2)$
考虑用范德蒙恒等式展开Bad

$$Bad(n)=\sum_{i=0}^{n+1}\sum_{j=0}^{6}C_{6n-6}^{i-j}C_6^j$$
$$=\sum_{j=0}^{6}C_6^j\sum_{i=0}^{n+1-j}C_{6(n-1)}^{i}$$
当j=1的时候后边的部分就是$Bad(n-1)$，j!=1的时候照样可以加一些组合数来完成。
这样我们就解决了第一部分。

然而这道题真正恶心的是第二部分，无良出题人竟然卡常！
首先，我们可以把$\omega$全部放到外面预处理，然后，我们其实有一个可以两次DFT做FFT的方法。
设C[i]=complex(A[i]+B[i],A[i]-B[i])，对C做DFT，然后自乘做IDFT，可以发现这时候的实部是我们要求的四倍。
然后你也不一定能卡过去~

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;
typedef long double db;

const int N=6*1e4+5,Mo=1e5+3;

const db pi=M_PI;

struct complex{
    db r,i;
    complex (db _r=0,db _i=0) {r=_r;i=_i;}
    friend complex operator + (complex x,complex y) {return complex(x.r+y.r,x.i+y.i);}
    friend complex operator - (complex x ,complex y) {return complex(x.r-y.r,x.i-y.i);}
    friend complex operator * (complex x,complex y) {return complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}A[N<<1],B[N<<1],T[N<<1],W[N<<1],t[N<<1];


int n,k,Len,len,lg,Id[18][N<<1],Dp[N],fact[Mo],inv[Mo],Bad[N];

int pwr(int x,int y) {
    int z=1;
    for(;y;y>>=1,x=(ll)x*x%Mo)
        if (y&1) z=(ll)z*x%Mo;
    return z;
}

int C(int m,int n) {
    if (n>m||m<0||n<0) return 0;
    if (m<Mo&&n<Mo) return (ll)fact[m]*inv[n]%Mo*inv[m-n]%Mo;
    return (ll)C(m/Mo,n/Mo)*C(m%Mo,n%Mo)%Mo;
}

void prepare() {
    fact[0]=1;fo(i,1,Mo-1) fact[i]=(ll)fact[i-1]*i%Mo;
    inv[Mo-1]=pwr(fact[Mo-1],Mo-2);fd(i,Mo-2,0) inv[i]=(ll)inv[i+1]*(i+1)%Mo;

    Bad[0]=1;
    fo(i,1,N) {
        int now=(Bad[i-1]+C(6*i-6,i+1))%Mo;
        fo(j,0,6) {
            (Bad[i]+=now*C(6,j)%Mo)%=Mo;
            (now+=Mo-C(6*i-6,i+1-j))%=Mo;
        }
    }

    Dp[1]=9;
    fo(i,2,N) {
        int sum=(pwr(2,6*i-12)-Bad[i-2]+Mo)%Mo;
        fo(k,0,4) {
            if (k) (sum+=C(6*i-12,i-k))%=Mo;
            (Dp[i]+=(ll)sum*C(4,k)*pwr(3,k)%Mo)%=Mo;
        }
    }

    for(Len=1,lg=0;Len<=N;) {
        Len<<=1;lg++;
        fo(i,0,Len-1) {
            int p=0;
            for(int j=i,k=0;k<lg;k++,j>>=1) p=(p<<1)+(j&1);
            Id[lg][i]=p;
        }
    }
    W[0]=complex(1,0);
    fo(i,1,Len) W[i]=complex(cos(2*i*pi/Len),sin(2*i*pi/Len));
}

int F[N<<1],G[N<<1],H[N<<1],Ans[N<<1];

void pre(int n) {
    for(len=1,lg=0;len<=n;len<<=1) lg++;
}

void DFT(complex *a,int flag) {
    fo(i,0,len-1) t[Id[lg][i]]=a[i];
    for(int m=2;m<=len;m<<=1) {
        int half=m/2,times=Len/m;
        fo(i,0,half-1) {
            complex w=(flag==1)?W[i*times]:W[Len-i*times];
            for(int j=i;j<len;j+=m) {
                complex u=t[j],v=t[j+half]*w;
                t[j]=u+v;t[j+half]=u-v;
            }
        }
    }
    fo(i,0,len-1) a[i]=t[i];
    if (flag==-1) fo(i,0,len-1) a[i].r/=len;
}

void FFT() {
    if (len<=32) {
        fo(i,0,len-1)
            fo(j,0,len-1)
                (Ans[i+j]+=(ll)G[i]*H[j]%Mo)%=Mo;
        return;
    }
    fo(i,0,len-1) T[i]=complex(G[i]+H[i],G[i]-H[i]);
    DFT(T,1);
    fo(i,0,len-1) T[i]=T[i]*T[i];
    DFT(T,-1);
    fo(i,0,len-1) Ans[i]=(ll)(T[i].r/4+0.5)%Mo;
}

void solve(int l,int r) {
    if (l>=r) return;
    int mid=l+r>>1;
    solve(l,mid);solve(mid+1,r);
    pre(r-l+1);
    G[0]=1;fo(i,l,mid) G[i-l+1]=F[i];
    H[0]=1;fo(i,mid+1,r) H[i-mid]=F[i];
    FFT();
    fo(i,l,r) F[i]=Ans[i-l+1];
    fo(i,0,len) G[i]=H[i]=Ans[i]=0;
}

int ty;
int main() {
    prepare();
    for(scanf("%d",&ty);ty;ty--) {
        scanf("%d%d",&n,&k);
        fo(i,1,n) F[i]=Dp[i];
        solve(1,n);
        int ans=0;
        fo(i,k,n) (ans+=F[i])%=Mo;
        printf("%d\n",(ans+Mo)%Mo);
    }
    return 0;
}