Description：

给出 $m$ 和质数 $p$ ，求 $C_m^0,C_m^1,C_m^2,\dots,C_m^m$ 中有多少是 $p^0$ 的倍数不是 $p^1$ 的倍数，是 $p^1$ 的倍数不是 $p^2$ 的倍数 $\dots\dots$

$1\le t \le 5000,2\le m,p\le 10^{18}$

Solution：

库默尔定理 $Kummer\ Theorem $ ：组合数 $C^m_{n+m}$ 中所含的质数 $p$ 的幂次为：

$\begin{align}\sum_{i=1}^\infty(\lfloor\frac{n+m}{p^i}\rfloor-\lfloor\frac n {p^i}\rfloor-\lfloor\frac m{p^i}\rfloor)\end{align}$

如果 $m+n$ 在 $p^i$ 上进位，那么 $\lfloor\frac{n+m}{p^i}\rfloor-\lfloor\frac n {p^i}\rfloor-\lfloor\frac m{p^i}\rfloor=1$ ，否则 $\lfloor\frac{n+m}{p^i}\rfloor-\lfloor\frac n {p^i}\rfloor-\lfloor\frac m{p^i}\rfloor=0$ ，也就是说组合数 $C^m_{n+m}$ 中所含 $p$ 的幂次为 $n+m$ 在 $p$ 进制下的进位次数。

$C_m^k$ 中所含 $p$ 的幂次为 $\begin{align}\sum_{i=1}^\infty(\lfloor\frac{m}{p^i}\rfloor-\lfloor\frac k {p^i}\rfloor-\lfloor\frac {m-k}{p^i}\rfloor)\end{align}$ ，当时 $\lfloor\frac{m}{p^i}\rfloor-\lfloor\frac k {p^i}\rfloor-\lfloor\frac {m-k}{p^i}\rfloor\ne0$ 会有 $1$ 的贡献。

再考虑什么时候 $\lfloor\frac{m}{p^i}\rfloor-\lfloor\frac k {p^i}\rfloor-\lfloor\frac {m-k}{p^i}\rfloor\ne1$ ，发现只能是原来一个 $p^i$ 被分到了后两部份中，也就是说 $m$ 中 $p^i+k(k<p^i)$ 被分成了两份，这两份都不大于 $p^i$ ，那么一定有 $k\%p^i>m\%p^i$ 。问题就是有多少数贡献为 $0,1,\dots,m$ ，于是就可以 $DP$ 了，设 $f[i][j][0/1]$ 表示已经处理到 $m$ 在 $p$ 进制下从低到高的第 $i$ 位，产生了 $j$ 的贡献，当前后缀比 $m$ 小 $/$ 大，如果当前位为 $x$ ，状态为 $f[i][j][0]$ ，后缀相当于取模，如果后缀较小，那么也就相当于取模小于 $m$ ，所以不会产生贡献，那么如果这一位的值 $b\in[0,x]$ ，那么后面 $0$ 的状态都可以转移，如果 $b\in[0,x)$ ，那么后面 $1$ 的状态可以转移。如果状态为 $f[i][j][1]$ 那么在这一位上会产生一的贡献，所以应从 $j-1$ 里转移，如果 $b\in[x,p-1]$ ，那么后面一的可以转移，如果 $b\in[x+1,p-1]$ ，那么后面 $0$ 的可以转移到这个状态。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

typedef long long ll;

ll m,p;

#define MAXL 100

ll bit[MAXL],tot = 0;

ll f[MAXL][MAXL][2];

void work()

{

    scanf("%lld%lld",&m,&p);

    tot = 0;

    while(m > 0)

    {

        bit[++tot] = m % p;

        m = m / p;

    }

    memset(f,0,sizeof(f));

    f[1][0][0] = bit[1] + 1;f[1][1][1] = p - bit[1] - 1;

    for(int i = 2;i <= tot;++i)

    {

        int x = bit[i];

        for(int j = 0;j <= i;++j)

        {

            f[i][j][0] += (x + 1) * f[i - 1][j][0] + x * f[i - 1][j][1];

            f[i][j + 1][1] += (p - x) * f[i - 1][j][1] + (p - x - 1) * f[i - 1][j][0];

        }

    }

    for(int i = 0;f[tot][i][0] != 0;++i)printf("%lld ",f[tot][i][0]);

    puts("");

    return;

}

int main()

{

    int testcases;

    scanf("%d",&testcases);

    while(testcases--)work();

    return 0;

}