[集训队作业2018]复读机——生成函数+单位根反演

题面

　　uoj#450

解析

　　本文中用$m$表示原题中的$k$

　　设第$i$个复读机复读$t_i$次，最后答案等于：$$sum_{sum_{i=1}^mt_i=n}frac{n!}{prod t_i!}prod[d|t_i]\ =n!sum_{sum_{i=1}^mt_i=n}prod frac{[d|t_i]}{t_i!}$$

　　构造生成函数$A(x)=sum_{i=0}^{infty}[d|i]frac{1}{i!}x^i$，那么答案等于$A^m$的$n$次项系数乘以$n!$

　　用单位根反演化简$A(x)$：$$egin{align*}A(x)&=sum_{i=0}^{infty}[d|i]frac{1}{i!}x^i\&=sum_{i=0}^{infty}frac{x^i}{i!}*frac{1}{d}sum_{j=0}^{d-1}w_d^{ij}\&=frac{1}{d}sum_{j=0}^{d-1}sum_{i=0}^{infty}frac{(w_d^j)^i}{i!}x^i\&=frac{1}{d}sum_{j=0}^{d-1}e^{w_d^jx}end{align*}$$

　　$$A^m=frac{1}{d^m}(sum_{j=0}^{d-1}e^{w_d^jx})^m$$

　　$d=1$时，答案为$m^n$

　　$d=2$时，$$egin{align*}A^m&=frac{1}{d^m}(sum_{j=0}^{d-1}e^{w_d^jx})^m\&=frac{1}{2^m}(e^{w_2^0x}+e^{w_2^1x})^m\&=frac{1}{2^m}sum_{i=0}^{m}inom{m}{i}e^{(w_2^0i+w_2^1(m-i))x}end{align*}$$

　　　　答案为：$$egin{align*}Ans&=n!*frac{1}{2^m}sum_{i=0}^minom{m}{i}frac{(w_2^0i+w_2^1(m-i))^n}{n!}\&=frac{1}{2^m}sum_{i=0}^minom{m}{i}(w_2^0i+w_2^1(m-i))^nend{align*}$$

　　$d=3$时，$$egin{align*}A^m&=frac{1}{d^m}(sum_{j=0}^{d-1}e^{w_d^jx})^m\&=frac{1}{3^m}(e^{w_3^0x}+e^{w_3^1x}+e^{w_3^2x})^m\&=frac{1}{3^m}sum_{i=0}^{m}sum_{j=0}^{m-i}inom{m}{i}inom{m-i}{j}e^{(w_3^0i+w_3^1j+w_3^2(m-i-j))x}end{align*}$$

　　　　答案为：$$egin{align*}Ans&=n!*frac{1}{3^m}sum_{i=0}^msum_{j=0}^{m-i}inom{m}{i}inom{m-i}{j}frac{(w_3^0i+w_3^1j+w_3^2(m-i-j))^n}{n!}\&=frac{1}{3^m}sum_{i=0}^msum_{j=0}^{m-i}inom{m}{i}inom{m-i}{j}(w_3^0i+w_3^1j+w_3^2(m-i-j))^nend{align*}$$

　　单位根$w_d^1=g^{frac{mod-1}{d}}$，$g$为$mod$的原根，本题中等于$7$

　　$O(M^{d-1}log N)$

　代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
const int mod = 19491001;

int add(int x, int y)
{
    return x + y < mod? x + y: x + y - mod;
}

int rdc(int x, int y)
{
    return x - y < 0? x - y + mod: x - y;
}

ll qpow(ll x, int y)
{
    ll ret = 1;
    while(y)
    {
        if(y&1)
            ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}

int n, m, d;

namespace d2{
    void work()
    {
        ll ans = 0, mul = 1;
        for(int i = 0; i <= m; ++i)
        {
            ans = add(ans, qpow(rdc(2 * i, m), n) * mul % mod);
            mul = (mul * (m - i) % mod) * qpow(i + 1, mod - 2) % mod;
        }
        printf("%lld", ans * qpow(qpow(2, m), mod - 2) % mod);
    }
}

namespace d3{
    ll c[1005][1005];

    void init()
    {
        for(int i = 0; i <= m; ++i)
        {
            c[i][0] = 1;
            for(int j = 1; j <= i; ++j)
                c[i][j] = add(c[i-1][j-1], c[i-1][j]);
        }
    }

    void work()
    {
        init();
        ll w1 = qpow(7, (mod - 1) / 3), w2 = w1 * w1 % mod, ans = 0;
        for(int i = 0; i <= m; ++i)
            for(int j = 0; j <= m - i; ++j)
                ans = add(ans, (c[m][i] * c[m-i][j] % mod) * qpow(add(i, add(w1 * j % mod, w2 * (m - i - j) % mod)), n) % mod);
        printf("%lld", ans * qpow(qpow(3, m), mod - 2) % mod);
    }
}

int main()
{
    scanf("%d%d%d", &n, &m, &d);
    if(d == 1)
    {
        printf("%lld", qpow(m, n));
        return 0;
    }
    if(d == 2)
        d2::work();
    else
        d3::work();
    return 0;
}