diff --git a/.verify-helper/timestamps.remote.json b/.verify-helper/timestamps.remote.json
index 5192d41..95d1a5a 100644
--- a/.verify-helper/timestamps.remote.json
+++ b/.verify-helper/timestamps.remote.json
@@ -32,18 +32,26 @@
 "test/yosupo-enumerate-palindromes.test.cpp": "2023-05-31 16:27:00 +0900",
 "test/yosupo-enumerate-quotients.test.cpp": "2023-05-03 12:22:21 +0900",
 "test/yosupo-enumerate-quotients.test.py": "2023-05-23 13:25:17 +0900",
+"test/yosupo-exp-of-fps.test.cpp": "2024-01-20 00:58:37 +0900",
+"test/yosupo-inv-of-fps.test.cpp": "2024-01-20 00:58:37 +0900",
 "test/yosupo-lca.1.test.cpp": "2023-06-16 15:41:49 +0900",
 "test/yosupo-lca.2.test.cpp": "2023-08-01 18:34:30 +0900",
 "test/yosupo-line-add-get-min.test.cpp": "2023-06-19 07:11:52 +0900",
+"test/yosupo-log-of-fps.test.cpp": "2024-01-20 00:58:37 +0900",
+"test/yosupo-partition-function_1.test.cpp": "2024-01-20 00:58:37 +0900",
+"test/yosupo-partition-function_2.test.cpp": "2024-01-20 00:58:37 +0900",
 "test/yosupo-point-add-rectangle-sum.test.cpp": "2023-06-25 14:14:02 +0900",
 "test/yosupo-point-set-range-composite.1.test.cpp": "2023-08-29 23:08:45 +0900",
 "test/yosupo-point-set-range-composite.2.test.cpp": "2023-08-29 23:08:45 +0900",
 "test/yosupo-point-set-range-composite.test.py": "2023-09-07 14:26:13 +0900",
+"test/yosupo-polynomial-taylor-shift.test.cpp": "2024-01-20 00:58:37 +0900",
+"test/yosupo-pow-of-fps.test.cpp": "2024-01-20 00:58:37 +0900",
 "test/yosupo-range-affine-range-sum.1.test.cpp": "2023-08-29 23:08:45 +0900",
 "test/yosupo-range-affine-range-sum.2.test.cpp": "2023-08-29 23:08:45 +0900",
 "test/yosupo-range-chmin-chmax-add-range-sum.1.test.cpp": "2023-08-01 17:59:54 +0900",
 "test/yosupo-range-chmin-chmax-add-range-sum.2.test.cpp": "2023-08-01 17:59:54 +0900",
 "test/yosupo-segment-add-get-min.test.cpp": "2023-06-19 08:48:44 +0900",
+"test/yosupo-sharp-p-subset-sum.test.cpp": "2024-01-20 00:58:37 +0900",
 "test/yosupo-shortest-path.test.cpp": "2023-06-06 14:52:29 +0900",
 "test/yosupo-shortest-path.test.py": "2023-05-26 10:38:14 +0900",
 "test/yosupo-unionfind.test.cpp": "2023-04-28 12:54:27 +0900",
diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 0000000..e268551
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,13 @@
+{
+    "files.associations": {
+        "array": "cpp",
+        "deque": "cpp",
+        "list": "cpp",
+        "string": "cpp",
+        "unordered_map": "cpp",
+        "vector": "cpp",
+        "string_view": "cpp",
+        "initializer_list": "cpp",
+        "valarray": "cpp"
+    }
+}
\ No newline at end of file
diff --git a/cpp/fps.hpp b/cpp/fps.hpp
new file mode 100644
index 0000000..77449d8
--- /dev/null
+++ b/cpp/fps.hpp
@@ -0,0 +1,224 @@
+#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;
+}
\ No newline at end of file
diff --git a/test/yosupo-exp-of-fps.test.cpp b/test/yosupo-exp-of-fps.test.cpp
new file mode 100644
index 0000000..7a07aab
--- /dev/null
+++ b/test/yosupo-exp-of-fps.test.cpp
@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/exp_of_formal_power_series"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n;
+    std::cin >> n;
+    FPS<mint> a(n);
+    for(auto& x : a) std::cin >> x;
+    auto f = exp(a);
+    for(auto e : f) std::cout << e << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-inv-of-fps.test.cpp b/test/yosupo-inv-of-fps.test.cpp
new file mode 100644
index 0000000..98c8570
--- /dev/null
+++ b/test/yosupo-inv-of-fps.test.cpp
@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/inv_of_formal_power_series"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n;
+    std::cin >> n;
+    FPS<mint> a(n);
+    for(auto& x : a) std::cin >> x;
+    auto f = inv(a);
+    for(auto e : f) std::cout << e << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-log-of-fps.test.cpp b/test/yosupo-log-of-fps.test.cpp
new file mode 100644
index 0000000..758c80a
--- /dev/null
+++ b/test/yosupo-log-of-fps.test.cpp
@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/log_of_formal_power_series"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n;
+    std::cin >> n;
+    FPS<mint> a(n);
+    for(auto& x : a) std::cin >> x;
+    auto f = log(a);
+    for(auto e : f) std::cout << e << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-partition-function_1.test.cpp b/test/yosupo-partition-function_1.test.cpp
new file mode 100644
index 0000000..3749e06
--- /dev/null
+++ b/test/yosupo-partition-function_1.test.cpp
@@ -0,0 +1,19 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/partition_function"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n;
+    std::cin >> n;
+    FPS<mint> f(n + 1, 0);
+    for(int i = 1; i <= n; i ++) {
+        mint inv = mint(i).inv();
+        for(int j = i; j <= n; j += i) f[j] += inv;
+    }
+    f = exp(f);
+    for(auto x : f) std::cout << x << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-partition-function_2.test.cpp b/test/yosupo-partition-function_2.test.cpp
new file mode 100644
index 0000000..4cef1dd
--- /dev/null
+++ b/test/yosupo-partition-function_2.test.cpp
@@ -0,0 +1,23 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/partition_function"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n;
+    std::cin >> n;
+    FPS<mint> f(n + 1, 0);
+    for(int i = 0;; i ++) {
+        if(i * (3 * i - 1) > n * 2) break;
+        f[i * (3 * i - 1) / 2] = (i % 2 ? -1 : 1);
+    }
+    for(int i = -1;; i --) {
+        if(i * (3 * i - 1) > n * 2) break;
+        f[i * (3 * i - 1) / 2] = (i % 2 ? -1 : 1);
+    }
+    f = inv(f);
+    for(auto x : f) std::cout << x << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-polynomial-taylor-shift.test.cpp b/test/yosupo-polynomial-taylor-shift.test.cpp
new file mode 100644
index 0000000..ab7d56a
--- /dev/null
+++ b/test/yosupo-polynomial-taylor-shift.test.cpp
@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/polynomial_taylor_shift"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n, c;
+    std::cin >> n >> c;
+    FPS<mint> a(n);
+    for(auto& x : a) std::cin >> x;
+    auto f = taylor_shift(a, c);
+    for(auto e : f) std::cout << e << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-pow-of-fps.test.cpp b/test/yosupo-pow-of-fps.test.cpp
new file mode 100644
index 0000000..6b0daf8
--- /dev/null
+++ b/test/yosupo-pow-of-fps.test.cpp
@@ -0,0 +1,16 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/pow_of_formal_power_series"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+
+int main() {
+    long long n, m;
+    std::cin >> n >> m;
+    FPS<mint> a(n);
+    for(auto& x : a) std::cin >> x;
+    auto f = pow(a, m);
+    for(auto e : f) std::cout << e << " ";
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/test/yosupo-sharp-p-subset-sum.test.cpp b/test/yosupo-sharp-p-subset-sum.test.cpp
new file mode 100644
index 0000000..3a29a5e
--- /dev/null
+++ b/test/yosupo-sharp-p-subset-sum.test.cpp
@@ -0,0 +1,26 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/sharp_p_subset_sum"
+
+#include "../cpp/fps.hpp"
+
+using mint = modint998244353;
+const int mod = 998244353;
+
+int main() {
+    long long n, t;
+    std::cin >> n >> t;
+    std::vector<mint> inv(t + 1, 1);
+    for(int i = 2; i <= t; i ++) inv[i] = -inv[mod % i] * (mod / i);
+    std::vector<int> cnt(t + 1, 0);
+    for(int i = 0; i < n; i ++) {
+        int x; std::cin >> x;
+        cnt[x] ++;
+    }
+    FPS<mint> f(t + 1);
+    for(int i = 1; i <= t; i ++) for(int j = 1; i * j <= t; j ++) {
+        f[i * j] += mint(cnt[i]) * inv[j] * (j & 1 ? 1 : -1);
+    }
+    f = exp(f);
+    for(int i = 1; i <= t; i ++) std::cout << f[i] << " ";
+    std::cout << std::endl;
+    return 0;
+}