根据Lucas定理,等价于在$P$进制下每一位分别求组合数最后乘积模$P$。
因为答案为$0$的并不好算,所以可以考虑用$n+1$减去其它所有的答案。
那么每一位的组合数都不能是$0$,那么这就保证了$k$的每一位都不大于$n$,所以无需考虑$k\leq n$这个限制。
求出模$P$下每个数的指标,考虑数位DP,设$f[i][j]$表示考虑了最后$i$位,组合数的指标之和模$\varphi(P)$为$j$的方案数。
那么转移就是卷积的形式,FFT加速即可。
时间复杂度$O(P\log P\log n)$。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int P=51061,F=P-1,O=29,N=131100;
const double pi=acos(-1.0);
char s[P];
int flag,n=131072,len,i,j,po[P],ind[P],fac[P],inv[P],a[P],ans[P],pos[N],f[N],g[N];
struct comp{
double r,i;comp(double _r=0,double _i=0){r=_r,i=_i;}
comp operator+(const comp&x){return comp(r+x.r,i+x.i);}
comp operator-(const comp&x){return comp(r-x.r,i-x.i);}
comp operator*(const comp&x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
comp conj(){return comp(r,-i);}
}A[N],B[N];
inline int C(int n,int m){return 1LL*fac[n]*inv[m]*inv[n-m]%P;}
void FFT(comp*a,int t){
for(int i=1;i<n;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
for(int d=0;(1<<d)<n;d++){
int m=1<<d,m2=m<<1;double o=pi*2/m2*t;comp _w(cos(o),sin(o));
for(int i=0;i<n;i+=m2){
comp w(1,0);
for(int j=0;j<m;j++){
comp&A=a[i+j+m],&B=a[i+j],t=w*A;
A=B-t;B=B+t;w=w*_w;
}
}
}
if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
}
void mul(int k){
if(!flag){
flag=1;
for(i=0;i<=k;i++)f[ind[C(k,i)]]++;
for(i=0;i<F;i++)f[i]%=O;
return;
}
for(i=0;i<F;i++)g[i]=0;
for(i=0;i<=k;i++)g[ind[C(k,i)]]++;
for(i=0;i<F;i++)g[i]%=O;
for(i=0;i<F;i++)A[i]=comp(f[i],g[i]);
for(i=F;i<n;i++)A[i]=comp(0,0);
FFT(A,1);
for(i=0;i<n;i++){
j=(n-i)&(n-1);
B[i]=(A[i]*A[i]-(A[j]*A[j]).conj())*comp(0,-0.25);
}
FFT(B,-1);
for(i=0;i<F;i++)f[i]=0;
for(i=0;i<n;i++)f[i%F]=(f[i%F]+int(B[i].r+0.5))%O;
}
int main(){
for(po[0]=i=1;i<P-1;i++)po[i]=po[i-1]*2%P;
for(i=0;i<P-1;i++)ind[po[i]]=i;
for(fac[0]=fac[1]=inv[0]=inv[1]=1,i=2;i<P;i++)inv[i]=1LL*(P-inv[P%i])*(P/i)%P;
for(i=2;i<P;i++)fac[i]=1LL*fac[i-1]*i%P,inv[i]=1LL*inv[i-1]*inv[i]%P;
scanf("%s",s+1);
len=strlen(s+1);
for(i=1;i<=len;i++)a[i]=s[len-i+1]-'0',ans[0]=(ans[0]*10+s[i]-'0')%O;
ans[0]++;
ans[0]%=O;
j=__builtin_ctz(n)-1;
for(i=0;i<n;i++)pos[i]=pos[i>>1]>>1|((i&1)<<j);
while(len){
for(a[0]=0,i=len;i;i--)a[i-1]+=a[i]%P*10,a[i]/=P;
while(len&&!a[len])len--;
mul(a[0]/10);
}
for(i=0;i<F;i++)ans[po[i]]=f[i];
for(i=1;i<P;i++)ans[0]=(ans[0]+O-ans[i])%O;
for(i=0;i<P;i++)if(ans[i]<=9)putchar(ans[i]+'0');else putchar(ans[i]-10+'A');
return 0;
}