几何

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)个木棒的方案数为：

[inom 4 acdot3^a ]

​　　然后是棱（指的是。。不与顶点相连的木棒），对应的这里取走的数量就应该是至少(k-a)个：去掉之前的与顶点相连的木棒，不与顶点相连的木棒总共有(6k-12)个，为了方便，我们定义一个(s(6k-12,k-a))表示从(6*k-12)个木棒中取至少(k-a)个木棒的方案数，那么有：

[s(6k-12,k-a)=sumlimits_{j=k-a}^{6k-12}inom {6k-12}{j} ]

​　　所以我们可以得到：

[d[k]=sumlimits_{a=0}^4inom 4 acdot 3^asumlimits_{j=k-a}^{6k-12}inom{6k-12}j ]

​　　现在的问题就是怎么快速求后面的那个玩意，也就是(s)

​ 　　

​　　注意到这个其实是一个。。杨辉三角上面一行的一个后缀和，而根据组合数的递推式子(inom i j=inom {i-1}{j-1}+inom {i-1}{j})，如果说我们知道(sum_i=sumlimits_{k=j}^{i}inom i k)，那么(sum_{i+1}=sumlimits_{k=j}^{i+1}inom {i+1} k)是可以直接得到的，乘个(2)然后再加上(inom i {j-1})就好了，具体的话其实就相当于变成了(sumlimits_{k=j}^iinom i k+inom i {k-1})，也就是(sumlimits_{k=j}^{i}inom {i+1} k)，然后再把后面缺的东西补上就好了

​　　也就是说现在我们可以得到(sumlimits_{j=k}^{6k-12}inom{6k-12}{j})，那么对于上面的(j)从(k-a)开始枚举就直接加一下就好了（反正。。(a)很小）

​　　这样我们就可以预处理出(d)数组了

　　

​　　最后是答案的求解：记(g(n,k))表示在前(n)个四面体中保留(k)个，那么有：

[g(n,k)=g(n-1,k-1)+g(n-1,k)*d[n] ]

​　　然后最终的答案就是：

[Ans=sumlimits_{i=K}^Ng(N,i) ]

​　　那现在就是怎么求(g(n,k))了，这里有一个比较套路的处理方式：我们可以构造这样一个多项式(P_n(x))：

[P_n(x)=(x+d[n])P_{n-1}(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("
");
}
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("
");
}
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: 
",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: 
",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
",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
",s[i]);
	for (int o=1;o<=T;++o){
		scanf("%d%d",&n,&K);
		Solve(n,K);
	}
}