cpp/fps.hpp

• View this file on GitHub
• Last update: 2024-01-20 00:58:37+09:00
• Include: #include "cpp/fps.hpp"

Depends on

• 四則演算において自動で mod を取るクラス (cpp/modint.hpp)
• a^(-1) mod MODを返す (cpp/number-theory.hpp)

Verified with

• test/yosupo-exp-of-fps.test.cpp
• test/yosupo-inv-of-fps.test.cpp
• test/yosupo-log-of-fps.test.cpp
• test/yosupo-partition-function_1.test.cpp
• test/yosupo-partition-function_2.test.cpp
• test/yosupo-polynomial-taylor-shift.test.cpp
• test/yosupo-pow-of-fps.test.cpp
• test/yosupo-sharp-p-subset-sum.test.cpp

Code

#pragma once

#include <algorithm>
#include "number-theory.hpp"

template <typename mint> struct FPS : std::vector<mint> {
    using std::vector<mint>::vector;
 
    // constructor
    FPS(const std::vector<mint>& r) : std::vector<mint>(r) {}
 
    // core operator
    inline FPS pre(int siz) const {
        return FPS(begin(*this), begin(*this) + std::min((int)this->size(), siz));
    }
    inline FPS rev() const {
        FPS res = *this;
        reverse(begin(res), end(res));
        return res;
    }
    inline FPS& normalize() {
        while (!this->empty() && this->back() == 0) this->pop_back();
        return *this;
    }
 
    // basic operator
    inline FPS operator - () const noexcept {
        FPS res = (*this);
        for (int i = 0; i < (int)res.size(); ++i) res[i] = -res[i];
        return res;
    }
    inline FPS operator + (const mint& v) const { return FPS(*this) += v; }
    inline FPS operator + (const FPS& r) const { return FPS(*this) += r; }
    inline FPS operator - (const mint& v) const { return FPS(*this) -= v; }
    inline FPS operator - (const FPS& r) const { return FPS(*this) -= r; }
    inline FPS operator * (const mint& v) const { return FPS(*this) *= v; }
    inline FPS operator * (const FPS& r) const { return FPS(*this) *= r; }
    inline FPS operator / (const mint& v) const { return FPS(*this) /= v; }
    inline FPS operator << (int x) const { return FPS(*this) <<= x; }
    inline FPS operator >> (int x) const { return FPS(*this) >>= x; }
    inline FPS& operator += (const mint& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] += v;
        return *this;
    }
    inline FPS& operator += (const FPS& r) {
        if (r.size() > this->size()) this->resize(r.size());
        for (int i = 0; i < (int)r.size(); ++i) (*this)[i] += r[i];
        return this->normalize();
    }
    inline FPS& operator -= (const mint& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] -= v;
        return *this;
    }
    inline FPS& operator -= (const FPS& r) {
        if (r.size() > this->size()) this->resize(r.size());
        for (int i = 0; i < (int)r.size(); ++i) (*this)[i] -= r[i];
        return this->normalize();
    }
    inline FPS& operator *= (const mint& v) {
        for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= v;
        return *this;
    }
    inline FPS& operator *= (const FPS& r) {
        return *this = convolution((*this), r);
    }
    inline FPS& operator /= (const mint& v) {
        assert(v != 0);
        mint iv = v.inv();
        for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= iv;
        return *this;
    }
    inline FPS& operator <<= (int x) {
        FPS res(x, 0);
        res.insert(res.end(), begin(*this), end(*this));
        return *this = res;
    }
    inline FPS& operator >>= (int x) {
        if((int) this->size() <= x) return *this = FPS<mint> (1, 0);
        FPS res;
        res.insert(res.end(), begin(*this) + x, end(*this));
        return *this = res;
    }
    inline mint eval(const mint& v){
        mint res = 0;
        for (int i = (int)this->size()-1; i >= 0; --i) {
            res *= v;
            res += (*this)[i];
        }
        return res;
    }
    inline friend FPS gcd(const FPS& f, const FPS& g) {
        if (g.empty()) return f;
        return gcd(g, f % g);
    }
 
    // advanced operation
    // df/dx
    inline friend FPS diff(const FPS& f) {
        int n = (int)f.size();
        FPS res(n-1);
        for (int i = 1; i < n; ++i) res[i-1] = f[i] * i;
        return res;
    }
 
    // \int f dx
    inline friend FPS intg(const FPS& f) {
        int n = (int)f.size();
        FPS res(n+1, 0);
        for (int i = 0; i < n; ++i) res[i+1] = f[i] / (i+1);
        return res;
    }
 
    // inv(f), f[0] must not be 0
    inline friend FPS inv(const FPS& f, int deg) {
        assert(f[0] != 0);
        if (deg < 0) deg = (int)f.size();
        FPS res({mint(1) / f[0]});
        for (int i = 1; i < deg; i <<= 1) {
            res = (res + res - res * res * f.pre(i << 1)).pre(i << 1);
        }
        res.resize(deg);
        return res;
    }
    inline friend FPS inv(const FPS& f) {
        return inv(f, f.size());
    }
 
    // division, r must be normalized (r.back() must not be 0)
    inline FPS& operator /= (const FPS& r) {
        assert(!r.empty());
        assert(r.back() != 0);
        this->normalize();
        if (this->size() < r.size()) {
            this->clear();
            return *this;
        }
        int need = (int)this->size() - (int)r.size() + 1;
        *this = ((*this).rev().pre(need) * inv(r.rev(), need)).pre(need).rev();
        return *this;
    }
    inline FPS& operator %= (const FPS &r) {
        assert(!r.empty());
        assert(r.back() != 0);
        this->normalize();
        FPS q = (*this) / r;
        return *this -= q * r;
    }
    inline FPS operator / (const FPS& r) const { return FPS(*this) /= r; }
    inline FPS operator % (const FPS& r) const { return FPS(*this) %= r; }
 
    // log(f) = \int f'/f dx, f[0] must be 1
    inline friend FPS log(const FPS& f, int deg) {
        assert(f[0] == 1);
        FPS res = intg(diff(f) * inv(f, deg));
        res.resize(deg);
        return res;
    }
    inline friend FPS log(const FPS& f) {
        return log(f, f.size());
    }
 
    // exp(f), f[0] must be 0
    inline friend FPS exp(const FPS& f, int deg) {
        assert(f[0] == 0);
        FPS res(1, 1);
        for (int i = 1; i < deg; i <<= 1) {
            res = res * (f.pre(i<<1) - log(res, i<<1) + 1).pre(i<<1);
        }
        res.resize(deg);
        return res;
    }
    inline friend FPS exp(const FPS& f) {
        return exp(f, f.size());
    }
 
    // pow(f) = exp(e * log f)
    inline friend FPS pow(const FPS& f, long long e, int deg) {
        if(e == 0) {
            auto ret = FPS(deg, 0);
            ret[0] = 1;
            return ret;
        }
        long long i = 0;
        while (i < (int)f.size() && f[i] == 0) ++i;
        if (i == (int)f.size()) return FPS(deg, 0);
        if ((i >= 1 and e >= deg) or i * e >= deg) return FPS(deg, 0);
        mint k = f[i];
        FPS res = exp(log((f >> i) / k, deg) * e, deg) * k.pow(e) << (e * i);
        res.resize(deg);
        return res;
    }
    inline friend FPS pow(const FPS& f, long long e) {
        return pow(f, e, f.size());
    }

    inline friend FPS taylor_shift(FPS f, mint a) {
        int n = f.size();
        std::vector<mint> fac(n, 1), inv(n, 1), finv(n, 1);
        int mod = mint::mod();
        for(int i = 2; i < n; i ++) {
            fac[i] = fac[i - 1] * i;
            inv[i] = -inv[mod % i] * (mod / i);
            finv[i] = finv[i - 1] * inv[i];
        }
        for(int i = 0; i < n; i ++) f[i] *= fac[i];
        std::reverse(f.begin(), f.end());
        FPS<mint> g(n, 1);
        for(int i = 1; i < n; i ++) g[i] = g[i - 1] * a * inv[i];
        f = (f * g).pre(n);
        std::reverse(f.begin(), f.end());
        for(int i = 0; i < n; i ++) f[i] *= finv[i];
        return f;
    }
};

template <typename mint> FPS<mint> modpow(const FPS<mint> &f, long long n, const FPS<mint> &m) {
    if (n == 0) return FPS<mint>(1, 1);
    auto t = modpow(f, n / 2, m);
    t = (t * t) % m;
    if (n & 1) t = (t * f) % m;
    return t;
}

#line 2 "cpp/fps.hpp"

#include <algorithm>
#line 2 "cpp/number-theory.hpp"

#include <numeric>
#include <vector>
#line 2 "cpp/modint.hpp"

/**
 * @file modint.hpp
 * @brief 四則演算において自動で mod を取るクラス
 */

#include <iostream>
#include <utility>
#include <limits>
#include <type_traits>
#include <cstdint>
#include <cassert>

namespace detail {
    static constexpr std::uint16_t prime32_bases[] {
        15591,  2018,  166, 7429,  8064, 16045, 10503,  4399,  1949,  1295, 2776,  3620,   560,  3128,  5212,  2657,
         2300,  2021, 4652, 1471,  9336,  4018,  2398, 20462, 10277,  8028, 2213,  6219,   620,  3763,  4852,  5012,
         3185,  1333, 6227, 5298,  1074,  2391,  5113,  7061,   803,  1269, 3875,   422,   751,   580,  4729, 10239,
          746,  2951,  556, 2206,  3778,   481,  1522,  3476,   481,  2487, 3266,  5633,   488,  3373,  6441,  3344,
           17, 15105, 1490, 4154,  2036,  1882,  1813,   467,  3307, 14042, 6371,   658,  1005,   903,   737,  1887,
         7447,  1888, 2848, 1784,  7559,  3400,   951, 13969,  4304,   177,   41, 19875,  3110, 13221,  8726,   571,
         7043,  6943, 1199,  352,  6435,   165,  1169,  3315,   978,   233, 3003,  2562,  2994, 10587, 10030,  2377,
         1902,  5354, 4447, 1555,   263, 27027,  2283,   305,   669,  1912,  601,  6186,   429,  1930, 14873,  1784,
         1661,   524, 3577,  236,  2360,  6146,  2850, 55637,  1753,  4178, 8466,   222,  2579,  2743,  2031,  2226,
         2276,   374, 2132,  813, 23788,  1610,  4422,  5159,  1725,  3597, 3366, 14336,   579,   165,  1375, 10018,
        12616,  9816, 1371,  536,  1867, 10864,   857,  2206,  5788,   434, 8085, 17618,   727,  3639,  1595,  4944,
         2129,  2029, 8195, 8344,  6232,  9183,  8126,  1870,  3296,  7455, 8947, 25017,   541, 19115,   368,   566,
         5674,   411,  522, 1027,  8215,  2050,  6544, 10049,   614,   774, 2333,  3007, 35201,  4706,  1152,  1785,
         1028,  1540, 3743,  493,  4474,  2521, 26845,  8354,   864, 18915, 5465,  2447,    42,  4511,  1660,   166,
         1249,  6259, 2553,  304,   272,  7286,    73,  6554,   899,  2816, 5197, 13330,  7054,  2818,  3199,   811,
          922,   350, 7514, 4452,  3449,  2663,  4708,   418,  1621,  1171, 3471,    88, 11345,   412,  1559,   194,
    };

    static constexpr bool is_SPRP(std::uint32_t n, std::uint32_t a) noexcept {
        std::uint32_t d = n - 1;
        std::uint32_t s = 0;
        while ((d & 1) == 0) {
            ++s;
            d >>= 1;
        }
        std::uint64_t cur = 1;
        std::uint64_t pw = d;
        while (pw) {
            if (pw & 1) cur = (cur * a) % n;
            a = (static_cast<std::uint64_t>(a) * a) % n;
            pw >>= 1;
        }
        if (cur == 1) return true;
        for (std::uint32_t r = 0; r < s; ++r) {
            if (cur == n - 1) return true;
            cur = (cur * cur) % n;
        }
        return false;
    }

    // 32ビット符号なし整数の素数判定
    // 参考: M. Forisek and J. Jancina, “Fast Primality Testing for Integers That Fit into a Machine Word,” presented at the Conference on Current Trends in Theory and Practice of Informatics, 2015.
    [[nodiscard]]
    static constexpr bool is_prime32(std::uint32_t x) noexcept {
        if (x == 2 || x == 3 || x == 5 || x == 7) return true;
        if (x % 2 == 0 || x % 3 == 0 || x % 5 == 0 || x % 7 == 0) return false;
        if (x < 121) return (x > 1);
        std::uint64_t h = x;
        h = ((h >> 16) ^ h) * 0x45d9f3b;
        h = ((h >> 16) ^ h) * 0x45d9f3b;
        h = ((h >> 16) ^ h) & 0xff;
        return is_SPRP(x, prime32_bases[h]);
    }
}

/// @brief static_modint と dynamic_modint の実装を CRTP によって行うためのクラステンプレート
/// @tparam Modint このクラステンプレートを継承するクラス
template <class Modint>
class modint_base {
public:
    /// @brief 保持する値の型
    using value_type = std::uint32_t;

    /// @brief 0 で初期化します。
    constexpr modint_base() noexcept
        : m_value{ 0 } {}

    /// @brief @c value の剰余で初期化します。
    /// @param value 初期化に使う値
    template <class SignedIntegral, std::enable_if_t<std::is_integral_v<SignedIntegral> && std::is_signed_v<SignedIntegral>>* = nullptr>
    constexpr modint_base(SignedIntegral value) noexcept
        : m_value{ static_cast<value_type>((static_cast<long long>(value) % Modint::mod() + Modint::mod()) % Modint::mod()) } {}

    /// @brief @c value の剰余で初期化します。
    /// @param value 初期化に使う値
    template <class UnsignedIntegral, std::enable_if_t<std::is_integral_v<UnsignedIntegral> && std::is_unsigned_v<UnsignedIntegral>>* = nullptr>
    constexpr modint_base(UnsignedIntegral value) noexcept
        : m_value{ static_cast<value_type>(value % Modint::mod()) } {}

    /// @brief 保持している値を取得します。
    /// @return 保持している値
    [[nodiscard]]
    constexpr value_type value() const noexcept {
        return m_value;
    }

    /// @brief 保持している値をインクリメントして、剰余を取ります。
    /// @return @c *this
    constexpr Modint& operator++() noexcept {
        ++m_value;
        if (m_value == Modint::mod()) {
            m_value = 0;
        }
        return static_cast<Modint&>(*this);
    }

    /// @brief 保持している値をインクリメントして、剰余を取ります。
    /// @return @c *this
    constexpr Modint operator++(int) noexcept {
        auto x = static_cast<const Modint&>(*this);
        ++*this;
        return x;
    }

    /// @brief 保持している値をデクリメントして、剰余を取ります。
    /// @return @c *this
    constexpr Modint& operator--() noexcept {
        if (m_value == 0) {
            m_value = Modint::mod();
        }
        --m_value;
        return static_cast<Modint&>(*this);
    }

    /// @brief 保持している値をデクリメントして、剰余を取ります。
    /// @return @c *this
    constexpr Modint operator--(int) noexcept {
        auto x = static_cast<const Modint&>(*this);
        --*this;
        return x;
    }

    /// @brief 保持している値に @c x の持つ値を足して、剰余を取ります。
    /// @param x 足す数
    /// @return @c *this
    constexpr Modint& operator+=(const Modint& x) noexcept {
        m_value += x.m_value;
        if (m_value >= Modint::mod()) {
            m_value -= Modint::mod();
        }
        return static_cast<Modint&>(*this);
    }

    /// @brief 保持している値から @c x の持つ値を引いて、剰余を取ります。
    /// @param x 引く数
    /// @return @c *this
    constexpr Modint& operator-=(const Modint& x) noexcept {
        m_value -= x.m_value;
        if (m_value >= Modint::mod()) {
            m_value += Modint::mod();
        }
        return static_cast<Modint&>(*this);
    }

    /// @brief 保持している値に @c x の持つ値を掛けて、剰余を取ります。
    /// @param x 掛ける数
    /// @return @c *this
    constexpr Modint& operator*=(const Modint& x) noexcept {
        m_value = static_cast<value_type>(static_cast<std::uint64_t>(m_value) * x.m_value % Modint::mod());
        return static_cast<Modint&>(*this);
    }

    /// @brief 保持している値を @c x の持つ値で割って、剰余を取ります。
    /// @remark 時間計算量： @f$O(\log x)@f$
    /// @param x 割る数
    /// @return @c *this
    constexpr Modint& operator/=(const Modint& x) noexcept {
        return *this *= x.inv();
    }

    /// @brief 自身のコピーを返します。
    /// @return @c *this
    [[nodiscard]]
    constexpr Modint operator+() const noexcept {
        return static_cast<const Modint&>(*this);
    }

    /// @brief 自身の反数を返します。
    /// @return 自身の反数
    [[nodiscard]]
    constexpr Modint operator-() const noexcept {
        return 0 - static_cast<const Modint&>(*this);
    }

    /// @brief 自身の @c n 乗を返します。
    /// @remark 時間計算量： @f$O(\log n)@f$
    /// @param n 指数
    /// @return 自身の @c n 乗
    [[nodiscard]]
    constexpr Modint pow(unsigned long long n) const noexcept {
        Modint x = 1;
        Modint y = static_cast<const Modint&>(*this);
        while (n) {
            if (n & 1) {
                x *= y;
            }
            y *= y;
            n >>= 1;
        }
        return x;
    }

    /// @brief 自身の逆数を返します。
    /// @remark 時間計算量： @f$O(\log value)@f$
    /// @return 自身の逆数
    [[nodiscard]]
    constexpr Modint inv() const noexcept {
        long long a = Modint::mod();
        long long b = m_value;
        long long x = 0;
        long long y = 1;
        while (b) {
            auto t = a / b;
            auto u = a - t * b;
            a = b;
            b = u;
            u = x - t * y;
            x = y;
            y = u;
        }
        assert(a == 1 && "The inverse element does not exist.");
        x %= Modint::mod();
        if (x < 0) {
            x += Modint::mod();
        }
        return x;
    }

    /// @brief @c x に @c y を足したオブジェクトを返します。
    /// @param x 足される数
    /// @param y 足す数
    /// @return @c x に @c y を足したオブジェクト
    [[nodiscard]]
    friend constexpr Modint operator+(const Modint& x, const Modint& y) noexcept {
        return std::move(Modint{ x } += y);
    }

    /// @brief @c x から @c y を引いたオブジェクトを返します。
    /// @param x 引かれる数
    /// @param y 引く数
    /// @return @c x から @c y を引いたオブジェクト
    [[nodiscard]]
    friend constexpr Modint operator-(const Modint& x, const Modint& y) noexcept {
        return std::move(Modint{ x } -= y);
    }

    /// @brief @c x に @c y を掛けたオブジェクトを返します。
    /// @param x 掛けられる数
    /// @param y 掛ける数
    /// @return @c x に @c y を掛けたオブジェクト
    [[nodiscard]]
    friend constexpr Modint operator*(const Modint& x, const Modint& y) noexcept {
        return std::move(Modint{ x } *= y);
    }

    /// @brief @c x を @c y で割ったオブジェクトを返します。
    /// @param x 割られる数
    /// @param y 割る数
    /// @return @c x を @c y で割ったオブジェクト
    [[nodiscard]]
    friend constexpr Modint operator/(const Modint& x, const Modint& y) noexcept {
        return std::move(Modint{ x } /= y);
    }

    /// @brief @c x と @c y の保持する値が等しいかどうかを調べます。
    /// @return @c x と @c y の保持する値が等しければ @c true 、そうでなければ @c false
    [[nodiscard]]
    friend constexpr bool operator==(const Modint& x, const Modint& y) noexcept {
        return x.m_value == y.m_value;
    }

    /// @brief @c x と @c y の保持する値が等しくないかどうかを調べます。
    /// @return @c x と @c y の保持する値が等しければ @c false 、そうでなければ @c true
    [[nodiscard]]
    friend constexpr bool operator!=(const Modint& x, const Modint& y) noexcept {
        return not (x == y);
    }

    /// @brief 入力ストリームから符号付き整数を読み取り、 @c x に格納します。
    /// @tparam CharT 入力ストリームの文字型
    /// @tparam Traits 入力ストリームの文字トレイト
    /// @param is 入力ストリーム
    /// @param x 入力を受け取るオブジェクト
    /// @return @c is
    template <class CharT, class Traits>
    friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& is, Modint& x) {
        long long tmp;
        is >> tmp;
        x = tmp;
        return is;
    }

    /// @brief 出力ストリームに @c x の保持する値を出力します。
    /// @tparam CharT 出力ストリームの文字型
    /// @tparam Traits 出力ストリームの文字トレイト
    /// @param os 出力ストリーム
    /// @param x 出力するオブジェクト
    /// @return @c os
    template <class CharT, class Traits>
    friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const Modint& x) {
        os << x.value();
        return os;
    }

protected:
    value_type m_value;
};

/// @brief コンパイル時に法が決まるとき、四則演算において自動で mod を取るクラス
/// @tparam Mod 法
template <std::uint32_t Mod>
class static_modint : public modint_base<static_modint<Mod>> {
    static_assert(Mod > 0 && Mod <= std::numeric_limits<std::uint32_t>::max() / 2);

private:
    using base_type = modint_base<static_modint<Mod>>;

public:
    using typename base_type::value_type;

    /// @brief 法を取得します。
    /// @return 法
    [[nodiscard]]
    static constexpr value_type mod() noexcept {
        return Mod;
    }

    /// @brief 0 で初期化します。
    constexpr static_modint() noexcept
        : base_type{} {}

    /// @brief @c value の剰余で初期化します。
    /// @param value 初期化に使う値
    template <class SignedIntegral, std::enable_if_t<std::is_integral_v<SignedIntegral>>* = nullptr>
    constexpr static_modint(SignedIntegral value) noexcept
        : base_type{value} {}

    /// @brief 自身の逆数を返します。
    /// @remark 時間計算量： @f$O(\log value)@f$
    /// @return 自身の逆数
    [[nodiscard]]
    constexpr static_modint inv() const noexcept {
        if constexpr (detail::is_prime32(Mod)) {
            assert(this->m_value != 0 && "The inverse element of zero does not exist.");
            return this->pow(Mod - 2);
        }
        else {
            return base_type::inv();
        }
    }
};

/// @brief 実行時に法が決まるとき、四則演算において自動で mod を取るクラス
/// @tparam ID このIDごとに法を設定することができます
template <int ID>
class dynamic_modint : public modint_base<dynamic_modint<ID>> {
private:
    using base_type = modint_base<dynamic_modint<ID>>;

public:
    using typename base_type::value_type;

    /// @brief 法を取得します。
    /// @return 法
    [[nodiscard]]
    static value_type mod() noexcept {
        return modulus;
    }

    /// @brief 法を設定します。
    /// @param m 新しい法
    static void set_mod(value_type m) noexcept {
        assert(m > 0 && m <= std::numeric_limits<value_type>::max() / 2);
        modulus = m;
    }

    /// @brief 0 で初期化します。
    constexpr dynamic_modint() noexcept
        : base_type{} {}

    /// @brief @c value の剰余で初期化します。
    /// @param value 初期化に使う値
    template <class SignedIntegral, std::enable_if_t<std::is_integral_v<SignedIntegral>>* = nullptr>
    constexpr dynamic_modint(SignedIntegral value) noexcept
        : base_type{value} {}

private:
    inline static value_type modulus = 998244353;
};

using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
using modint = dynamic_modint<-1>;
#line 6 "cpp/number-theory.hpp"

/**
 * @brief a^(-1) mod MODを返す
 * @param a long long
 * @param MOD long long
 * @return long long
 */
long long modinv(long long a, long long MOD) {
    long long b = MOD, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; std::swap(a, b);
        u -= t * v; std::swap(u, v);
    }
    u %= MOD; 
    if (u < 0) u += MOD;
    return u;
}

/**
 * @brief a^n mod MODを返す
 * @param a long long
 * @param n long long
 * @param MOD long long
 * @return long long
 */
long long modpow(long long a, long long n, long long MOD) {
    long long res = 1;
    a %= MOD;
    if(n < 0) {
        n = -n;
        a = modinv(a, MOD);
    }
    while (n > 0) {
        if (n & 1) res = res * a % MOD;
        a = a * a % MOD;
        n >>= 1;
    }
    return res;
}

/**
 * @brief 2式の連立合同式を、mが互いに素になるように変形する
 * @param r1 long long
 * @param m1 long long
 * @param r2 long long
 * @param m2 long long
 * @note 矛盾する場合、r1 = r2 = m1 = m2 = -1となる
 */
void coprimize_simulaneous_congruence_equation(long long& r1, long long& m1, long long& r2, long long& m2) {
    long long g = std::gcd(m1, m2);
    if((r2 - r1) % g != 0) {
        r1 = r2 = m1 = m2 = -1;
        return;
    }
    m1 /= g, m2 /= g;
    long long gi = std::gcd(g, m1);
    long long gj = g / gi;
    do {
        g = std::gcd(gi, gj);
        gi *= g, gj /= g;
    } while(g != 1);
    m1 *= gi, m2 *= gj;
    r1 %= m1, r2 %= m2;
}

/**
 * @brief 連立合同式を解く
 * @param r vector<long long> 余りの配列
 * @param m vector<long long> modの配列
 * @return std::pair<long long, long long> (解, LCM) 解なしのときは{-1, -1}
 */
std::pair<long long, long long> crt(const std::vector<long long>& r, const std::vector<long long>& m) {
    assert(r.size() == m.size());
    if(r.size() == 0) return {0, 1};
    int n = (int)r.size();
    long long m_lcm = m[0];
    long long ans = r[0] % m[0];
    for (int i = 1; i < n; i++) {
        long long rr = r[i] % m[i], mm = m[i];
        coprimize_simulaneous_congruence_equation(ans, m_lcm, rr, mm);
        if(m_lcm == -1) return {-1, -1};
        long long t = ((rr - ans) * modinv(m_lcm, mm)) % mm;
        if(t < 0) t += mm;
        ans += t * m_lcm;
        m_lcm *= mm;
    }
    return {ans, m_lcm};
}

/**
 * @brief 連立合同式の最小の非負整数解 % MODを求める
 * @param r vector<long long> 余りの配列
 * @param m vector<long long> modの配列
 * @param MOD long long
 * @return std::pair<long long, long long> (最小解 % MOD, LCM % MOD) 解なしのときは{-1, -1}
 */
std::pair<long long, long long> crt(const std::vector<long long>& r, const std::vector<long long>& m, long long MOD) {
    assert(r.size() == m.size());
    if(r.size() == 0) return {0, 1};
    int n = (int)r.size();
    std::vector<long long> r2 = r, m2 = m;
    // mを互いに素にする
    for(int i = 1; i < n; i++) {
        for(int j = 0; j < i; j++) {
            coprimize_simulaneous_congruence_equation(r2[i], m2[i], r2[j], m2[j]);
            if(m2[i] == -1) return {-1, -1};
        }
    }

    m2.push_back(MOD);
    std::vector<long long> prod(n+1, 1); // m2[0] * ... * m2[i - 1] mod m2[i]
    std::vector<long long> x(n+1, 0); // i番目までの解 mod m2[i]
    for(int i = 0; i < n; i++) {
        long long t = (r2[i] - x[i]) * modinv(prod[i], m2[i]) % m2[i];
        if(t < 0) t += m2[i];
        for(int j = i + 1; j <= n; j++) {
            (x[j] += t * prod[j]) %= m2[j];
            (prod[j] *= m2[i]) %= m2[j];
        }
    }
    return {x[n], prod[n]};
}

/**
 * @brief 畳み込み
 */
namespace NTT {
    /**
     * @brief 原子根
     * @param MOD int
     * @return int
     */
    int calc_primitive_root(int MOD) {
        if (MOD == 2) return 1;
        if (MOD == 167772161) return 3;
        if (MOD == 469762049) return 3;
        if (MOD == 754974721) return 11;
        if (MOD == 998244353) return 3;
        int divs[20] = {};
        divs[0] = 2;
        int cnt = 1;
        long long x = (MOD - 1) >> 1;
        while (x % 2 == 0) x >>= 1;
        for (long long i = 3; i * i <= x; i += 2) {
            if (x % i == 0) {
                divs[cnt ++] = i;
                while (x % i == 0) x /= i;
            }
        }
        if (x > 1) divs[cnt++] = x;
        for (int g = 2;; ++ g) {
            bool ok = true;
            for (int i = 0; i < cnt; i++) {
                if (modpow(g, (MOD - 1) / divs[i], MOD) == 1) {
                    ok = false;
                    break;
                }
            }
            if (ok) return g;
        }
    }

    /**
     * @brief 畳み込みのサイズを2のべき乗にする
     */
    int get_fft_size(int N, int M) {
        int size_a = 1, size_b = 1;
        while (size_a < N) size_a <<= 1;
        while (size_b < M) size_b <<= 1;
        return std::max(size_a, size_b) << 1;
    }

    /**
     * @brief NTT
     */
    template<class mint> void trans(std::vector<mint>& v, bool inv = false) {
        if (v.empty()) return;
        int N = (int) v.size();
        int MOD = v[0].mod();
        int PR = calc_primitive_root(MOD);
        static bool first = true;
        static std::vector<long long> vbw(30), vibw(30);
        if (first) {
            first = false;
            for (int k = 0; k < 30; ++ k) {
                vbw[k] = modpow(PR, (MOD - 1) >> (k + 1), MOD);
                vibw[k] = modinv(vbw[k], MOD);
            }
        }
        for (int i = 0, j = 1; j < N - 1; ++ j) {
            for (int k = N >> 1; k > (i ^= k); k >>= 1);
            if (i > j) std::swap(v[i], v[j]);
        }
        for (int k = 0, t = 2; t <= N; ++ k, t <<= 1) {
            long long bw = vbw[k];
            if (inv) bw = vibw[k];
            for (int i = 0; i < N; i += t) {
                mint w = 1;
                for (int j = 0; j < (t >> 1); ++ j) {
                    int j1 = i + j, j2 = i + j + (t >> 1);
                    mint c1 = v[j1], c2 = v[j2] * w;
                    v[j1] = c1 + c2;
                    v[j2] = c1 - c2;
                    w *= bw;
                }
            }
        }
        if (inv) {
            long long invN = modinv(N, MOD);
            for (int i = 0; i < N; ++ i) v[i] = v[i] * invN;
        }
    }

    static constexpr int MOD0 = 754974721;
    static constexpr int MOD1 = 167772161;
    static constexpr int MOD2 = 469762049;
    using mint0 = static_modint<MOD0>;
    using mint1 = static_modint<MOD1>;
    using mint2 = static_modint<MOD2>;
    static const mint1 imod0 = 95869806; // modinv(MOD0, MOD1);
    static const mint2 imod1 = 104391568; // modinv(MOD1, MOD2);
    static const mint2 imod01 = 187290749; // imod1 / MOD0;

    /**
     * @brief 配列のサイズが小さいときの畳み込み
     * @param T mint, long long
     * @param A vector<T>
     * @param B vector<T>
     * @return vector<T>
     */
    template<class T> std::vector<T> naive
    (const std::vector<T>& A, const std::vector<T>& B) {
        if (A.empty() || B.empty()) return {};
        int N = (int) A.size(), M = (int) B.size();
        std::vector<T> res(N + M - 1);
        for (int i = 0; i < N; ++ i)
            for (int j = 0; j < M; ++ j)
                res[i + j] += A[i] * B[j];
        return res;
    }
};

/**
 * @brief modintの畳み込み
 * @param A vector<mint>
 * @param B vector<mint>
 * @return vector<mint>
 */
template<class mint> std::vector<mint> convolution
(const std::vector<mint>& A, const std::vector<mint>& B) {
    if (A.empty() || B.empty()) return {};
    int N = (int) A.size(), M = (int) B.size();
    if (std::min(N, M) < 30) return NTT::naive(A, B);
    int MOD = A[0].mod();
    int size_fft = NTT::get_fft_size(N, M);
    if (MOD == 998244353) {
        std::vector<mint> a(size_fft), b(size_fft), c(size_fft);
        for (int i = 0; i < N; ++i) a[i] = A[i];
        for (int i = 0; i < M; ++i) b[i] = B[i];
        NTT::trans(a), NTT::trans(b);
        std::vector<mint> res(size_fft);
        for (int i = 0; i < size_fft; ++i) res[i] = a[i] * b[i];
        NTT::trans(res, true);
        res.resize(N + M - 1);
        return res;
    }
    std::vector<NTT::mint0> a0(size_fft, 0), b0(size_fft, 0), c0(size_fft, 0);
    std::vector<NTT::mint1> a1(size_fft, 0), b1(size_fft, 0), c1(size_fft, 0);
    std::vector<NTT::mint2> a2(size_fft, 0), b2(size_fft, 0), c2(size_fft, 0);
    for (int i = 0; i < N; ++ i) {
        a0[i] = A[i].value();
        a1[i] = A[i].value();
        a2[i] = A[i].value();
    }
    for (int i = 0; i < M; ++ i) {
        b0[i] = B[i].value();
        b1[i] = B[i].value();
        b2[i] = B[i].value();
    }
    NTT::trans(a0), NTT::trans(a1), NTT::trans(a2), 
    NTT::trans(b0), NTT::trans(b1), NTT::trans(b2);
    for (int i = 0; i < size_fft; ++i) {
        c0[i] = a0[i] * b0[i];
        c1[i] = a1[i] * b1[i];
        c2[i] = a2[i] * b2[i];
    }
    NTT::trans(c0, true), NTT::trans(c1, true), NTT::trans(c2, true);
    static const mint mod0 = NTT::MOD0, mod01 = mod0 * NTT::MOD1;
    std::vector<mint> res(N + M - 1);
    for (int i = 0; i < N + M - 1; ++ i) {
        int y0 = c0[i].value();
        int y1 = (NTT::imod0 * (c1[i] - y0)).value();
        int y2 = (NTT::imod01 * (c2[i] - y0) - NTT::imod1 * y1).value();
        res[i] = mod01 * y2 + mod0 * y1 + y0;
    }
    return res;
}

/**
 * @brief long longの畳み込み
 * @param A vector<long long>
 * @param B vector<long long>
 * @return vector<long long>
 */
std::vector<long long> convolution_ll
(const std::vector<long long>& A, const std::vector<long long>& B) {
    if (A.empty() || B.empty()) return {};
    int N = (int) A.size(), M = (int) B.size();
    if (std::min(N, M) < 30) return NTT::naive(A, B);
    int size_fft = NTT::get_fft_size(N, M);
    std::vector<NTT::mint0> a0(size_fft, 0), b0(size_fft, 0), c0(size_fft, 0);
    std::vector<NTT::mint1> a1(size_fft, 0), b1(size_fft, 0), c1(size_fft, 0);
    std::vector<NTT::mint2> a2(size_fft, 0), b2(size_fft, 0), c2(size_fft, 0);
    for (int i = 0; i < N; ++ i) {
        a0[i] = A[i];
        a1[i] = A[i];
        a2[i] = A[i];
    }
    for (int i = 0; i < M; ++ i) {
        b0[i] = B[i];
        b1[i] = B[i];
        b2[i] = B[i];
    }
    NTT::trans(a0), NTT::trans(a1), NTT::trans(a2), 
    NTT::trans(b0), NTT::trans(b1), NTT::trans(b2);
    for (int i = 0; i < size_fft; ++ i) {
        c0[i] = a0[i] * b0[i];
        c1[i] = a1[i] * b1[i];
        c2[i] = a2[i] * b2[i];
    }
    NTT::trans(c0, true), NTT::trans(c1, true), NTT::trans(c2, true);
    static const long long mod0 = NTT::MOD0, mod01 = mod0 * NTT::MOD1;
    static const __int128_t mod012 = (__int128_t)mod01 * NTT::MOD2;
    std::vector<long long> res(N + M - 1);
    for (int i = 0; i < N + M - 1; ++ i) {
        int y0 = c0[i].value();
        int y1 = (NTT::imod0 * (c1[i] - y0)).value();
        int y2 = (NTT::imod01 * (c2[i] - y0) - NTT::imod1 * y1).value();
        __int128_t tmp = (__int128_t)mod01 * y2 + (__int128_t)mod0 * y1 + y0;
        if(tmp < (mod012 >> 1)) res[i] = tmp;
        else res[i] = tmp - mod012;
    }
    return res;
}
#line 5 "cpp/fps.hpp"

template <typename mint> struct FPS : std::vector<mint> {
    using std::vector<mint>::vector;
 
    // constructor
    FPS(const std::vector<mint>& r) : std::vector<mint>(r) {}
 
    // core operator
    inline FPS pre(int siz) const {
        return FPS(begin(*this), begin(*this) + std::min((int)this->size(), siz));
    }
    inline FPS rev() const {
        FPS res = *this;
        reverse(begin(res), end(res));
        return res;
    }
    inline FPS& normalize() {
        while (!this->empty() && this->back() == 0) this->pop_back();
        return *this;
    }
 
    // basic operator
    inline FPS operator - () const noexcept {
        FPS res = (*this);
        for (int i = 0; i < (int)res.size(); ++i) res[i] = -res[i];
        return res;
    }
    inline FPS operator + (const mint& v) const { return FPS(*this) += v; }
    inline FPS operator + (const FPS& r) const { return FPS(*this) += r; }
    inline FPS operator - (const mint& v) const { return FPS(*this) -= v; }
    inline FPS operator - (const FPS& r) const { return FPS(*this) -= r; }
    inline FPS operator * (const mint& v) const { return FPS(*this) *= v; }
    inline FPS operator * (const FPS& r) const { return FPS(*this) *= r; }
    inline FPS operator / (const mint& v) const { return FPS(*this) /= v; }
    inline FPS operator << (int x) const { return FPS(*this) <<= x; }
    inline FPS operator >> (int x) const { return FPS(*this) >>= x; }
    inline FPS& operator += (const mint& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] += v;
        return *this;
    }
    inline FPS& operator += (const FPS& r) {
        if (r.size() > this->size()) this->resize(r.size());
        for (int i = 0; i < (int)r.size(); ++i) (*this)[i] += r[i];
        return this->normalize();
    }
    inline FPS& operator -= (const mint& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] -= v;
        return *this;
    }
    inline FPS& operator -= (const FPS& r) {
        if (r.size() > this->size()) this->resize(r.size());
        for (int i = 0; i < (int)r.size(); ++i) (*this)[i] -= r[i];
        return this->normalize();
    }
    inline FPS& operator *= (const mint& v) {
        for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= v;
        return *this;
    }
    inline FPS& operator *= (const FPS& r) {
        return *this = convolution((*this), r);
    }
    inline FPS& operator /= (const mint& v) {
        assert(v != 0);
        mint iv = v.inv();
        for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= iv;
        return *this;
    }
    inline FPS& operator <<= (int x) {
        FPS res(x, 0);
        res.insert(res.end(), begin(*this), end(*this));
        return *this = res;
    }
    inline FPS& operator >>= (int x) {
        if((int) this->size() <= x) return *this = FPS<mint> (1, 0);
        FPS res;
        res.insert(res.end(), begin(*this) + x, end(*this));
        return *this = res;
    }
    inline mint eval(const mint& v){
        mint res = 0;
        for (int i = (int)this->size()-1; i >= 0; --i) {
            res *= v;
            res += (*this)[i];
        }
        return res;
    }
    inline friend FPS gcd(const FPS& f, const FPS& g) {
        if (g.empty()) return f;
        return gcd(g, f % g);
    }
 
    // advanced operation
    // df/dx
    inline friend FPS diff(const FPS& f) {
        int n = (int)f.size();
        FPS res(n-1);
        for (int i = 1; i < n; ++i) res[i-1] = f[i] * i;
        return res;
    }
 
    // \int f dx
    inline friend FPS intg(const FPS& f) {
        int n = (int)f.size();
        FPS res(n+1, 0);
        for (int i = 0; i < n; ++i) res[i+1] = f[i] / (i+1);
        return res;
    }
 
    // inv(f), f[0] must not be 0
    inline friend FPS inv(const FPS& f, int deg) {
        assert(f[0] != 0);
        if (deg < 0) deg = (int)f.size();
        FPS res({mint(1) / f[0]});
        for (int i = 1; i < deg; i <<= 1) {
            res = (res + res - res * res * f.pre(i << 1)).pre(i << 1);
        }
        res.resize(deg);
        return res;
    }
    inline friend FPS inv(const FPS& f) {
        return inv(f, f.size());
    }
 
    // division, r must be normalized (r.back() must not be 0)
    inline FPS& operator /= (const FPS& r) {
        assert(!r.empty());
        assert(r.back() != 0);
        this->normalize();
        if (this->size() < r.size()) {
            this->clear();
            return *this;
        }
        int need = (int)this->size() - (int)r.size() + 1;
        *this = ((*this).rev().pre(need) * inv(r.rev(), need)).pre(need).rev();
        return *this;
    }
    inline FPS& operator %= (const FPS &r) {
        assert(!r.empty());
        assert(r.back() != 0);
        this->normalize();
        FPS q = (*this) / r;
        return *this -= q * r;
    }
    inline FPS operator / (const FPS& r) const { return FPS(*this) /= r; }
    inline FPS operator % (const FPS& r) const { return FPS(*this) %= r; }
 
    // log(f) = \int f'/f dx, f[0] must be 1
    inline friend FPS log(const FPS& f, int deg) {
        assert(f[0] == 1);
        FPS res = intg(diff(f) * inv(f, deg));
        res.resize(deg);
        return res;
    }
    inline friend FPS log(const FPS& f) {
        return log(f, f.size());
    }
 
    // exp(f), f[0] must be 0
    inline friend FPS exp(const FPS& f, int deg) {
        assert(f[0] == 0);
        FPS res(1, 1);
        for (int i = 1; i < deg; i <<= 1) {
            res = res * (f.pre(i<<1) - log(res, i<<1) + 1).pre(i<<1);
        }
        res.resize(deg);
        return res;
    }
    inline friend FPS exp(const FPS& f) {
        return exp(f, f.size());
    }
 
    // pow(f) = exp(e * log f)
    inline friend FPS pow(const FPS& f, long long e, int deg) {
        if(e == 0) {
            auto ret = FPS(deg, 0);
            ret[0] = 1;
            return ret;
        }
        long long i = 0;
        while (i < (int)f.size() && f[i] == 0) ++i;
        if (i == (int)f.size()) return FPS(deg, 0);
        if ((i >= 1 and e >= deg) or i * e >= deg) return FPS(deg, 0);
        mint k = f[i];
        FPS res = exp(log((f >> i) / k, deg) * e, deg) * k.pow(e) << (e * i);
        res.resize(deg);
        return res;
    }
    inline friend FPS pow(const FPS& f, long long e) {
        return pow(f, e, f.size());
    }

    inline friend FPS taylor_shift(FPS f, mint a) {
        int n = f.size();
        std::vector<mint> fac(n, 1), inv(n, 1), finv(n, 1);
        int mod = mint::mod();
        for(int i = 2; i < n; i ++) {
            fac[i] = fac[i - 1] * i;
            inv[i] = -inv[mod % i] * (mod / i);
            finv[i] = finv[i - 1] * inv[i];
        }
        for(int i = 0; i < n; i ++) f[i] *= fac[i];
        std::reverse(f.begin(), f.end());
        FPS<mint> g(n, 1);
        for(int i = 1; i < n; i ++) g[i] = g[i - 1] * a * inv[i];
        f = (f * g).pre(n);
        std::reverse(f.begin(), f.end());
        for(int i = 0; i < n; i ++) f[i] *= finv[i];
        return f;
    }
};

template <typename mint> FPS<mint> modpow(const FPS<mint> &f, long long n, const FPS<mint> &m) {
    if (n == 0) return FPS<mint>(1, 1);
    auto t = modpow(f, n / 2, m);
    t = (t * t) % m;
    if (n & 1) t = (t * f) % m;
    return t;
}

Back to top page