math/FPSbase.hpp

• View this file on GitHub
• Last update: 2025-10-21 09:42:15+09:00
• Include: #include "math/FPSbase.hpp"

Depends on

• math/NTT.hpp
• math/convolution.hpp
• math/modint.hpp

Code

#pragma once
#include "convolution"
#include <bits/stdc++.h>
using namespace std;

template <long long FPS_MOD>
struct FPS : vector<modint<FPS_MOD>> {
    using mint = modint<FPS_MOD>;
    using vector<mint>::size;
    using vector<mint>::empty;
    using vector<mint>::resize;
    using vector<mint>::operator[];
    FPS() : vector<mint>() {}
    FPS(size_t n) : vector<mint>(n) {}
    FPS(size_t n, const mint& x) : vector<mint>(n, x) {}
    FPS(const initializer_list<mint>& il) : vector<mint>(il) {}
    template<class It>
    FPS(It first, It last) : vector<mint>(first, last) {}
    FPS(const vector<mint>& v) : vector<mint>(v) {}
    FPS(vector<mint>&& v) : vector<mint>(move(v)) {}
    static vector<mint> to_vec(const FPS& f) {
        return vector<mint>(f.begin(), f.end());
    }
    static FPS from_vec(const vector<mint>& v) { return FPS(v); }
    static FPS from_vec(vector<mint>&& v) { return FPS(move(v)); }
    FPS low(int n) const {
        int take = min(n, (int)this->size());
        return FPS(this->begin(), this->begin() + take);
    }
    FPS& shrink() {
        while (!this->empty() && this->back() == mint(0)) this->pop_back();
        return *this;
    }
    FPS operator-() const {
        FPS res = *this;
        for (auto& x : res) x = -x;
        return res;
    }
    FPS& operator+=(const FPS& x) {
        if (size() < x.size()) resize(x.size());
        for (int i = 0; i < (int)x.size(); i++) (*this)[i] += x[i];
        return *this;
    }
    FPS& operator-=(const FPS& x) {
        if (size() < x.size()) resize(x.size());
        for (int i = 0; i < (int)x.size(); i++) (*this)[i] -= x[i];
        return *this;
    }
    FPS& operator*=(const FPS& x) {
        auto va = to_vec(*this);
        auto vb = to_vec(x);
        *this = from_vec(convolution(va, vb));
        return *this;
    }
    FPS operator*(const FPS& x) const { FPS t = *this; t *= x; return t; }
    FPS& operator*=(const mint& c) {
        for (auto& y : *this) y *= c;
        return *this;
    }
    FPS operator*(const mint& x) const { FPS t = *this; t *= x; return t; }
    FPS& operator/=(const mint& c) {
        mint invc = c.inv();
        for (auto& y : *this) y *= invc;
        return *this;
    }
    FPS operator/(const mint& x) const { FPS t = *this; t /= x; return t; }
    FPS operator+(const FPS& x) const { return FPS(*this) += x; }
    FPS operator-(const FPS& x) const { return FPS(*this) -= x; }
    FPS diff() const {
        if (empty()) return FPS();
        FPS g(size() - 1);
        for (int i = 1; i < (int)size(); i++) g[i - 1] = (*this)[i] * mint(i);
        return g;
    }
    FPS integ() const {
        FPS g(size() + 1);
        static vector<mint> invs = {mint(0), mint(1)};
        int n = g.size();
        for (int i=(int)invs.size(); i<n; i++) invs.push_back(-invs[FPS_MOD % i] * mint(FPS_MOD / i));
        for (int i=0; i<(int)size(); i++) g[i+1] = (*this)[i] * invs[i+1];
        return g;
    }
    FPS inv(int deg = -1) const {
        assert(!this->empty() && (*this)[0] != 0);
        int n = this->size();
        if (deg == -1) deg = n;
        FPS res; res.resize(1); res[0] = (*this)[0].inv();
        for (int m=1; m<deg; m<<=1) {
            int cut = 2 * m;
            FPS f_low = this->low(cut);
            auto vres = to_vec(res);
            auto vres2 = vres;
            auto v_f_low = to_vec(f_low);
            auto t1 = convolution(vres, vres2);
            auto t2 = convolution(t1, v_f_low);
            FPS t = from_vec(move(t2));
            t.resize(cut);
            res.resize(cut);
            for (int i=0; i<cut; i++) {
                mint ti = (i < (int)t.size() ? t[i] : mint(0));
                res[i] = res[i] * mint(2) - ti;
            }
            res = res.low(cut);
        }
        res.resize(deg);
        return res;
    }
    pair<FPS, FPS> divmod(const FPS& g) const {
        FPS f = *this;
        if (f.size() < g.size()) return {FPS{mint(0)}, f};
        int n = f.size(), m = g.size();
        FPS rf = f, rg = g;
        reverse(rf.begin(), rf.end());
        reverse(rg.begin(), rg.end());
        auto vrf = to_vec(rf);
        auto vrg = to_vec(rg);
        FPS tmp = from_vec(convolution(vrf, to_vec(rg.inv(n - m + 1))));
        FPS q = tmp.low(n - m + 1);
        reverse(q.begin(), q.end());
        FPS r = f - g * q;
        r.resize(m - 1);
        return {q, r};
    }
    FPS log(int deg = -1) const {
        assert(!this->empty() && (*this)[0] == mint(1));
        int n = this->size();
        if (deg == -1) deg = n;
        FPS f = diff();
        FPS invf = this->inv(deg);
        auto vf = to_vec(f), vinvf = to_vec(invf);
        FPS g = from_vec(convolution(vf, vinvf)).low(deg - 1).integ();
        g.resize(deg);
        return g;
    }
    FPS exp(int deg = -1) const {
        assert(!this->empty() && (*this)[0] == mint(0));
        int n = this->size();
        if (deg == -1) deg = n;
        FPS g; g.resize(1); g[0] = mint(1);
        for (int m=1; m<deg; m<<=1) {
            int cut = 2 * m;
            FPS g_log = g.log(cut);
            FPS q = this->low(cut) - g_log;
            q[0] += mint(1);
            auto vg = to_vec(g), vq = to_vec(q);
            FPS ng = from_vec(convolution(vg, vq)).low(cut);
            g = move(ng);
        }
        g.resize(deg);
        return g;
    }
    FPS pow(long long k, int n=-1) const {
        int v = 0;
        while(v < (int)this->size() && (*this)[v].val == 0) v++;
        if (v == (int)this->size()) return FPS(n, mint(0));
        if (n == -1) n = this->size();
    }
};

#line 2 "math/modint.hpp"
#include <bits/stdc++.h>
using namespace std;

template <long long modint_MOD>
struct modint {
    long long val;
    constexpr modint(long long x=0) noexcept {
        x %= modint_MOD;
        if (x < 0) x += modint_MOD;
        val = x;
    }
    explicit constexpr operator long long() const noexcept { return val; }
    template<typename Int>
    static constexpr long long adjust(Int x) noexcept {
        static_assert(std::is_integral<Int>::value, "Int must be integral");
        __int128 t = static_cast<__int128>(x);
        __int128 m = static_cast<__int128>(modint_MOD);
        t %= m;
        if (t < 0) t += m;
        return static_cast<long long>(t);
    }
    constexpr modint operator+(const modint& x) const noexcept { return modint(*this) += x; }
    constexpr modint operator-(const modint& x) const noexcept { return modint(*this) -= x; }
    constexpr modint operator*(const modint& x) const noexcept { return modint(*this) *= x; }
    constexpr modint operator/(const modint& x) const noexcept { return modint(*this) /= x; }
    constexpr modint& operator+=(const modint& x) noexcept {
        if (modint_MOD <= (val += x.val)) val -= modint_MOD;
        return *this;
    }
    constexpr modint& operator-=(const modint& x) noexcept {
        if ((val -= x.val) < 0) val += modint_MOD;
        return *this;
    }
    constexpr modint& operator*=(const modint& x) noexcept {
        __int128 t = (__int128)val * x.val;
        val = (long long)(t % modint_MOD);
        if (val < 0) val += modint_MOD;
        return *this;
    }
    constexpr modint& operator/=(modint x) noexcept {
        *this *= x.inv();
        return *this;
    }
    constexpr modint operator-() const noexcept { return modint() - *this; }
    constexpr modint& operator++() noexcept { return *this += 1; }
    constexpr modint& operator--() noexcept { return *this -= 1; }
    constexpr modint operator++(int) noexcept { modint tmp = *this; ++*this; return tmp; }
    constexpr modint operator--(int) noexcept { modint tmp = *this; --*this; return tmp; }
    constexpr bool operator==(const modint& x) const noexcept { return val == x.val; }
    constexpr bool operator!=(const modint& x) const noexcept { return val != x.val; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    constexpr modint& operator+=(Int x) noexcept {
        long long y = adjust<Int>(x);
        if (modint_MOD <= (val += y)) val -= modint_MOD;
        return *this;
    }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    constexpr modint& operator-=(Int x) noexcept {
        long long y = adjust<Int>(x);
        if ((val -= y) < 0) val += modint_MOD;
        return *this;
    }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    constexpr modint& operator*=(Int x) noexcept {
        long long y = adjust<Int>(x);
        __int128 t = (__int128)val * y;
        val = (long long)(t % modint_MOD);
        if (val < 0) val += modint_MOD;
        return *this;
    }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    constexpr modint& operator/=(Int x) noexcept {
        modint t = modint(adjust<Int>(x));
        *this /= t;
        return *this;
    }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator+(const modint& a, Int b) noexcept { modint t = a; t += b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator+(Int a, const modint& b) noexcept { modint t = b; t += a; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator-(const modint& a, Int b) noexcept { modint t = a; t -= b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator-(Int a, const modint& b) noexcept { modint t = modint(adjust<Int>(a)); t -= b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator*(const modint& a, Int b) noexcept { modint t = a; t *= b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator*(Int a, const modint& b) noexcept { modint t = b; t *= a; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator/(const modint& a, Int b) noexcept { modint t = a; t /= b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr modint operator/(Int a, const modint& b) noexcept { modint t = modint(adjust<Int>(a)); t /= b; return t; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr bool operator==(const modint& a, Int b) noexcept { return a.val == adjust<Int>(b); }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr bool operator==(Int a, const modint& b) noexcept { return adjust<Int>(a) == b.val; }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr bool operator!=(const modint& a, Int b) noexcept { return !(a == b); }
    template<typename Int, typename = std::enable_if_t<std::is_integral<Int>::value>>
    friend constexpr bool operator!=(Int a, const modint& b) noexcept { return !(a == b); }
    friend constexpr istream& operator>> (istream& i, modint& x) { long long y; i >> y; x = modint(y); return i; }
    friend constexpr ostream& operator<< (ostream& o, const modint& x) { o << x.val; return o; }
    constexpr modint pow(long long n) const noexcept {
        if (n < 0) return inv().pow(-n);
        modint x = *this, re = 1;
        while (0 < n){
            if (n & 1) re *= x;
            x *= x;
            n >>= 1;
        }
        return re;
    }
    constexpr modint inv() const noexcept { return pow(modint_MOD - 2); }
};

template <long long modint_MOD>
struct factor {
    vector<modint<modint_MOD>> fac, ifac;
    factor(int n) : fac(n+1), ifac(n+1) {
        fac[0] = 1;
        for (int i=1; i<=n; i++) fac[i] = fac[i-1] * i;
        ifac[n] = fac[n].inv();
        for (int i=n; 1<=i; i--) ifac[i-1] = ifac[i] * i;
    }
    modint<modint_MOD> P(int n, int r) const noexcept { return (r<0 || n<r) ? 0 : fac[n] * ifac[n-r]; }
    modint<modint_MOD> C(int n, int r) const noexcept { return (r<0 || n<r) ? 0 : P(n, r) * ifac[r]; }
    modint<modint_MOD> H(int n, int r) const noexcept { return (n==0 && r==0) ? 1 : C(n+r-1, r); }
};
#line 4 "math/NTT.hpp"
using namespace std;

template <long long NTT_MOD>
void NTT(vector<modint<NTT_MOD>>& a, bool inv=false, int g_=3) {
    modint<NTT_MOD> g = g_;
    int n = (int)a.size();
    for (int i=1, j=0; i<n; i++) {
        int bit = n >> 1;
        for (; j&bit; bit>>=1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int len=2; len<=n; len<<=1) {
        modint<NTT_MOD> wlen = g.pow((NTT_MOD - 1) / len);
        if (inv) wlen = wlen.inv();
        for (int i=0; i<n; i+=len) {
            modint<NTT_MOD> w = 1;
            for (int j=0; j<len/2; j++) {
                modint<NTT_MOD> u = a[i+j];
                modint<NTT_MOD> v = a[i+j+len/2]*w;
                a[i+j] = u + v;
                a[i+j+len/2] = u - v;
                w *= wlen;
            }
        }
    }
    if (inv) {
        modint<NTT_MOD> inv_n = modint<NTT_MOD>(n).inv();
        for (auto& x : a) x *= inv_n;
    }
}
#line 4 "math/convolution.hpp"
using namespace std;

template<long long convolution_MOD>
vector<modint<convolution_MOD>> convolution(vector<modint<convolution_MOD>> a,vector<modint<convolution_MOD>> b){
    int n = 1;
    int s = (int)(a.size()+b.size()-1);
    while (n < s) n<<=1;
    a.resize(n);
    b.resize(n);
    NTT(a);
    NTT(b);
    for (int i=0; i<n; i++) a[i] *= b[i];
    NTT(a,true);
    a.resize(s);
    return a;
}
#line 4 "math/FPSbase.hpp"
using namespace std;

template <long long FPS_MOD>
struct FPS : vector<modint<FPS_MOD>> {
    using mint = modint<FPS_MOD>;
    using vector<mint>::size;
    using vector<mint>::empty;
    using vector<mint>::resize;
    using vector<mint>::operator[];
    FPS() : vector<mint>() {}
    FPS(size_t n) : vector<mint>(n) {}
    FPS(size_t n, const mint& x) : vector<mint>(n, x) {}
    FPS(const initializer_list<mint>& il) : vector<mint>(il) {}
    template<class It>
    FPS(It first, It last) : vector<mint>(first, last) {}
    FPS(const vector<mint>& v) : vector<mint>(v) {}
    FPS(vector<mint>&& v) : vector<mint>(move(v)) {}
    static vector<mint> to_vec(const FPS& f) {
        return vector<mint>(f.begin(), f.end());
    }
    static FPS from_vec(const vector<mint>& v) { return FPS(v); }
    static FPS from_vec(vector<mint>&& v) { return FPS(move(v)); }
    FPS low(int n) const {
        int take = min(n, (int)this->size());
        return FPS(this->begin(), this->begin() + take);
    }
    FPS& shrink() {
        while (!this->empty() && this->back() == mint(0)) this->pop_back();
        return *this;
    }
    FPS operator-() const {
        FPS res = *this;
        for (auto& x : res) x = -x;
        return res;
    }
    FPS& operator+=(const FPS& x) {
        if (size() < x.size()) resize(x.size());
        for (int i = 0; i < (int)x.size(); i++) (*this)[i] += x[i];
        return *this;
    }
    FPS& operator-=(const FPS& x) {
        if (size() < x.size()) resize(x.size());
        for (int i = 0; i < (int)x.size(); i++) (*this)[i] -= x[i];
        return *this;
    }
    FPS& operator*=(const FPS& x) {
        auto va = to_vec(*this);
        auto vb = to_vec(x);
        *this = from_vec(convolution(va, vb));
        return *this;
    }
    FPS operator*(const FPS& x) const { FPS t = *this; t *= x; return t; }
    FPS& operator*=(const mint& c) {
        for (auto& y : *this) y *= c;
        return *this;
    }
    FPS operator*(const mint& x) const { FPS t = *this; t *= x; return t; }
    FPS& operator/=(const mint& c) {
        mint invc = c.inv();
        for (auto& y : *this) y *= invc;
        return *this;
    }
    FPS operator/(const mint& x) const { FPS t = *this; t /= x; return t; }
    FPS operator+(const FPS& x) const { return FPS(*this) += x; }
    FPS operator-(const FPS& x) const { return FPS(*this) -= x; }
    FPS diff() const {
        if (empty()) return FPS();
        FPS g(size() - 1);
        for (int i = 1; i < (int)size(); i++) g[i - 1] = (*this)[i] * mint(i);
        return g;
    }
    FPS integ() const {
        FPS g(size() + 1);
        static vector<mint> invs = {mint(0), mint(1)};
        int n = g.size();
        for (int i=(int)invs.size(); i<n; i++) invs.push_back(-invs[FPS_MOD % i] * mint(FPS_MOD / i));
        for (int i=0; i<(int)size(); i++) g[i+1] = (*this)[i] * invs[i+1];
        return g;
    }
    FPS inv(int deg = -1) const {
        assert(!this->empty() && (*this)[0] != 0);
        int n = this->size();
        if (deg == -1) deg = n;
        FPS res; res.resize(1); res[0] = (*this)[0].inv();
        for (int m=1; m<deg; m<<=1) {
            int cut = 2 * m;
            FPS f_low = this->low(cut);
            auto vres = to_vec(res);
            auto vres2 = vres;
            auto v_f_low = to_vec(f_low);
            auto t1 = convolution(vres, vres2);
            auto t2 = convolution(t1, v_f_low);
            FPS t = from_vec(move(t2));
            t.resize(cut);
            res.resize(cut);
            for (int i=0; i<cut; i++) {
                mint ti = (i < (int)t.size() ? t[i] : mint(0));
                res[i] = res[i] * mint(2) - ti;
            }
            res = res.low(cut);
        }
        res.resize(deg);
        return res;
    }
    pair<FPS, FPS> divmod(const FPS& g) const {
        FPS f = *this;
        if (f.size() < g.size()) return {FPS{mint(0)}, f};
        int n = f.size(), m = g.size();
        FPS rf = f, rg = g;
        reverse(rf.begin(), rf.end());
        reverse(rg.begin(), rg.end());
        auto vrf = to_vec(rf);
        auto vrg = to_vec(rg);
        FPS tmp = from_vec(convolution(vrf, to_vec(rg.inv(n - m + 1))));
        FPS q = tmp.low(n - m + 1);
        reverse(q.begin(), q.end());
        FPS r = f - g * q;
        r.resize(m - 1);
        return {q, r};
    }
    FPS log(int deg = -1) const {
        assert(!this->empty() && (*this)[0] == mint(1));
        int n = this->size();
        if (deg == -1) deg = n;
        FPS f = diff();
        FPS invf = this->inv(deg);
        auto vf = to_vec(f), vinvf = to_vec(invf);
        FPS g = from_vec(convolution(vf, vinvf)).low(deg - 1).integ();
        g.resize(deg);
        return g;
    }
    FPS exp(int deg = -1) const {
        assert(!this->empty() && (*this)[0] == mint(0));
        int n = this->size();
        if (deg == -1) deg = n;
        FPS g; g.resize(1); g[0] = mint(1);
        for (int m=1; m<deg; m<<=1) {
            int cut = 2 * m;
            FPS g_log = g.log(cut);
            FPS q = this->low(cut) - g_log;
            q[0] += mint(1);
            auto vg = to_vec(g), vq = to_vec(q);
            FPS ng = from_vec(convolution(vg, vq)).low(cut);
            g = move(ng);
        }
        g.resize(deg);
        return g;
    }
    FPS pow(long long k, int n=-1) const {
        int v = 0;
        while(v < (int)this->size() && (*this)[v].val == 0) v++;
        if (v == (int)this->size()) return FPS(n, mint(0));
        if (n == -1) n = this->size();
    }
};

Back to top page