diff --git a/math/composition_of_formal_power_series_large/checker.cpp b/math/composition_of_formal_power_series_large/checker.cpp
new file mode 100644
index 000000000..6a66d5330
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/checker.cpp
@@ -0,0 +1,62 @@
+// https://github.com/MikeMirzayanov/testlib/blob/master/checkers/wcmp.cpp
+
+// The MIT License (MIT)
+
+// Copyright (c) 2015 Mike Mirzayanov
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#include "testlib.h"
+
+using namespace std;
+
+int main(int argc, char * argv[])
+{
+    setName("compare sequences of tokens");
+    registerTestlibCmd(argc, argv);
+
+    int n = 0;
+    string j, p;
+
+    while (!ans.seekEof() && !ouf.seekEof()) 
+    {
+        n++;
+
+        ans.readWordTo(j);
+        ouf.readWordTo(p);
+        
+        if (j != p)
+            quitf(_wa, "%d%s words differ - expected: '%s', found: '%s'", n, englishEnding(n).c_str(), compress(j).c_str(), compress(p).c_str());
+    }
+
+    if (ans.seekEof() && ouf.seekEof())
+    {
+        if (n == 1)
+            quitf(_ok, "\"%s\"", compress(j).c_str());
+        else
+            quitf(_ok, "%d tokens", n);
+    }
+    else
+    {
+        if (ans.seekEof())
+            quitf(_wa, "Participant output contains extra tokens");
+        else
+            quitf(_wa, "Unexpected EOF in the participants output");
+    }
+}
diff --git a/math/composition_of_formal_power_series_large/gen/example_00.in b/math/composition_of_formal_power_series_large/gen/example_00.in
new file mode 100644
index 000000000..cfe630312
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/example_00.in
@@ -0,0 +1,3 @@
+5
+5 4 3 2 1
+0 1 2 3 4
diff --git a/math/composition_of_formal_power_series_large/gen/hack.cpp b/math/composition_of_formal_power_series_large/gen/hack.cpp
new file mode 100644
index 000000000..1a2518d69
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/hack.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int n = N_MAX;
+    vector<int> f(n),g(n);
+    for (int i = 0; i < n; i++){
+        f[i] = gen.uniform(0ll,MOD-1);
+    }
+    for (int i = gen.uniform(0,n-1); i < n; i++) {
+        if(i)g[i] = gen.uniform(0ll,MOD-1);
+        else g[i]=0;
+    }
+
+    printf("%d\n",n);
+
+    for (int i = 0; i < n; i++) {
+        printf("%d", f[i]);
+        if (i != n - 1) printf(" ");
+    }
+    printf("\n");
+    for (int i = 0; i < n; i++) {
+        printf("%d", g[i]);
+        if (i != n - 1) printf(" ");
+    }
+    printf("\n");
+    return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/gen/hack2.cpp b/math/composition_of_formal_power_series_large/gen/hack2.cpp
new file mode 100644
index 000000000..59a981c81
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/hack2.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int n = N_MAX;
+    vector<int> f(n),g(n);
+    for (int i = 0; i < n; i++){
+        f[i] = gen.uniform(0ll,MOD-1);
+    }
+    for (int i = gen.uniform(0,10); i < n; i++) {
+        if(i)g[i] = gen.uniform(0ll,MOD-1);
+        else g[i]=0;
+    }
+
+    printf("%d\n",n);
+
+    for (int i = 0; i < n; i++) {
+        printf("%d", f[i]);
+        if (i != n - 1) printf(" ");
+    }
+    printf("\n");
+    for (int i = 0; i < n; i++) {
+        printf("%d", g[i]);
+        if (i != n - 1) printf(" ");
+    }
+    printf("\n");
+    return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/gen/max_random.cpp b/math/composition_of_formal_power_series_large/gen/max_random.cpp
new file mode 100644
index 000000000..dfc562419
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/max_random.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = N_MAX;
+  vector<int> f(n), g(n);
+  for (int i = 0; i < n; i++) { f[i] = gen.uniform(0ll, MOD - 1); }
+  for (int i = 0; i < n; i++) {
+    if (i)
+      g[i] = gen.uniform(0ll, MOD - 1);
+    else
+      g[i] = 0;
+  }
+
+  printf("%d\n", n);
+
+  for (int i = 0; i < n; i++) {
+    printf("%d", f[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  for (int i = 0; i < n; i++) {
+    printf("%d", g[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/gen/mid.cpp b/math/composition_of_formal_power_series_large/gen/mid.cpp
new file mode 100644
index 000000000..54142f866
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/mid.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform(1ll, N_MID_MAX);
+  vector<int> f(n), g(n);
+  for (int i = 0; i < n; i++) { f[i] = gen.uniform(0ll, MOD - 1); }
+  for (int i = 0; i < n; i++) {
+    if (i)
+      g[i] = gen.uniform(0ll, MOD - 1);
+    else
+      g[i] = 0;
+  }
+
+  printf("%d\n", n);
+
+  for (int i = 0; i < n; i++) {
+    printf("%d", f[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  for (int i = 0; i < n; i++) {
+    printf("%d", g[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/gen/random.cpp b/math/composition_of_formal_power_series_large/gen/random.cpp
new file mode 100644
index 000000000..f1377b4d9
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/random.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform(1ll, N_MAX);
+  vector<int> f(n), g(n);
+  for (int i = 0; i < n; i++) { f[i] = gen.uniform(0ll, MOD - 1); }
+  for (int i = 0; i < n; i++) {
+    if (i)
+      g[i] = gen.uniform(0ll, MOD - 1);
+    else
+      g[i] = 0;
+  }
+
+  printf("%d\n", n);
+
+  for (int i = 0; i < n; i++) {
+    printf("%d", f[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  for (int i = 0; i < n; i++) {
+    printf("%d", g[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/gen/small.cpp b/math/composition_of_formal_power_series_large/gen/small.cpp
new file mode 100644
index 000000000..343f7f0dd
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/gen/small.cpp
@@ -0,0 +1,34 @@
+#include <cstdio>
+#include "../params.h"
+#include "random.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform(1ll, N_SMALL_MAX);
+  vector<int> f(n), g(n);
+  for (int i = 0; i < n; i++) { f[i] = gen.uniform(0ll, MOD - 1); }
+  for (int i = 0; i < n; i++) {
+    if (i)
+      g[i] = gen.uniform(0ll, MOD - 1);
+    else
+      g[i] = 0;
+  }
+
+  printf("%d\n", n);
+
+  for (int i = 0; i < n; i++) {
+    printf("%d", f[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  for (int i = 0; i < n; i++) {
+    printf("%d", g[i]);
+    if (i != n - 1) printf(" ");
+  }
+  printf("\n");
+  return 0;
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/hash.json b/math/composition_of_formal_power_series_large/hash.json
new file mode 100644
index 000000000..250be0bae
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/hash.json
@@ -0,0 +1,66 @@
+{
+  "example_00.in": "2eb7de1f18b179172d224b4f373ecbb661dd167cb16322b550c32ef9a972fcdd",
+  "example_00.out": "4b3202ac36cd9062d26063aa07945dbc320715575042bb3ecbcf0077fb40d734",
+  "hack2_00.in": "ccadf8a5d7c4e7a6e01dfa588c8e2998632c9c96b35e01edf9e84006e3a39991",
+  "hack2_00.out": "41677c75d3a6f9d68755cabfa6cbbf53459ad86602227a229f26f3ab2ebfc6f7",
+  "hack2_01.in": "88d13fe1aaf20bae12e26e8782a935860b68b1d8b112bc93e2e15f9d3ccd9810",
+  "hack2_01.out": "d4ad030cde84df469cefdfe7d76984ac664eb0119e2a7708f361acc109ab8659",
+  "hack2_02.in": "4a174035d278904d5ed519c534fd5977f3109d36c8d13d8341212d5c046c9855",
+  "hack2_02.out": "3251b1ae4d93c57b3a220bfded4539e2a8f744c8bdee1a03bd74d7e7c16e2fde",
+  "hack_00.in": "786b6a019d97a8c13d41e6b7d495451bafc8dcf65450db3b931256b45370f301",
+  "hack_00.out": "5655df70fdb49a3e426a6ed21b5d6e3ed312c8525afa0ca94d438a71c35e5799",
+  "hack_01.in": "66a9244b8941adafd268416854e8e202deaebfa3e8a4ff9b69dc6d6d61343608",
+  "hack_01.out": "78069962c5b9e7dfcfa0b8babcbe45dabe49cf20d12c4ffb419845921b45f09c",
+  "hack_02.in": "2b274c3a33d2c764cfea13c7042a63fba0f45f73e6be4b412f6edbcdff982648",
+  "hack_02.out": "8cffad7dac2ff4ef81d435558330feaeb482fe35d783e53e7f7f6e276b35b62e",
+  "max_random_00.in": "9404320bd03fe869c53cbe50511a45d262f037388dd5d2e41074db1c8a6fffa6",
+  "max_random_00.out": "2e032a16bc09a8e289e7a2470ec8f70b71904ec433c60df30da260a2fb6fd7dd",
+  "max_random_01.in": "6329c2e34e70c85d0d57734d742937361bbcb5a3cb36d1f85b919c8544877593",
+  "max_random_01.out": "361d7ebd0178a95f236eb33f27fbcc192deb7daba88779ed646454876b12f909",
+  "max_random_02.in": "a50c6b629a04643cc4001ab6b33792975097316f395791fa13576a25e36f8275",
+  "max_random_02.out": "211146427e5f5c43f315e3bb21cd9d586d539b9df164ec48fe56ed77e985896c",
+  "max_random_03.in": "d1e9d03fcbec7ac084f01b94c12051536f0d89f148edfca0059850e9925568b6",
+  "max_random_03.out": "4598aedf372145999617e25a845fe6bfd6895e8fe1be6ac2e42dd4be660c730f",
+  "max_random_04.in": "cfe035208770b5b22b248d272872896b9fde8edec65646e974ee3cbaf7418ae6",
+  "max_random_04.out": "927add67f5047b2d214f84df7847d5001bd15b643790a747df9da3fb196ebe44",
+  "mid_00.in": "5974d2d92af1e6fb8772a352b9ac77a3be932d8210af541d3d98a5435e5e4eec",
+  "mid_00.out": "d941b589a815e52383ccf9d20328404603ce8b3928ccdbe9ffbe60ceaf3cadb9",
+  "mid_01.in": "ff53fc06fa315f3862b9d16e576f619a53d0ad357eee42df883e7e58d6799b83",
+  "mid_01.out": "bd5f65f5ae79bbd6da7af00597551fdfa396d3fc2c201fb0016c48e13aff16e6",
+  "mid_02.in": "3992571daa5bbd1464533d931006b69a9f68594d9777953f2380b9659dc50a43",
+  "mid_02.out": "4d3e09d05b44abe8f8b0a740a129de0ce884ace84c9e91feaaa83292c024c86e",
+  "mid_03.in": "9e1de46c00747ef9bdee5eef1f4edbb6bfca56f37c00463b74c417e5f15a92f4",
+  "mid_03.out": "7d71c7428994fb3ef14de2971c5156babe499734ef8cbd6bcc657de791826a48",
+  "mid_04.in": "eb66072ce88bc0a19793fdceeb956622ae5728c9da2ea5fb3b37c4bfe9d7d046",
+  "mid_04.out": "3c720b376814e917da2dd09513b6e02f04ce234d1e45fbaf1baaff2f02c4f041",
+  "random_00.in": "8148b655cb81c1114e8f9abb1a5cc488d326ff82ca53774da1b660353de535e5",
+  "random_00.out": "fd8127c466a06e10ca62583b27be261f4bdf408fd99e5cfb7102378e4bc74cc2",
+  "random_01.in": "a64521ec14f32af5c1ce0da9fcc9cc59ecfb1c4bef671e9e8505faa10b339900",
+  "random_01.out": "9f4a7b13b0de679a7bc6c3f3a379c498825f59968524c9366342f1cf3b45496e",
+  "random_02.in": "383a2eeb79e2d048e78352911b5b69b1e53b4d066e5570ad47d80a4ff66c6957",
+  "random_02.out": "0f4c3d0e18e72a3d9f9cfeab51df4d9d7d2f7ec709eea8c3cef1c6a8023fd7cb",
+  "random_03.in": "51d5c90f2eb26b7cce557646722a117b68367bf26d9779411ef6a0e840acd5d0",
+  "random_03.out": "ffb5f15ab1d014ad1d69e32c271c49d7eadfd6949a195e47fd3fee9987e1f26b",
+  "random_04.in": "3030c0e8acdf602b27bb87d2d080a8a3b06b2581e2e44574cc4285bdda0e4a9e",
+  "random_04.out": "27fc3b93d669f832c43ca0de5dfc8b20cdd016d7bb5408acec3565012210d1ce",
+  "small_00.in": "d018ae901616ea66ba42457bf205efa6106d897ae24862db07ec34272709474c",
+  "small_00.out": "55cac1e5b3088c32ae2210231f585795e621c439d578a3f622854116901905dc",
+  "small_01.in": "4b243ec8ac49dd0b54af3035339a470361ed98189dba1bf50f91761b851eedf9",
+  "small_01.out": "2dd6019379733130e3c0cb7c1d2aa60318e5b6ade39aab875ff3033d0468cb2a",
+  "small_02.in": "8e1746e2b59b22c1a0f878b1bb6145d49dd583f92d10c56bfa3d7b775d42def6",
+  "small_02.out": "4598a8610e26a7af95861e836b141048e4529a7064beb373beadb2e603652779",
+  "small_03.in": "cb91ccacdfa0ef9baa38a29a0b1dfcd7eb2a77188449366f07296b63cff53fea",
+  "small_03.out": "e11527191c91d06e888ceda99d7257b6f10e6676e0c18bae5418f67a33007ec4",
+  "small_04.in": "3cb79fa759f604728767c1985bfcd862ff12561cc27699f7fa739ab47e0b13bd",
+  "small_04.out": "355b7847714ae67e59677ade79b408280f7cd2551bd515d5e78e34a58ff0bf7f",
+  "small_05.in": "5af9f59362f6a142f3ea4dc09317a20e425ad076d81cba0b184aa243b1337b0f",
+  "small_05.out": "a4074b02f113d3e411bea15bf6c485faf85fcc5392c98a292aa95d2790708c11",
+  "small_06.in": "d73255d167f40b2836eb7582a4918d5c606ec152a22d38267700ce2a594c3f13",
+  "small_06.out": "119c2b97c93980a811107009df1747735629b479518c3b119dc6d11c38b24dcc",
+  "small_07.in": "8de4fb5f474c0d2ee3c7387159e55f739e4ef3e77e759e30a8b96f5fd4730244",
+  "small_07.out": "6a631b8c7ec0638305d8b6b073a7a80fecdeff78bf52e4007ee7490cec5014c7",
+  "small_08.in": "f9c438cdafa47aeffed4f54ccdbaa302ebc46a57950b4855cb010167b2fcdfa5",
+  "small_08.out": "9a9daa342941a9de5a31da6820f7e5ae92cead46f6231d5fb4b227401a769084",
+  "small_09.in": "067ad2b7b0ef92d46bf9c6f582414f5435ab63ce2aa57ef268c25e0a178df05e",
+  "small_09.out": "9393c6c1840811e58610a95083f69ed438a3ef856a91009a3eaa72e5ee2e4e5f"
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/info.toml b/math/composition_of_formal_power_series_large/info.toml
new file mode 100644
index 000000000..5a38b75f9
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/info.toml
@@ -0,0 +1,40 @@
+title = 'Composition of Formal Power Series (Large)'
+timelimit = 10.0
+forum = "https://github.com/yosupo06/library-checker-problems/issues/1112"
+
+[[tests]]
+    name = "example.in"
+    number = 1
+[[tests]]
+    name = "random.cpp"
+    number = 5
+[[tests]]
+    name = "mid.cpp"
+    number = 5
+[[tests]]
+    name = "max_random.cpp"
+    number = 5
+[[tests]]
+    name = "small.cpp"
+    number = 10
+[[tests]]
+    name = "hack.cpp"
+    number = 3
+[[tests]]
+    name = "hack2.cpp"
+    number = 3
+[[solutions]]
+    name = "small_correct.cpp"
+    allow_re = true
+[[solutions]]
+    name = "naive.cpp"
+    allow_re = true
+[[solutions]]
+    name = "wa.cpp"
+    wrong = true
+    allow_tle = true
+[params]
+    N_MAX = 131_072
+    N_MID_MAX = 8_000
+    N_SMALL_MAX = 10
+    MOD =998_244_353
diff --git a/math/composition_of_formal_power_series_large/sol/correct.cpp b/math/composition_of_formal_power_series_large/sol/correct.cpp
new file mode 100644
index 000000000..c0e94966e
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/sol/correct.cpp
@@ -0,0 +1,575 @@
+#include <vector>
+#include <cassert>
+#include <cstdio>
+#include <algorithm>
+#include <array>
+
+using namespace std;
+
+using ll = long long;
+using u32 = unsigned int;
+using u64 = unsigned long long;
+
+template <class T>
+using vc = vector<T>;
+template <class T>
+using vvc = vector<vc<T>>;
+
+#define vv(type, name, h, ...) \
+  vector<vector<type>> name(h, vector<type>(__VA_ARGS__))
+
+// https://trap.jp/post/1224/
+#define FOR1(a) for (ll _ = 0; _ < ll(a); ++_)
+#define FOR2(i, a) for (ll i = 0; i < ll(a); ++i)
+#define FOR3(i, a, b) for (ll i = a; i < ll(b); ++i)
+#define FOR4(i, a, b, c) for (ll i = a; i < ll(b); i += (c))
+#define FOR1_R(a) for (ll i = (a)-1; i >= ll(0); --i)
+#define FOR2_R(i, a) for (ll i = (a)-1; i >= ll(0); --i)
+#define FOR3_R(i, a, b) for (ll i = (b)-1; i >= ll(a); --i)
+#define overload4(a, b, c, d, e, ...) e
+#define overload3(a, b, c, d, ...) d
+#define FOR(...) overload4(__VA_ARGS__, FOR4, FOR3, FOR2, FOR1)(__VA_ARGS__)
+#define FOR_R(...) overload3(__VA_ARGS__, FOR3_R, FOR2_R, FOR1_R)(__VA_ARGS__)
+
+#define all(x) x.begin(), x.end()
+#define len(x) ll(x.size())
+#define elif else if
+
+#define eb emplace_back
+#define fi first
+#define se second
+
+// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2)
+int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
+int topbit(u32 x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
+int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
+int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
+
+template <typename T>
+T floor(T a, T b) {
+  return a / b - (a % b && (a ^ b) < 0);
+}
+template <typename T>
+T ceil(T x, T y) {
+  return floor(x + y - 1, y);
+}
+
+template <typename mint>
+mint inv(int n) {
+  static const int mod = mint::get_mod();
+  static vector<mint> dat = {0, 1};
+  assert(0 <= n);
+  if (n >= mod) n %= mod;
+  while (len(dat) <= n) {
+    int k = len(dat);
+    int q = (mod + k - 1) / k;
+    dat.eb(dat[k * q - mod] * mint::raw(q));
+  }
+  return dat[n];
+}
+
+template <typename mint>
+mint fact(int n) {
+  static const int mod = mint::get_mod();
+  assert(0 <= n && n < mod);
+  static vector<mint> dat = {1, 1};
+  while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * mint::raw(len(dat)));
+  return dat[n];
+}
+
+template <typename mint>
+mint fact_inv(int n) {
+  static vector<mint> dat = {1, 1};
+  if (n < 0) return mint(0);
+  while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * inv<mint>(len(dat)));
+  return dat[n];
+}
+
+template <int mod>
+struct modint {
+  static constexpr u32 umod = u32(mod);
+  static_assert(umod < u32(1) << 31);
+  u32 val;
+
+  static modint raw(u32 v) {
+    modint x;
+    x.val = v;
+    return x;
+  }
+  constexpr modint() : val(0) {}
+  constexpr modint(u32 x) : val(x % umod) {}
+  constexpr modint(u64 x) : val(x % umod) {}
+  constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){};
+  constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){};
+  bool operator<(const modint &other) const { return val < other.val; }
+  modint &operator+=(const modint &p) {
+    if ((val += p.val) >= umod) val -= umod;
+    return *this;
+  }
+  modint &operator-=(const modint &p) {
+    if ((val += umod - p.val) >= umod) val -= umod;
+    return *this;
+  }
+  modint &operator*=(const modint &p) {
+    val = u64(val) * p.val % umod;
+    return *this;
+  }
+  modint &operator/=(const modint &p) {
+    *this *= p.inverse();
+    return *this;
+  }
+  modint operator-() const { return modint::raw(val ? mod - val : u32(0)); }
+  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; }
+  bool operator==(const modint &p) const { return val == p.val; }
+  bool operator!=(const modint &p) const { return val != p.val; }
+  modint inverse() const {
+    int a = val, b = mod, u = 1, v = 0, t;
+    while (b > 0) {
+      t = a / b;
+      swap(a -= t * b, b), swap(u -= t * v, v);
+    }
+    return modint(u);
+  }
+  modint pow(ll n) const {
+    assert(n >= 0);
+    modint ret(1), mul(val);
+    while (n > 0) {
+      if (n & 1) ret *= mul;
+      mul *= mul;
+      n >>= 1;
+    }
+    return ret;
+  }
+  static constexpr int get_mod() { return mod; }
+  // (n, r), r は 1 の 2^n 乗根
+  static constexpr pair<int, int> ntt_info() {
+    if (mod == 998244353) return {23, 31};
+    return {-1, -1};
+  }
+  static constexpr bool can_ntt() { return ntt_info().fi != -1; }
+};
+
+using modint998 = modint<998244353>;
+
+template <class T>
+vc<T> convolution_naive(const vc<T> &a, const vc<T> &b) {
+  int n = int(a.size()), m = int(b.size());
+  if (n > m) return convolution_naive<T>(b, a);
+  if (n == 0) return {};
+  vc<T> ans(n + m - 1);
+  if (n <= 16 && (T::get_mod() < (1 << 30))) {
+    for (int k = 0; k < n + m - 1; ++k) {
+      int s = max(0, k - m + 1);
+      int t = min(n, k + 1);
+      u64 sm = 0;
+      for (int i = s; i < t; ++i) { sm += u64(a[i].val) * (b[k - i].val); }
+      ans[k] = sm;
+    }
+  } else {
+    FOR(i, n) FOR(j, m) ans[i + j] += a[i] * b[j];
+  }
+  return ans;
+}
+
+// 任意の環でできる
+template <typename T>
+vc<T> convolution_karatsuba(const vc<T> &f, const vc<T> &g) {
+  const int thresh = 30;
+  if (min(len(f), len(g)) <= thresh) return convolution_naive(f, g);
+  int n = max(len(f), len(g));
+  int m = ceil(n, 2);
+  vc<T> f1, f2, g1, g2;
+  if (len(f) < m) f1 = f;
+  if (len(f) >= m) f1 = {f.begin(), f.begin() + m};
+  if (len(f) >= m) f2 = {f.begin() + m, f.end()};
+  if (len(g) < m) g1 = g;
+  if (len(g) >= m) g1 = {g.begin(), g.begin() + m};
+  if (len(g) >= m) g2 = {g.begin() + m, g.end()};
+  vc<T> a = convolution_karatsuba(f1, g1);
+  vc<T> b = convolution_karatsuba(f2, g2);
+  FOR(i, len(f2)) f1[i] += f2[i];
+  FOR(i, len(g2)) g1[i] += g2[i];
+  vc<T> c = convolution_karatsuba(f1, g1);
+  vc<T> F(len(f) + len(g) - 1);
+  assert(2 * m + len(b) <= len(F));
+  FOR(i, len(a)) F[i] += a[i], c[i] -= a[i];
+  FOR(i, len(b)) F[2 * m + i] += b[i], c[i] -= b[i];
+  if (c.back() == T(0)) c.pop_back();
+  FOR(i, len(c)) if (c[i] != T(0)) F[m + i] += c[i];
+  return F;
+}
+
+template <class mint>
+void ntt(vector<mint> &a, bool inverse) {
+  assert(mint::can_ntt());
+  const int rank2 = mint::ntt_info().fi;
+  const int mod = mint::get_mod();
+  static array<mint, 30> root, iroot;
+  static array<mint, 30> rate2, irate2;
+  static array<mint, 30> rate3, irate3;
+
+  assert(rank2 != -1 && len(a) <= (1 << max(0, rank2)));
+
+  static bool prepared = 0;
+  if (!prepared) {
+    prepared = 1;
+    root[rank2] = mint::ntt_info().se;
+    iroot[rank2] = mint(1) / root[rank2];
+    FOR_R(i, rank2) {
+      root[i] = root[i + 1] * root[i + 1];
+      iroot[i] = iroot[i + 1] * iroot[i + 1];
+    }
+    mint prod = 1, iprod = 1;
+    for (int i = 0; i <= rank2 - 2; i++) {
+      rate2[i] = root[i + 2] * prod;
+      irate2[i] = iroot[i + 2] * iprod;
+      prod *= iroot[i + 2];
+      iprod *= root[i + 2];
+    }
+    prod = 1, iprod = 1;
+    for (int i = 0; i <= rank2 - 3; i++) {
+      rate3[i] = root[i + 3] * prod;
+      irate3[i] = iroot[i + 3] * iprod;
+      prod *= iroot[i + 3];
+      iprod *= root[i + 3];
+    }
+  }
+
+  int n = int(a.size());
+  int h = topbit(n);
+  assert(n == 1 << h);
+  if (!inverse) {
+    int len = 0;
+    while (len < h) {
+      if (h - len == 1) {
+        int p = 1 << (h - len - 1);
+        mint rot = 1;
+        FOR(s, 1 << len) {
+          int offset = s << (h - len);
+          FOR(i, p) {
+            auto l = a[i + offset];
+            auto r = a[i + offset + p] * rot;
+            a[i + offset] = l + r;
+            a[i + offset + p] = l - r;
+          }
+          rot *= rate2[topbit(~s & -~s)];
+        }
+        len++;
+      } else {
+        int p = 1 << (h - len - 2);
+        mint rot = 1, imag = root[2];
+        for (int s = 0; s < (1 << len); s++) {
+          mint rot2 = rot * rot;
+          mint rot3 = rot2 * rot;
+          int offset = s << (h - len);
+          for (int i = 0; i < p; i++) {
+            u64 mod2 = u64(mod) * mod;
+            u64 a0 = a[i + offset].val;
+            u64 a1 = u64(a[i + offset + p].val) * rot.val;
+            u64 a2 = u64(a[i + offset + 2 * p].val) * rot2.val;
+            u64 a3 = u64(a[i + offset + 3 * p].val) * rot3.val;
+            u64 a1na3imag = (a1 + mod2 - a3) % mod * imag.val;
+            u64 na2 = mod2 - a2;
+            a[i + offset] = a0 + a2 + a1 + a3;
+            a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3));
+            a[i + offset + 2 * p] = a0 + na2 + a1na3imag;
+            a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
+          }
+          rot *= rate3[topbit(~s & -~s)];
+        }
+        len += 2;
+      }
+    }
+  } else {
+    mint coef = mint(1) / mint(len(a));
+    FOR(i, len(a)) a[i] *= coef;
+    int len = h;
+    while (len) {
+      if (len == 1) {
+        int p = 1 << (h - len);
+        mint irot = 1;
+        FOR(s, 1 << (len - 1)) {
+          int offset = s << (h - len + 1);
+          FOR(i, p) {
+            u64 l = a[i + offset].val;
+            u64 r = a[i + offset + p].val;
+            a[i + offset] = l + r;
+            a[i + offset + p] = (mod + l - r) * irot.val;
+          }
+          irot *= irate2[topbit(~s & -~s)];
+        }
+        len--;
+      } else {
+        int p = 1 << (h - len);
+        mint irot = 1, iimag = iroot[2];
+        FOR(s, (1 << (len - 2))) {
+          mint irot2 = irot * irot;
+          mint irot3 = irot2 * irot;
+          int offset = s << (h - len + 2);
+          for (int i = 0; i < p; i++) {
+            u64 a0 = a[i + offset + 0 * p].val;
+            u64 a1 = a[i + offset + 1 * p].val;
+            u64 a2 = a[i + offset + 2 * p].val;
+            u64 a3 = a[i + offset + 3 * p].val;
+            u64 x = (mod + a2 - a3) * iimag.val % mod;
+            a[i + offset] = a0 + a1 + a2 + a3;
+            a[i + offset + 1 * p] = (a0 + mod - a1 + x) * irot.val;
+            a[i + offset + 2 * p] = (a0 + a1 + 2 * mod - a2 - a3) * irot2.val;
+            a[i + offset + 3 * p] = (a0 + 2 * mod - a1 - x) * irot3.val;
+          }
+          irot *= irate3[topbit(~s & -~s)];
+        }
+        len -= 2;
+      }
+    }
+  }
+}
+
+template <class mint>
+vector<mint> convolution_ntt(vector<mint> a, vector<mint> b) {
+  if (a.empty() || b.empty()) return {};
+  int n = int(a.size()), m = int(b.size());
+  int sz = 1;
+  while (sz < n + m - 1) sz *= 2;
+
+  // sz = 2^k のときの高速化。分割統治的なやつで損しまくるので。
+  if ((n + m - 3) <= sz / 2) {
+    auto a_last = a.back(), b_last = b.back();
+    a.pop_back(), b.pop_back();
+    auto c = convolution(a, b);
+    c.resize(n + m - 1);
+    c[n + m - 2] = a_last * b_last;
+    FOR(i, len(a)) c[i + len(b)] += a[i] * b_last;
+    FOR(i, len(b)) c[i + len(a)] += b[i] * a_last;
+    return c;
+  }
+
+  a.resize(sz), b.resize(sz);
+  bool same = a == b;
+  ntt(a, 0);
+  if (same) {
+    b = a;
+  } else {
+    ntt(b, 0);
+  }
+  FOR(i, sz) a[i] *= b[i];
+  ntt(a, 1);
+  a.resize(n + m - 1);
+  return a;
+}
+
+template <typename mint>
+vc<mint> convolution(const vc<mint> &a, const vc<mint> &b) {
+  int n = len(a), m = len(b);
+  if (!n || !m) return {};
+  if (min(n, m) <= 50) return convolution_karatsuba<mint>(a, b);
+  return convolution_ntt(a, b);
+}
+
+template <typename T>
+vc<vc<T>> convolution2d(vc<vc<T>> &f, vc<vc<T>> &g) {
+  auto shape = [&](vc<vc<T>> &f) -> pair<int, int> {
+    ll H = len(f);
+    ll W = (H == 0 ? 0 : len(f[0]));
+    return {H, W};
+  };
+  auto [H1, W1] = shape(f);
+  auto [H2, W2] = shape(g);
+  ll H = H1 + H2 - 1;
+  ll W = W1 + W2 - 1;
+
+  vc<T> ff(H1 * W);
+  vc<T> gg(H2 * W);
+  FOR(x, H1) FOR(y, W1) ff[W * x + y] = f[x][y];
+  FOR(x, H2) FOR(y, W2) gg[W * x + y] = g[x][y];
+  auto hh = convolution(ff, gg);
+  vc<vc<T>> h(H, vc<T>(W));
+  FOR(x, H) FOR(y, W) h[x][y] = hh[W * x + y];
+  return h;
+}
+
+// a^0, ..., a^N
+template <typename mint>
+vc<mint> powertable_1(mint a, ll N) {
+  // table of a^i
+  vc<mint> f(N + 1, 1);
+  FOR(i, N) f[i + 1] = a * f[i];
+  return f;
+}
+
+// f(x) -> f(x+c)
+template <typename mint>
+vc<mint> poly_taylor_shift(vc<mint> f, mint c) {
+  ll N = len(f);
+  FOR(i, N) f[i] *= fact<mint>(i);
+  auto b = powertable_1<mint>(c, N);
+  FOR(i, N) b[i] *= fact_inv<mint>(i);
+  reverse(all(f));
+  f = convolution(f, b);
+  f.resize(N);
+  reverse(all(f));
+  FOR(i, N) f[i] *= fact_inv<mint>(i);
+  return f;
+}
+
+// n, m 次多項式 (n>=m) a, b → n-m 次多項式 c
+// c[i] = sum_j b[j]a[i+j]
+template <typename mint>
+vc<mint> middle_product(vc<mint> &a, vc<mint> &b) {
+  assert(len(a) >= len(b));
+  if (b.empty()) return vc<mint>(len(a) - len(b) + 1);
+  if (min(len(b), len(a) - len(b) + 1) <= 60) {
+    return middle_product_naive(a, b);
+  }
+  int n = 1 << topbit(2 * len(a) - 1);
+  vc<mint> fa(n), fb(n);
+  copy(a.begin(), a.end(), fa.begin());
+  copy(b.rbegin(), b.rend(), fb.begin());
+  ntt(fa, 0), ntt(fb, 0);
+  FOR(i, n) fa[i] *= fb[i];
+  ntt(fa, 1);
+  fa.resize(len(a));
+  fa.erase(fa.begin(), fa.begin() + len(b) - 1);
+  return fa;
+}
+
+template <typename mint>
+vc<mint> middle_product_naive(vc<mint> &a, vc<mint> &b) {
+  vc<mint> res(len(a) - len(b) + 1);
+  FOR(i, len(res)) FOR(j, len(b)) res[i] += b[j] * a[i + j];
+  return res;
+}
+
+// https://noshi91.hatenablog.com/entry/2024/03/16/224034
+// O(Nlog^2N)
+template <typename mint>
+vc<mint> composition(vc<mint> f, vc<mint> g) {
+  // 必要ではないが簡略化のために，[x^0]g=0 に帰着しておく.
+  const int N = len(f) - 1;
+  assert(len(f) == N + 1 && len(g) == N + 1);
+  if (N == -1) return {};
+  mint c = g[0];
+  g[0] = 0;
+  f = poly_taylor_shift<mint>(f, c);
+
+  auto degree = [&](vvc<mint> &F) -> pair<int, int> {
+    return {len(F) - 1, len(F[0]) - 1};
+  };
+
+  /*
+  f, g は N 次とする.
+  合成は (f_i) -> sum_i f_ig(x)^i である. これの転置を考える.
+  (h_j) に対して [x^0]h(x^-1)sum_if_ig(x)^i = sum_if_i[x^0]h(x^{-1})g(x)^i
+  であるから [x^N]rev_h(x)g(x)^i を計算するというのが転置である．
+  したがって合成は [x^N]rev_h(x)g(x)^i 列挙の転置をとれば計算できる.
+
+  転置問題は, 次のような反復により解ける：
+  ・入力を reverse する.
+  ・次を k=0,1,2,...,K で繰り返す（K=lg n 程度）.
+  　・g だけから定まる何らかの G_k(x,y) と畳み込む
+  　・適切な部分を抜き出す
+
+  よって合成は次の計算による
+  ・次を k=m,...,2,1,0 で繰り返す
+  　・適切な部分に入れる
+  　・G_k(x,y) との畳み込みの転置でうつす
+  ・reverse して出力する.
+
+  G_k(x,y) は k の昇順に計算できて，これを逆順に処理するため，
+  ここでは再帰関数により実装する．
+  */
+  auto dfs = [&](auto &dfs, vvc<mint> G) -> vvc<mint> {
+    auto [n, m] = degree(G);
+    if (n == 0) {
+      // [x^0]g=0 としていたので G=1 である.
+      // 転置問題ではこの時点での多項式がそのまま出力となる.
+      // よって入力をそのまま多項式に入れて返す.
+      vv(mint, F, n + 1, m + 1);
+      FOR(j, len(f)) F[0][j] = f[j];
+      return F;
+    }
+    vvc<mint> G1 = G;
+    FOR(i, n + 1) {
+      if (i % 2 == 0) continue;
+      FOR(j, m + 1) G1[i][j] = -G[i][j];
+    }
+    vvc<mint> G2 = convolution2d(G, G1);
+    int n2 = n / 2;
+    FOR(i, n2 + 1) G2[i] = G2[2 * i];
+    G2.resize(n2 + 1);
+    vvc<mint> F = dfs(dfs, G2);
+    assert(degree(F) == degree(G2));
+    // 転置問題では G1 と畳み込んだあと n の parity に応じた部分を抜き出す.
+    // よってまず dfs の出力を抜き出す前の位置に入れて, 畳み込みの転置をとる.
+
+    // 畳み込みの転置をとるとは何なのか.
+    // (n,m) 次の多項式積 H(x,y) G(x,y)
+    // H, G を適切に配置した (n+1)(2m+1)-1 次多項式を用意して
+    // 畳み込んで 2(n+1)(2m+1)-2 次多項式を得る.
+    // この適切な (2n,2m) 次多項式を抜き出す.
+    // 転置をとると次のようになる.
+    // F を適切に配置した (2n,2m) 次多項式を作る.
+    // ここから 2(n+1)(2m+1)-2 次多項式を作る.
+    // これと適切な多項式の middle product をとって
+    // (n+1)(2m+1)-1 次多項式を得る.
+    vc<mint> A(2 * (n + 1) * (2 * m + 1) - 1);
+    FOR(i, len(F)) FOR(j, len(F[0])) {
+      int i1 = 2 * i + (n % 2);
+      int idx = i1 * (2 * m + 1) + j;
+      A[idx] = F[i][j];
+    }
+    vc<mint> B((n + 1) * (2 * m + 1));
+    FOR(i, n + 1) FOR(j, m + 1) {
+      int idx = i * (2 * m + 1) + j;
+      B[idx] = G1[i][j];
+    }
+
+    A = middle_product<mint>(A, B);
+    assert(len(A) == (n + 1) * (2 * m + 1));
+    vv(mint, res, n + 1, m + 1);
+    FOR(i, n + 1) FOR(j, m + 1) {
+      int idx = (2 * m + 1) * i + j;
+      res[i][j] = A[idx];
+    }
+    return res;
+  };
+  vv(mint, G, N + 1, 2);
+  G[0][0] = 1;
+  FOR(i, N + 1) G[i][1] = -g[i];
+  vvc<mint> F = dfs(dfs, G);
+  assert(degree(F) == (pair<int, int>{N, 1}));
+  vc<mint> ANS(N + 1);
+  FOR(i, N + 1) ANS[i] = F[i][0];
+  // 最後に reverse する
+  reverse(all(ANS));
+  return ANS;
+}
+
+using mint = modint998;
+
+void solve() {
+  int N;
+  scanf("%d", &N);
+  vc<mint> F(N), G(N);
+  FOR(i, N) {
+    int x;
+    scanf("%d", &x);
+    F[i] = x;
+  }
+  FOR(i, N) {
+    int x;
+    scanf("%d", &x);
+    G[i] = x;
+  }
+  F = composition<mint>(F, G);
+  FOR(i, N) {
+    if (i) printf(" ");
+    printf("%d", F[i].val);
+  }
+  printf("\n");
+}
+
+signed main() { solve(); }
diff --git a/math/composition_of_formal_power_series_large/sol/naive.cpp b/math/composition_of_formal_power_series_large/sol/naive.cpp
new file mode 100644
index 000000000..9940498fd
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/sol/naive.cpp
@@ -0,0 +1,433 @@
+#include <iostream>
+#include <vector>
+#include <math.h>
+#include <assert.h>
+#include <algorithm>
+#include <cstdio>
+#include <cstdint>
+
+#define MOD 998244353LL
+using namespace ::std;
+
+template <typename T, typename F>
+struct FPS_BASE : vector<T> {
+  using vector<T>::vector;
+  using P = FPS_BASE<T, F>;
+  F fft;
+  FPS_BASE() {}
+  inline P operator+(T x) const noexcept { return P(*this) += x; }
+  inline P operator-(T x) const noexcept { return P(*this) -= x; }
+  inline P operator*(T x) const noexcept { return P(*this) *= x; }
+  inline P operator/(T x) const noexcept { return P(*this) /= x; }
+  inline P operator<<(int x) noexcept { return P(*this) <<= x; }
+  inline P operator>>(int x) noexcept { return P(*this) >>= x; }
+  inline P operator+(const P& x) const noexcept { return P(*this) += x; }
+  inline P operator-(const P& x) const noexcept { return P(*this) -= x; }
+  inline P operator-() const noexcept { return P(1, T(0)) -= P(*this); }
+  inline P operator*(const P& x) const noexcept { return P(*this) *= x; }
+  inline P operator/(const P& x) const noexcept { return P(*this) /= x; }
+  inline P operator%(const P& x) const noexcept { return P(*this) %= x; }
+  bool operator==(P x) {
+    for (int i = 0; i < (int)max((*this).size(), x.size()); ++i) {
+      if (i >= (int)(*this).size() && x[i] != T()) return 0;
+      if (i >= (int)x.size() && (*this)[i] != T()) return 0;
+      if (i < (int)min((*this).size(), x.size()))
+        if ((*this)[i] != x[i]) return 0;
+    }
+    return 1;
+  }
+  P& operator+=(T x) {
+    if (this->size() == 0) this->resize(1, T(0));
+    (*this)[0] += x;
+    return (*this);
+  }
+  P& operator-=(T x) {
+    if (this->size() == 0) this->resize(1, T(0));
+    (*this)[0] -= x;
+    return (*this);
+  }
+  P& operator*=(T x) {
+    for (int i = 0; i < (int)this->size(); ++i) { (*this)[i] *= x; }
+    return (*this);
+  }
+  P& operator/=(T x) { return (*this) *= (T(1) / x); }
+  P& operator<<=(int x) {
+    P ret(x, T(0));
+    ret.insert(ret.end(), begin(*this), end(*this));
+    return (*this) = ret;
+  }
+  P& operator>>=(int x) {
+    P ret;
+    ret.insert(ret.end(), begin(*this) + x, end(*this));
+    return (*this) = ret;
+  }
+  P& operator+=(const P& x) {
+    if (this->size() < x.size()) this->resize(x.size(), T(0));
+    for (int i = 0; i < (int)x.size(); ++i) { (*this)[i] += x[i]; }
+    return (*this);
+  }
+  P& operator-=(const P& x) {
+    if (this->size() < x.size()) this->resize(x.size(), T(0));
+    for (int i = 0; i < (int)x.size(); ++i) { (*this)[i] -= x[i]; }
+    return (*this);
+  }
+  P& operator*=(const P& x) { return (*this) = F()(*this, x); }
+  P& operator/=(P x) {
+    if (this->size() < x.size()) {
+      this->clear();
+      return (*this);
+    }
+    const int n = this->size() - x.size() + 1;
+    return (*this) = (rev().pre(n) * x.rev().inv(n)).pre(n).rev(n);
+  }
+  P& operator%=(const P& x) { return ((*this) -= *this / x * x); }
+  inline void print() {
+    for (int i = 0; i < (int)(*this).size(); ++i)
+      cerr << (*this)[i] << " \n"[i == (int)(*this).size() - 1];
+    if ((int)(*this).size() == 0) cerr << endl;
+  }
+  inline P& shrink() {
+    while ((*this).back() == 0) (*this).pop_back();
+    return (*this);
+  }
+  inline P pre(int sz) const {
+    return P(begin(*this), begin(*this) + min((int)this->size(), sz));
+  }
+  P rev(int deg = -1) {
+    P ret(*this);
+    if (deg != -1) ret.resize(deg, T(0));
+    reverse(begin(ret), end(ret));
+    return ret;
+  }
+  P inv(int deg = -1) {
+    assert((*this)[0] != T(0));
+    const int n = deg == -1 ? this->size() : deg;
+    P ret({T(1) / (*this)[0]});
+    for (int i = 1; i < n; i <<= 1) {
+      ret *= (-ret * pre(i << 1) + 2).pre(i << 1);
+    }
+    return ret.pre(n);
+  }
+  inline P dot(const P& x) {
+    P ret(*this);
+    for (int i = 0; i < int(min(this->size(), x.size())); ++i) {
+      ret[i] *= x[i];
+    }
+    return ret;
+  }
+  P diff() {
+    if ((int)(*this).size() <= 1) return P();
+    P ret(*this);
+    for (int i = 0; i < (int)ret.size(); i++) { ret[i] *= i; }
+    return ret >> 1;
+  }
+  P integral() {
+    P ret(*this);
+    for (int i = 0; i < (int)ret.size(); i++) { ret[i] /= i + 1; }
+    return ret << 1;
+  }
+  P log(int deg = -1) {
+    assert((*this)[0] == T(1));
+    const int n = deg == -1 ? this->size() : deg;
+    return (diff() * inv(n)).pre(n - 1).integral();
+  }
+  P exp(int deg = -1) {
+    assert((*this)[0] == T(0));
+    const int n = deg == -1 ? this->size() : deg;
+    P ret({T(1)});
+    for (int i = 1; i < n; i <<= 1) {
+      ret = ret * (pre(i << 1) + 1 - ret.log(i << 1)).pre(i << 1);
+    }
+    return ret.pre(n);
+  }
+  P pow(int c, int deg = -1) {
+    // const int n=deg==-1?this->size():deg;
+    // long long i=0;
+    // P ret(*this);
+    // while(i!=(int)this->size()&&ret[i]==0)i++;
+    // if(i==(int)this->size())return P(n,0);
+    // if(i*c>=n)return P(n,0);
+    // T k=ret[i];
+    // return ((((ret>>i)/k).log()*c).exp()*(k.pow(c))<<(i*c)).pre(n);
+    P x(*this);
+    P ret(1, 1);
+    while (c) {
+      if (c & 1) {
+        ret *= x;
+        if (~deg) ret = ret.pre(deg);
+      }
+      x *= x;
+      if (~deg) x = x.pre(deg);
+      c >>= 1;
+    }
+    return ret;
+  }
+  P sqrt(int deg = -1) {
+    const int n = deg == -1 ? this->size() : deg;
+    if ((*this)[0] == T(0)) {
+      for (int i = 1; i < (int)this->size(); i++) {
+        if ((*this)[i] != T(0)) {
+          if (i & 1) return {};
+          if (n - i / 2 <= 0) break;
+          auto ret = (*this >> i).sqrt(n - i / 2) << (i / 2);
+          if ((int)ret.size() < n) ret.resize(n, T(0));
+          return ret;
+        }
+      }
+      return P(n, 0);
+    }
+    P ret({T(1)});
+    for (int i = 1; i < n; i <<= 1) {
+      ret = (ret + pre(i << 1) * ret.inv(i << 1)).pre(i << 1) / T(2);
+    }
+    return ret.pre(n);
+  }
+  P shift(int c) {
+    const int n = this->size();
+    P f(*this), g(n, 0);
+    for (int i = 0; i < n; ++i) f[i] *= T(i).fact();
+    for (int i = 0; i < n; ++i) g[i] = T(c).pow(i) / T(i).fact();
+    g = g.rev();
+    f *= g;
+    f >>= n - 1;
+    for (int i = 0; i < n; ++i) f[i] /= T(i).fact();
+    return f;
+  }
+  T eval(T x) {
+    T res = 0;
+    for (int i = (int)this->size() - 1; i >= 0; --i) {
+      res *= x;
+      res += (*this)[i];
+    }
+    return res;
+  }
+  vector<T> multipoint_eval(const vector<T>& x) {
+    const int n = x.size();
+    P* v = new P[2 * n - 1];
+    for (int i = 0; i < n; i++) v[i + n - 1] = {T() - x[i], T(1)};
+    for (int i = n - 2; i >= 0; i--) { v[i] = v[i * 2 + 1] * v[i * 2 + 2]; }
+    v[0] = P(*this) % v[0];
+    v[0].shrink();
+    for (int i = 1; i < n * 2 - 1; i++) {
+      v[i] = v[(i - 1) / 2] % v[i];
+      v[i].shrink();
+    }
+    vector<T> res(n);
+    for (int i = 0; i < n; i++) res[i] = v[i + n - 1][0];
+    return res;
+  }
+  P slice(int s, int e, int k) {
+    P res;
+    for (int i = s; i < e; i += k) res.push_back((*this)[i]);
+    return res;
+  }
+  T nth_term(P q, int64_t x) {
+    if (x == 0) return (*this)[0] / q[0];
+    P p(*this);
+    P q2 = q;
+    for (int i = 1; i < (int)q2.size(); i += 2) q2[i] *= -1;
+    q *= q2;
+    p *= q2;
+    return p.slice(x % 2, p.size(), 2).nth_term(q.slice(0, q.size(), 2), x / 2);
+  }
+
+  //(*this)(t(x))
+  P manipulate(P t, int deg) {
+    P s = P(*this);
+    if (deg == 0) return P();
+    if ((int)t.size() == 1) return P{s.eval(t[0])};
+    int k = min((int)::sqrt(deg / (::log2(deg) + 1)) + 1, (int)t.size());
+    int b = deg / k + 1;
+    P t2 = t.pre(k);
+    vector<P> table(s.size() / 2 + 1, P{1});
+    for (int i = 1; i < (int)table.size(); i++) {
+      table[i] = ((table[i - 1]) * t2).pre(deg);
+    }
+    auto f = [&](auto f, auto l, auto r, int deg) -> P {
+      if (r - l == 1) return P{*l};
+      auto m = l + (r - l) / 2;
+      return f(f, l, m, deg) + (table[m - l] * f(f, m, r, deg)).pre(deg);
+    };
+    P ans = P();
+    P tmp = f(f, s.begin(), s.end(), deg);
+    P tmp2 = P{1};
+    T tmp3 = T(1);
+    P tmp4 = t2.diff().inv(deg);
+    for (int i = 0; i < b; ++i) {
+      ans += (tmp * tmp2).pre(deg) / tmp3;
+      tmp = (tmp.diff() * tmp4).pre(deg);
+      tmp2 = (tmp2 * (t - t2)).pre(deg);
+      tmp3 *= T(i + 1);
+    }
+    return ans;
+  }
+  //(*this)(t(x))
+  P manipulate2(P t, int deg) {
+    P ans = P();
+    P s = (*this).rev();
+    for (int i = 0; i < (int)s.size(); ++i) { ans = (ans * t + s[i]).pre(deg); }
+    return ans;
+  }
+};
+
+template <typename Mint>
+struct _FPS9 {
+  template <typename T>
+  T operator()(T s, T t) {
+    if (s == decltype(s)()) return decltype(s)();
+    if (t == decltype(t)()) return decltype(t)();
+    auto ntt = [](auto v, const bool& inv) {
+      const int n = v.size();
+      assert(Mint::get_mod() >= 3 && Mint::get_mod() % 2 == 1);
+      decltype(v) tmp(n, Mint());
+      Mint root = inv ? Mint(Mint::root()).inv() : Mint::root();
+      for (int b = n >> 1; b >= 1; b >>= 1, v.swap(tmp)) {
+        Mint w = root.pow((Mint::get_mod() - 1) / (n / b)), p = 1;
+        for (int i = 0; i < n; i += b * 2, p *= w)
+          for (int j = 0; j < b; ++j) {
+            v[i + j + b] *= p;
+            tmp[i / 2 + j] = v[i + j] + v[i + b + j];
+            tmp[n / 2 + i / 2 + j] = v[i + j] - v[i + b + j];
+          }
+      }
+      if (inv) v /= n;
+      return v;
+    };
+    const int n = s.size() + t.size() - 1;
+    int h = 1;
+    while ((1 << h) < n) h++;
+    s.resize((1 << h), Mint(0));
+    t.resize((1 << h), Mint(0));
+    return ntt(ntt(s, 0).dot(ntt(t, 0)), 1).pre(n);
+  }
+};
+template <typename Mint>
+using fps = FPS_BASE<Mint, _FPS9<Mint>>;
+
+class mint {
+  using u64 = std::uint_fast64_t;
+
+public:
+  u64 a;
+  constexpr mint(const long long x = 0) noexcept
+      : a(x >= 0 ? x % get_mod() : get_mod() - (-x) % get_mod()) {}
+  constexpr u64& value() noexcept { return a; }
+  constexpr const u64& value() const noexcept { return a; }
+  constexpr mint operator+(const mint rhs) const noexcept {
+    return mint(*this) += rhs;
+  }
+  constexpr mint operator-(const mint rhs) const noexcept {
+    return mint(*this) -= rhs;
+  }
+  constexpr mint operator*(const mint rhs) const noexcept {
+    return mint(*this) *= rhs;
+  }
+  constexpr mint operator/(const mint rhs) const noexcept {
+    return mint(*this) /= rhs;
+  }
+  constexpr mint& operator+=(const mint rhs) noexcept {
+    a += rhs.a;
+    if (a >= get_mod()) a -= get_mod();
+    return *this;
+  }
+  constexpr mint& operator-=(const mint rhs) noexcept {
+    if (a < rhs.a) a += get_mod();
+    a -= rhs.a;
+    return *this;
+  }
+  constexpr mint& operator*=(const mint rhs) noexcept {
+    a = a * rhs.a % get_mod();
+    return *this;
+  }
+  constexpr mint operator++(int) noexcept {
+    a += 1;
+    if (a >= get_mod()) a -= get_mod();
+    return *this;
+  }
+  constexpr mint operator--(int) noexcept {
+    if (a < 1) a += get_mod();
+    a -= 1;
+    return *this;
+  }
+  constexpr mint& operator/=(mint rhs) noexcept {
+    u64 exp = get_mod() - 2;
+    while (exp) {
+      if (exp % 2) { *this *= rhs; }
+      rhs *= rhs;
+      exp /= 2;
+    }
+    return *this;
+  }
+  constexpr bool operator==(mint x) noexcept { return a == x.a; }
+  constexpr bool operator!=(mint x) noexcept { return a != x.a; }
+  constexpr bool operator<(mint x) noexcept { return a < x.a; }
+  constexpr bool operator>(mint x) noexcept { return a > x.a; }
+  constexpr bool operator<=(mint x) noexcept { return a <= x.a; }
+  constexpr bool operator>=(mint x) noexcept { return a >= x.a; }
+  constexpr static int root() {
+    mint root = 2;
+    while (root.pow((get_mod() - 1) >> 1).a == 1) root++;
+    return root.a;
+  }
+  constexpr mint pow(long long n) {
+    long long x = a;
+    mint ret = 1;
+    while (n > 0) {
+      if (n & 1) (ret *= x);
+      (x *= x) %= get_mod();
+      n >>= 1;
+    }
+    return ret;
+  }
+  constexpr mint inv() { return pow(get_mod() - 2); }
+  static vector<mint> fac, ifac;
+  static bool init;
+  constexpr static int mx = 10000001;
+  void build() {
+    init = 0;
+    fac.resize(mx);
+    ifac.resize(mx);
+    fac[0] = 1, ifac[0] = 1;
+    for (int i = 1; i < mx; i++) fac[i] = fac[i - 1] * i;
+    ifac[mx - 1] = fac[mx - 1].inv();
+    for (int i = mx - 2; i >= 0; i--) ifac[i] = ifac[i + 1] * (i + 1);
+  }
+  mint comb(long long b) {
+    if (init) build();
+    if (a == 0 && b == 0) return 1;
+    if ((long long)a < b) return 0;
+    return fac[a] * ifac[a - b] * ifac[b];
+  }
+  mint fact() {
+    if (init) build();
+    return fac[a];
+  }
+  mint fact_inv() {
+    if (init) build();
+    return ifac[a];
+  }
+  friend ostream& operator<<(ostream& lhs, const mint& rhs) noexcept {
+    lhs << rhs.a;
+    return lhs;
+  }
+  friend istream& operator>>(istream& lhs, mint& rhs) noexcept {
+    lhs >> rhs.a;
+    return lhs;
+  }
+  constexpr static u64 get_mod() { return MOD; }
+};
+vector<mint> mint::fac;
+vector<mint> mint::ifac;
+bool mint::init = 1;
+
+int main() {
+  int n;
+  cin >> n;
+  assert(n <= 100);
+  fps<mint> f(n), g(n);
+  for (int i = 0; i < n; ++i) cin >> f[i];
+  for (int i = 0; i < n; ++i) cin >> g[i];
+  auto h = f.manipulate2(g, n);
+  h.resize(n, 0);
+  for (int i = 0; i < n; ++i) cout << h[i] << " \n"[i == n - 1];
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/sol/small_correct.cpp b/math/composition_of_formal_power_series_large/sol/small_correct.cpp
new file mode 100644
index 000000000..270710dd3
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/sol/small_correct.cpp
@@ -0,0 +1,448 @@
+#include <iostream>
+#include <vector>
+#include <math.h>
+#include <assert.h>
+#include <cstdint>
+#include <algorithm>
+#define MOD 998244353LL
+using namespace ::std;
+
+template <typename T, typename F>
+struct FPS_BASE : vector<T> {
+  using vector<T>::vector;
+  using P = FPS_BASE<T, F>;
+  F fft;
+  FPS_BASE() {}
+  inline P operator+(T x) const noexcept { return P(*this) += x; }
+  inline P operator-(T x) const noexcept { return P(*this) -= x; }
+  inline P operator*(T x) const noexcept { return P(*this) *= x; }
+  inline P operator/(T x) const noexcept { return P(*this) /= x; }
+  inline P operator<<(int x) noexcept { return P(*this) <<= x; }
+  inline P operator>>(int x) noexcept { return P(*this) >>= x; }
+  inline P operator+(const P& x) const noexcept { return P(*this) += x; }
+  inline P operator-(const P& x) const noexcept { return P(*this) -= x; }
+  inline P operator-() const noexcept { return P(1, T(0)) -= P(*this); }
+  inline P operator*(const P& x) const noexcept { return P(*this) *= x; }
+  inline P operator/(const P& x) const noexcept { return P(*this) /= x; }
+  inline P operator%(const P& x) const noexcept { return P(*this) %= x; }
+  bool operator==(P x) {
+    for (int i = 0; i < (int)max((*this).size(), x.size()); ++i) {
+      if (i >= (int)(*this).size() && x[i] != T()) return 0;
+      if (i >= (int)x.size() && (*this)[i] != T()) return 0;
+      if (i < (int)min((*this).size(), x.size()))
+        if ((*this)[i] != x[i]) return 0;
+    }
+    return 1;
+  }
+  P& operator+=(T x) {
+    if (this->size() == 0) this->resize(1, T(0));
+    (*this)[0] += x;
+    return (*this);
+  }
+  P& operator-=(T x) {
+    if (this->size() == 0) this->resize(1, T(0));
+    (*this)[0] -= x;
+    return (*this);
+  }
+  P& operator*=(T x) {
+    for (int i = 0; i < (int)this->size(); ++i) { (*this)[i] *= x; }
+    return (*this);
+  }
+  P& operator/=(T x) { return (*this) *= (T(1) / x); }
+  P& operator<<=(int x) {
+    P ret(x, T(0));
+    ret.insert(ret.end(), begin(*this), end(*this));
+    return (*this) = ret;
+  }
+  P& operator>>=(int x) {
+    if ((int)(*this).size() <= x) return (*this) = P();
+    P ret;
+    ret.insert(ret.end(), begin(*this) + x, end(*this));
+    return (*this) = ret;
+  }
+  P& operator+=(const P& x) {
+    if (this->size() < x.size()) this->resize(x.size(), T(0));
+    for (int i = 0; i < (int)x.size(); ++i) { (*this)[i] += x[i]; }
+    return (*this);
+  }
+  P& operator-=(const P& x) {
+    if (this->size() < x.size()) this->resize(x.size(), T(0));
+    for (int i = 0; i < (int)x.size(); ++i) { (*this)[i] -= x[i]; }
+    return (*this);
+  }
+  P& operator*=(const P& x) { return (*this) = F()(*this, x); }
+  P& operator/=(P x) {
+    if (this->size() < x.size()) {
+      this->clear();
+      return (*this);
+    }
+    const int n = this->size() - x.size() + 1;
+    return (*this) = (rev().pre(n) * x.rev().inv(n)).pre(n).rev(n);
+  }
+  P& operator%=(const P& x) { return ((*this) -= *this / x * x); }
+  inline void print() {
+    for (int i = 0; i < (int)(*this).size(); ++i)
+      cerr << (*this)[i] << " \n"[i == (int)(*this).size() - 1];
+    if ((int)(*this).size() == 0) cerr << endl;
+  }
+  inline P& shrink() {
+    while ((*this).back() == 0) (*this).pop_back();
+    return (*this);
+  }
+  inline P pre(int sz) const {
+    return P(begin(*this), begin(*this) + min((int)this->size(), sz));
+  }
+  P rev(int deg = -1) {
+    P ret(*this);
+    if (deg != -1) ret.resize(deg, T(0));
+    reverse(begin(ret), end(ret));
+    return ret;
+  }
+  P inv(int deg = -1) {
+    assert((*this)[0] != T(0));
+    const int n = deg == -1 ? this->size() : deg;
+    P ret({T(1) / (*this)[0]});
+    for (int i = 1; i < n; i <<= 1) {
+      ret *= (-ret * pre(i << 1) + 2).pre(i << 1);
+    }
+    return ret.pre(n);
+  }
+  inline P dot(const P& x) {
+    P ret(*this);
+    for (int i = 0; i < int(min(this->size(), x.size())); ++i) {
+      ret[i] *= x[i];
+    }
+    return ret;
+  }
+  P diff() {
+    if ((int)(*this).size() <= 1) return P();
+    P ret(*this);
+    for (int i = 0; i < (int)ret.size(); i++) { ret[i] *= i; }
+    return ret >> 1;
+  }
+  P integral() {
+    P ret(*this);
+    for (int i = 0; i < (int)ret.size(); i++) { ret[i] /= i + 1; }
+    return ret << 1;
+  }
+  P log(int deg = -1) {
+    assert((*this)[0] == T(1));
+    const int n = deg == -1 ? this->size() : deg;
+    return (diff() * inv(n)).pre(n - 1).integral();
+  }
+  P exp(int deg = -1) {
+    assert((*this)[0] == T(0));
+    const int n = deg == -1 ? this->size() : deg;
+    P ret({T(1)});
+    for (int i = 1; i < n; i <<= 1) {
+      ret = ret * (pre(i << 1) + 1 - ret.log(i << 1)).pre(i << 1);
+    }
+    return ret.pre(n);
+  }
+  P pow(int c, int deg = -1) {
+    // const int n=deg==-1?this->size():deg;
+    // long long i=0;
+    // P ret(*this);
+    // while(i!=(int)this->size()&&ret[i]==0)i++;
+    // if(i==(int)this->size())return P(n,0);
+    // if(i*c>=n)return P(n,0);
+    // T k=ret[i];
+    // return ((((ret>>i)/k).log()*c).exp()*(k.pow(c))<<(i*c)).pre(n);
+    P x(*this);
+    P ret(1, 1);
+    while (c) {
+      if (c & 1) {
+        ret *= x;
+        if (~deg) ret = ret.pre(deg);
+      }
+      x *= x;
+      if (~deg) x = x.pre(deg);
+      c >>= 1;
+    }
+    return ret;
+  }
+  P sqrt(int deg = -1) {
+    const int n = deg == -1 ? this->size() : deg;
+    if ((*this)[0] == T(0)) {
+      for (int i = 1; i < (int)this->size(); i++) {
+        if ((*this)[i] != T(0)) {
+          if (i & 1) return {};
+          if (n - i / 2 <= 0) break;
+          auto ret = (*this >> i).sqrt(n - i / 2) << (i / 2);
+          if ((int)ret.size() < n) ret.resize(n, T(0));
+          return ret;
+        }
+      }
+      return P(n, 0);
+    }
+    P ret({T(1)});
+    for (int i = 1; i < n; i <<= 1) {
+      ret = (ret + pre(i << 1) * ret.inv(i << 1)).pre(i << 1) / T(2);
+    }
+    return ret.pre(n);
+  }
+  P shift(int c) {
+    const int n = this->size();
+    P f(*this), g(n, 0);
+    for (int i = 0; i < n; ++i) f[i] *= T(i).fact();
+    for (int i = 0; i < n; ++i) g[i] = T(c).pow(i) / T(i).fact();
+    g = g.rev();
+    f *= g;
+    f >>= n - 1;
+    for (int i = 0; i < n; ++i) f[i] /= T(i).fact();
+    return f;
+  }
+  T eval(T x) {
+    T res = 0;
+    for (int i = (int)this->size() - 1; i >= 0; --i) {
+      res *= x;
+      res += (*this)[i];
+    }
+    return res;
+  }
+  vector<T> multipoint_eval(const vector<T>& x) {
+    const int n = x.size();
+    P* v = new P[2 * n - 1];
+    for (int i = 0; i < n; i++) v[i + n - 1] = {T() - x[i], T(1)};
+    for (int i = n - 2; i >= 0; i--) { v[i] = v[i * 2 + 1] * v[i * 2 + 2]; }
+    v[0] = P(*this) % v[0];
+    v[0].shrink();
+    for (int i = 1; i < n * 2 - 1; i++) {
+      v[i] = v[(i - 1) / 2] % v[i];
+      v[i].shrink();
+    }
+    vector<T> res(n);
+    for (int i = 0; i < n; i++) res[i] = v[i + n - 1][0];
+    return res;
+  }
+  P slice(int s, int e, int k) {
+    P res;
+    for (int i = s; i < e; i += k) res.push_back((*this)[i]);
+    return res;
+  }
+  T nth_term(P q, int64_t x) {
+    if (x == 0) return (*this)[0] / q[0];
+    P p(*this);
+    P q2 = q;
+    for (int i = 1; i < (int)q2.size(); i += 2) q2[i] *= -1;
+    q *= q2;
+    p *= q2;
+    return p.slice(x % 2, p.size(), 2).nth_term(q.slice(0, q.size(), 2), x / 2);
+  }
+
+  //(*this)(t(x))
+  P manipulate(P t, int deg) {
+    P s = P(*this);
+    if (deg == 0) return P();
+    if ((int)t.size() == 1) return P{s.eval(t[0])};
+    int k = min((int)::sqrt(deg / (::log2(deg) + 1)) + 1, (int)t.size());
+    int b = deg / k + 1;
+    P t2 = t.pre(k);
+    vector<P> table(s.size() / 2 + 1, P{1});
+    for (int i = 1; i < (int)table.size(); i++) {
+      table[i] = ((table[i - 1]) * t2).pre(deg);
+    }
+    auto f = [&](auto f, auto l, auto r, int deg) -> P {
+      if (r - l == 1) return P{*l};
+      auto m = l + (r - l) / 2;
+      return f(f, l, m, deg) + (table[m - l] * f(f, m, r, deg)).pre(deg);
+    };
+    P ans = P();
+    P tmp = f(f, s.begin(), s.end(), deg);
+    P tmp2 = P{1};
+    T tmp3 = T(1);
+    int tmp5 = -1;
+    P tmp6 = t2.diff();
+    if (tmp6 == P()) {
+      for (int i = 0; i < b; ++i) {
+        if (i >= (int)s.size()) break;
+        ans += (tmp2 * s[i]).pre(deg);
+        tmp2 = (tmp2 * (t - t2)).pre(deg);
+      }
+    } else {
+      while (t2[++tmp5] == T())
+        ;
+      P tmp4 = (tmp6 >> (tmp5 - 1)).inv(deg);
+      for (int i = 0; i < b; ++i) {
+        ans += (tmp * tmp2).pre(deg) / tmp3;
+        tmp = ((tmp.diff() >> (tmp5 - 1)) * tmp4).pre(deg);
+        tmp2 = (tmp2 * (t - t2)).pre(deg);
+        tmp3 *= T(i + 1);
+      }
+    }
+    return ans;
+  }
+  //(*this)(t(x))
+  P manipulate2(P t, int deg) {
+    P ans = P();
+    P s = (*this).rev();
+    for (int i = 0; i < (int)s.size(); ++i) { ans = (ans * t + s[i]).pre(deg); }
+    return ans;
+  }
+  void debug() {
+    for (int i = 0; i < (int)(*this).size(); ++i)
+      cerr << (*this)[i] << " \n"[i == (int)(*this).size() - 1];
+  }
+};
+
+template <typename Mint>
+struct _FPS9 {
+  template <typename T>
+  T operator()(T s, T t) {
+    if (s == decltype(s)()) return decltype(s)();
+    if (t == decltype(t)()) return decltype(t)();
+    auto ntt = [](auto v, const bool& inv) {
+      const int n = v.size();
+      assert(Mint::get_mod() >= 3 && Mint::get_mod() % 2 == 1);
+      decltype(v) tmp(n, Mint());
+      Mint root = inv ? Mint(Mint::root()).inv() : Mint::root();
+      for (int b = n >> 1; b >= 1; b >>= 1, v.swap(tmp)) {
+        Mint w = root.pow((Mint::get_mod() - 1) / (n / b)), p = 1;
+        for (int i = 0; i < n; i += b * 2, p *= w)
+          for (int j = 0; j < b; ++j) {
+            v[i + j + b] *= p;
+            tmp[i / 2 + j] = v[i + j] + v[i + b + j];
+            tmp[n / 2 + i / 2 + j] = v[i + j] - v[i + b + j];
+          }
+      }
+      if (inv) v /= n;
+      return v;
+    };
+    const int n = s.size() + t.size() - 1;
+    int h = 1;
+    while ((1 << h) < n) h++;
+    s.resize((1 << h), Mint(0));
+    t.resize((1 << h), Mint(0));
+    return ntt(ntt(s, 0).dot(ntt(t, 0)), 1).pre(n);
+  }
+};
+template <typename Mint>
+using fps = FPS_BASE<Mint, _FPS9<Mint>>;
+
+class mint {
+  using u64 = std::uint_fast64_t;
+
+public:
+  u64 a;
+  constexpr mint(const long long x = 0) noexcept
+      : a(x >= 0 ? x % get_mod() : get_mod() - (-x) % get_mod()) {}
+  constexpr u64& value() noexcept { return a; }
+  constexpr const u64& value() const noexcept { return a; }
+  constexpr mint operator+(const mint rhs) const noexcept {
+    return mint(*this) += rhs;
+  }
+  constexpr mint operator-(const mint rhs) const noexcept {
+    return mint(*this) -= rhs;
+  }
+  constexpr mint operator*(const mint rhs) const noexcept {
+    return mint(*this) *= rhs;
+  }
+  constexpr mint operator/(const mint rhs) const noexcept {
+    return mint(*this) /= rhs;
+  }
+  constexpr mint& operator+=(const mint rhs) noexcept {
+    a += rhs.a;
+    if (a >= get_mod()) a -= get_mod();
+    return *this;
+  }
+  constexpr mint& operator-=(const mint rhs) noexcept {
+    if (a < rhs.a) a += get_mod();
+    a -= rhs.a;
+    return *this;
+  }
+  constexpr mint& operator*=(const mint rhs) noexcept {
+    a = a * rhs.a % get_mod();
+    return *this;
+  }
+  constexpr mint operator++(int) noexcept {
+    a += 1;
+    if (a >= get_mod()) a -= get_mod();
+    return *this;
+  }
+  constexpr mint operator--(int) noexcept {
+    if (a < 1) a += get_mod();
+    a -= 1;
+    return *this;
+  }
+  constexpr mint& operator/=(mint rhs) noexcept {
+    u64 exp = get_mod() - 2;
+    while (exp) {
+      if (exp % 2) { *this *= rhs; }
+      rhs *= rhs;
+      exp /= 2;
+    }
+    return *this;
+  }
+  constexpr bool operator==(mint x) noexcept { return a == x.a; }
+  constexpr bool operator!=(mint x) noexcept { return a != x.a; }
+  constexpr bool operator<(mint x) noexcept { return a < x.a; }
+  constexpr bool operator>(mint x) noexcept { return a > x.a; }
+  constexpr bool operator<=(mint x) noexcept { return a <= x.a; }
+  constexpr bool operator>=(mint x) noexcept { return a >= x.a; }
+  constexpr static int root() {
+    mint root = 2;
+    while (root.pow((get_mod() - 1) >> 1).a == 1) root++;
+    return root.a;
+  }
+  constexpr mint pow(long long n) {
+    long long x = a;
+    mint ret = 1;
+    while (n > 0) {
+      if (n & 1) (ret *= x);
+      (x *= x) %= get_mod();
+      n >>= 1;
+    }
+    return ret;
+  }
+  constexpr mint inv() { return pow(get_mod() - 2); }
+  static vector<mint> fac, ifac;
+  static bool init;
+  constexpr static int mx = 10000001;
+  void build() {
+    init = 0;
+    fac.resize(mx);
+    ifac.resize(mx);
+    fac[0] = 1, ifac[0] = 1;
+    for (int i = 1; i < mx; i++) fac[i] = fac[i - 1] * i;
+    ifac[mx - 1] = fac[mx - 1].inv();
+    for (int i = mx - 2; i >= 0; i--) ifac[i] = ifac[i + 1] * (i + 1);
+  }
+  mint comb(long long b) {
+    if (init) build();
+    if (a == 0 && b == 0) return 1;
+    if ((long long)a < b) return 0;
+    return fac[a] * ifac[a - b] * ifac[b];
+  }
+  mint fact() {
+    if (init) build();
+    return fac[a];
+  }
+  mint fact_inv() {
+    if (init) build();
+    return ifac[a];
+  }
+  friend ostream& operator<<(ostream& lhs, const mint& rhs) noexcept {
+    lhs << rhs.a;
+    return lhs;
+  }
+  friend istream& operator>>(istream& lhs, mint& rhs) noexcept {
+    lhs >> rhs.a;
+    return lhs;
+  }
+  constexpr static u64 get_mod() { return MOD; }
+};
+vector<mint> mint::fac;
+vector<mint> mint::ifac;
+bool mint::init = 1;
+
+int main() {
+  int n;
+  cin >> n;
+  assert(n <= 8000);
+  fps<mint> f(n), g(n);
+  for (int i = 0; i < n; ++i) cin >> f[i];
+  for (int i = 0; i < n; ++i) cin >> g[i];
+  auto p = f.manipulate(g, n);
+  p.resize(n, 0);
+  for (int i = 0; i < n; ++i) cout << p[i] << " \n"[i == n - 1];
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/sol/wa.cpp b/math/composition_of_formal_power_series_large/sol/wa.cpp
new file mode 100644
index 000000000..bcea419ab
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/sol/wa.cpp
@@ -0,0 +1,447 @@
+#include <iostream>
+#include <vector>
+#include <math.h>
+#include <assert.h>
+#include <algorithm>
+#include <cstdint>
+
+#define MOD 998244353LL
+using namespace::std;
+
+template<typename T,typename F>
+struct FPS_BASE:vector<T>{
+    using vector<T>::vector;
+    using P=FPS_BASE<T,F>;
+    F fft;
+    FPS_BASE(){}
+    inline P operator +(T x)const noexcept{return P(*this)+=x;}
+    inline P operator -(T x)const noexcept{return P(*this)-=x;}
+    inline P operator *(T x)const noexcept{return P(*this)*=x;}
+    inline P operator /(T x)const noexcept{return P(*this)/=x;}
+    inline P operator <<(int x)noexcept{return P(*this)<<=x;}
+    inline P operator >>(int x)noexcept{return P(*this)>>=x;}
+    inline P operator +(const P& x)const noexcept{return P(*this)+=x;}
+    inline P operator -(const P& x)const noexcept{return P(*this)-=x;}
+    inline P operator -()const noexcept{return P(1,T(0))-=P(*this);}
+    inline P operator *(const P& x)const noexcept{return P(*this)*=x;}
+    inline P operator /(const P& x)const noexcept{return P(*this)/=x;}
+    inline P operator %(const P& x)const noexcept{return P(*this)%=x;}
+    bool operator ==(P x){
+        for(int i=0;i<(int)max((*this).size(),x.size());++i){
+            if(i>=(int)(*this).size()&&x[i]!=T())return 0;
+            if(i>=(int)x.size()&&(*this)[i]!=T())return 0;
+            if(i<(int)min((*this).size(),x.size()))if((*this)[i]!=x[i])return 0;
+        }
+        return 1;
+    }
+    P &operator +=(T x){
+        if(this->size()==0)this->resize(1,T(0));
+        (*this)[0]+=x;
+        return (*this);
+    }
+    P &operator -=(T x){
+        if(this->size()==0)this->resize(1,T(0));
+        (*this)[0]-=x;
+        return (*this);
+    }
+    P &operator *=(T x){
+        for(int i=0;i<(int)this->size();++i){
+            (*this)[i]*=x;
+        }
+        return (*this);
+    }
+    P &operator /=(T x){
+        return (*this)*=(T(1)/x);
+    }
+    P &operator <<=(int x){
+        P ret(x,T(0));
+        ret.insert(ret.end(),begin(*this),end(*this));
+        return (*this)=ret;
+    }
+    P &operator >>=(int x){
+        P ret;
+        ret.insert(ret.end(),begin(*this)+x,end(*this));
+        return (*this)=ret;
+    }
+    P &operator +=(const P& x){
+        if(this->size()<x.size())this->resize(x.size(),T(0));
+        for(int i=0;i<(int)x.size();++i){
+            (*this)[i]+=x[i];
+        }
+        return (*this);
+    }
+    P &operator -=(const P& x){
+        if(this->size()<x.size())this->resize(x.size(),T(0));
+        for(int i=0;i<(int)x.size();++i){
+            (*this)[i]-=x[i];
+        }
+        return (*this);
+    }
+    P &operator *=(const P& x){
+        return (*this)=F()(*this,x);
+    }
+    P &operator /=(P x){
+        if(this->size()<x.size()) {
+            this->clear();
+            return (*this);
+        }
+        const int n=this->size()-x.size()+1;
+        return (*this) = (rev().pre(n)*x.rev().inv(n)).pre(n).rev(n);
+    }
+    P &operator %=(const P& x){
+        return ((*this)-=*this/x*x);
+    }
+    inline void print(){
+        for(int i=0;i<(int)(*this).size();++i)cerr<<(*this)[i]<<" \n"[i==(int)(*this).size()-1];
+        if((int)(*this).size()==0)cerr<<endl;
+    }
+    inline P& shrink(){while((*this).back()==0)(*this).pop_back();return (*this);}
+    inline P pre(int sz)const{
+        return P(begin(*this),begin(*this)+min((int)this->size(),sz));
+    }
+    P rev(int deg=-1){
+        P ret(*this);
+        if(deg!=-1)ret.resize(deg,T(0));
+        reverse(begin(ret),end(ret));
+        return ret;
+    }
+    P inv(int deg=-1){
+        assert((*this)[0]!=T(0));
+        const int n=deg==-1?this->size():deg;
+        P ret({T(1)/(*this)[0]});
+        for(int i=1;i<n;i<<=1){
+            ret*=(-ret*pre(i<<1)+2).pre(i<<1);
+        }
+        return ret.pre(n);
+    }
+    inline P dot(const P& x){
+        P ret(*this);
+        for(int i=0;i<int(min(this->size(),x.size()));++i){
+            ret[i]*=x[i];
+        }
+        return ret;
+    }
+    P diff(){
+        if((int)(*this).size()<=1)return P();
+        P ret(*this);
+        for(int i=0;i<(int)ret.size();i++){
+            ret[i]*=i;
+        }
+        return ret>>1;
+    }
+    P integral(){
+        P ret(*this);
+        for(int i=0;i<(int)ret.size();i++){
+            ret[i]/=i+1;
+        }
+        return ret<<1;
+    }
+    P log(int deg=-1){
+        assert((*this)[0]==T(1));
+        const int n=deg==-1?this->size():deg;
+        return (diff()*inv(n)).pre(n-1).integral();
+    }
+    P exp(int deg=-1){
+        assert((*this)[0]==T(0));
+        const int n=deg==-1?this->size():deg;
+        P ret({T(1)});
+        for(int i=1;i<n;i<<=1){
+            ret=ret*(pre(i<<1)+1-ret.log(i<<1)).pre(i<<1);
+        }
+        return ret.pre(n);
+    }
+    P pow(int c,int deg=-1){
+        // const int n=deg==-1?this->size():deg;
+        // long long i=0;
+        // P ret(*this);
+        // while(i!=(int)this->size()&&ret[i]==0)i++;
+        // if(i==(int)this->size())return P(n,0);
+        // if(i*c>=n)return P(n,0);
+        // T k=ret[i];
+        // return ((((ret>>i)/k).log()*c).exp()*(k.pow(c))<<(i*c)).pre(n);
+        P x(*this);
+        P ret(1,1);
+        while(c) {
+            if(c&1){
+                ret*=x;
+                if(~deg)ret=ret.pre(deg);
+            }
+            x*=x;
+            if(~deg)x=x.pre(deg);
+            c>>=1;
+        }
+        return ret;
+    }
+    P sqrt(int deg=-1){
+        const int n=deg==-1?this->size():deg;
+        if((*this)[0]==T(0)) {
+            for(int i=1;i<(int)this->size();i++) {
+                if((*this)[i]!=T(0)) {
+                    if(i&1)return{};
+                    if(n-i/2<=0)break;
+                    auto ret=(*this>>i).sqrt(n-i/2)<<(i/2);
+                    if((int)ret.size()<n)ret.resize(n,T(0));
+                    return ret;
+                }
+            }
+            return P(n,0);
+        }
+        P ret({T(1)});
+        for(int i=1;i<n;i<<=1){
+            ret=(ret+pre(i<<1)*ret.inv(i<<1)).pre(i<<1)/T(2);
+        }
+        return ret.pre(n);
+    }
+    P shift(int c){
+        const int n=this->size();
+        P f(*this),g(n,0);
+        for(int i=0;i<n;++i)f[i]*=T(i).fact();
+        for(int i=0;i<n;++i)g[i]=T(c).pow(i)/T(i).fact();
+        g=g.rev();
+        f*=g;
+        f>>=n-1;
+        for(int i=0;i<n;++i)f[i]/=T(i).fact();
+        return f;
+    }
+    T eval(T x){
+        T res=0;
+        for(int i=(int)this->size()-1;i>=0;--i){
+            res*=x;
+            res+=(*this)[i];
+        }
+        return res;
+    }
+    vector<T> multipoint_eval(const vector<T>&x){
+        const int n=x.size();
+        P* v=new P[2*n-1];
+        for(int i=0;i<n;i++)v[i+n-1]={T()-x[i],T(1)};
+        for(int i=n-2;i>=0;i--){v[i]=v[i*2+1]*v[i*2+2];}
+        v[0]=P(*this)%v[0];v[0].shrink();
+        for(int i=1;i<n*2-1;i++){
+            v[i]=v[(i-1)/2]%v[i];
+            v[i].shrink();
+        }
+        vector<T>res(n);
+        for(int i=0;i<n;i++)res[i]=v[i+n-1][0];
+        return res;
+    }
+    P slice(int s,int e,int k){
+        P res;
+        for(int i=s;i<e;i+=k)res.push_back((*this)[i]);
+        return res;
+    }
+    T nth_term(P q,int64_t x){
+        if(x==0)return (*this)[0]/q[0];
+        P p(*this);
+        P q2=q;
+        for(int i=1;i<(int)q2.size();i+=2)q2[i]*=-1;
+        q*=q2;
+        p*=q2;
+        return p.slice(x%2,p.size(),2).nth_term(q.slice(0,q.size(),2),x/2);
+    }
+    
+    //(*this)(t(x))
+    P manipulate(P t,int deg){
+        P s=P(*this);
+        if(deg==0)return P();
+        if((int)t.size()==1)return P{s.eval(t[0])};
+        int k=min((int)::sqrt(deg/(::log2(deg)+1))+1,(int)t.size());
+        int b=deg/k+1;
+        P t2=t.pre(k);
+        vector<P>table(s.size()/2+1,P{1});
+        for(int i=1;i<(int)table.size();i++){
+            table[i]=((table[i-1])*t2).pre(deg);
+        }
+        auto f=[&](auto f,auto l,auto r,int deg)->P{
+            if(r-l==1)return P{*l};
+            auto m=l+(r-l)/2;
+            return f(f,l,m,deg)+(table[m-l]*f(f,m,r,deg)).pre(deg);
+        };
+        P ans=P();
+        P tmp=f(f,s.begin(),s.end(),deg);
+        P tmp2=P{1};
+        T tmp3=T(1);
+        P tmp4=t2.diff().inv(deg);
+        for(int i=0;i<b;++i){
+            ans+=(tmp*tmp2).pre(deg)/tmp3;
+            tmp=(tmp.diff()*tmp4).pre(deg);
+            tmp2=(tmp2*(t-t2)).pre(deg);
+            tmp3*=T(i+1);
+        }
+        return ans;
+    }
+    //(*this)(t(x))
+    P manipulate2(P t,int deg){
+        P ans=P();
+        P s=(*this).rev();
+        for(int i=0;i<(int)s.size();++i){
+            ans=(ans*t+s[i]).pre(deg);
+        }
+        return ans;
+    }
+};
+
+template<typename Mint>
+struct _FPS9{
+    template<typename T>
+    T operator()(T s,T t){
+        if(s==decltype(s)())return decltype(s)();
+        if(t==decltype(t)())return decltype(t)();
+        auto ntt=[](auto v,const bool& inv){
+            const int n=v.size();
+            assert(Mint::get_mod()>=3&&Mint::get_mod()%2==1);
+            decltype(v) tmp(n,Mint());
+            Mint root=inv?Mint(Mint::root()).inv():Mint::root();
+            for(int b=n>>1;b>=1;b>>=1,v.swap(tmp)){
+                Mint w=root.pow((Mint::get_mod()-1)/(n/b)),p=1;
+                for(int i=0;i<n;i+=b*2,p*=w)for(int j=0;j<b;++j){
+                    v[i+j+b]*=p;
+                    tmp[i/2+j]=v[i+j]+v[i+b+j];
+                    tmp[n/2+i/2+j]=v[i+j]-v[i+b+j];
+                }
+            }
+            if(inv)v/=n;
+            return v;
+        };
+        const int n=s.size()+t.size()-1;
+        int h=1;
+        while((1<<h)<n)h++;
+        s.resize((1<<h),Mint(0));
+        t.resize((1<<h),Mint(0));
+        return ntt(ntt(s,0).dot(ntt(t,0)),1).pre(n);
+    }
+};
+template<typename Mint>using fps=FPS_BASE<Mint,_FPS9<Mint>>;
+
+class mint {
+  using u64 = std::uint_fast64_t;
+    public:
+    u64 a;
+    constexpr mint(const long long x = 0)noexcept:a(x>=0?x%get_mod():get_mod()-(-x)%get_mod()){}
+    constexpr u64 &value()noexcept{return a;}
+    constexpr const u64 &value() const noexcept {return a;}
+    constexpr mint operator+(const mint rhs)const noexcept{return mint(*this) += rhs;}
+    constexpr mint operator-(const mint rhs)const noexcept{return mint(*this)-=rhs;}
+    constexpr mint operator*(const mint rhs) const noexcept {return mint(*this) *= rhs;}
+    constexpr mint operator/(const mint rhs) const noexcept {return mint(*this) /= rhs;}
+    constexpr mint &operator+=(const mint rhs) noexcept {
+        a += rhs.a;
+        if (a >= get_mod())a -= get_mod();
+        return *this;
+    }
+    constexpr mint &operator-=(const mint rhs) noexcept {
+        if (a<rhs.a)a += get_mod();
+        a -= rhs.a;
+        return *this;
+    }
+    constexpr mint &operator*=(const mint rhs) noexcept {
+        a = a * rhs.a % get_mod();
+        return *this;
+    }
+    constexpr mint operator++(int) noexcept {
+        a += 1;
+        if (a >= get_mod())a -= get_mod();
+        return *this;
+    }
+    constexpr mint operator--(int) noexcept {
+        if (a<1)a += get_mod();
+        a -= 1;
+        return *this;
+    }
+    constexpr mint &operator/=(mint rhs) noexcept {
+        u64 exp=get_mod()-2;
+        while (exp) {
+            if (exp % 2) {
+                *this *= rhs;
+            }
+            rhs *= rhs;
+            exp /= 2;
+        }
+        return *this;
+    }
+    constexpr bool operator==(mint x) noexcept {
+        return a==x.a;
+    }
+    constexpr bool operator!=(mint x) noexcept {
+        return a!=x.a;
+    }
+    constexpr bool operator<(mint x) noexcept {
+        return a<x.a;
+    }
+    constexpr bool operator>(mint x) noexcept {
+        return a>x.a;
+    }
+    constexpr bool operator<=(mint x) noexcept {
+        return a<=x.a;
+    }
+    constexpr bool operator>=(mint x) noexcept {
+        return a>=x.a;
+    }
+    constexpr static int root(){
+        mint root = 2;
+        while(root.pow((get_mod()-1)>>1).a==1)root++;
+        return root.a;
+    }
+    constexpr mint pow(long long n){
+        long long x=a;
+        mint ret = 1;
+        while(n>0) {
+            if(n&1)(ret*=x);
+            (x*=x)%=get_mod();
+            n>>=1;
+        }
+        return ret;
+    }
+    constexpr mint inv(){
+        return pow(get_mod()-2);
+    }
+    static vector<mint> fac,ifac;
+    static bool init;
+    constexpr static int mx=10000001;
+    void build(){
+        init=0;
+        fac.resize(mx);
+        ifac.resize(mx);
+        fac[0]=1,ifac[0]=1;
+        for(int i=1;i<mx;i++)fac[i]=fac[i-1]*i;
+        ifac[mx-1]=fac[mx-1].inv();
+        for(int i=mx-2;i>=0;i--)ifac[i]=ifac[i+1]*(i+1);
+    }
+    mint comb(long long b){
+        if(init)build();
+        if(a==0&&b==0)return 1;
+        if((long long)a<b)return 0;
+        return fac[a]*ifac[a-b]*ifac[b];
+    }
+    mint fact(){
+        if(init)build();
+        return fac[a];
+    }
+    mint fact_inv(){
+        if(init)build();
+        return ifac[a];
+    }
+    friend ostream& operator<<(ostream& lhs, const mint& rhs) noexcept {
+        lhs << rhs.a;
+        return lhs;
+    }
+    friend istream& operator>>(istream& lhs,mint& rhs) noexcept {
+        lhs >> rhs.a;
+        return lhs;
+    }
+    constexpr static u64 get_mod(){return MOD;}
+};
+vector<mint> mint::fac;
+vector<mint> mint::ifac;
+bool mint::init=1;
+
+int main(){
+    int n;
+    cin>>n;
+    fps<mint>f(n),g(n);
+    for(int i=0;i<n;++i)cin>>f[i];
+    for(int i=0;i<n;++i)cin>>g[i];
+    auto h=f.manipulate(g,n);
+    h.resize(n,0);
+    for(int i=0;i<n;++i)cout<<h[i]<<" \n"[i==n-1];
+}
\ No newline at end of file
diff --git a/math/composition_of_formal_power_series_large/task.md b/math/composition_of_formal_power_series_large/task.md
new file mode 100644
index 000000000..4a6c856d9
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/task.md
@@ -0,0 +1,40 @@
+## @{keyword.statement}
+
+@{lang.en}
+You are given two formal power series $f(x) = \sum_{i=0}^{N-1} a_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ and $g(x) = \sum_{i=0}^{N-1} b_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ with $b_0 \ne 0$.
+Calculate the first $N$ terms of $f(g(x)) = \sum_{i=0}^{\infty} c_i x^i$.
+In other words, find $h(x) = \sum_{i=0}^{N-1} c_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ such that
+
+$$h(x) \equiv \sum_{i=0}^{N-1} a_i g(x)^i \pmod{x^N}.$$
+
+@{lang.ja}
+形式的冪級数 $f(x) = \sum_{i=0}^{N-1} a_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ と $g(x) = \sum_{i=0}^{N-1} b_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ が与えられます ($b_0 = 0$)。$f(g(x)) = \sum_{i=0}^{\infty} c_i x^i$ の先頭 $N$ 項を求めてください。つまり、
+
+$$h(x) \equiv \sum_{i=0}^{N-1} a_i g(x)^i \pmod{x^N}$$
+
+となる $h(x) = \sum_{i=0}^{N-1} c_i x^i \in \mathbb{F}\_{@{param.MOD}}[[x]]$ を求めてください。
+
+@{lang.end}
+
+## @{keyword.constraints}
+
+- $1 \leq N \leq @{param.N_MAX}$
+- $0 \leq a_i, b_i < @{param.MOD}$
+- $b_0 = 0$
+## @{keyword.input}
+
+```
+$N$
+$a_0$ $a_1$ $\cdots$ $a_{N - 1}$
+$b_0$ $b_1$ $\cdots$ $b_{N - 1}$
+```
+
+## @{keyword.output}
+
+```
+$c_0$ $c_1$ $\cdots$ $c_{N - 1}$
+```
+
+## @{keyword.sample}
+
+@{example.example_00}
diff --git a/math/composition_of_formal_power_series_large/verifier.cpp b/math/composition_of_formal_power_series_large/verifier.cpp
new file mode 100644
index 000000000..fcbc8ea97
--- /dev/null
+++ b/math/composition_of_formal_power_series_large/verifier.cpp
@@ -0,0 +1,25 @@
+#include <iostream>
+
+#include "params.h"
+#include "testlib.h"
+
+int main() {
+  registerValidation();
+  int n = inf.readInt(1, N_MAX);
+  inf.readChar('\n');
+  for (int i = 0; i < n; i++) {
+    inf.readInt(0, MOD - 1);
+    if (i != n - 1) inf.readSpace();
+  }
+  inf.readChar('\n');
+  for (int i = 0; i < n; i++) {
+    if (i)
+      inf.readInt(0, MOD - 1);
+    else
+      inf.readInt(0, 0);
+    if (i != n - 1) inf.readSpace();
+  }
+  inf.readChar('\n');
+  inf.readEof();
+  return 0;
+}