#316949 (C++17) No.215 素数サイコロと合成数サイコロ (3-Hard)

• 問題: No.215 素数サイコロと合成数サイコロ (3-Hard)
• ユーザー: yosupot
• 提出日時: 2019-02-16 23:05:00
• 言語: C++17 (gcc 15.2.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 3,113 ms / 4,000 ms
• コード長: 10,910 bytes
• コンパイル時間: 3,569 ms
• コンパイル使用メモリ: 224,952 KB
• 最終ジャッジ日時: 2025-01-06 21:12:40
• ジャッジサーバーID （参考情報）: judge5 / judge3

ファイルパターン	結果
other	AC * 2

ソースコード

#include <bits/stdc++.h>

using namespace std;
using uint = unsigned int;
using ll = long long;
using ull = unsigned long long;
constexpr ll TEN(int n) { return (n==0) ? 1 : 10*TEN(n-1); }
template<class T> using V = vector<T>;
template<class T> using VV = V<V<T>>;

template <class T> ostream& operator<<(ostream& os, const V<T>& v) {
    os << "[";
    for (auto d : v) os << d << ", ";
    return os << "]";
}

template <class T, class U>
ostream& operator<<(ostream& os, const pair<T, U>& p) {
    return os << "P(" << p.first << ", " << p.second << ")";
}

template <uint MD> struct ModInt {
    using M = ModInt;
    const static M G;
    uint v;
    ModInt() : v(0) {}
    ModInt(ll _v) { set_v(_v % MD + MD); }
    M& set_v(uint _v) {
        v = (_v < MD) ? _v : _v - MD;
        return *this;
    }
    explicit operator bool() const { return v != 0; }
    M operator-() const { return M() - *this; }
    M operator+(const M& r) const { return M().set_v(v + r.v); }
    M operator-(const M& r) const { return M().set_v(v + MD - r.v); }
    M operator*(const M& r) const { return M().set_v(ull(v) * r.v % MD); }
    M operator/(const M& r) const { return *this * r.inv(); }
    M& operator+=(const M& r) { return *this = *this + r; }
    M& operator-=(const M& r) { return *this = *this - r; }
    M& operator*=(const M& r) { return *this = *this * r; }
    M& operator/=(const M& r) { return *this = *this / r; }
    M pow(ll n) const {
        M x = *this, r = 1;
        while (n) {
            if (n & 1) r *= x;
            x *= x;
            n >>= 1;
        }
        return r;
    }
    M inv() const { return pow(MD - 2); }
    friend ostream& operator<<(ostream& os, const M& r) { return os << r.v; }
};
using Mint = ModInt<TEN(9) + 7>;

using D = double;
const D PI = acos(D(-1));
using Pc = complex<D>;

void fft(bool type, V<Pc>& a) {
    int n = int(a.size()), s = 0;
    while ((1 << s) < n) s++;
    assert(1 << s == n);

    static V<Pc> ep[30];
    if (!ep[s].size()) {
        for (int i = 0; i < n; i++) {
            ep[s].push_back(polar<D>(1, i * 2 * PI / n));
        }
    }
    V<Pc> b(n);
    for (int i = 1; i <= s; i++) {
        int w = 1 << (s - i);
        for (int y = 0; y < n / 2; y += w) {
            Pc now = ep[s][y];
            if (type) now = conj(now);
            for (int x = 0; x < w; x++) {
                auto l = a[y << 1 | x];
                auto u = now, v = a[y << 1 | x | w];
                auto r = Pc(u.real() * v.real() - u.imag() * v.imag(),
                            u.real() * v.imag() + u.imag() * v.real());
                b[y | x] = l + r;
                b[y | x | n >> 1] = l - r;
            }
        }
        swap(a, b);
    }
}

template <class Mint, int K = 3, int SHIFT = 11>
V<Mint> multiply(const V<Mint>& a, const V<Mint>& b) {
    int A = int(a.size()), B = int(b.size());
    if (!A || !B) return {};
    int lg = 0;
    while ((1 << lg) < A + B - 1) lg++;
    int N = 1 << lg;

    VV<Pc> x(K, V<Pc>(N)), y(K, V<Pc>(N));
    for (int ph = 0; ph < K; ph++) {
        V<Pc> z(N);
        for (int i = 0; i < N; i++) {
            D nx = 0, ny = 0;
            if (i < A) nx = (a[i].v >> (ph * SHIFT)) & ((1 << SHIFT) - 1);
            if (i < B) ny = (b[i].v >> (ph * SHIFT)) & ((1 << SHIFT) - 1);
            z[i] = Pc(nx, ny);
        }
        fft(false, z);
        for (int i = 0; i < N; i++) {
            z[i] *= 0.5;
        }
        for (int i = 0; i < N; i++) {
            int j = (i) ? N - i : 0;
            x[ph][i] = Pc(z[i].real() + z[j].real(), z[i].imag() - z[j].imag());
            y[ph][i] =
                    Pc(z[i].imag() + z[j].imag(), -z[i].real() + z[j].real());
        }
    }
    VV<Pc> z(K, V<Pc>(N));
    for (int xp = 0; xp < K; xp++) {
        for (int yp = 0; yp < K; yp++) {
            int zp = (xp + yp) % K;
            for (int i = 0; i < N; i++) {
                if (xp + yp < K) {
                    z[zp][i] += x[xp][i] * y[yp][i];
                } else {
                    z[zp][i] += x[xp][i] * y[yp][i] * Pc(0, 1);
                }
            }
        }
    }
    for (int ph = 0; ph < K; ph++) {
        fft(true, z[ph]);
    }
    V<Mint> c(A + B - 1);
    Mint base = 1;
    for (int ph = 0; ph < 2 * K - 1; ph++) {
        for (int i = 0; i < A + B - 1; i++) {
            if (ph < K) {
                c[i] += Mint(ll(round(z[ph][i].real() / N))) * base;
            } else {
                c[i] += Mint(ll(round(z[ph - K][i].imag() / N))) * base;
            }
        }
        base *= 1 << SHIFT;
    }
    return c;
}

template <class D> struct Poly {
    V<D> v;
    Poly(const V<D>& _v = {}) : v(_v) { shrink(); }
    void shrink() { while (v.size() && !v.back()) v.pop_back(); }

    int size() const { return int(v.size()); }
    D freq(int p) const { return (p < size()) ? v[p] : D(0); }

    Poly operator+(const Poly& r) const {
        auto n = max(size(), r.size());
        V<D> res(n);
        for (int i = 0; i < n; i++) res[i] = freq(i) + r.freq(i);
        return res;
    }
    Poly operator-(const Poly& r) const {
        int n = max(size(), r.size());
        V<D> res(n);
        for (int i = 0; i < n; i++) res[i] = freq(i) - r.freq(i);
        return res;
    }
    Poly operator*(const Poly& r) const { return {multiply(v, r.v)}; }
    Poly operator*(const D& r) const {
        int n = size();
        V<D> res(n);
        for (int i = 0; i < n; i++) res[i] = v[i] * r;
        return res;
    }
    Poly operator/(const Poly& r) const {
        if (size() < r.size()) return {{}};
        int n = size() - r.size() + 1;
        return (rev().pre(n) * r.rev().inv(n)).pre(n).rev();
    }
    Poly operator%(const Poly& r) const { return *this - (*this) / r * r; }
    Poly operator<<(int s) const {
        V<D> res(size() + s);
        for (int i = 0; i < size(); i++) res[i + s] = v[i];
        return res;
    }
    Poly operator>>(int s) const {
        if (size() <= s) return Poly();
        V<D> res(size() - s);
        for (int i = 0; i < size() - s; i++) res[i] = v[i + s];
        return res;
    }
    Poly& operator+=(const Poly& r) { return *this = *this + r; }
    Poly& operator-=(const Poly& r) { return *this = *this - r; }
    Poly& operator*=(const Poly& r) { return *this = *this * r; }
    Poly& operator*=(const D& r) { return *this = *this * r; }
    Poly& operator/=(const Poly& r) { return *this = *this / r; }
    Poly& operator%=(const Poly& r) { return *this = *this % r; }
    Poly& operator<<=(const size_t& n) { return *this = *this << n; }
    Poly& operator>>=(const size_t& n) { return *this = *this >> n; }

    Poly pre(int le) const { return {{v.begin(), v.begin() + min(size(), le)}}; }
    Poly rev(int n = -1) const {
        V<D> res = v;
        if (n != -1) res.resize(n);
        reverse(begin(res), end(res));
        return Poly(res);
    }
    // f * f.inv() = 1 + g(x)x^m
    Poly inv(int m) const {
        Poly res = Poly({D(1) / freq(0)});
        for (int i = 1; i < m; i *= 2) {
            res = (res * D(2) - res * res * pre(2 * i)).pre(2 * i);
        }
        return res;
    }
    // TODO: reuse inv
    Poly pow_mod(ll n, const Poly& mod) {
        Poly x = *this, r = {{1}};
        while (n) {
            if (n & 1) r = r * x % mod;
            x = x * x % mod;
            n >>= 1;
        }
        return r;
    }

    friend ostream& operator<<(ostream& os, const Poly& p) {
        if (p.size() == 0) return os << "0";
        for (auto i = 0; i < p.size(); i++) {
            if (p.v[i]) {
                os << p.v[i] << "x^" << i;
                if (i != p.size() - 1) os << "+";
            }
        }
        return os;
    }
};

template <class Mint> Poly<Mint> berlekamp_massey(const V<Mint>& s) {
    int n = int(s.size());
    V<Mint> b = {Mint(-1)}, c = {Mint(-1)};
    Mint y = Mint(1);
    for (int ed = 1; ed <= n; ed++) {
        int l = int(c.size()), m = int(b.size());
        Mint x = 0;
        for (int i = 0; i < l; i++) {
            x += c[i] * s[ed - l + i];
        }
        b.push_back(0);
        m++;
        if (!x) continue;
        Mint freq = x / y;
        if (l < m) {
            // use b
            auto tmp = c;
            c.insert(begin(c), m - l, Mint(0));
            for (int i = 0; i < m; i++) {
                c[m - 1 - i] -= freq * b[m - 1 - i];
            }
            b = tmp;
            y = x;
        } else {
            // use c
            for (int i = 0; i < m; i++) {
                c[l - 1 - i] -= freq * b[m - 1 - i];
            }
        }
    }
    return c;
}

const int MD = 610 * 13;
const int MC = 310;
const V<int> pd = {2, 3, 5, 7, 11, 13};
const V<int> cd = {4, 6, 8, 9, 10, 12};

int main() {
    ll n; int x, y;
    cin >> n >> x >> y;
    V<Mint> co(MD+1); co[0] = 1;
    {
        int c = x;
        V<Mint> buf[6];
        for (int i = 0; i < 6; i++) {
            buf[i] = V<Mint>(MD);
            buf[i][0] = 1;
        }
        for (int nc = 0; nc < c; nc++) {
            for (int i = 0; i < 6; i++) {
                for (int nw = MD-1; nw >= 0; nw--) {
                    buf[i][nw] = 0;
                    if (i) buf[i][nw] += buf[i-1][nw];
                    int d = pd[i];
                    if (nw < d) continue;
                    buf[i][nw] += buf[i][nw-d];
                }
            }
        }
        auto co2 = multiply(co, buf[5]);
        co2.resize(co.size());
        co = co2;
    }
    {
        int c = y;
        V<Mint> buf[6];
        for (int i = 0; i < 6; i++) {
            buf[i] = V<Mint>(MD);
            buf[i][0] = 1;
        }
        for (int nc = 0; nc < c; nc++) {
            for (int i = 0; i < 6; i++) {
                for (int nw = MD-1; nw >= 0; nw--) {
                    buf[i][nw] = 0;
                    if (i) buf[i][nw] += buf[i-1][nw];
                    int d = cd[i];
                    if (nw < d) continue;
                    buf[i][nw] += buf[i][nw-d];
                }
            }
        }
        auto co2 = multiply(co, buf[5]);
        co2.resize(co.size());
        co = co2;
    }
    co[0] = -1;
    V<Mint> dp(2*MD);
    dp[0] = 1;
    for (int i = 0; i < 2*MD; i++) {
        for (int j = 1; j < int(co.size()); j++) {
            if (i+j >= 2*MD) break;
            dp[i+j] += dp[i]*co[j];
        }
    }
    auto pol = Poly<Mint>({0, 1}).pow_mod(n, berlekamp_massey<Mint>(dp));
    V<Mint> buf(MD); buf[0] = 1;
    V<Mint> sm(MD);
    for (int i = 0; i < MD; i++) {
        //v[i] * f
        for (int j = 0; j < MD; j++) {
            sm[j] += pol.freq(i)*buf[j];
        }
        for (int j = 1; j < MD; j++) {
            buf[j] += buf[0]*co[j];
        }
        for (int j = 0; j < MD-1; j++) {
            buf[j] = buf[j+1];
        }
        buf[MD-1] = 0;
    }
    Mint ans = 0;
    for (int i = 0; i < MD; i++) {
        ans += sm[i];
    }
    cout << ans << endl;
    return 0;
}