[BZOJ2629]binomial (高精度+Lucas定理+离散对数+FFT)

题意：对于给定的n和p，求对于所有的0<=i<p,满足C(n,k)%p=i的k的个数。

注：p虽然要输入，但是题目标注了所以测试点的p是固定的。

首先需要用正确的姿势理解lucas定理，比如求C(n,r)%p，就是将n和r分别转换为p进制，然后依次算组合数乘起来。n是一个高精度数，求C(n,r)的过程中，n不断模p得到的数（即n的p进制表示）是固定的。

就是长这样：n0!/((n0-k0)!*k0!) * n1!/((n1-k1)!*k1!)。。。其中n0,n1...都是定值，k0,k1...都是枚举的值，最终的k是由这些k0,k1...决定。这样容易让人想到生成函数。由于是在模意义下做，每个数都在0到p-1之间。而分母是乘起来而不是加起来，我们可以将每个数取离散对数，将乘法变成加法，但是注意是不能对0取对数的，好在0我们可以特殊处理，即Lucas定理的乘起来的若干个组合数中如果有ki>ni则为0。

然后对分母FFT一下就可以了。注意答案是模29，基本操作是在模51061下做，两个模数不要搞混了，答案模的29比较小可以直接用double来FFT，卷积起来不会溢出的。。（因为如果模998244353就太明显了）

#include<iostream>#include<algorithm>#include<cstdio>#include<cstring>#include<cmath>#include<vector>#define rep(i,a,b) for(int i=a;i<=b;++i)#define erp(i,a,b) for(int i=a;i>=b;--i)#define LL long long#define DB doubleusing namespace std;const DB PI = acos(-1);const int mo = 51061; // primitive root of which is 7.const int mn = 65540;const int lim = 100000000ll;struct bign{LL a[50];int len;bign() { memset(a, 0, sizeof a); len = 1; } void operator /= (LL b){erp(i, len, 2) a[i-1] += a[i]%b*lim, a[i] /= b;for (a[1]/=b; len>1&&!a[len]; --len); }bign operator / (LL b){bign t = *this; t /= b; return t;}LL operator % (LL b){LL res = 0, tl = lim%b;erp(i, len, 1) res = (res*tl%b + a[i]%b) % b;return res;}bool zero() { return len==1&&!a[1]; }    void readstr(const char*s){len = strlen(s);for (int i = 0, j; j=(len-i-1)/8+1, i<len; ++i)a[j] = a[j]*10 + s[i]-'0';len = (len+7) / 8;}};int jc[mn], rjc[mn], lnd[mn], pw[mn];void math_table(){jc[0] = jc[1] = rjc[0] = rjc[1] = 1;rep(i, 2, mo) jc[i]=1ll*jc[i-1]*i%mo, rjc[i]=1ll*(mo-mo/i)*rjc[mo%i]%mo;rep(i, 2, mo) rjc[i] = 1ll*rjc[i-1]*rjc[i]%mo;pw[0]=1; rep(i,1,mo-1) pw[i] = 7ll*pw[i-1]%mo;rep(i,0,mo-1) lnd[pw[i]]=i;}int fenmu(int n, int k){return lnd[1ll*rjc[k]*rjc[n-k]%mo];}struct cpx{DB r, i;cpx() {}cpx(DB a, DB b):r(a), i(b){};inline cpx operator + (cpx b) { return cpx(r+b.r, i+b.i); }inline cpx operator - (cpx b) { return cpx(r-b.r, i-b.i); }inline cpx operator * (cpx b) { return cpx(r*b.r - i*b.i, r*b.i+i*b.r); }};namespace FFT{int r[mn*2];cpx a[mn*2], b[mn*2];void fft(cpx*a, int flag, int N){rep(i, 0, N-1) if (i<r[i]) swap(a[i], a[r[i]]);for (int i = 1; i<N; i<<=1){cpx wn(cos(PI/i), flag*sin(PI/i));for (int j = 0; j<N; j+=i<<1){cpx w(1, 0);for (int k = 0; k<i; ++k, w=w*wn){cpx x = a[j+k], y = w*a[j+k+i];a[j+k] = x+y, a[j+k+i] = x-y;}}}if (flag<0) rep(i, 0, N-1) a[i].r /= N;}void polymul(int*p1, int*p2, int m) //卷积后保存在p1中{rep(i, 0, mo-1) a[i]=cpx(p1[i], 0), b[i]=cpx(p2[i], 0);rep(i, mo, m) a[i] = b[i] = cpx(0, 0);fft(a, 1, m); fft(b, 1, m);rep(i, 0, m-1) a[i] = a[i]*b[i];fft(a, -1, m);rep(i, 0, m-1) p1[i] = a[i].r+0.5;}}int res[mn*2], a[mn*2], ans[mn];void calc(vector<int> num, int tot){int m = 131072, tmp = 1;sort(num.begin(), num.end());for (int i = 0; i<(int)num.size(); ++i)tmp = tmp*(num[i]+1)%29;ans[0] = (tot+29-tmp)%29;rep(j, 0, num[0])(++res[fenmu(num[0], j)]) %= 29;rep(i, 0, m) FFT::r[i] = (FFT::r[i>>1]>>1)|((i&1)<<(17-1));for (int i = 1; i<(int)num.size(); ++i){memset(a, 0, sizeof a);rep(j, 0, num[i]) (++a[fenmu(num[i], j)]) %= 29;FFT::polymul(res, a, m);for (int j = mo-1; j<m; res[j++]=0)res[j%(mo-1)] += res[j]%29;for (int j = 0; j<mo-1; res[j++]%=29);}int fz = 1;for (int i = 0; i<(int)num.size(); ++i)fz = 1ll*fz*jc[num[i]]%mo;fz = lnd[fz];for (int i = 0; i<mo-1; ++i)ans[pw[(i+fz)%(mo-1)]] = res[i];}char mystr[2333];char to29(int x) { x%=29; return x<10?'0'+x:'A'+x-10; }int main(){scanf("%s%*d", mystr);math_table();vector<int> num;bign tmp; tmp.readstr(mystr);int tot = tmp%29 + 1;while (!tmp.zero()) num.push_back(tmp%mo), tmp /= mo;calc(num, tot);for (int i = 0; i<mo; ++i) putchar(to29(ans[i]));puts("");return 0;}