几何
Portal --> broken qwq
Decsription
Solution
首先当然是要膜lyy啦%%%
咳咳
（为了避免混淆，这里强调一下接下来的内容中小\(k\)是。。各种枚举什么的时候用的，大\(K\)才是题目中给的\(K\)，\(n\)和\(N\)同理！）
如果说我们知道了破坏一个k-四面体的方案数\(d[k]\)，那么统计答案的时候我们就可以考虑用递推之类的方式处理了（至于细节先不管，放到后面考虑），那所以首先来考虑破坏一个k-四面体的情况，这里我们需要分两大类讨论：
（1）\(k=1\)：这个时候比较特殊，因为每条棱上面只有一根木棍，连接的两个顶点会相互影响，所以我们将这种情况单独拿出来考虑，从样例手算一下可以得到\(d[1]=9\)
（2）\(k>1\)：这个时候每个顶点不会像（1）中那样相互影响（因为每条棱至少有两根木棒嘛），那么这个时候我们可以把统计的答案分成两个部分：顶点和棱
首先是顶点（指的是。。与顶点相连的木棒），假设这类木棒中取走了\(a\)个：总共有\(4\)个点，每个点有\(3\)个可取的木棒，根据题目的限制，显然每个顶点的\(3\)个木棒中最多只有一个能被取走，所以我们可以得到取走\(a\)个木棒的方案数为：
然后是棱（指的是。。不与顶点相连的木棒），对应的这里取走的数量就应该是至少\(k-a\)个：去掉之前的与顶点相连的木棒，不与顶点相连的木棒总共有\(6k-12\)个，为了方便，我们定义一个\(s(6k-12,k-a)\)表示从\(6*k-12\)个木棒中取至少\(k-a\)个木棒的方案数，那么有：
所以我们可以得到：
现在的问题就是怎么快速求后面的那个玩意，也就是\(s\)
注意到这个其实是一个。。杨辉三角上面一行的一个后缀和，而根据组合数的递推式子\(\binom i j=\binom {i-1}{j-1}+\binom {i-1}{j}\)，如果说我们知道\(sum_i=\sum\limits_{k=j}^{i}\binom i k\)，那么\(sum_{i+1}=\sum\limits_{k=j}^{i+1}\binom {i+1} k\)是可以直接得到的，乘个\(2\)然后再加上\(\binom i {j-1}\)就好了，具体的话其实就相当于变成了\(\sum\limits_{k=j}^i\binom i k+\binom i {k-1}\)，也就是\(\sum\limits_{k=j}^{i}\binom {i+1} k\)，然后再把后面缺的东西补上就好了
也就是说现在我们可以得到\(\sum\limits_{j=k}^{6k-12}\binom{6k-12}{j}\)，那么对于上面的\(j\)从\(k-a\)开始枚举就直接加一下就好了（反正。。\(a\)很小）
这样我们就可以预处理出\(d\)数组了
最后是答案的求解：记\(g(n,k)\)表示在前\(n\)个四面体中保留\(k\)个，那么有：
然后最终的答案就是：
那现在就是怎么求\(g(n,k)\)了，这里有一个比较套路的处理方式：我们可以构造这样一个多项式\(P_n(x)\)：
会发现其实\(P_N(x)\)中的\(x^i\)的系数就是\(g(N,i)\)，那所以我们分治FFT一下把这个东西求出来就好了
mark：为了保证精度。。预处理一下单位根（虽然说好像是常规操作qwq不过平时写FFT没有这个习惯。。导致误差很大qwq所以。。以后这个习惯要改一下）
以及组合数的话因为带进去的数可能会大于模数所以需要使用Lucas定理
代码大概长这个样子
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#define ll long long
#define vct vector<ll>
using namespace std;
const int N=6*(1e4)+10,MOD=1e5+3;
const double pi=acos(-1);
int mul(int x,int y){return 1LL*x*y%MOD;}
int plu(int x,int y){return (1LL*x+y)%MOD;}
int ksm(int x,int y){
int ret=1,base=x;
for (;y;y>>=1,base=mul(base,base))
if (y&1) ret=mul(ret,base);
return ret;
}
struct cmplx{/*{{{*/
double a,b;
cmplx(){}
cmplx(double x,double y){a=x; b=y;}
void init(){a=0;b=0;}
friend cmplx operator + (cmplx x,cmplx y)
{return cmplx(x.a+y.a,x.b+y.b);}
friend cmplx operator - (cmplx x,cmplx y)
{return cmplx(x.a-y.a,x.b-y.b);}
friend cmplx operator * (cmplx x,cmplx y)
{return cmplx(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);}
};/*}}}*/
namespace FFT{/*{{{*/
const int N=::N*4;
cmplx A[N],B[N],w[20][N];
int rev[N],vis[20];
int len;
void get_len(int n){
for (int i=0;i<len;++i) A[i].init(),B[i].init();
int bit=0;
for (len=1;len<=n;++bit,len<<=1);
rev[0]=0;
for (int i=1;i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
for (int i=2,k=0;i<=len;i<<=1,++k){
if (vis[k]) continue;
for (int j=0;j<(i>>1);++j)
w[k][j]=cmplx(cos(j*pi/(i/2)),sin(j*pi/(i/2)));
vis[k]=1;
}
}
void fft(cmplx *a,int op){
cmplx W,u,v;
for (int i=0;i<len;++i)
if (rev[i]>i) swap(a[rev[i]],a[i]);
for (int step=2,k=0;step<=len;step<<=1,++k){
for (int st=0;st<len;st+=step){
for (int i=0;i<(step>>1);++i){
W=w[k][i]; W.b*=op;
v=a[st+i+(step>>1)]*W;
u=a[st+i];
a[st+i]=u+v;
a[st+i+(step>>1)]=u-v;
}
}
}
if (op==-1)
for (int i=0;i<len;++i) a[i].a/=len;
}
void work(){
fft(A,1);
fft(B,1);
for (int i=0;i<len;++i) A[i]=A[i]*B[i];
fft(A,-1);
}
void calc(vct &a,vct &b){
int n=a.size(),m=b.size();
get_len(n+m);
for (int i=0;i<n;++i) A[i]=cmplx(a[i],0);
for (int i=0;i<m;++i) B[i]=cmplx(b[i],0);
work();
}
}/*}}}*/
int bad[N],d[N],s[N],fac[MOD+10],invfac[MOD+10];
vct rec[N*20];
int n,m,T,K,tot;
ll C(int n,int m){
if (n<0||m<0||n<m) return 0;
if (n<MOD&&m<MOD) return mul(fac[n],mul(invfac[m],invfac[n-m]));
return mul(C(n/MOD,m/MOD),C(n%MOD,m%MOD));
}
void prework(int n){
fac[0]=1;
for (int i=1;i<MOD;++i) fac[i]=mul(fac[i-1],i);
invfac[MOD-1]=ksm(fac[MOD-1],MOD-2);
for (int i=MOD-2;i>=0;--i) invfac[i]=mul(invfac[i+1],i+1);
int tmp;
s[3]=0;
for (int j=3;j<=6*3-12;++j)
s[3]=plu(s[3],C(6*3-12,j));
for (int i=4;i<=n;++i){
tmp=s[i-1];
for (int j=6*(i-1)-12;j<6*i-12;++j) {
tmp=plu(mul(tmp,2),C(j,i-2));
}
tmp=plu(tmp,MOD-C(6*i-12,i-1));
s[i]=tmp;
}
d[1]=9;
for (int i=2;i<=n;++i){
d[i]=0;
tmp=plu(s[i],MOD-C(6*i-12,i));
for (int a=0;a<=4;++a){
tmp=plu(tmp,C(6*i-12,i-a));
d[i]=plu(d[i],mul(tmp,mul(C(4,a),ksm(3,a))));
}
}
//for (int i=1;i<=10;++i) printf("%d ",d[i]); printf("\n");
}
ll get_val(double x){return ((ll)round(x))%MOD;}
void print(int x){
for (int i=0;i<rec[x].size();++i) printf("%I64d ",rec[x][i]);
printf("\n");
}
int solve(int l,int r){
int mid=l+r>>1,lc,rc,now,len=r-l+1;
now=++tot;
if (l==r){
rec[now].resize(2);
rec[now][0]=1;
rec[now][1]=d[l];
/*printf("#%d: \n",l);
print(now);*/
return now;
}
lc=solve(l,mid);
rc=solve(mid+1,r);
FFT::calc(rec[lc],rec[rc]);
rec[now].resize(len+1);
for (int i=0;i<=len;++i) rec[now][i]=get_val(FFT::A[i].a);
//printf("#%d %d: \n",l,r);
//print(now);
return now;
}
void Solve(int n,int K){
tot=0;
int id=solve(1,n);
int ans=0;
for (int i=K;i<=n;++i) ans=plu(ans,rec[id][i]%MOD);
printf("%d\n",ans);
}
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
scanf("%d",&T);
prework(60000);
/*FFT::get_len(5);
FFT::A[0]=cmplx(1,0); FFT::A[1]=cmplx(1,0); FFT::B[0]=cmplx(1,0); FFT::B[1]=cmplx(1,0);
FFT::work();*/
//for (int i=1;i<=600000;++i) printf("%d\n",s[i]);
for (int o=1;o<=T;++o){
scanf("%d%d",&n,&K);
Solve(n,K);
}
}
