HDU-6036 Division Game（ntt模板）

Posted 2020-11-25 artoriax

题面

Description

There are (k) piles of stones in a circle, numbered from (0) to (k - 1), where the number of the stones in each pile is (n) initially. You can do some round operations, where the initial round is numbered as the (1)-st round.

The operation of the (i)-th round is to modify the pile of stones numbered ((i - 1) mod k). In each round, you should remove from this pile some stones (at least one stone), satisfying that the number of stones in this pile before this operation is a multiple of the number of stones in this pile after operation, which means that you ought to remain at least one stone in this pile.

The game is ended if there exists at least one pile containing only one stone. Given two positive integers (n) and (k), your task is to calculate for each pile the number of the possible operation plans that it is the last operated pile before the game is ended.

The integer (n) may be very large, so the prime-factor decomposition of (n) will be given, in other words, if (n = prod_{i = 1}^{m}{p_i^{e_i}}), then the integers (m) and ((p_i, e_i)) ((1 leq i leq m)) will be given, but the integer (n) will not.

The answer may be very large, so you only need to give the value of the answer modulo (985661441).

Input

The input contains multiple test cases.

For each test case:

The first line contains two positive integers (m) and (k), satisfying that (1 leq m, k leq 10).

In next (m) lines, the (i)-th line contains two positive integers (p_i) and (e_i), satisfying that (2 leq p_i leq 10^9,) (e_i geq 1,) (sum_{i = 1}^{m}{e_i} leq 10^5).

It is guaranteed that (p_1, p_2, cdots, p_m) are distinct.

About (200) test cases in total, where no more than (5) cases satisfy (sum_{i = 1}^{m}{e_i} geq 10^4).

Output

For each test case, output " Case #(x): (y_0) (y_1) (cdots) (y_{k - 1})" in one line (without quotes), where (x) indicates the case number starting from (1) and (y_i) ((0 leq i < k)) denotes the number of the possible operation plans modulo (985661441) for the pile numbered (i) of corresponding case.

Sample Input

1 1
2 2
2 1
3 1
5 1
1 2
2 3
2 2
2 4
5 4

Sample Output

Case #1: 2
Case #2: 3
Case #3: 6 4
Case #4: 1499980 1281085

题意

给定k堆石子，每堆有n个，n很大，以(n = prod_{i = 1}^{m}{p_i^{e_i}})形式给出p和e

按照堆的编号顺序取石子，取之前的石子数必须是取之后的石子数的倍数，至少取1个石子，当某一堆石子个数为1，游戏结束。

对于每堆石子，问使它成为第一个石子数取为1的方案数

题解

参考博客:

https://www.cnblogs.com/flipped/p/7617961.html

https://blog.csdn.net/qq_36553623/article/details/78154095

设(f(x))是每堆石子x次取为1的方案数，那么x次不取为1的方案数就是(f(x+1))，同理，x-1次不取为1的方案数就是$f(x) $

那么我们停止在第i堆的方案数就是
[ f(x+1)^{i-1}*f(x)^{k-i+1} ]
考虑x的取值范围，设(w=sum{e_i})

则(1 le x le w)，最少1次，最多w次全部取完

但考虑到第i堆前面的，w次必然取完，所以，第i堆最多取到w-1次。

如果i=1，则可以取到w次，因为i之前没有其他石子堆了，这样就需要分类讨论了

然后就是如何计算f(x)

我们设g(x)为x次取完所有质因子的方案数；则g(x)可以用挡板法计算：

对于每个质因子，在x次内取完，即(e_i)分成x组，允许空组，即
[ inom{e_i+x-1}{x-1} ]
那么根据乘法原理，g(x)即
[ g(x)=prod_{i=1}^m inom{e_i+x-1}{x-1} ]
考虑g(x)和f(x)的关系

由于每次至少取1个石子，所以不能有一次什么数都不取，而g(x)显然可能出现这种错误的状况

我们根据容斥定理
[ f(x)=随意取的方案数-某一次什么都没有取的方案数 \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+某两次什么都没有取的方案数 \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-某三次什么都没有取的方案数 \\~~~~~~~~~~~~~~~~~~~~~~~~~+.... ]
某一次什么都没有取的，即从x中选出1次什么都不取，剩下的x-1次正常取，即
[ inom{x}{1}*g(x-1) ]
同理计算出某两次，某三次没有取的，合并即得到
[ f(x)=sum_{y=0}^x (-1)^{x-y}g(y)inom{x}{y} ]
这个可以化成卷积的形式，如下
[ frac {f(x)} {x!}=sum_{y=0}^x frac {(-1)^{x-y}}{(x-y)!} frac {g(y)} {y!} ]
左边是x-y，右边是y，相加为x，可以放到两个数组里卷积。

由于模数(985661441=235×2^{22}+1)，符合ntt形式，可以用ntt优化，其原根为3

代码

//不加inline TLE >5000ms 加inline 1994ms
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 985661441;
const int N = 3e5 + 50;
namespace comb {
    const int N = 2e5 + 50;
    ll fac[N] = {1, 1}, inv[N] = {1, 1}, f[N] = {1, 1};
    void init() {
        for (int i = 2; i < N; i++) {
            fac[i] = fac[i - 1] * i % mod;
            f[i] = (mod - mod / i) * f[mod % i] % mod;//线性递推逆元
            inv[i] = inv[i - 1] * f[i] % mod;//阶乘的逆元
        }
    }
    inline ll C(int a, int b) {
        if (b > a || b < 0) return 0;
        return fac[a] * inv[b] % mod * inv[a - b] % mod;
    }
};

inline ll qpow(ll a, ll b) {
    ll ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ans;
}

namespace ntt {
    ll w[N];
    int G = 3;
    inline void ntt(ll a[], int n, int op) {//ntt
        for (int i = 0, j = 0; i < n; i++) {
            if (i > j) swap(a[i], a[j]);
            for (int l = n >> 1; (j ^= l) < l; l >>= 1);
        }
        for (int i = 2; i <= n; i <<= 1) {
            for (int j = 0, m = i >> 1; j < n; j += i) {
                for (int k = 0; k < m; k++) {
                    ll b = w[n / i * k] * a[j + m + k] % mod;
                    a[j + m + k] = (a[j + k] - b + mod) % mod;
                    a[j + k] = (a[j + k] + b) % mod;
                }
            }
        }
        if (op == -1) {
            reverse(a + 1, a + n);
            ll invN = qpow(n, mod - 2);
            for (int i = 0; i < n; i++) a[i] = a[i] * invN % mod;
        }
    }
    void getroot() {//求原根
        static int prime[1000], cnt;
        int n = mod - 1;
        for (int i = 2; i * i <= n; i++) {
            if (n % i == 0) {
                prime[cnt++] = i;
                while (n % i == 0) n /= i;
            }
        }
        if (n > 1) prime[cnt++] = n;
        for (int i = 1; i < mod; i++) {
            if (qpow(i, mod - 1) == 1) {
                bool flag = true;
                for (int j = 0; j < cnt; j++) {
                    if (qpow(i, (mod - 1) / prime[j]) == 1) {
                        flag = false;
                        break;
                    }
                }
                if (flag) {
                    G = i;
                    break;
                }
            }
        }
    }
    inline void conv(int n, ll x[], ll y[]) {//卷积
        ll g = qpow(G, (mod - 1) / n);
        w[0] = 1;
        for (int i = 1; i < n; i++) w[i] = w[i - 1] * g % mod;
        ntt(x, n, 1); ntt(y, n, 1);
        for (int i = 0; i < n; i++) x[i] = x[i] * y[i] % mod;
        ntt(x, n, -1);
    }
};
ll p[20], e[20];
ll g[N];
ll x[N], y[N];
ll f[N];
int main() {
    comb::init();
    int cse = 0;
    int m, k;
    while (~scanf("%d%d", &m, &k)) {
        ll w = 0;
        for (int i = 1; i <= m; i++) {
            scanf("%lld%lld", &p[i], &e[i]);
            w += e[i];
        }
        for (int i = 0; i <= w; i++) g[i] = 1;
        for (int i = 1; i <= m; i++) {
            for (int j = 0; j <= w; j++) {
                g[j] = g[j] * comb::C(e[i] + j - 1, j - 1) % mod;//挡板法
            }
        }
        int n = 1;
        while (n < w + 1) n <<= 1;
        n <<= 1;
        for (int i = 0; i < n; i++) x[i] = 0, y[i] = 0;
        for (int i = 0; i <= w; i++) {
            x[i] = comb::inv[i];
            if (i & 1) x[i] = -x[i];
            y[i] = g[i] * comb::inv[i] % mod;
        }
        ntt::conv(n, x, y);
        for (int i = 0; i < n; i++) f[i] = 0;
        for (int i = 1; i <= w; i++) {
            f[i] = x[i] * comb::fac[i] % mod;
        }
        printf("Case #%d:", ++cse);
        for (int i = 1; i <= k; i++) {
            ll tot = 0;
            for (int j = 1; j <= w; j++) {//循环还是从1到w，但是由于f[w+1]为0，所以除了i=1的时候，都取不到f[w]
                tot = (tot + qpow(f[j + 1], i - 1) * qpow(f[j], k - i + 1) % mod) % mod;
            }
            printf(" %lld", tot);
        }
        puts("");
    }
    return 0;
}