codeforces 582D - Number of Binominal Coefficients

emmmmmmmmmmmmmm
这题模拟赛的时候做过
当时back的时候推柿子推得很痛苦
现在再做还是不会….不过推柿子变得熟练了很多…

题解：
因为组合数可以写成阶乘的形式$C_n^m=\dfrac{n!}{m!(n-m)!}$
又因为$p$是质数，组合数能否被$p^k$整除取决于他含多少个$p$作为因子
所以我们不妨先考虑阶乘$n!$含多少个$p$
假设$n!$含有$k$个$p$作为因子，有$k=\sum_{i=1}^{∞}\lfloor \frac{n}{p^i}\rfloor$ （易证）
那么$C_n^m$含有的$p$的个数就是$\sum_{i=1}^{∞}\lfloor \frac{n}{p^i}\rfloor -\lfloor \frac{m}{p^i}\rfloor -\lfloor \frac{n-m}{p^i}\rfloor$
写在$p$进制下，容易发现对于每个i，那个计算式只会是0或者1，且取决于m和n-m加起来在第i位是否产生了进位
所以问$n,m<=A$时有多少个$C_n^m$被$p^k$整除，就是问有多少对数$x,y$，满足$x+y<=A$且在$p$进制下产生了至少k次进位

我们将$A$转成$p$进制，做数位dp
$dp[i][x][0/1][0/1]$表示dp到第i位，已经产生了x次进位，前面是否顶格，下一位是否会进位的方案数

code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 3500;
const ll Mod = 1e9+7;

int n;
ll p,m;
char str[maxn];
int t[maxn],t2[maxn],L,len,A[maxn],tp;

ll f[2][maxn][2][2];
ll cal(const ll x){ return x*(x+1)/2%Mod; }
ll dp()
{
    int now=0; f[now][0][1][0]=1;
    for(int j=tp;j>=1;j--)
    {
        now=!now;
        ll k=A[j];
        for(int i=0;i<tp;i++)
        {
            ll tmp=f[!now][i][0][0]%Mod; f[!now][i][0][0]=0;
            if(tmp)
            {
                f[now][i+1][0][1]+=cal(p-1)*tmp%Mod;
                f[now][i][0][0]+=cal(p)*tmp%Mod;
            }
            tmp=f[!now][i][0][1]%Mod; f[!now][i][0][1]=0;
            if(tmp)
            {
                f[now][i+1][0][1]+=cal(p)*tmp%Mod;
                f[now][i][0][0]+=cal(p-1)*tmp%Mod;
            }
            tmp=f[!now][i][1][0]%Mod; f[!now][i][1][0]=0;
            if(tmp)
            {
                f[now][i][1][0]+=(k+1)*tmp%Mod;
                f[now][i+1][1][1]+=k*tmp%Mod;
                f[now][i][0][0]+=cal(k)*tmp%Mod;
                f[now][i+1][0][1]+=cal(k-1)*tmp%Mod;
            }
            tmp=f[!now][i][1][1]%Mod; f[!now][i][1][1]=0;
            if(tmp)
            {
                f[now][i][1][0]+=(p-k-1)*tmp%Mod;
                f[now][i+1][1][1]+=(p-k)*tmp%Mod;
                f[now][i][0][0]+=(cal(p-1)-cal(p-k-1))*tmp%Mod;
                f[now][i+1][0][1]+=(cal(p)-cal(p-k))*tmp%Mod;
            }
        }
    }
    ll re=0;
    for(int x=m;x<tp;x++)
    {
        for(int i=0;i<2;i++)
            (re+=f[now][x][i][0]%Mod)%=Mod;
    }
    if(re<0) re+=Mod;
    return re;
}

ll divide()
{
    ll now=0;
    for(int i=L;i<=len;i++)
    {
        now=now*10+t[i];
        t2[i]=now/p; now%=p;
    }
    for(int i=1;i<=len;i++) t[i]=t2[i],t2[i]=0;
    while(!t[L]&&L<=len) L++;
    return now;
}

int main()
{
    scanf("%I64d%I64d",&p,&m);
    scanf("%s",str); len=strlen(str);
    for(int i=1;i<=len;i++) t[i]=str[i-1]-'0';
    tp=0; while(L<=len) A[++tp]=divide();

    if(tp<=m) puts("0");
    else printf("%I64d\n",dp());

    return 0;
}