CodeChef IMGOD Apun Hi Bhagwan Hai !

問題文

https://www.codechef.com/problems/IMGOD

問題概要

整数N (1 <= N <= 1e5) のみが与えられる.
分割する対象が区別できるものでの N の分割数 S(N, i) mod 163577857 (1 <= i <= N) をすべて求めよ.

解法

これは第二種スターリング数というものらしい.

S(n, k) = (1 / k!) * Σ(-1)^(k - m) * kCm * m^n = Σ ( (-1)^(k - m) / (k - m)! ) * ( m^n / m! )

という風に変形すると畳み込みの形で表せることがわかります.
なのでFFTをしたい気持ちになるけど値が大きいので誤差が怖い.

実はこの問題の mod の値 163577857 = 39 * 2^22 + 1 は NTT (高速剰余変換) が使える mod であることがググるとわかります. ちなみに最小の原始根は 23
なので整数のまま誤差なく計算できるので安心.

ソースコード

https://www.codechef.com/viewsolution/20267065

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

int mod_inv(ll a, ll m) {
    ll b = m, u = 1, v = 0;
    while (b > 0) {
        ll t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    return (u % m + m) % m;
}

ll mod_pow(ll x, ll y, ll md) {
    ll res = 1;
    while (y) {
        if (y & 1) res = res * x % md;
        x = x * x % md;
        y >>= 1;
    }
    return res;
}

template <int Mod, int PrimitiveRoot>
class fast_modulo_transform {
public:
    static vector<int> fft(vector<int> v, bool inv) {
        const int N = v.size();
        assert((N ^ (N & -N)) == 0);
        int ww = mod_pow(PrimitiveRoot, (Mod - 1) / N, Mod);
        if (inv) ww = mod_inv(ww, Mod);
        for (int m = N; m >= 2; m >>= 1) {
            const int mh = m >> 1;
            int w = 1;
            for (int i = 0; i < mh; ++i) {
                for (int j = i; j < N; j += m) {
                    const int k = j + mh;
                    int x = v[j] - v[k];
                    if (x < 0) x += Mod;
                    v[j] += -Mod + v[k];
                    if (v[j] < 0) v[j] += Mod;
                    v[k] = (1LL * w * x) % Mod;
                }
                w = (1LL * w * ww) % Mod;
            }
            ww = (1LL * ww * ww) % Mod;
        }

        int i = 0;
        for (int j = 1; j < N - 1; ++j) {
            for (int k = N >> 1; k > (i ^= k); k >>= 1);
            if (j < i) swap(v[i], v[j]);
        }
        if (inv) {
            const int inv_n = mod_inv(N, Mod);
            for (auto& x : v) {
                x = (1LL * x * inv_n) % Mod;
                assert(0 <= x && x < Mod);
            }
        }
        return v;
    }
    static vector<int> convolution(vector<int> f, vector<int> g) {
        int sz = 1;
        const int m = f.size() + g.size() - 1;
        while (sz < m) sz *= 2;
        f.resize(sz), g.resize(sz);
        f = fast_modulo_transform::fft(move(f), false); g = fast_modulo_transform::fft(move(g), false);
        for (int i = 0; i < sz; ++i) {
            f[i] = (1LL * f[i] * g[i]) % Mod;
        }

        return fast_modulo_transform::fft(move(f), true);
    }

    static int get_mod() {
        return Mod;
    }
};

const ll mod = 163577857;

using fmt = fast_modulo_transform<163577857, 23>;

int main()
{
    int N;
    cin >> N;
    vector<int> tmp(N + 1);
    tmp[0] = 1;
    for (int i = 1; i <= N; i++) {
        tmp[i] = (ll)tmp[i - 1] * i % mod;
    }
    tmp[N] = mod_inv(tmp[N], mod);
    for (int i = N; i >= 1; i--) {
        tmp[i - 1] = (ll)tmp[i] * i % mod;
    }
    vector<int> v1 = tmp, v2 = tmp;
    for (int i = 0, val = 1; i <= N; i++) {
        v1[i] = (ll)v1[i] * val % mod;
        val = (mod - val) % mod;
    }
    for (int i = 0; i <= N; i++) {
        v2[i] = (ll)v2[i] * mod_pow(i, N, mod) % mod;
    }
    auto res = fmt::convolution(v1, v2);
    for (int i = 1; i <= N; i++) {
        printf("%d%c", res[i], " \n"[i == N]);
    }
    return 0;
}

感想

これコンテスト中に通したかったなあ (第二種スターリング数であることまではわかってたので)
それはそうと NTT 速すぎてビビった.