题目大意
定义一个$S-$四面体表示六条边由$S$根不同的木棍组成，定义一种染色方法合法当且仅当至少有$S$根木棍被染色且与每个顶点相邻的三根木棍中至多有一根被染色，求有$N$个$S=1,2...N$四面体，求至少染$K$个的方案数。
题解
单独考虑$S=1$四面体，染它有$9$中方案，否则将与每个顶点相邻的$12$条边单独拿出来计算。
考虑其余四面体被染色的方案数：他们中有$k$个顶点$(k=0,1,2,3,4)$相邻的木棍三条中有一条被染色，方案是$C_4^k imes 3^k$，剩余的边则有$sumlimits_{i=S-k}^{6S-12}C_{6S-12}^{i}$，用乘法原理乘一下即可。
由于$k$很小，可以只求$sumlimits_{i=S}^{6S-12}C_{6S-12}^{i}$，其余加上组合数即可。而这个式子是杨辉三角上的某一行的一个后缀，所以是可以递推的，考虑上一行答案对下一行答案的贡献，只需要乘二再加上上一行后缀左侧的那个组合数即可。
我们现在已经解决了每一个$S$的方案数，接下来就是求染出至少$K$个的方案数，也只需要递推，设$G_x$表示染$x$个的方案数$V$表示当染前四边形的方案数，$G_x=G_x+V imes G_{x-1}$。这是$N$个简单最高次项次数是一次的多项式的卷积，用分治$FFT$解决即可。
有两个细节，由于模数只有$10^5+3$，所以算组合数要用$lucas$定理。并且，$FFT$中每一位可能达到$10^{15}$级别，精度可能会爆炸，所以要优化或者使用$longspace double$。
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define LL long long
#define LD long double
#define mod 100003
#define M 100020
using namespace std;
int read(){
int nm=0,fh=1; int cw=getchar();
for(;!isdigit(cw);cw=getchar()) if(cw=='-') fh=-fh;
for(;isdigit(cw);cw=getchar()) nm=nm*10+(cw-'0');
return nm*fh;
}
const int P[]={1,12,54,108,81};
const LD PI=3.141592653589793238462643383;
int mul(int x,int y){return (LL)x*(LL)y%mod;}
int add(int x,int y){return (x+y)>=mod?x+y-mod:x+y;}
int mus(int x,int y){return (x-y)<0?x-y+mod:x-y;}
void upd(int &x,int y){x=add(x,y);}
int qpow(int x,int sq){
int res=1;
for(;sq;sq>>=1,x=mul(x,x)) if(sq&1) res=mul(res,x);
return res;
}
int n,m,fac[mod],ifac[mod],G[M<<2],lg[M<<2];
int C(int tot,int tk){if(tot<0||tk<0||tot<tk)return 0;return mul(fac[tot],mul(ifac[tot-tk],ifac[tk]));}
int lucas(int tot,int tk){if(!tk) return 1;return mul(C(tot%mod,tk%mod),lucas(tot/mod,tk/mod));}
int rev[M<<2],F[M<<2],S[M];
struct comp{
LD r,d; comp(){r=d=0.0;}
comp(LD _r,LD _d){r=_r,d=_d;}
comp operator +(const comp&k)const{return comp(r+k.r,d+k.d);}
comp operator -(const comp&k)const{return comp(r-k.r,d-k.d);}
comp operator *(const comp&k)const{return comp(r*k.r-d*k.d,r*k.d+d*k.r);}
}A[M<<2],B[M<<2];
void FFT(comp *x,int len,LD kd){
for(int i=1;i<len;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
for(int tt=1;tt<len;tt<<=1){
comp unit(cos(PI*kd/(tt*1.0)),sin(PI*kd/(tt*1.0)));
for(int st=0;st<len;st+=(tt<<1)){
comp now(1.0,0.0);
for(int pos=st;pos<st+tt;pos++,now=now*unit){
comp t1=x[pos],t2=x[pos+tt]*now;
x[pos]=t1+t2,x[pos+tt]=t1-t2;
}
}
}
if(kd<0.0){for(int i=0;i<len;i++) x[i].r/=(len*1.0);}
}
void tms(int *x,int *x1,int *x2,int n1,int n2){
if(n1+n2<=120){
memset(S,0,sizeof(int)*(n1+n2+1));
for(int i=0;i<=n1;i++){
for(int j=0;j<=n2;j++) S[i+j]+=mul(x1[i],x2[j]);
}
for(int i=0;i<=n1+n2;i++) x[i]=S[i]%mod;
}
else{
int len=1,nw=-1;for(;len<=n1+n2+1;len<<=1,nw++);
for(int i=1;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<nw);
for(int i=0;i<=n1;i++) A[i]=comp(x1[i]*1.0,0.0);
for(int i=0;i<=n2;i++) B[i]=comp(x2[i]*1.0,0.0);
for(int i=n1+1;i<=len;i++) A[i]=comp(0.0,0.0);
for(int i=n2+1;i<=len;i++) B[i]=comp(0.0,0.0);
FFT(A,len,1.0),FFT(B,len,1.0);
for(int i=0;i<len;i++) A[i]=A[i]*B[i]; FFT(A,len,-1.0);
for(int i=0;i<=n1+n2;i++) x[i]=llround(A[i].r)%mod;
}
}
void solve(int *x,int L,int R){
if(L==R){x[0]=1,x[1]=G[L];return;}
int mid=((L+R)>>1),ls,rs; ls=mid-L+1,rs=R-mid;
solve(x,L,mid),solve(x+ls+1,mid+1,R);
tms(x,x,x+ls+1,ls,rs);
}
int main(){
fac[0]=ifac[0]=1,G[1]=9,G[2]=243,G[3]=16224,G[4]=46489;
for(int i=1;i<mod;i++) fac[i]=mul(fac[i-1],i),ifac[i]=qpow(fac[i],mod-2);
for(int i=5,K,pre,last=3797,rem;i<M;i++){
for(K=(i-2)*6,pre=K-6;pre<K;pre++) last=add(add(last,last),lucas(pre,i-2));
last=mus(last,lucas(K,i-1)),rem=last;
for(int k=0;k<=4;k++) upd(G[i],mul(P[k],rem)),upd(rem,lucas(K,i-k-1));
}
for(int T=read(),ans=0;T;--T,ans=0){
n=read(),m=read(),solve(F,1,n);
for(int i=n;i>=m;i--) upd(ans,F[i]);
printf("%d
",ans);
}
return 0;
}