#1137968 (C++23) No.3404 形式群法則

• 問題: No.3404 形式群法則
• ユーザー: sclara
• 提出日時: 2025-12-02 16:14:40
• 言語: C++23 (gcc 13.3.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 104 ms / 2,000 ms
• コード長: 17,211 bytes
• コンパイル時間: 3,347 ms
• コンパイル使用メモリ: 301,220 KB
• 実行使用メモリ: 7,848 KB
• 最終ジャッジ日時: 2025-12-10 23:30:48
• 合計ジャッジ時間: 4,821 ms
• ジャッジサーバーID （参考情報）: judge5 / judge1

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

ソースコード

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef std::pair<long long, long long> P;
typedef std::priority_queue<P, std::vector<P>, std::greater<P>> PQ;
typedef std::complex<double> cd;

struct P3 {
    long long first, second, third;
};

struct P3P {
    P first, second, third;
};

struct compP3{
    bool operator()(const P3 &p1,const P3 &p2) const {
        if (p1.first != p2.first) return p1.first < p2.first;
        if (p1.second != p2.second) return p1.second < p2.second;
        else return p1.third < p2.third;
    }
};

struct gcompP3{
    bool operator()(const P3 &p1,const P3 &p2) const {
        if (p1.first != p2.first) return p1.first > p2.first;
        if (p1.second != p2.second) return p1.second > p2.second;
        else return p1.third > p2.third;
    }
};

const double PI = acos(-1.0);

bool ckran(int a, int n) {
    return (a >= 0 && a < n);
}

void yn (bool f) {
    if (f) std::cout << "Yes" << '\n';
    else std::cout << "No" << '\n';
}

long long pplus(P a) {
    return a.first + a.second;
}

long long pminus(P a) {
    return a.first - a.second;
}

long long ptime(P a) {
    return a.first * a.second;
}

long long pdiv(P a) {
    return a.first / a.second;
}

template<typename T, typename U>
bool chmax(T& a, U b) {
    if (a < b) {
        a = b;
        return true;
    } else {
        return false;
    }
}

template<typename T, typename U>
bool chmin(T& a, U b) {
    if (a > b) {
        a = b;
        return true;
    } else {
        return false;
    }
}

template<typename T>
void outspace(std::vector<T> a) {
    int n = a.size();
    for (int i = 0; i < n; ++i) {
        std::cout << a[i];
        if (i != n - 1) std::cout << " ";
        else std::cout << std::endl;
    }
}

void outspace(P a) {
    std::cout << a.first << ' ' << a.second << '\n';
}

template<typename T>
void outendl(std::vector<T> a) {
    int n = a.size();
    for (int i = 0; i < n; ++i) {
        std::cout << a[i] << '\n';
    }
}

void outdouble(long double a) {
    std::cout << std::fixed << std::setprecision(10) << a << std::endl;
}

std::vector<long long> lltovec(long long n, long long base = 10, long long minsize = -1) {
    std::vector<long long> res;
    while (minsize-- > 0 || n > 0) {
        res.push_back(n % base);
        n /= base;
    }
    // std::reverse(res.begin(), res.end());
    return res;
}

long long vectoll(std::vector<long long> vec, long long base = 10) {
    long long res = 0;
    std::reverse(vec.begin(), vec.end());
    for (auto i : vec) {
        res *= base;
        res += i;
    }
    std::reverse(vec.begin(), vec.end());
    return res;
}

class uftree {

    private:
    std::vector<int> union_find_tree;
    std::vector<int> rank;
    std::vector<int> nodes_count;

    public:

    uftree (int uftree_size) {
        union_find_tree.resize(uftree_size);
        rank.resize(uftree_size);
        nodes_count.resize(uftree_size);

        for (int i = 0; i < uftree_size; ++i) {
            union_find_tree[i] = i;
            rank[i] = 0;
            nodes_count[i] = 1;
        }
    }

    int find (int n) {
        if (union_find_tree[n] == n) return n;
        else {
            union_find_tree[n] = find(union_find_tree[n]);
            return union_find_tree[n];
        }
    }

    void unite (int x, int y) {
        x = find(x);
        y = find(y);

        if (x == y) return;

        if (rank[x] > rank[y]) {
            union_find_tree[y] = union_find_tree[x];
            nodes_count[x] += nodes_count[y];
        } else {
            union_find_tree[x] = union_find_tree[y];
            nodes_count[y] += nodes_count[x];
            if (rank[x] == rank[y]) rank[y]++;
        }
    }

    bool connected (int x, int y) {
        x = find(x);
        y = find(y);
        return x == y;
    }
    
    long long groupCount (int x) {
        x = find(x);
        return nodes_count[x];
    }

};


class maxSegmentTree {
    int n;
    std::vector<long long> tree, lazy;
    const long long MINI = -4e18;

    public:
    maxSegmentTree(int size) {
        n = size;
        tree.assign(4 * n, MINI);
        lazy.assign(4 * n, MINI);
    }

    void push(int node, int start, int end) {
        if (lazy[node] != MINI) {
            // 遅延値を現在のノードに適用
            tree[node] = std::max(tree[node], lazy[node]);

            // 子ノードに遅延値を伝播
            if (start != end) {
                lazy[node * 2] = std::max(lazy[node * 2], lazy[node]);
                lazy[node * 2 + 1] = std::max(lazy[node * 2 + 1], lazy[node]);
            }

            // 現在のノードの遅延値をクリア
            lazy[node] = MINI;
        }
    }

    void updateRange(int l, int r, long long value, int node = 1, int start = 0, int end = -1) {
        if (end == -1) end = n - 1;

        push(node, start, end);

        if (start > r || end < l) {
            // 完全に範囲外
            return;
        }

        if (start >= l && end <= r) {
            // 完全に範囲内
            lazy[node] = value;
            push(node, start, end);
            return;
        }

        // 部分的に範囲が重なる場合
        int mid = (start + end) / 2;
        updateRange(l, r, value, node * 2, start, mid);
        updateRange(l, r, value, node * 2 + 1, mid + 1, end);
        tree[node] = std::max(tree[node * 2], tree[node * 2 + 1]);
    }

    long long queryRange(int l, int r, int node = 1, int start = 0, int end = -1) {
        if (end == -1) end = n - 1;

        push(node, start, end);

        if (start > r || end < l) {
            // 完全に範囲外
            return MINI;
        }

        if (start >= l && end <= r) {
            // 完全に範囲内
            return tree[node];
        }

        // 部分的に範囲が重なる場合
        int mid = (start + end) / 2;
        long long leftQuery = queryRange(l, r, node * 2, start, mid);
        long long rightQuery = queryRange(l, r, node * 2 + 1, mid + 1, end);
        return std::max(leftQuery, rightQuery);
    }
};

class sumSegmentTree {
    int n;
    std::vector<long long> tree, lazy;

    public:
    sumSegmentTree(int size) {
        n = size;
        tree.assign(4 * n, 0);
        lazy.assign(4 * n, 0);
    }

    void push(int node, int start, int end) {
        if (lazy[node] != 0) {
            // 遅延値を現在のノードに適用
            tree[node] += (end - start + 1) * lazy[node];

            // 子ノードに遅延値を伝播
            if (start != end) {
                lazy[node * 2] += lazy[node];
                lazy[node * 2 + 1] += lazy[node];
            }

            // 現在のノードの遅延値をクリア
            lazy[node] = 0;
        }
    }

    void updateRange(int l, int r, long long value, int node = 1, int start = 0, int end = -1) {
        if (end == -1) end = n - 1;

        push(node, start, end);

        if (start > r || end < l) {
            // 完全に範囲外
            return;
        }

        if (start >= l && end <= r) {
            // 完全に範囲内
            lazy[node] += value;
            push(node, start, end);
            return;
        }

        // 部分的に範囲が重なる場合
        int mid = (start + end) / 2;
        updateRange(l, r, value, node * 2, start, mid);
        updateRange(l, r, value, node * 2 + 1, mid + 1, end);
        tree[node] = tree[node * 2] + tree[node * 2 + 1];
    }

    long long queryRange(int l, int r, int node = 1, int start = 0, int end = -1) {
        if (end == -1) end = n - 1;

        push(node, start, end);

        if (start > r || end < l) {
            // 完全に範囲外
            return 0;
        }

        if (start >= l && end <= r) {
            // 完全に範囲内
            return tree[node];
        }

        // 部分的に範囲が重なる場合
        int mid = (start + end) / 2;
        long long leftQuery = queryRange(l, r, node * 2, start, mid);
        long long rightQuery = queryRange(l, r, node * 2 + 1, mid + 1, end);
        return leftQuery + rightQuery;
    }

    long long lower_bound(int l, long long value, int node = 1, int start = 0, int end = -1, long long curval = 0) {
        if (end == -1) {
            end = n - 1;
            push(node, start, end);
            curval += queryRange(0, l - 1);
            if (tree[node] < value + curval) return n;
        }

        if (start == end) {
            if (tree[node] + curval < value) return start + 1;
            else return start;
        }

        int mid = (start + end) / 2;
        push(node * 2, start, mid);
        push(node * 2 + 1, mid + 1, end);
        if (tree[node * 2] + curval >= value) {
            return lower_bound(l, value, node * 2, start, mid, curval);
        } else {
            curval += tree[node * 2];
            return lower_bound(l, value, node * 2 + 1, mid + 1, end, curval);
        }
    }
};

static const long long MOD = 998244353;
static const int MAXN = 1000000;  // 必要に応じて大きく設定

// 階乗・階乗逆元の配列
static long long fact[MAXN+1], invFact[MAXN+1];

// 繰り返し二乗法 (a^b mod M)
long long modpow(long long a, long long b, long long M) {
    long long ret = 1 % M;
    a %= M;
    while (b > 0) {
        if (b & 1) ret = (ret * a) % M;
        a = (a * a) % M;
        b >>= 1;
    }
    return ret;
}

// 前処理: 階乗と逆元のテーブルを作る
void initFactorials() {
    // 階乗テーブル
    fact[0] = 1;
    for (int i = 1; i <= MAXN; i++) {
        fact[i] = fact[i-1] * i % MOD;
    }
    // 階乗逆元テーブル
    invFact[MAXN] = modpow(fact[MAXN], MOD - 2, MOD);  // フェルマーの小定理で逆元を計算
    for (int i = MAXN; i >= 1; i--) {
        invFact[i-1] = invFact[i] * i % MOD;
    }
}

// 組み合わせ数 C(n, r)
long long comb(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] * invFact[r] % MOD * invFact[n - r] % MOD;
}

vector<pair<long long, long long>> rlell(vector<long long> a) {
    vector<pair<long long, long long>> res;
    if (a.empty()) return res;
    res.push_back({a[0], 1});
    long long n = a.size();
    for (int i = 1; i < n; ++i) {
        if (a[i] == a[i - 1]) {
            res.back().second++;
        } else {
            res.push_back({a[i], 1});
        }
    }
    return res;
}

vector<pair<char, long long>> rlest(string a) {
    vector<pair<char, long long>> res;
    if (a.empty()) return res;
    res.push_back({a[0], 1});
    long long n = a.size();
    for (int i = 1; i < n; ++i) {
        if (a[i] == a[i - 1]) {
            res.back().second++;
        } else {
            res.push_back({a[i], 1});
        }
    }
    return res;
}

// --------------------- 多項式ユーティリティ ---------------------

// a(x) * b(x) を計算して、次数 n-1 までで切り捨てる
vector<long long> poly_mul(const vector<long long>& a,
                       const vector<long long>& b,
                       int n,
                       long long mod) {
    if (a.empty() || b.empty()) return {};
    int na = (int)a.size();
    int nb = (int)b.size();
    int sz = min(na + nb - 1, n);
    vector<long long> res(sz, 0);

    for (int i = 0; i < na; i++) {
        if (a[i] == 0) continue;
        for (int j = 0; j < nb && i + j < n; j++) {
            if (b[j] == 0) continue;
            long long val = res[i + j] + a[i] * b[j];
            if (val >= (1LL << 62)) val %= mod; // 念のためデカすぎるのを抑える
            res[i + j] = val % mod;
        }
    }
    return res;
}

// res += a (mod mod), 長さは min(res.size(), a.size())
void poly_add_inplace(vector<long long>& res,
                      const vector<long long>& a,
                      long long mod) {
    int m = min(res.size(), a.size());
    for (int i = 0; i < m; i++) {
        res[i] += a[i];
        if (res[i] >= mod) res[i] -= mod;
    }
}

// res += c * a (mod mod), 長さは min(res.size(), a.size())
void poly_add_scaled(vector<long long>& res,
                     const vector<long long>& a,
                     long long c,
                     long long mod) {
    if (c % mod == 0) return;
    c %= mod;
    if (c < 0) c += mod;
    int m = min(res.size(), a.size());
    for (int i = 0; i < m; i++) {
        long long val = res[i] + c * a[i];
        if (val >= (1LL << 62)) val %= mod; // 念のため
        res[i] = val % mod;
    }
}

// --------------------- Brent–Kung 風平方分割合成 ---------------------
//
// h(x) = f(g(x)) mod x^n を計算する
// ・f, g: 係数ベクトル (f[i] は x^i の係数)
// ・g(0) == 0 を想定（形式的べき級数の合成）
// ・mod は素数（最大 1e9 程度）
// ・返り値は長さ n のベクトル
//
vector<long long> compose_brent_kung(
    const vector<long long>& f_in,
    const vector<long long>& g_in,
    int n,
    long long mod
) {
    if (n == 0) return {};

    vector<long long> f = f_in;
    vector<long long> g = g_in;

    if ((int)f.size() > n) f.resize(n);
    if ((int)g.size() > n) g.resize(n);
    f.resize(n, 0);
    g.resize(n, 0);

    // g(0) == 0 のチェック
    if (g[0] % mod != 0) {
        cerr << "Error: this implementation assumes g(0) == 0 "
             << "for formal power series composition.\n";
        // 必要ならここで素直な Horner 法にフォールバックしてもよい
        // が、今回はエラー表示だけにしておく
    }

    if (n == 1) {
        // h(0) = f(g(0)) = f(0) （g(0) = 0 を仮定）
        vector<long long> h(1);
        h[0] = (f[0] % mod + mod) % mod;
        return h;
    }

    // m ≒ sqrt(n)
    int m = (int)ceil(sqrt((long double)n));
    int k = (n + m - 1) / m;  // ブロック数

    // 1) baby steps: g^0, g^1, ..., g^{m-1}
    vector<vector<long long>> pow_g(m);
    pow_g[0] = vector<long long>(1, 1); // g^0 = 1
    if (m > 1) pow_g[1] = g;
    for (int i = 2; i < m; i++) {
        pow_g[i] = poly_mul(pow_g[i - 1], g, n, mod);
    }

    // 2) G = g^m
    vector<long long> G;
    if (m == 1) {
        G = g;
    } else {
        G = poly_mul(pow_g[m - 1], g, n, mod); // g^{m-1} * g = g^m
    }

    // 3) giant steps: G^0, G^1, ..., G^{k-1}
    vector<vector<long long>> pow_G(k);
    pow_G[0] = vector<long long>(1, 1); // G^0 = 1
    for (int j = 1; j < k; j++) {
        pow_G[j] = poly_mul(pow_G[j - 1], G, n, mod);
    }

    // 4) f をブロックに分割しつつ、
    //    H_j(x) = F_j(g(x)) を計算して、
    //    h(x) += H_j(x) * G(x)^j を足していく
    vector<long long> h(n, 0);

    for (int j = 0; j < k; j++) {
        int start = j * m;
        if (start >= n) break;

        // F_j(x) = sum_{i=0}^{m-1} f[start + i] x^i
        // H_j(x) = F_j(g(x)) = sum_{i=0}^{m-1} f[start + i] * g(x)^i
        vector<long long> Hj(n, 0);

        for (int i = 0; i < m; i++) {
            int idx = start + i;
            if (idx >= n) break;
            long long coef = f[idx] % mod;
            if (coef == 0) continue;
            poly_add_scaled(Hj, pow_g[i], coef, mod);
        }

        // T = Hj(x) * G(x)^j (mod x^n)
        vector<long long> T = poly_mul(Hj, pow_G[j], n, mod);
        poly_add_inplace(h, T, mod);
    }

    for (int i = 0; i < n; i++) {
        h[i] %= mod;
        if (h[i] < 0) h[i] += mod;
    }
    return h;
}

// vector<int> cl = {-1, 0, 1, 0, -1};

// vector<ll> cl = {-1, 0, 1, 1, 1, 0, -1, -1, -1, 0};

int main()
{
    ll n, m, p;
    cin >> n >> m >> p;
    vector<ll> a(n);
    for (int i = 0; i < n; ++i) {
        cin >> a[i];
        a[i] %= p;
    }
    vector<ll> sin(n + 1, 0), cos(n + 1, 0);
    cos[0] = 1;
    ll cur = 1;
    for (int i = 1; i <= n; ++i) {
        cur *= i;
        cur %= p;
        if (i % 2 == 1) {
            sin[i] = modpow(cur, p - 2, p);
            if (i % 4 == 3) {
                sin[i] *= -1;
                sin[i] += p;
            }
        } else {
            cos[i] = modpow(cur, p - 2, p);
            if (i % 4 == 2) {
                cos[i] *= -1;
                cos[i] += p;
            }
        }
    }
    vector<ll> sec(n + 1, 0);
    sec[0] = 1;
    for (int i = 1; i <= n; ++i) {
        for (int j = 0; j < i; ++j) {
            sec[i] -= cos[i - j] * sec[j];
            sec[i] %= p;
        }
        sec[i] += p;
        sec[i] %= p;
    }
    vector<ll> tan(n + 1, 0);
    for (int i = 1; i <= n; ++i) {
        for (int j = 0; j <= i; ++j) {
            tan[i] += sin[j] * sec[i - j];
            tan[i] %= p;
        }
        tan[i] += p;
        tan[i] %= p;
    }
    vector<ll> arctan(n + 1, 0);
    for (int i = 1; i <= n; i += 2) {
        arctan[i] = modpow(i, p - 2, p) * (m % p);
        arctan[i] %= p;
        if (i % 4 == 3) {
            arctan[i] *= -1;
            arctan[i] += p;
        }
    }
    reverse(a.begin(), a.end());
    a.push_back(0);
    reverse(a.begin(), a.end());
    vector<ll> gx = compose_brent_kung(arctan, a, n + 1, p);
    vector<ll> ans = compose_brent_kung(tan, gx, n + 1, p);
    reverse(ans.begin(), ans.end());
    ans.pop_back();
    reverse(ans.begin(), ans.end());
    outspace(ans);
}