[JZOJ5330]密码

题目大意

给定n,p,k，求∑i<=j<=n[pk|Cji]
其中，n<=101000，p<=109且p为质数，k<=109

分析

首先我们看一个组合数，怎么样才能被pk整除呢？我们有Cnm=m!n!(m−n)!，我们可以分别算出上下阶乘因子含有多少个p，记为cntp()再相减。对于一个阶乘n!我们有cntp(n)=∑∞i=1[npi]（下取整）。那么如果cnt(m)-cnt(n)-cnt(m-n)>=k，该组合数能被pk整除。
这样不能看出什么，我们把cnt都合起来，原式变成∑∞i=1[mpi]−[npi]−[m−npi]，可以发现，对于一个i，多项式的值只有可能是0或者1，1就意味着有因子p了嘛。有lucas定理，我们容易想到把n，m，m-n看成p进制数，那么多项式的值为1的时候，就意味着二进制下n,m−n的第i−1位发生了进位，导致m=n+(m−n)中，第i位m多了1。
这个其实是库默尔定理。
那么我们尝试构造n和m-n，问题就变成：有多少对(x,y)，满足十进制时x,y<=n，p进制时x+y有K次或更多次进位。
然后就可以用数位DP做了，P进制下的数位DP。考虑设计状态。
状态[i][j][0/1][0/1]分别表示从高位往低位考虑到第i位，一共有j位进位了，x+y是否贴着上界N，i+1位是否进位。f[i][j][0,1][0,1]就表示此时的方案数。然后算的时候，对于每个状态我们计算x和y在i+1位中的满足条件方案数，乘上当前的f，转移到新的f。转移十分恶心，我推了好久，系数是一堆等差数列。

代码

#include<cstdio>#include<algorithm>#include<cmath>#include<cstring>#include<bitset> using namespace std;typedef long long ll;typedef double db;#define fo(i,j,k) for(i=j;i<=k;i++)#define fd(i,j,k) for(i=j;i>=k;i--)const int N=1005,mo=1e9+7;ll f[N][N][2][2],a[N],b[N],a00,a01,a10,a11,ax,ay,bx,by,siz,o00,o10,i,p,K,n,j,ans;char s[N];void predo(){    while (a[0])    {        fd(i,a[0],2)        {            a[i-1]+=a[i]%p*10;            a[i]/=p;        }        b[++b[0]]=a[1]%p;        a[1]/=p;        while (a[0]&&a[a[0]]==0) a[0]--;    }}int main(){    freopen("t3.in","r",stdin);    //freopen("password.out","w",stdout);    scanf("%s %d %d",s+1,&p,&K);    a[0]=strlen(s+1);    fo(i,1,a[0]) a[a[0]-i+1]=s[i]-'0';    predo();    n=b[0];    fo(i,1,n/2) swap(b[i],b[n-i+1]);    f[1][0][0][0]=(b[1]+1)*b[1]/2%mo;    f[1][0][1][0]=b[1]+1;    f[1][1][0][1]=(b[1])*(b[1]-1)/2%mo;//等下重新推     f[1][1][1][1]=b[1];    o00=a00=(1+p)*p/2;    a01=p*p-a00;    o10=a10=p*(p-1)/2;    a11=p*p-a10;    a00%=mo;a01%=mo;a10%=mo;a11%=mo;    // 第i+1位还没有考虑方案，只考虑第i位。     fo(i,1,n)        fo(j,0,i)        {                //f[i+1][j+1][0][1], 0~p-1,                //考虑i+1位放置情况                 ax=(1+b[i+1])*b[i+1]/2;                siz=(p-b[i+1]-1)+1;                ay=(b[i+1]+1+p)*siz/2+(p-siz)*p-o00;                     ax%=mo;ay%=mo;                bx=b[i+1]*(b[i+1]-1)/2;                siz=(p-b[i+1])+1;                by=(b[i+1]+p)*siz/2+(p-siz)*p-o10;                bx%=mo;by%=mo;                (f[i+1][j][0][0]+=f[i][j][0][0]*a00+f[i][j][0][1]*a01+f[i][j][1][0]*ax+f[i][j][1][1]*ay)%=mo;                (f[i+1][j+1][0][1]+=f[i][j][0][0]*a10+f[i][j][0][1]*a11+f[i][j][1][0]*bx+f[i][j][1][1]*by)%=mo;                (f[i+1][j][1][0]+=f[i][j][1][0]*(b[i+1]+1)+f[i][j][1][1]*(p-1-(b[i+1]+1)+1))%=mo;                (f[i+1][j+1][1][1]+=f[i][j][1][0]*(b[i+1]-1+1)+f[i][j][1][1]*(p-1-b[i+1]+1))%=mo;        }    fo(j,K,n)    {        ans=(ans+f[n][j][0][0]+f[n][j][1][0])%mo;    }    printf("%lld\n",ans);}