[JZOJ5330]密码

题目大意

给定n,p,k，求$\sum_{i<=j<=n} [ p^k|C_i^j ]$
其中，$n<=10^{1000}，p<=10^9且p为质数，k<=10^9$

分析

首先我们看一个组合数，怎么样才能被$p^k$整除呢？我们有$C_m^n= \dfrac{m!}{n!(m-n)!}$，我们可以分别算出上下阶乘因子含有多少个p，记为$cnt_p()$再相减。对于一个阶乘$n!我们有cnt_p(n)=\sum_{i=1}^{∞}[\dfrac{n}{p^i}]（下取整）$。那么如果cnt(m)-cnt(n)-cnt(m-n)>=k，该组合数能被$p^k$整除。
这样不能看出什么，我们把cnt都合起来，原式变成$\sum_{i=1}^{∞}[\dfrac{m}{p^i}]-[\dfrac{n}{p^i}]-[\dfrac{m-n}{p^i}]$，可以发现，对于一个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);
}