[XSY 1501][组合数＋分治FFT] 几何

要是我们能预处理出i－四面体的方案数，那么我们查询时分治FFT即可
注意FFT时要预处理单位复数根来保证精度
考虑如何预处理出i－四面体的方案数
我们可以打表得到一个num数组
num[0]=1
num[1]=12
num[2]=54
num[3]=108
num[4]=81
numi代表2-四面体去掉i条边的方案数
那么我们可以得到k-四面体的方案数(k>=3)：
${\sum }_{i=0}^{4}nu{m}_i\ast {\sum }_{j=max(0,k-i)}^{6\ast k-12}\left(\genfrac{}{}{0ex}{}{6\ast k-12}{j}\right)$
k=1,2时，方案数分别为9,243
我们的问题转化成了如何对每个k求
${\sum }_{i=0}^k\left(\genfrac{}{}{0ex}{}{6\ast k-12}{i}\right)$
可以轻松地$o(n)$求出
时间复杂度:$O(n+nlo{g}_{2}nlo{g}_{2}n)$

#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
const int Mod=100003;
int n,K;
#define Maxn 60010
#define pi acos(-1)
int f[Maxn];
int val[Maxn];
int fact[Maxn*2],inv[Maxn*2];
int num[5]={1,12,54,108,81};
inline int c(int i,int j){
    if(i<j)return 0;
    return 1ll*fact[i]*inv[i-j]%Mod*inv[j]%Mod;
}
int C(int i,int j){
    if(!i||!j)return 1;
    return 1ll*c(i%Mod,j%Mod)*C(i/Mod,j/Mod)%Mod;
}
inline int FP(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;
}

struct cp{
    double a,b;
    friend cp operator +(cp &num1,cp &num2){return (cp){num1.a+num2.a,num1.b+num2.b};}
    friend cp operator -(cp &num1,cp &num2){return (cp){num1.a-num2.a,num1.b-num2.b};}
    friend cp operator *(cp &num1,cp &num2){return (cp){num1.a*num2.a-num1.b*num2.b,num1.a*num2.b+num1.b*num2.a};}
    friend cp operator /(cp &num,int t){return (cp){num.a/t,num.b/t};}
}t1[Maxn*2],t2[Maxn*2],*wr[16];
int rev[Maxn*2],len,bit;
inline void wpre(){
    for(int i=1,k=0;i<len;i<<=1,++k)
        if(wr[k]==0){
            wr[k]=new cp[i];
            for(int j=0;j<i;++j)wr[k][j]=(cp){cos(j*pi/i),sin(j*pi/i)};
        }
}
inline void FFT(cp *A,int t){
    for(int i=0;i<len;++i)
        if(i<rev[i])swap(A[i],A[rev[i]]);
    wpre();
    for(int i=1,f=0;i<len;i<<=1,++f){
        for(int j=0;j<len;j+=i<<1){
            for(int k=0;k<i;++k){
                cp x=A[j+k];
                cp w=wr[f][k];
                w.b*=t;
                cp y=w*A[j+k+i];
                A[j+k]=x+y;
                A[j+k+i]=x-y;
            }
        }
    }
    if(t==-1)
        for(int i=0;i<len;++i)A[i]=A[i]/len;
}
inline void Mul(cp *A,cp *B){
    FFT(A,1);FFT(B,1);
    for(int i=0;i<len;++i)A[i]=A[i]*B[i];
    FFT(A,-1);
}

vector<int> g[Maxn<<2];
void solve(int k,int l,int r){
    g[k].clear();
    if(l==r){
        g[k].push_back(1);
        g[k].push_back(val[l]);
        return;
    }
    int mid=(l+r)>>1;
    solve(k<<1,l,mid);
    solve(k<<1|1,mid+1,r);

    int l1=g[k<<1].size()-1,l2=g[k<<1|1].size()-1;
    len=1;bit=0;
    while(len<=l1+l2)len<<=1,bit++;
    for(int i=0;i<len;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    for(int i=0;i<len;++i)t1[i]=t2[i]=(cp){0,0};
    for(int i=0;i<=l1;++i)t1[i]=(cp){1.0*g[k<<1][i],0};
    for(int i=0;i<=l2;++i)t2[i]=(cp){1.0*g[k<<1|1][i],0};
    Mul(t1,t2);
    for(int i=0;i<=l1+l2;++i){
        ll now=(ll)(t1[i].a+0.5);
        g[k].push_back(now%Mod);
    }
}


int main(){
    fact[0]=1;
    for(register int i=1;i<Mod;++i)fact[i]=1ll*fact[i-1]*i%Mod;
    inv[0]=inv[1]=1;
    for(register int i=2;i<Mod;++i)inv[i]=1ll*(Mod-Mod/i)*inv[Mod%i]%Mod;
    for(register int i=2;i<Mod;++i)inv[i]=1ll*inv[i-1]*inv[i]%Mod;
    val[1]=9;val[2]=243;
    val[3]=42;
    for(register int i=4;i<=60000;++i){
        int now=1,tmp=val[i-1];
        for(register int j=i-1;j>=max(i-5,0);--j){
            now=now+C(6,i-j);
            if(now>=Mod)now-=Mod;
            val[i]=(val[i]+1ll*C(6*(i-1)-12,j)*now)%Mod;
            tmp=(tmp-C(6*(i-1)-12,j)+Mod)%Mod;
        }
        val[i]=(64ll*tmp+val[i])%Mod;
        val[i]=(val[i]+C(6*(i-1)-12,i))%Mod;
    }
    for(register int i=3;i<=60000;++i){
        int now=(FP(2,6*i-12)-val[i]+Mod)%Mod;
        val[i]=0;
        for(register int j=0;j<5;++j){
            if(i>=j)now=(now+C(6*i-12,i-j));
            if(now>=Mod)now-=Mod;
            val[i]=(val[i]+1ll*now*num[j])%Mod;
        }
    }
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&K);
        solve(1,1,n);
        int Ans=0;
        for(int i=K;i<=n;++i)Ans=(Ans+g[1][i])%Mod;
        printf("%d\n",Ans);
    }
    return 0;
}