#565839 (C++17) No.1145 Sums of Powers

提出ソース
結果

問題	No.1145 Sums of Powers
ユーザー	Chanyuh
提出日時	2020-10-11 13:09:10
言語	C++17 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	265 ms / 2,000 ms
コード長	17,313 bytes
コンパイル時間	2,139 ms
コンパイル使用メモリ	146,696 KB
最終ジャッジ日時	2025-01-15 06:35:25
ジャッジサーバーID （参考情報）	judge1 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
other	AC * 6
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include<iostream>
#include<string>
#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
#include<functional>
#include<iomanip>
#include<queue>
#include<ciso646>
#include<random>
#include<map>
#include<set>
#include<complex>
#include<bitset>
#include<stack>
#include<unordered_map>
#include<utility>
#include<tuple>
#include<cassert>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
const ll mod = 998244353;
const ll INF = (ll)1000000007 * 1000000007;
typedef pair<int, int> P;
#define stop char nyaa;cin>>nyaa;
#define rep(i,n) for(int i=0;i<n;i++)
#define per(i,n) for(int i=n-1;i>=0;i--)
#define Rep(i,sta,n) for(int i=sta;i<n;i++)
#define Per(i,sta,n) for(int i=n-1;i>=sta;i--)
#define rep1(i,n) for(int i=1;i<=n;i++)
#define per1(i,n) for(int i=n;i>=1;i--)
#define Rep1(i,sta,n) for(int i=sta;i<=n;i++)
typedef long double ld;
const ld eps = 1e-8;
const ld pi = acos(-1.0);
typedef pair<ll, ll> LP;
int dx[4]={1,-1,0,0};
int dy[4]={0,0,1,-1};
template<class T>bool chmax(T &a, const T &b) {if(a<b){a=b;return 1;}return 0;}
template<class T>bool chmin(T &a, const T &b) {if(b<a){a=b;return 1;}return 0;}

template<int mod>
struct ModInt {
    long long x;
    static constexpr int MOD = mod;
 
    ModInt() : x(0) {}
    ModInt(long long y) : x(y >= 0 ? y % mod : (mod - (-y) % mod) % mod) {}

    explicit operator int() const {return x;}
 
    ModInt &operator+=(const ModInt &p) {
        if((x += p.x) >= mod) x -= mod;
        return *this;
    }
    ModInt &operator-=(const ModInt &p) {
        if((x += mod - p.x) >= mod) x -= mod;
        return *this;
    }
    ModInt &operator*=(const ModInt &p) {
        x = (int)(1LL * x * p.x % mod);
        return *this;
    }
    ModInt &operator/=(const ModInt &p) {
        *this *= p.inverse();
        return *this;
    }
 
    ModInt operator-() const { return ModInt(-x); }
    ModInt operator+(const ModInt &p) const { return ModInt(*this) += p; }
    ModInt operator-(const ModInt &p) const { return ModInt(*this) -= p; }
    ModInt operator*(const ModInt &p) const { return ModInt(*this) *= p; }
    ModInt operator/(const ModInt &p) const { return ModInt(*this) /= p; }
    ModInt operator%(const ModInt &p) const { return ModInt(0); }         
 
    bool operator==(const ModInt &p) const { return x == p.x; }
    bool operator!=(const ModInt &p) const { return x != p.x; }
 
    ModInt inverse() const{
        int a = x, b = mod, u = 1, v = 0, t;
        while(b > 0) {
            t = a / b;
            a -= t * b;
            swap(a, b);
            u -= t * v;
            swap(u, v);
        }
        return ModInt(u);
    }

    ModInt power(long long n) const {
        ModInt ret(1), mul(x);
        while(n > 0) {
            if(n & 1)
            ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }

    ModInt power(const ModInt p) const{
        return ((ModInt)x).power(p.x);
    }

    friend ostream &operator<<(ostream &os, const ModInt<mod> &p) {
        return os << p.x;
    }
    friend istream &operator>>(istream &is, ModInt<mod> &a) {
        long long x;
        is >> x;
        a = ModInt<mod>(x);
        return (is);
    }
};

template <typename Mint>
struct NumberTheoreticTransformFriendlyModInt {

    vector<Mint> dw, idw;
    int max_base;
    Mint root;

    NumberTheoreticTransformFriendlyModInt() {
        const unsigned mod = Mint::MOD;
        assert(mod >= 3 && mod % 2 == 1);
        auto tmp = mod - 1;
        max_base = 0;
        while(tmp % 2 == 0)
            tmp >>= 1, max_base++;
        root = 2;
        while(root.power((mod - 1) >> 1) == 1)
            root += 1;
        assert(root.power(mod - 1) == 1);
        dw.resize(max_base);
        idw.resize(max_base);
        for(int i = 0; i < max_base; i++) {
            dw[i] = -root.power((mod - 1) >> (i + 2));
            idw[i] = Mint(1) / dw[i];
        }
    }

    void ntt(vector<Mint> &a) {
        const int n = (int)a.size();
        assert((n & (n - 1)) == 0);
        assert(__builtin_ctz(n) <= max_base);
        for(int m = n; m >>= 1;) {
            Mint w = 1;
            for(int s = 0, k = 0; s < n; s += 2 * m) {
                for(int i = s, j = s + m; i < s + m; ++i, ++j) {
                    auto x = a[i], y = a[j] * w;
                    a[i] = x + y, a[j] = x - y;
                }
                w *= dw[__builtin_ctz(++k)];
            }
        }
    }

    void intt(vector<Mint> &a, bool f = true) {
        const int n = (int)a.size();
        assert((n & (n - 1)) == 0);
        assert(__builtin_ctz(n) <= max_base);
        for(int m = 1; m < n; m *= 2) {
            Mint w = 1;
            for(int s = 0, k = 0; s < n; s += 2 * m) {
                for(int i = s, j = s + m; i < s + m; ++i, ++j) {
                    auto x = a[i], y = a[j];
                    a[i] = x + y, a[j] = (x - y) * w;
                }
                w *= idw[__builtin_ctz(++k)];
            }
        }
        if(f) {
            Mint inv_sz = Mint(1) / n;
            for(int i = 0; i < n; i++)
                a[i] *= inv_sz;
        }
    }

    vector<Mint> multiply(vector<Mint> a, vector<Mint> b) {
        int need = a.size() + b.size() - 1;
        int nbase = 1;
        while((1 << nbase) < need)
            nbase++;
        int sz = 1 << nbase;
        a.resize(sz, 0);
        b.resize(sz, 0);
        ntt(a);
        ntt(b);
        Mint inv_sz = Mint(1) / sz;
        for(int i = 0; i < sz; i++)
            a[i] *= b[i] * inv_sz;
        intt(a, false);
        a.resize(need);
        return a;
    }
};


template <typename T>
struct FormalPowerSeries : vector<T> {
    using vector<T>::vector;
    using P = FormalPowerSeries;
 
    using MULT = function<P(P, P)>;
    using FFT = function<void(P &)>;
    using SQRT = function<T(T)>;
 
    static MULT &get_mult() {
        static MULT mult = nullptr;
        return mult;
    }
 
    static void set_mult(MULT f) {
        get_mult() = f;
    }
 
    static FFT &get_fft() {
        static FFT fft = nullptr;
        return fft;
    }
 
    static FFT &get_ifft() {
        static FFT ifft = nullptr;
        return ifft;
    }
 
    static void set_fft(FFT f, FFT g) {
        get_fft() = f;
        get_ifft() = g;
    }
 
    static SQRT &get_sqrt() {
        static SQRT sqr = nullptr;
        return sqr;
    }
 
    static void set_sqrt(SQRT sqr) {
        get_sqrt() = sqr;
    }
 
    void shrink() {
        while(this->size() && this->back() == T(0))
            this->pop_back();
    }
 
    P operator+(const P &r) const { return P(*this) += r; }
 
    P operator+(const T &v) const { return P(*this) += v; }
 
    P operator-(const P &r) const { return P(*this) -= r; }
 
    P operator-(const T &v) const { return P(*this) -= v; }
 
    P operator*(const P &r) const { return P(*this) *= r; }
 
    P operator*(const T &v) const { return P(*this) *= v; }
 
    P operator/(const P &r) const { return P(*this) /= r; }
 
    P operator%(const P &r) const { return P(*this) %= r; }
 
    P &operator+=(const P &r) {
        if(r.size() > this->size())
            this->resize(r.size());
        for(int i = 0; i < r.size(); i++)
            (*this)[i] += r[i];
        return *this;
    }
 
    P &operator+=(const T &r) {
        if(this->empty())
            this->resize(1);
        (*this)[0] += r;
        return *this;
    }
 
    P &operator-=(const P &r) {
        if(r.size() > this->size())
            this->resize(r.size());
        for(int i = 0; i < r.size(); i++)
            (*this)[i] -= r[i];
        shrink();
        return *this;
    }
 
    P &operator-=(const T &r) {
        if(this->empty())
            this->resize(1);
        (*this)[0] -= r;
        shrink();
        return *this;
    }
 
    P &operator*=(const T &v) {
        const int n = (int)this->size();
        for(int k = 0; k < n; k++)
            (*this)[k] *= v;
        return *this;
    }
 
    P &operator*=(const P &r) {
        if(this->empty() || r.empty()) {
            this->clear();
            return *this;
        }
        assert(get_mult() != nullptr);
        return *this = get_mult()(*this, r);
    }
 
    P &operator%=(const P &r) { return *this -= *this / r * r; }
 
    P operator-() const {
        P ret(this->size());
        for(int i = 0; i < this->size(); i++)
            ret[i] = -(*this)[i];
        return ret;
    }
 
    P &operator/=(const P &r) {
        if(this->size() < r.size()) {
            this->clear();
            return *this;
        }
        int n = this->size() - r.size() + 1;
        return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n);
    }
 
    P dot(P r) const {
        P ret(min(this->size(), r.size()));
        for(int i = 0; i < ret.size(); i++)
            ret[i] = (*this)[i] * r[i];
        return ret;
    }
 
    P pre(int sz) const { return P(begin(*this), begin(*this) + min((int)this->size(), sz)); }
 
    P operator>>(int sz) const {
        if(this->size() <= sz)
            return {};
        P ret(*this);
        ret.erase(ret.begin(), ret.begin() + sz);
        return ret;
    }
 
    P operator<<(int sz) const {
        P ret(*this);
        ret.insert(ret.begin(), sz, T(0));
        return ret;
    }
 
    P rev(int deg = -1) const {
        P ret(*this);
        if(deg != -1)
            ret.resize(deg, T(0));
        reverse(begin(ret), end(ret));
        return ret;
    }
 
    P diff() const {
        const int n = (int)this->size();
        P ret(max(0, n - 1));
        for(int i = 1; i < n; i++)
            ret[i - 1] = (*this)[i] * T(i);
        return ret;
    }
 
    P integral() const {
        const int n = (int)this->size();
        P ret(n + 1);
        ret[0] = T(0);
        for(int i = 0; i < n; i++)
            ret[i + 1] = (*this)[i] / T(i + 1);
        return ret;
    }
 
    // F(0) must not be 0
    P inv(int deg = -1) const {
        assert(((*this)[0]) != T(0));
        const int n = (int)this->size();
        if(deg == -1)
            deg = n;
        if(get_fft() != nullptr) {
            P ret(*this);
            ret.resize(deg, T(0));
            return ret.inv_fast();
        }
        P ret({T(1) / (*this)[0]});
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1);
        }
        return ret.pre(deg);
    }
 
    // F(0) must be 1
    P log(int deg = -1) const {
        assert((*this)[0] == 1);
        const int n = (int)this->size();
        if(deg == -1)
            deg = n;
        return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
    }
 
    P sqrt(int deg = -1) const {
        const int n = (int)this->size();
        if(deg == -1)
            deg = n;
        if((*this)[0] == T(0)) {
            for(int i = 1; i < n; i++) {
                if((*this)[i] != T(0)) {
                    if(i & 1)
                        return {};
                    if(deg - i / 2 <= 0)
                        break;
                    auto ret = (*this >> i).sqrt(deg - i / 2);
                    if(ret.empty())
                        return {};
                    ret = ret << (i / 2);
                    if(ret.size() < deg)
                        ret.resize(deg, T(0));
                    return ret;
                }
            }
            return P(deg, 0);
        }
 
        P ret;
        if(get_sqrt() == nullptr) {
            assert((*this)[0] == T(1));
            ret = {T(1)};
        } else {
            auto sqr = get_sqrt()((*this)[0]);
            if(sqr * sqr != (*this)[0])
                return {};
            ret = {T(sqr)};
        }
 
        T inv2 = T(1) / T(2);
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
        }
        return ret.pre(deg);
    }
 
    // F(0) must be 0
    P exp(int deg = -1) const {
        assert((*this)[0] == T(0));
        const int n = (int)this->size();
        if(deg == -1)
            deg = n;
        if(get_fft() != nullptr) {
            P ret(*this);
            ret.resize(deg, T(0));
            return ret.exp_rec();
        }
        P ret({T(1)});
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1);
        }
        return ret.pre(deg);
    }
 
    P online_convolution_exp(const P &conv_coeff) const {
        const int n = (int)conv_coeff.size();
        assert((n & (n - 1)) == 0);
        vector<P> conv_ntt_coeff;
        auto &fft = get_fft();
        auto &ifft = get_ifft();
        for(int i = n; i >= 1; i >>= 1) {
            P g(conv_coeff.pre(i));
            fft(g);
            conv_ntt_coeff.emplace_back(g);
        }
        P conv_arg(n), conv_ret(n);
        auto rec = [&](auto rec, int l, int r, int d) -> void {
            if(r - l <= 16) {
                for(int i = l; i < r; i++) {
                    T sum = 0;
                    for(int j = l; j < i; j++)
                        sum += conv_arg[j] * conv_coeff[i - j];
                    conv_ret[i] += sum;
                    conv_arg[i] = i == 0 ? T(1) : conv_ret[i] / i;
                }
            } else {
                int m = (l + r) / 2;
                rec(rec, l, m, d + 1);
                P pre(r - l);
                for(int i = 0; i < m - l; i++)
                    pre[i] = conv_arg[l + i];
                fft(pre);
                for(int i = 0; i < r - l; i++)
                    pre[i] *= conv_ntt_coeff[d][i];
                ifft(pre);
                for(int i = 0; i < r - m; i++)
                    conv_ret[m + i] += pre[m + i - l];
                rec(rec, m, r, d + 1);
            }
        };
        rec(rec, 0, n, 0);
        return conv_arg;
    }
 
    P exp_rec() const {
        assert((*this)[0] == T(0));
        const int n = (int)this->size();
        int m = 1;
        while(m < n)
            m *= 2;
        P conv_coeff(m);
        for(int i = 1; i < n; i++)
            conv_coeff[i] = (*this)[i] * i;
        return online_convolution_exp(conv_coeff).pre(n);
    }
 
    P inv_fast() const {
        assert(((*this)[0]) != T(0));
 
        const int n = (int)this->size();
        P res{T(1) / (*this)[0]};
 
        for(int d = 1; d < n; d <<= 1) {
            P f(2 * d), g(2 * d);
            for(int j = 0; j < min(n, 2 * d); j++)
                f[j] = (*this)[j];
            for(int j = 0; j < d; j++)
                g[j] = res[j];
            get_fft()(f);
            get_fft()(g);
            for(int j = 0; j < 2 * d; j++)
                f[j] *= g[j];
            get_ifft()(f);
            for(int j = 0; j < d; j++) {
                f[j] = 0;
                f[j + d] = -f[j + d];
            }
            get_fft()(f);
            for(int j = 0; j < 2 * d; j++)
                f[j] *= g[j];
            get_ifft()(f);
            for(int j = 0; j < d; j++)
                f[j] = res[j];
            res = f;
        }
        return res.pre(n);
    }
 
    P pow(int64_t k, int deg = -1) const {
        const int n = (int)this->size();
        if(deg == -1)
            deg = n;
        for(int i = 0; i < n; i++) {
            if((*this)[i] != T(0)) {
                T rev = T(1) / (*this)[i];
                P ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].power(k));
                if(i * k > deg)
                    return P(deg, T(0));
                ret = (ret << (i * k)).pre(deg);
                if(ret.size() < deg)
                    ret.resize(deg, T(0));
                return ret;
            }
        }
        return *this;
    }
 
    T eval(T x) const {
        T r = 0, w = 1;
        for(auto &v : *this) {
            r += w * v;
            w *= x;
        }
        return r;
    }
 
    P pow_mod(int64_t n, P mod) const {
        P modinv = mod.rev().inv();
        auto get_div = [&](P base) {
            if(base.size() < mod.size()) {
                base.clear();
                return base;
            }
            int n = base.size() - mod.size() + 1;
            return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n);
        };
        P x(*this), ret{1};
        while(n > 0) {
            if(n & 1) {
                ret *= x;
                ret -= get_div(ret) * mod;
            }
            x *= x;
            x -= get_div(x) * mod;
            n >>= 1;
        }
        return ret;
    }
};
using modint = ModInt<mod>;
using FPS = FormalPowerSeries<modint>;

int n,m;
modint a[100010];
NumberTheoreticTransformFriendlyModInt<modint> ntt;

FPS::P calc(int l,int r){
    if(r-l<1) return {1};
    if(r-l==1) return {1,-a[l]};
    int mid=(l+r)/2;
    auto res= ntt.multiply(calc(l,mid),calc(mid,r));
    return FPS::P(begin(res),end(res));
}


void solve(){
    auto mult=[&](const FPS::P &a,const FPS::P &b){
        auto f=ntt.multiply(a,b);
        return FPS::P(begin(f),end(f));
    };
    FPS::set_mult(mult);
    FPS::set_fft([&](FPS::P &a){return ntt.ntt(a);},[&](FPS::P &a){return ntt.intt(a);});

    cin >> n >> m;
    rep(i,n) cin >> a[i];
    FPS F=calc(0,n);
    F=-F.log(m+2);
    F=F.diff();
    rep(i,m){
        cout << F[i] << " ";
    }
    cout << "" << endl;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout << fixed << setprecision(50);
    solve();
}
yukicoder

結果

ソースコード
