JZOJ 4349 几何

Description

定义一个$S$-四面体为一个每条边由S根依次首尾相连的木棍组成的四面体。
对一个$S$-四面体的一种合法的破坏方法为破坏不少于$S$根木棍且每个节点至少与两根木棍相连。
现在给出$n$个四面体，分别为$1$-四面体，$2$-四面体……$n$-四面体。
一种破坏这套四面体的合法方法为一次至少破坏掉$k$个四面体，问破坏者套四面体的方法有多少种。（不考虑空间同构因素，不考虑顺序因素，对某个四面体破坏了不同的棍视为不同的方案）

Data Constraint

Solution

先求出破坏一个四面体的方法。
首先可以暴力算出$1$-四面体有$9$种破坏方法。
对于$k$-四面体($k$>$1$)，由于要保证每个节点至少要与$2$根木棍相连，也就意味着一个节点连接的三根木棍要么全部完好，要么只能有一根被破坏。
于是我们可以枚举有多少个节点会被破坏掉$1$根木棍，若有$a$个节点会被破坏掉木棍，那么方案数为$C_{4}^a*3^a$*$left_k$，$left_{k,a}$为剩下的$6k$-$12$根木棍的合法破坏方法。

由于剩下的木棍中至少要破坏掉$k-a$根木棍，显然可得

$$left_{k,a}=\sum_{i=k-a}^{6k-12}C_{6k-12}^i$$
设 $$good(k)=\sum_{i=k}^{6k-12}C_{6k-12}^i$$
那显然 $left_{k,a}$显然可以由 $good(k)$加加减减几个组合数得到。

那我们怎么快速求出$good(k)$呢？
设

$$bad(k)=\sum_{i=0}^{k+1}C_{6k}^i$$
那么有
$$good(k)=\sum_{i=k}^{6k-12}C_{6k-12}^i=2^{6k-12}-\sum_{i=0}^{(k-2)+1}C_{6(k-2)}^i=2^{6k-12}-bad(k-2)$$
那我们的问题也就转换成了如何快速地求 $bad(k)$。

将$bad(k)$表示出来，并用范德蒙恒等式展开

$$bad(k)=\sum_{i=0}^{k+1}C_{6k}^i=\sum_{b=0}^6 \sum_{j=0}^{k+1-b}C_6^b*C_{6k-6}^j=\sum_{b=0}^6 C_6^b*(\sum_{j=0}^{k+1-b}*C_{6k-6}^j)$$
可以发现 $(\sum_{j=0}^{k+1-b}*C_{6k-6}^j)$与 $bad(k-1)$不过几个组合数的差别，那我们就可以从 $bad(k-1)$用 $O$($1$)的时间求出 $(\sum_{j=0}^{k+1-b}*C_{6k-6}^j)$。

那$bad(k)$就可以用$O$($1$)的时间通过$bad(k-1)$求出了。

那我们破坏一个$k$-四面体的合法方法也就可以顺势求出来了，记作$B_k$。
考虑构造多项式$A(x)$=$\Pi_{i=1}^{n} (B_ix+1)$
那么从中选择$k$个多面体的方案数显然为$[x^k]$$A(x)$。
现在构造出多项式$A(x)$，就能解决问题了，构造都是个老套路了，用分治$FFT$便可解决，时间复杂度$O$($n log^2 n$)

$PS$:垃圾出题人还卡常，害得我分治FFT超时，还要打两次DFT版本，分治区间过小时还要直接预处理，还要预处理蝴蝶变换，还要手打复数，恶xin……

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath> 

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)
#define eps 0.1 

using namespace std;
typedef long long ll;
typedef double db;
const ll N=62e3,M=4*N,K=6e4+10,mo=1e5+3;
const db pi=acos(-1);
ll xs[N],cs[N],jc[M],ny[M],bad[N],S[7],W[5],ycl[N];
ll f[1000];
int n,m,j,k,l,i,ttt;

struct Z{
    db x,y;
    Z(db _x=0,db _y=0){x=_x;y=_y;}
}a[M],b[M],c[M],t[M],tri[M];

Z operator +(Z a,Z b){return Z(a.x+b.x,a.y+b.y);}
Z operator -(Z a,Z b){return Z(a.x-b.x,a.y-b.y);}
Z operator *(Z a,Z b){return Z(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

void FFT(Z *a,int sig,int ws,int n)
{
    fo(i,0,n-1){
        int p=0;
        for(int tp=i,C=0;C<ws;++C,tp>>=1)p=(p<<1)+(tp&1);
        t[p]=a[i];
    }
    Z u;
    for(int m=2;m<=n;m<<=1){
        int half=m/2;
        fo(i,0,half-1){
            Z w(cos(i*pi*sig/half),sin(i*pi*sig/half));
            for(int j=i;j<n;j+=m){
                u.x=t[j+half].x*w.x-t[j+half].y*w.y;
                u.y=t[j+half].x*w.y+t[j+half].y*w.x;
                t[j+half].x=t[j].x-u.x; t[j+half].y=t[j].y-u.y;
                t[j].x=t[j].x+u.x; t[j].y=t[j].y+u.y;
            }
        }
    }
    fo(i,0,n-1)a[i]=t[i];
    if(sig==-1)fo(i,0,n-1)a[i].x=a[i].x/n;
}

void fz(int l,int r)
{
    if(l==r)return;
    if(l+750>r){
        f[0]=1;
        fo(i,1,r-l+1)f[i]=0;
        fo(i,l,r)
        fd(p,i-l+1,1)f[p]=(f[p-1]*xs[i]+f[p])%mo;
        fo(i,l,r)xs[i]=f[i-l+1];
        return;
    }
    int mid=(l+r)/2;
    fz(l,mid); fz(mid+1,r);
    int len=r-l+2;
    int y=0,cd=1;
    while(cd<len)cd<<=1,++y;
    a[0].x=a[0].y=1;
    fo(i,1,cd-1)a[i].x=a[i].y=0;
    fo(i,l,mid)a[i-l+1].x=xs[i];
    fo(i,mid+1,r)a[i-mid].y=xs[i];
    FFT(a,1,y,cd);
    db a1,a2,a3,a4;
    fo(i,0,cd-1){
        int j=(cd-i)&(cd-1);
        a1=(a[i].x+a[j].x)/2; a2=(a[i].y-a[j].y)/2;
        a4=(a[j].x-a[i].x)/2; a3=(a[i].y+a[j].y)/2;
        c[i].x=a1*a3-a2*a4;
        c[i].y=a1*a4+a2*a3;
    }
    FFT(c,-1,y,cd);
    fo(i,l,r)xs[i]=((ll)(c[i-l+1].x+eps))%mo;
}

ll ksm(ll o,ll t)
{
    ll y=1;
    for(;t;t>>=1,o=o*o%mo)
    if(t&1)y=y*o%mo;
    return y;
}

ll C(ll a,ll b)
{return b<0?0:(a<b?0:(a<mo?jc[a]*ny[b]%mo*ny[a-b]%mo:C(a%mo,b%mo)*C(a/mo,b/mo)%mo));}

int main()
{
    cin>>ttt;
    jc[0]=ny[0]=jc[1]=ny[1]=1;
    fo(i,2,mo-1)jc[i]=jc[i-1]*i%mo,ny[i]=ksm(jc[i],mo-2);
    S[0]=1; S[1]=6; S[2]=15; S[3]=20; S[4]=15; S[5]=6; S[6]=1;
    W[0]=1; W[1]=4; W[2]=6; W[3]=4; W[4]=1;
    bad[1]=22; 
    fo(i,2,K){
        ll op=(bad[i-1]+C(6*i-6,i+1))%mo;
        fo(b,0,6)
        if(i+1-b>=0)bad[i]=(bad[i]+op*S[b])%mo,op=(op-C(6*i-6,i+1-b)+mo)%mo;
    }
    fo(i,2,10){
        fo(j,0,i+1)ycl[i]=(ycl[i]+C(6*i,j))%mo;
    }
    fo(tt,1,ttt){
        cin>>n>>k;
        xs[1]=9; xs[2]=243; 
        fo(i,3,n){
            xs[i]=0;
            ll op=(ksm(2,6*i-12)-bad[i-2]+mo)%mo,fs=1;
            fo(l,0,4)xs[i]=(xs[i]+W[l]*op*fs)%mo,op=(op+C(6*i-12,i-l-1))%mo,fs=fs*3;
        }
        fz(1,n);
        ll ans=0;
        fo(i,k,n)ans=(ans+xs[i])%mo;
        printf("%lld\n",ans);
    }
}