bzoj 2629: binomial （FFT+DP+Lucas定理+短除法）

题目描述

传送门

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

　　注：C(n,k)=n!k!∗(n−k)! n<=p10,p=51061

题解

这道题因为n很大所以不能直接计算组合数。

但是我们可以利用Lucas定理进行变形，C(n,k)=C(n%p,k%p)∗C(n/p,k/p) 不断的除p，就可以得到一串C相乘的形式。分解成一串n%p其实就相当于是n的p进制数。因为C(n,k)在k>n的时候是等于0的，所以我们对于每一位来说只需要计算到p进制下这一位的值即可。这样组合出来的p进制数k一定是<=n的，对于剩余的情况都是余数为0的情况，最后用总数-∑p−1i=1ans[i]即可得到ans[0]

另f[i][j]表示计算到p进制下的第i位，余数为j的方案数。

对于每个i，预处理b[j] 表示C(ni,mi) 中余数为j的个数

f[i+1][j∗k%p]+=f[i][j]∗b[k]

如果我们求出p的原根，那么就可以把j*k变成num[j]+num[k]的形式，那么上面的式子就是一个卷积的形式，就可以用FFT来优化

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 150000
#define LL long long
#define pi acos(-1)
using namespace std;
char ch[40]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S'};
struct data{
double x,y;
data(double X=0,double Y=0) {
x=X,y=Y;
}
}b
,c
;
data operator -(data a,data b) {
return data(a.x-b.x,a.y-b.y);
}
data operator +(data a,data b){
return data(a.x+b.x,a.y+b.y);
}
data operator *(data a,data b){
return data(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);
}
int n,m,up[200],p,t[2000],k,l,a[2000],md
,ans
,f[20]
,R
,g,num
,L,mp
;
LL jc
,inv
;
char s[2000];
void calc()
{
for (l=0;k;) {
for (int i=k-1;i>0;--i) {
t[i-1]+=t[i]%p*10;
t[i]/=p;
}
a[++l]=t[0]%p;
t[0]/=p;
for (;k>0&&!t[k-1];--k);
}
}
int quickpow(LL num,int x,LL p)
{
LL ans=1; LL base=num%p;
while (x) {
if (x&1) ans=ans*base%p;
x>>=1;
base=base*base%p;
}
return ans;
}
LL C(int n,int m)
{
return (LL)jc
*inv[m]%p*inv[n-m]%p;
}
int get_g(int x)
{
if (x==2) return 1;
for (int i=2;i;i++) {
bool pd=1;
for (int j=2;j*j<x;j++)
if (quickpow(i,(x-1)/j,x)==1) {
pd=false;
break;
}
if (pd) return i;
}
}
void FFT(data a
,int opt)
{
for (int i=0;i<n;i++)
if (i<R[i]) swap(a[R[i]],a[i]);
for (int i=1;i<n;i<<=1) {
double s=(double)pi/i;
data wn=data(cos(s),opt*sin(s));
for (int p=i<<1,j=0;j<n;j+=p) {
data w=data(1,0);
for (int k=0;k<i;k++,w=w*wn) {
data x=a[j+k],y=w*a[j+k+i];
a[j+k]=x+y; a[j+k+i]=x-y;
}
}
}
}
int main()
{
freopen("tyrrell.in","r",stdin);
freopen("tyrrell.out","w",stdout);
scanf("%s",s);
k=strlen(s);
for (int i=k-1;i>=0;i--) t[k-i-1]=s[i]-'0';
scanf("%d",&p);
int cnt=0;
for (int i=k-1;i>=0;i--) {
cnt=cnt*10+t[i];
cnt%=29;
}
cnt=(cnt+1)%29;
calc();
g=get_g(p);
for (int i=0,j=1;i<p-1;i++,j=(j*g)%p)
num[j]=i,mp[i]=j;
m=2*p;
for (n=1;n<=m;n<<=1) L++;
for (int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
jc[0]=1;
for (int i=1;i<=p;i++) jc[i]=(jc[i-1]*i)%p;
for (int i=0;i<=p;i++) inv[i]=quickpow(jc[i],p-2,p)%p;
f[0][0]=1;
for (int i=1;i<=l;i++) {
memset(md,0,sizeof(md));
for (int j=0;j<=a[i];j++)
md[num[C(a[i],j)%p]]++;
for (int j=0;j<p;j++) md[j]%=29;
for (int j=0;j<p;j++) b[j].x=md[j],b[j].y=0;
for (int j=p;j<=n;j++) b[j].x=0,b[j].y=0;
FFT(b,1);
for (int j=0;j<p;j++) c[j].x=f[i-1][j],c[j].y=0;
for (int j=p;j<=n;j++) c[j].x=0,c[j].y=0;
FFT(c,1);
for (int j=0;j<=n;j++) b[j]=b[j]*c[j];
FFT(b,-1);
for (int j=0;j<n;j++) f[i][j%(p-1)]+=(int)(b[j].x/n+0.5);
for (int j=0;j<p;j++) f[i][j]%=29;
}
int sum=0;
for (int i=0;i<p;i++) ans[mp[i]]=f[l][i]%29,sum+=ans[mp[i]];
sum%=29;
ans[0]=(cnt-sum+29)%29;
for (int i=0;i<p;i++) putchar(ch[ans[i]%29]);
printf("\n");
}