結果
| 問題 | No.3770 Sum of Pow of Distinct Count |
| コンテスト | |
| ユーザー |
 2251799813685248
|
| 提出日時 | 2026-10-18 13:33:07 |
| 言語 | C++23 (gcc 15.3.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 425 ms / 4,000 ms |
| + 306µs | |
| コード長 | 127,287 bytes |
| 記録 | |
| コンパイル時間 | 5,763 ms |
| コンパイル使用メモリ | 424,140 KB |
| 実行使用メモリ | 142,704 KB |
| 最終ジャッジ日時 | 2026-10-18 13:33:19 |
| 合計ジャッジ時間 | 12,038 ms |
|
ジャッジサーバーID （参考情報） |
judge5_0 / judge1_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 15 |
ソースコード
// #define _GLIBCXX_DEBUG
#ifndef LIST__HPP_
#define LIST__HPP_
#include <algorithm>
#include <iostream>
#include <cassert>
#include <functional>
#include <utility>
using uint = unsigned int;
using ll = long long;
using namespace std;
template<typename T> concept has_identity = requires (){T();};
template<typename T> concept has_plus_operator = requires (T a, T b){a+b;};
template<typename T> concept has_compare_operator = requires (T a, T b){a < b;};
template<typename T> concept has_equal_operator = requires (T a, T b){a == b; a != b;};
template<typename T>
class list{
public:
using value_type = T;
constexpr list(): sz_(0), max_sz_(0), A(nullptr){}
constexpr list(uint sz): sz_(sz), max_sz_(sz){
if (sz == 0){A = nullptr; return;}
if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
A = new T[sz_]{};
if (A == nullptr){cerr << "MomoryError: failed to allocate memory" << endl; assert(false);}
}
constexpr list(uint sz, T elem): sz_(sz), max_sz_(sz){
if (sz == 0){A = nullptr; return;}
if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
A = new T[sz_]{};
if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
for (int i = 0; i < sz_; i++){
A[i] = elem;
}
}
constexpr list(initializer_list<T> init): sz_(init.size()), max_sz_(init.size()){
A = new T[sz_]{};
if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
for (int i = 0; i < sz_; i++){
A[i] = *(init.begin()+i);
}
};
// constexpr list<ll>(const range& r){
// if (r.step_ == 0){cerr << "ValueError: step argument must not be zero" << endl; assert(false);}
// if (r.step_ > 0){
// int imax = r.begin_ >= r.end_ ? 0 : (r.end_-r.begin_+r.step_-1)/r.step_;
// sz_ = imax;
// max_sz_ = imax;
// A = new T[sz_]{};
// if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
// for (int i = 0; i < imax; i++){
// int target = r.begin_+r.step_*i;
// A[i] = target;
// }
// }
// else{
// int imax = r.end_ >= r.begin_ ? 0 : (r.begin_-r.end_-r.step_-1)/(-r.step_);
// sz_ = imax;
// max_sz_ = imax;
// A = new T[sz_]{};
// if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
// for (int i = 0; i < imax; i++){
// int target = r.begin_+r.step_*i;
// A[i] = target;
// }
// }
// }
constexpr list(const list<T>& other): sz_(other.sz_), max_sz_(other.max_sz_){
A = new T[max_sz_]{};
for (int i = 0; i < sz_; i++){
A[i] = other.A[i];
}
}
constexpr list(list<T>&& other) noexcept : sz_(other.sz_), max_sz_(other.max_sz_){
A = other.A;
other.sz_ = 0;
other.max_sz_ = 0;
other.A = nullptr;
}
constexpr ~list() {
delete[] A;
}
class iterator{
public:
using difference_type = ll;
using value_type = T;
using pointer = T*;
using reference = T&;
using iterator_category = std::random_access_iterator_tag;
iterator() = delete;
constexpr iterator(const list* listptr_init, const int &idx_init): listptr(listptr_init), idx(idx_init){}
inline constexpr value_type& operator* () const {return listptr->A[this->check_idx()];}
inline constexpr value_type* operator->() const {return &listptr->A[this->check_idx()];}
inline constexpr bool operator==(const iterator &other) const {return idx == other.idx;}
inline constexpr bool operator!=(const iterator &other) const {return idx != other.idx;}
inline constexpr bool operator<(const iterator &other) const {return idx < other.idx;}
inline constexpr bool operator>(const iterator &other) const {return idx > other.idx;}
inline constexpr bool operator<=(const iterator &other) const {return idx <= other.idx;}
inline constexpr bool operator>=(const iterator &other) const {return idx >= other.idx;}
inline constexpr iterator operator++(){
this->idx++;
return *this;
}
inline constexpr iterator operator++(int){
iterator temp = *this;
this->idx++;
return temp;
}
inline constexpr iterator operator--(){
this->idx--;
return *this;
}
inline constexpr iterator operator--(int){
iterator temp = *this;
this->idx--;
return temp;
}
inline constexpr iterator& operator+=(int amount){
this->idx += amount;
return *this;
}
inline constexpr iterator operator+(int amount) const {
iterator temp = *this;
temp.idx += amount;
return temp;
}
inline constexpr iterator& operator-=(const int &amount){
this->idx -= amount;
return *this;
}
inline constexpr iterator operator-(const int &amount) const{
iterator temp = *this;
temp.idx -= amount;
return temp;
}
inline constexpr difference_type operator-(const iterator &other) const{
return idx - other.idx;
}
inline constexpr int check_idx() const {
if (idx < 0 || idx >= listptr->sz_){
cerr << "IteratorIndexError: index out of range. current size: " << listptr->sz_ << " target index: " << idx << endl; assert(false);
}
return idx;
}
private:
const list* listptr;
int idx;
};
inline constexpr iterator begin(){return iterator(this, 0);}
inline constexpr iterator begin() const {return iterator(this, 0);}
inline constexpr iterator end(){return iterator(this, sz_);}
inline constexpr iterator end() const {return iterator(this, sz_);}
inline constexpr list& operator=(const list& other){
if (this == &other){return *this;}
sz_ = other.sz_;
max_sz_ = other.max_sz_;
delete[] A;
A = new T[max_sz_]{};
for (int i = 0; i < sz_; i++){
A[i] = other.A[i];
}
return *this;
}
inline constexpr list& operator=(list&& other) noexcept {
if (this == &other){return *this;}
delete[] A;
sz_ = other.sz_;
max_sz_ = other.max_sz_;
A = other.A;
other.sz_ = 0;
other.max_sz_ = 0;
other.A = nullptr;
return *this;
}
inline constexpr bool operator==(const list& other) const {
if (sz_ != other.sz_){return false;}
bool res = true;
for (int i = 0; i < sz_; i++){
if (A[i] != other.A[i]){
res = false;
break;
}
}
return res;
}
inline constexpr bool operator!=(const list& other) const {
return !this->operator==(other);
}
inline constexpr bool operator<(const list& other) const {
int lcp_back = 0;
for (int i = 0, sz2 = min(sz_, other.sz_); i < sz2; i++){
if (A[i] != other.A[i]){break;}
lcp_back++;
}
if (lcp_back == sz_ && lcp_back == other.sz_){return false;}
else if (lcp_back == sz_){return true;}
else if (lcp_back == other.sz_){return false;}
else{return A[lcp_back] < other.A[lcp_back];}
}
inline constexpr bool operator<=(const list& other) const {
int lcp_back = 0;
for (int i = 0, sz2 = min(sz_, other.sz_); i < sz2; i++){
if (A[i] != other.A[i]){break;}
lcp_back++;
}
if (lcp_back == sz_ && lcp_back == other.sz_){return true;}
else if (lcp_back == sz_){return true;}
else if (lcp_back == other.sz_){return false;}
else{return A[lcp_back] < other.A[lcp_back];}
}
inline constexpr bool operator>(const list& other) const {return other < *this;}
inline constexpr bool operator>=(const list& other) const {return other <= *this;}
inline constexpr list& operator+=(const list& other){
for (int i = 0, sz2 = other.size(); i < sz2; i++){
append(other.A[i]);
}
return *this;
}
inline constexpr list operator+(const list& other) const {
return list(*this) += other;
}
inline constexpr list& operator+=(const T& val){
append(val);
return *this;
}
inline constexpr list operator+(const T& val) const {
return list(*this) += val;
}
inline constexpr void append(const T& val){
if (sz_ < max_sz_){
A[sz_] = val;
sz_++;
}
else{
if (max_sz_ < 2){max_sz_ = 2;}
else{max_sz_ += max_sz_>>1;}
auto B = new T[max_sz_]{};
B[sz_] = val;
for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
sz_++;
delete[] A;
A = B;
}
}
inline constexpr void append(T&& val){
if (sz_ < max_sz_){
A[sz_] = std::move(val);
sz_++;
}
else{
if (max_sz_ < 2){max_sz_ = 2;}
else{max_sz_ += max_sz_>>1;}
auto B = new T[max_sz_]{};
B[sz_] = std::move(val);
for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
sz_++;
delete[] A;
A = B;
}
}
inline constexpr T pop(){
if (sz_ == 0){cerr << "IndexError: pop from empty list" << endl; assert(false);}
sz_--;
return std::move(A[sz_]);
}
inline constexpr void reserve(uint new_sz){
if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
if ((int)new_sz <= max_sz_){return;}
max_sz_ = new_sz;
auto B = new T[new_sz]{};
if (B == nullptr){cerr << "MemoryError: failed to reallocate memory" << endl; assert(false);}
for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
delete[] A;
A = B;
}
inline constexpr void resize(uint new_sz, T elem = T()){
if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
if ((int)new_sz <= max_sz_){
while (sz_ > (int)new_sz){pop();}
for (int i = sz_, new_sz_int = new_sz; i < new_sz_int; i++){A[i] = elem;}
sz_ = new_sz;
return;
}
max_sz_ = new_sz + (new_sz >> 1);
auto B = new T[max_sz_]{};
if (B == nullptr){cerr << "MemoryError: failed to reallocate memory" << endl; assert(false);}
for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
for (int i = sz_, new_sz_int = new_sz; i < new_sz_int; i++){B[i] = elem;}
sz_ = new_sz;
delete[] A;
A = B;
}
template <class InputIterator> inline constexpr void assign(InputIterator first, InputIterator last){
reserve(last-first);
while (first != last){append(*first); first++;}
}
inline constexpr void assign(uint n, const T& val){
reserve(n);
sz_ = n;
for (uint i = 0; i < n; i++){
A[i] = val;
}
}
inline constexpr iterator insert(iterator pos, const T& val){
int idx = pos.check_idx();
append(val);
for (int i = idx; i < sz_-1; i++){swap(A[i], A[sz_-1]);}
return pos;
}
// inline constexpr void insert(iterator pos, iterator begin, iterator end)
inline constexpr iterator erase(iterator pos){
int idx = pos.check_idx();
for (int i = idx; i < sz_-1; i++){
swap(A[i], A[i+1]);
}
sz_ -= 1;
return pos;
}
inline constexpr iterator erase(iterator start, iterator end){
int startidx = start.check_idx(), endidx = end.check_idx();
if (startidx > endidx){cerr << "IndexError: invalid interval: [start, end) = []" << startidx << ", " << endidx << ")" << endl; assert(false);}
for (int i = startidx; i < endidx; i++){
A[i].~T();
}
if (startidx < endidx){
for (int i = endidx; i < sz_; i++){
swap(A[i+startidx-endidx], A[i]);
}
}
sz_ -= endidx-startidx;
return end;
}
inline constexpr void clear(){sz_ = 0; max_sz_ = 0; delete[] A;}
inline constexpr const T& operator[](int idx) const {
check_and_adjust_idx(idx);
return A[idx];
}
inline constexpr T& operator[](int idx){
check_and_adjust_idx(idx);
return A[idx];
}
inline constexpr T& front(){return (*this)[0];}
inline constexpr const T& front() const {return (*this)[0];}
inline constexpr T& back(){return (*this)[-1];}
inline constexpr const T& back() const {return (*this)[-1];}
inline constexpr list<T> operator[](int left, int right, int step = 1) const {
if (step == 0){cerr << "ValueError: step argument must not be zero" << endl; assert(false);}
if (step > 0){
int imax = left >= right ? 0 : (right-left+step-1)/step;
list<T> B(imax);
for (int i = 0; i < imax; i++){
int target = left+step*i;
if (target < 0 || target >= sz_){cerr << "IndexError: slice index out of range. current size: " << sz_ << " target index: " << target << endl; assert(false);}
B[i] = A[target];
}
return std::move(B);
}
else{
int imax = right >= left ? 0 : (left-right-step-1)/(-step);
list<T> B(imax);
for (int i = 0; i < imax; i++){
int target = left+step*i;
if (target < 0 || target >= sz_){cerr << "IndexError: slice index out of range. current size: " << sz_ << " target index: " << target << endl; assert(false);}
B[i] = A[target];
}
return std::move(B);
}
}
inline constexpr T* data() {return empty() ? nullptr : &A[0];}
inline constexpr const T* data() const {return empty() ? nullptr : &A[0];}
inline constexpr T sum() const requires has_plus_operator<T> && has_identity<T>{
auto res = T();
for (int i = 0; i < sz_; i++){
res += A[i];
}
return res;
}
template<typename U> inline constexpr U sum(U init) const requires has_plus_operator<T>{
for (int i = 0; i < sz_; i++){
init += A[i];
}
return init;
}
inline constexpr void sort(bool reverse = false) requires has_compare_operator<T> {
if (reverse){
std::sort(begin(), end(), std::move([](const T& a, const T& b){return b < a;}));
}
else{
std::sort(begin(), end());
}
}
inline constexpr void sort(const function<bool(T,T)>& f) requires has_compare_operator<T> {
std::sort(begin(), end(), f);
}
inline constexpr void reverse(){
std::reverse(begin(), end());
}
inline constexpr ll count(const T& val){
ll res = 0;
for (int i = 0; i < sz_; i++){res += A[i] == val;}
return res;
}
inline constexpr void apply(const function<T(T)>& f){for (int i = 0; i < sz_; i++){A[i] = f(A[i]);}}
inline constexpr list<pair<int, T>> RLE(){
list<pair<int, T>> res;
for (int l = 0, r = 0; l < sz_;){
while (r < sz_ && A[l] == A[r]){
r++;
}
res.append(make_pair(r-l, A[l]));
l = r;
}
return res;
}
inline constexpr iterator lower_bound(const T& val) const requires has_compare_operator<T> {return std::lower_bound(begin(), end(), val);}
inline constexpr iterator upper_bound(const T& val) const requires has_compare_operator<T> {return std::upper_bound(begin(), end(), val);}
inline constexpr ll size() const {return sz_;}
inline constexpr ll maxsize() const {return max_sz_;}
inline constexpr bool empty() const {return sz_ == 0;}
private:
inline constexpr void check_and_adjust_idx(int& idx) const {
if (idx >= sz_){cerr << "IndexError: index out of range. current size: " << sz_ << " target index: " << idx << endl; assert(false);}
if (idx < 0){
if (idx < -sz_){cerr << "IndexError: index out of range. current size: " << sz_ << " target index: " << idx << endl; assert(false);}
idx += sz_;
}
}
int sz_;
int max_sz_;
T* A;
};
#endif /* LIST__HPP_ */
#include <iostream>
#include <iterator>
#include <type_traits>
#include <utility>
#include <string>
#include <cmath>
#include <unordered_set>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <iomanip>
#include <cassert>
#include <functional>
#include <random>
#include <bitset>
// #include <boost/multiprecision/cpp_int.hpp>
using namespace std;
using ll = long long;
using lll = __int128_t;
using ull = unsigned long long;
#ifdef BOOST_VERSION
using bll = boost::multiprecision::cpp_int;
#endif
using ld = long double;
using pii = array<int,2>;
using pll = array<ll,2>;
using plll = array<lll,2>;
using pld = array<ld, 2>;
using tll = array<ll,3>;
#define vall(A) A.begin(), A.end()
template<typename T> inline bool btest(T K, int i){return K&(1ull<<i);}
#ifndef LEN_DEFINED
#define LEN_DEFINED
template<typename T> inline constexpr ll len(const T& A){return ssize(A);}
#endif
constexpr ll pow2ll[63] = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768,65536,131072,262144,524288,1048576,2097152,4194304,8388608,16777216,33554432,67108864,134217728,268435456,536870912,1073741824,2147483648,4294967296,8589934592,17179869184,34359738368,68719476736,137438953472,274877906944,549755813888,1099511627776,2199023255552,4398046511104,8796093022208,17592186044416,35184372088832,70368744177664,140737488355328,281474976710656,562949953421312,1125899906842624,2251799813685248,4503599627370496,9007199254740992,18014398509481984,36028797018963968,72057594037927936,144115188075855872,288230376151711744,576460752303423488,1152921504606846976,2305843009213693952,4611686018427387904};
constexpr ll pow10ll[19] = {1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000};
constexpr ll di[4] = {0,1,0,-1};
constexpr ll di8[8] = {0,1,1,1,0,-1,-1,-1};
constexpr ll dj[4] = {1,0,-1,0};
constexpr ll dj8[8] = {1,1,0,-1,-1,-1,0,1};
#ifndef FASTIO__HPP_
#define FASTIO__HPP_
#include <array>
#include <bit>
#include <charconv>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <iterator>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>
// C++20. Integer input must fit the destination type. For modint input, the
// magnitude must fit uint64_t; negative values are reduced modulo the modulus.
// FILE* is borrowed: it must outlive the reader/writer. Do not mix buffered
// readers on one stream. The default input buffer is intended for batch judges;
// interactive programs need BasicInput<1> and explicit output flush().
namespace fastio {
namespace detail {
template<class T> inline constexpr bool integer = std::is_integral_v<T>;
template<class T> struct unsigned_type { using type = std::make_unsigned_t<T>; };
template<> struct unsigned_type<bool> { using type = unsigned char; };
#if defined(__SIZEOF_INT128__)
template<> inline constexpr bool integer<__int128_t> = true;
template<> inline constexpr bool integer<__uint128_t> = true;
template<> struct unsigned_type<__int128_t> { using type = __uint128_t; };
template<> struct unsigned_type<__uint128_t> { using type = __uint128_t; };
#endif
template<class T> using unsigned_t = typename unsigned_type<T>::type;
template<class T> inline constexpr bool signed_integer = integer<T> && T(-1) < T(0);
template<class T> concept modular = requires(T x) { x.val; } &&
(requires { T::mod_value; } || requires { T::mod(); });
template<class T> constexpr auto modulus() {
if constexpr (requires { T::mod_value; }) return T::mod_value;
else return T::mod();
}
template<class T> concept indexed_range = requires(T& x) { std::size(x); x[0]; };
template<class T> concept tuple_like = requires { std::tuple_size<T>::value; };
template<class> inline constexpr bool unsupported = false;
inline constexpr auto digit_pairs = [] {
std::array<char, 200> result{};
for (int i = 0; i < 100; ++i) {
result[2*i] = char('0'+i/10);
result[2*i+1] = char('0'+i%10);
}
return result;
}();
inline constexpr auto digit_quads = [] {
std::array<std::array<char,4>,10000> result{};
for (unsigned i = 0; i < 10000; ++i) {
result[i] = {char('0'+i/1000),char('0'+i/100%10),
char('0'+i/10%10),char('0'+i%10)};
}
return result;
}();
inline char* leading_quad(char* p, unsigned value) {
if (value >= 1000) { std::memcpy(p,digit_quads[value].data(),4); return p+4; }
if (value >= 100) { *p++ = char('0'+value/100); value %= 100; }
else if (value < 10) { *p++ = char('0'+value); return p; }
std::memcpy(p,digit_pairs.data()+2*value,2);
return p+2;
}
inline char* decimal32(char* p, std::uint32_t value) {
if (value >= 100000000) {
p = leading_quad(p,value/100000000);
value %= 100000000;
std::memcpy(p,digit_quads[value/10000].data(),4);
std::memcpy(p+4,digit_quads[value%10000].data(),4);
return p+8;
}
if (value >= 10000) {
p = leading_quad(p,value/10000);
std::memcpy(p,digit_quads[value%10000].data(),4);
return p+4;
}
return leading_quad(p,value);
}
template<class T> inline constexpr bool output32 = [] {
if constexpr (integer<T>) return sizeof(T) <= 4 && !std::is_same_v<T,char>;
else if constexpr (modular<T>) return integer<decltype(T{}.val)> && sizeof(T{}.val) <= 4;
else return false;
}();
template<class T> inline char* decimal_small(char* p, const T& value) {
if constexpr (modular<T>) return decimal_small(p,value.val);
else {
std::uint32_t u = static_cast<std::uint32_t>(value);
if constexpr (signed_integer<T>) {
if (value < 0) { *p++ = '-'; u = 0-u; }
}
return decimal32(p,u);
}
}
} // namespace detail
template<std::size_t BufferSize = 1 << 18>
class BasicInput {
static_assert(BufferSize > 0);
std::FILE* file_;
char buffer_[BufferSize+1];
std::size_t position_ = 0, size_ = 0;
bool good_ = true;
int peek() {
if (position_ == size_) {
size_ = std::fread(buffer_, 1, BufferSize, file_);
buffer_[size_] = '\0';
position_ = 0;
if (size_ == 0) return EOF;
}
return static_cast<unsigned char>(buffer_[position_]);
}
int skip_space() {
int c;
while ((c = peek()) != EOF && c <= ' ') ++position_;
return c;
}
template<class U, bool SingleBlock = false> bool magnitude(U& value, bool& negative) {
int c = skip_space();
if (c == EOF) return false;
negative = c == '-';
if (c == '-' || c == '+') { ++position_; c = peek(); }
if (static_cast<unsigned>(c-'0') > 9) return false;
value = 0;
do {
const char* p = buffer_+position_;
if constexpr (sizeof(U) >= 4 && std::endian::native == std::endian::little) {
while (buffer_+size_-p >= 8) {
std::uint64_t block;
std::memcpy(&block,p,8);
// All eight bytes must be ASCII decimal digits.
if (((block+0x4646464646464646ULL) | (block-0x3030303030303030ULL)) & 0x8080808080808080ULL) break;
block -= 0x3030303030303030ULL;
block = (block*10+(block>>8)) & 0x00ff00ff00ff00ffULL;
block = (block*100+(block>>16)) & 0x0000ffff0000ffffULL;
block = (block*10000+(block>>32)) & 0xffffffffULL;
value = value*100000000+static_cast<std::uint32_t>(block);
p += 8;
if constexpr (SingleBlock || sizeof(U) == 4) break;
}
}
unsigned digit;
// The sentinel stops this loop without an EOF/buffer check for
// every digit. Only a token crossing a refill needs the outer loop.
while ((digit = static_cast<unsigned>(static_cast<unsigned char>(*p)-'0')) <= 9) {
value = value*10+digit;
++p;
}
position_ = static_cast<std::size_t>(p-buffer_);
if (position_ != size_) break;
c = peek();
} while (static_cast<unsigned>(c-'0') <= 9);
return true;
}
template<class T> bool read_one(T& value) {
if constexpr (std::is_same_v<T, char>) {
const int c = skip_space();
if (c == EOF) return false;
value = static_cast<char>(c);
++position_;
return true;
} else if constexpr (detail::integer<T>) {
detail::unsigned_t<T> u;
bool negative;
if (!magnitude(u,negative)) return false;
if constexpr (detail::signed_integer<T>) {
// Never negate the minimum signed value, even for int128.
value = negative && u != 0 ? -static_cast<T>(u-1)-1 : static_cast<T>(u);
} else {
if (negative || (std::is_same_v<T,bool> && u > 1)) return false;
value = static_cast<T>(u);
}
return true;
} else if constexpr (detail::modular<T>) {
std::uint64_t u;
bool negative;
// Most modular coefficients have at most ten digits. After the
// first eight, parse the short tail without another wide probe.
if (!magnitude<std::uint64_t,true>(u,negative)) return false;
const auto mod = detail::modulus<T>();
if (u >= mod) u %= mod;
value.val = negative && u != 0 ? mod-u : u;
return true;
} else if constexpr (std::is_same_v<T,std::string>) {
int c = skip_space();
if (c == EOF) return false;
value.clear();
do {
value.push_back(static_cast<char>(c));
++position_;
c = peek();
} while (c != EOF && c > ' ');
return true;
} else if constexpr (std::is_floating_point_v<T>) {
std::string token;
if (!read_one(token)) return false;
const char* first = token.data();
const char* last = first+token.size();
if (first != last && *first == '+') ++first;
const auto result = std::from_chars(first,last,value);
return result.ec == std::errc{} && result.ptr == last;
} else if constexpr (detail::indexed_range<T>) {
const auto n = std::size(value);
if constexpr (requires { std::data(value); }) {
auto* data = std::data(value);
for (decltype(std::size(value)) i = 0; i < n; ++i)
if (!read_one(data[i])) return false;
} else {
for (decltype(std::size(value)) i = 0; i < n; ++i)
if (!read_one(value[i])) return false;
}
return true;
} else if constexpr (detail::tuple_like<T>) {
return std::apply([&](auto&... x) { return (read_one(x) && ...); },value);
} else {
static_assert(detail::unsupported<T>, "Unsupported fastio input type");
}
}
public:
explicit BasicInput(std::FILE* file = stdin) : file_(file) {}
BasicInput(const BasicInput&) = delete;
BasicInput& operator=(const BasicInput&) = delete;
template<class... T> bool read(T&... values) {
if (!good_) return false;
return good_ = (read_one(values) && ...);
}
template<class T> BasicInput& operator>>(T& value) { read(value); return *this; }
explicit operator bool() const { return good_; }
};
template<std::size_t BufferSize = 1 << 18>
class BasicOutput {
static_assert(BufferSize > 0);
std::FILE* file_;
char buffer_[BufferSize];
std::size_t size_ = 0;
bool good_ = true;
bool drain() {
if (!good_) return false;
std::size_t done = 0;
while (done < size_) {
const auto n = std::fwrite(buffer_+done,1,size_-done,file_);
if (n == 0) { size_ = 0; return good_ = false; }
done += n;
}
size_ = 0;
return true;
}
bool put(char c) {
if (!good_ || (size_ == BufferSize && !drain())) return false;
buffer_[size_++] = c;
return true;
}
bool bytes(const char* data, std::size_t count) {
if (!good_) return false;
while (count != 0) {
if (size_ == BufferSize && !drain()) return false;
const auto room = BufferSize-size_;
const auto n = count < room ? count : room;
std::memcpy(buffer_+size_,data,n);
size_ += n; data += n; count -= n;
}
return true;
}
template<class U> bool unsigned_number(U value, bool negative) {
if constexpr (sizeof(U) <= 4 && BufferSize >= 12) {
if (!good_ || (BufferSize-size_ < 12 && !drain())) return false;
char* p = buffer_+size_;
if (negative) *p++ = '-';
p = detail::decimal32(p,value);
size_ = static_cast<std::size_t>(p-buffer_);
return true;
}
// 128-bit unsigned integers need at most 39 decimal digits.
char digits[3*sizeof(U)+2];
char* const end = digits+sizeof(digits);
char* p = end;
while (value >= 100) {
const U q = value/100;
const unsigned r = static_cast<unsigned>(value-q*100);
p -= 2;
std::memcpy(p,detail::digit_pairs.data()+2*r,2);
value = q;
}
if (value < 10) *--p = char('0'+static_cast<unsigned>(value));
else {
p -= 2;
std::memcpy(p,detail::digit_pairs.data()+2*static_cast<unsigned>(value),2);
}
if (negative) *--p = '-';
if constexpr (BufferSize >= sizeof(digits)) {
if (!good_ || (BufferSize-size_ < sizeof(digits) && !drain())) return false;
auto count = static_cast<std::size_t>(end-p);
char* dest = buffer_+size_;
size_ += count;
// Constant-size copies avoid a library call per integer. Each copy
// stays inside the initialized digits and the reserved buffer area.
if constexpr (sizeof(digits) >= 8)
while (count >= 8) {
std::memcpy(dest,p,8);
dest += 8; p += 8; count -= 8;
}
if (count >= 4) { std::memcpy(dest,p,4); dest += 4; p += 4; count -= 4; }
if (count >= 2) { std::memcpy(dest,p,2); dest += 2; p += 2; count -= 2; }
if (count) *dest = *p;
return true;
}
return bytes(p,static_cast<std::size_t>(end-p));
}
template<class T> bool write_one(const T& value) {
if constexpr (std::is_same_v<T,char>) {
return put(value);
} else if constexpr (detail::integer<T>) {
using U = detail::unsigned_t<T>;
U u = static_cast<U>(value);
bool negative = false;
if constexpr (detail::signed_integer<T>) {
negative = value < 0;
if (negative) u = U(0)-u;
}
return unsigned_number(u,negative);
} else if constexpr (detail::modular<T>) {
return write_one(value.val);
} else if constexpr (std::is_convertible_v<const T&,std::string_view>) {
const std::string_view text(value);
return bytes(text.data(),text.size());
} else if constexpr (std::is_floating_point_v<T>) {
char text[128];
const auto result = std::to_chars(text,text+sizeof(text),value);
if (result.ec != std::errc{}) return false;
return bytes(text,static_cast<std::size_t>(result.ptr-text));
} else if constexpr (detail::indexed_range<const T>) {
const auto n = std::size(value);
using Element = std::remove_cvref_t<decltype(value[0])>;
if constexpr (detail::indexed_range<Element> || detail::tuple_like<Element>){
for (decltype(std::size(value)) i = 0; i < n; ++i) {
if (i != 0 && !put('\n')) return false;
if (!write_one(value[i])) return false;
}
return true;
}
if constexpr (BufferSize >= 32 && detail::output32<Element> && requires { std::data(value); }) {
const auto* data = std::data(value);
char* p = buffer_+size_;
for (decltype(std::size(value)) i = 0; i < n; ++i) {
if (buffer_+BufferSize-p < 12) {
size_ = static_cast<std::size_t>(p-buffer_);
if (!drain()) return false;
p = buffer_;
}
if (i != 0) *p++ = ' ';
p = detail::decimal_small(p,data[i]);
}
size_ = static_cast<std::size_t>(p-buffer_);
return true;
}
for (decltype(std::size(value)) i = 0; i < n; ++i) {
if (i != 0 && !put(' ')) return false;
if (!write_one(value[i])) return false;
}
return true;
} else if constexpr (detail::tuple_like<T>) {
[[maybe_unused]] bool first = true;
return std::apply([&](const auto&... x) {
return (((std::exchange(first,false) || put(' ')) && write_one(x)) && ...);
},value);
} else {
static_assert(detail::unsupported<T>, "Unsupported fastio output type");
}
}
public:
explicit BasicOutput(std::FILE* file = stdout) : file_(file) {}
BasicOutput(const BasicOutput&) = delete;
BasicOutput& operator=(const BasicOutput&) = delete;
~BasicOutput() { flush(); }
// write concatenates its arguments; print separates them with spaces and
// appends a newline. Ranges and tuples use spaces between their elements.
template<class... T> bool write(const T&... values) {
if (!good_) return false;
return good_ = (write_one(values) && ...);
}
template<class... T> bool print(const T&... values) {
if (!good_) return false;
[[maybe_unused]] bool first = true;
good_ = (((std::exchange(first,false) || put(' ')) && write_one(values)) && ...);
return good_ && put('\n');
}
template<class... T> bool printflush(const T&... values) {
if (!good_) return false;
[[maybe_unused]] bool first = true;
good_ = (((std::exchange(first,false) || put(' ')) && write_one(values)) && ...);
return good_ && put('\n') && flush();
}
bool flush() { return drain() && (good_ = std::fflush(file_) == 0); }
template<class T> BasicOutput& operator<<(const T& value) { write(value); return *this; }
explicit operator bool() const { return good_; }
};
#ifdef LOCAL_a1073741824
using FastInput = BasicInput<1>;
#else
using FastInput = BasicInput<>;
#endif
using FastOutput = BasicOutput<>;
// template<class... T> bool read(T&... values) { return in.read(values...); }
// template<class... T> bool write(const T&... values) { return out.write(values...); }
// template<class... T> bool print(const T&... values) { return out.print(values...); }
// template<class... T> bool printflush(const T&... values) { return out.print(values...) && out.flush(); }
// inline bool flush() { return out.flush(); }
} // namespace fastio
#endif /* FASTIO__HPP_ */
#ifndef FPS_OPERATION_HPP_
#define FPS_OPERATION_HPP_
#include <array>
#include <algorithm>
#include <functional>
#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif
#ifndef DEQUE__HPP_
#define DEQUE__HPP_
#include <deque>
#include <iostream>
#include <iterator>
using namespace std;
using ll = long long;
using uint = unsigned int;
template<typename T>
class delist{
public:
inline constexpr delist(){A = deque<T>();}
inline constexpr delist(uint sz, const T& val = T()){A = deque<T>(sz, val);}
template <class InputIterator> inline constexpr delist(InputIterator first, InputIterator last){A = deque<T>(first, last);}
inline constexpr void append(const T& val){A.push_back(val);}
inline constexpr void appendleft(const T& val){A.push_front(val);}
inline constexpr T pop(){
if (A.empty()){cerr << "IndexError: pop from empty list" << endl;}
auto temp = std::move(A.back());
A.pop_back();
return std::move(temp);
}
inline constexpr T popleft(){
if (A.empty()){cerr << "IndexError: pop from empty list" << endl;}
auto temp = std::move(A.front());
A.pop_front();
return std::move(temp);
}
inline constexpr T& front(){return A.front();}
inline constexpr const T& front() const {return A.front();}
inline constexpr T& back(){return A.back();}
inline constexpr const T& back() const {return A.back();}
inline constexpr _Deque_iterator<T, T&, T*> begin(){return A.begin();}
inline constexpr _Deque_iterator<T, T&, T*> begin() const {return A.begin();}
inline constexpr _Deque_iterator<T, T&, T*> end(){return A.end();}
inline constexpr _Deque_iterator<T, T&, T*> end() const {return A.end();}
inline constexpr ll size(){return A.size();}
inline constexpr bool empty(){return A.empty();}
private:
deque<T> A;
};
#endif /* DEQUE__HPP_ */
#ifndef FPS__HPP_
#define FPS__HPP_
#include <iostream>
#include <iterator>
#include <cassert>
#ifndef MODINT__HPP_
#define MODINT__HPP_
#include <iostream>
#ifndef MATH_FUNCTION_HPP_
#define MATH_FUNCTION_HPP_
#include <array>
#include <cmath>
#include <cassert>
using namespace std;
using ll = long long;
using ld = long double;
using ull = unsigned long long;
using ulll = __uint128_t;
#ifdef BOOST_VERSION
using bll = boost::multiprecision::cpp_int;
#endif
constexpr ll getdigit(ull x){
if (x < 1000000000){
if (x < 100000){
return x == 0 ? 0 : x < 10 ? 1 : x < 100 ? 2 : x < 1000 ? 3 : x < 10000 ? 4 : 5;
}
else{
return x < 1000000 ? 6 : x < 10000000 ? 7 : x < 100000000 ? 8 : 9;
}
}
else{
if (x < 100000000000000){
return x < 10000000000 ? 10 : x < 100000000000 ? 11 : x < 1000000000000 ? 12 : x < 10000000000000 ? 13 : 14;
}
else{
return x < 1000000000000000 ? 15 : x < 10000000000000000 ? 16 : x < 100000000000000000 ? 17 : x < 1000000000000000000 ? 18 : x < 10000000000000000000u ? 19 : 20;
}
}
}
/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる
constexpr ll modpow(ll a, ull b, const ll m){
ll t = a%m;
ll ans = (m == 1 ? 0 : 1);
while (b > 0){
if (b&1){
ans = (ans*t)%m;
}
b >>= 1;
t = (t*t)%m;
}
return ans;
}
/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる。mはコンパイル時に決定している必要がある
template<ll m> constexpr ll modpow(ll a, ull b){
ll t = a%m;
ll ans = (m == 1 ? 0 : 1);
while (b > 0){
if (b&1){
ans = (ans*t)%m;
}
b >>= 1;
t = (t*t)%m;
}
return ans;
}
/// @brief a^nを返す。bに関して対数時間で計算できる。
template<typename T> constexpr T powll(T a, T n){
T t = a;
T ans = 1;
while (n > 0){
if (n%2){
ans *= t;
}
n >>= 1;
t *= t;
}
return ans;
}
/// @brief a^nを返す。bに関して対数時間で計算できる
template<> constexpr ll powll(ll a, ll n){
ll t = a;
ll ans = 1;
while (n > 0){
if (n&1){
ans *= t;
}
n >>= 1;
t *= t;
}
return ans;
}
/// @brief floor(sqrt(N))を返す。1.5×10^19まで対応
constexpr ll isqrt(ull N){
assert(N <= 15000000000000000000ull);
ull ret = sqrt(N);
while (ret*ret > N){
ret--;
}
while ((ret+1)*(ret+1) <= N){
ret++;
}
return ret;
}
#ifdef BOOST_VERSION
/// @brief floor(sqrt(N))を返す。多倍長整数に対応
constexpr bll isqrt_large(bll N){
bll ret = sqrt(N);
while (ret*ret > N){
ret--;
}
while ((ret+1)*(ret+1) <= N){
ret++;
}
return ret;
}
#endif
/// @brief floor(log_a(L))を返す
constexpr ll ilog(ll a, ll L){
__int128_t t = 1;
ll ans = 0;
while (t <= L){
ans++;
t *= a;
}
return ans-1;
}
/// @brief 有理数のfloorを求める。 floor(y/x)
template<typename T> constexpr inline T floor2(T y, T x){
if ((x^y) > 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return y/x;
}
else if ((x^y) < 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return -((y+x-1)/x);
}
else{
return y/x;
}
}
/// @brief 有理数のceilを求める。 ceil(y/x)
template<typename T> constexpr inline T ceil2(T y, T x){
if ((x^y) > 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return (y+x-1)/x;
}
else if ((x^y) < 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return -(y/x);
}
else{
return y/x;
}
}
/// @brief 一次不定方程式ax+by=gcd(a,b)の解を1つ見つける
/// @param a `a>=0`である必要がある
/// @param b `b>=0`である必要がある
/// @return {x,y,gcd(a,b)}
template<typename T>
constexpr array<T,3> axby1(T a, T b){
T x = 1, y = 0;
T z = 0, w = 1;
T tmp = 0;
while (b){
T p = a/b, q = a%b;
tmp = x - y * p; x = y; y = tmp;
tmp = z - w * p; z = w; w = tmp;
a = b; b = q;
}
return {x, z, a};
}
/// @brief 1/a mod Mを求める
template<typename T, typename U>
constexpr T inverse_mod(T a, U M){
auto temp = axby1(a,(T)M);
assert(temp[2] == 1);
return (M+temp[0])%M;
}
/// @brief sqrt(a) mod Mを求める。ないなら-1が返される。
template<ll M>
inline constexpr ll cipolla(ll a){
a %= M;
if (M == 2) return a;
if (a == 0) return 0;
ll z = (M-1)/2;
if (modpow<M>(a, z) != 1){return -1;}
int b = 0;
while (modpow<M>((b*b+M-a)%M, z) == 1){
b++;
}
array<ll,2> x{1,0};
array<ll,2> y{b, 1};
ll w = (b*b+M-a)%M;
z++;
while (z){
if (z&1){
ll temp = x[0];
x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
if (x[0] >= M){x[0] -= M;}
x[1] = temp*y[1]%M+x[1]*y[0]%M;
if (x[1] >= M){x[1] -= M;}
}
ll temp = y[0];
y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
if (y[0] >= M){y[0] -= M;}
y[1] = 2*temp*y[1]%M;
z >>= 1;
}
return x[0];
}
inline constexpr ll cipolla(ll a, const ll M){
a %= M;
if (M == 2) return a;
if (a == 0) return 0;
ll z = (M-1)/2;
if (modpow(a, z, M) != 1){return -1;}
int b = 0;
while (modpow((b*b+M-a)%M, z, M) == 1){
b++;
}
array<ll,2> x{1,0};
array<ll,2> y{b, 1};
ll w = (b*b+M-a)%M;
z++;
while (z){
if (z&1){
ll temp = x[0];
x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
if (x[0] >= M){x[0] -= M;}
x[1] = temp*y[1]%M+x[1]*y[0]%M;
if (x[1] >= M){x[1] -= M;}
}
ll temp = y[0];
y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
if (y[0] >= M){y[0] -= M;}
y[1] = 2*temp*y[1]%M;
z >>= 1;
}
return x[0];
}
/// @brief x以下の最大の2冪を返す。0は0が返る。
inline constexpr ull lowerpow2(ull x){
if (x == 0){return 0;}
return 1ull<<(63-__builtin_clzll(x));
}
/// @brief x以上の最小の2冪を返す。0は0が返る。
inline constexpr ull upperpow2(ull x){
if (x == 0){return 0;}
if (x == 1){return 1;}
return 1ull<<(64-__builtin_clzll(x-1));
}
/// @brief xのpopcountを求める
inline constexpr int popcount(ull x){
return __builtin_popcountll(x);
}
/// @brief xのbit lengthを求める。
inline constexpr int bit_length(ull x){
if (x == 0){return 0;}
return 64-__builtin_clzll(x);
}
#endif /* MATH_FUNCTION_HPP_ */
using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;
using namespace std;
template <uint M>
struct constant_modint {
static constexpr uint mod_value = M;
uint val;
constant_modint() : val(0) {}
template <class T> constant_modint(T x) {
ll y = (ll)(x % (ll)M);
if (y < 0) y += M;
val = y;
}
constant_modint& operator+=(const constant_modint& rhs) {
val += rhs.val;
if (val >= M) val -= M;
return *this;
}
constant_modint& operator-=(const constant_modint& rhs) {
if (val < rhs.val) val += M;
val -= rhs.val;
return *this;
}
constant_modint& operator*=(const constant_modint& rhs) {
val = (ull)val * rhs.val % M;
return *this;
}
constant_modint& operator/=(const constant_modint& rhs) {
return *this *= rhs.inv();
}
constant_modint operator+() const { return *this; }
constant_modint operator-() const { return constant_modint(0) - *this; }
friend constant_modint operator+(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) += rhs; }
friend constant_modint operator-(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) -= rhs; }
friend constant_modint operator*(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) *= rhs; }
friend constant_modint operator/(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) /= rhs; }
friend bool operator==(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val != rhs.val; }
constant_modint pow(ull n) const {
constant_modint res = 1, a = *this;
while (n) {
if (n&1) res *= a;
a *= a;
n >>= 1;
}
return res;
}
constant_modint inv() const {
return inverse_mod<ll,ll>(val, M);
}
friend std::ostream& operator<<(std::ostream& os, const constant_modint& m) {
return os << m.val;
}
friend std::istream& operator>>(std::istream& is, constant_modint& m) {
ll x;
is >> x;
m = constant_modint(x);
return is;
}
};
template <int id = -1>
struct dynamic_modint {
uint val;
static uint& mod() {
static uint M = 998244353;
return M;
}
static void set_mod(uint m) {
mod() = m;
}
dynamic_modint() : val(0) {}
template <class T> dynamic_modint(T x) {
ll y = x % (ll)mod();
if (y < 0) y += mod();
val = (uint)y;
}
dynamic_modint& operator+=(const dynamic_modint& rhs) {
val += rhs.val;
if (val >= mod()) val -= mod();
return *this;
}
dynamic_modint& operator-=(const dynamic_modint& rhs) {
if (val < rhs.val) val += mod();
val -= rhs.val;
return *this;
}
dynamic_modint& operator*=(const dynamic_modint& rhs) {
val = (ull)val * rhs.val % mod();
return *this;
}
dynamic_modint& operator/=(const dynamic_modint& rhs) {
return *this *= rhs.inv();
}
dynamic_modint operator+() const { return *this; }
dynamic_modint operator-() const { return dynamic_modint(0) - *this; }
friend dynamic_modint operator+(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) += rhs; }
friend dynamic_modint operator-(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) -= rhs; }
friend dynamic_modint operator*(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) *= rhs; }
friend dynamic_modint operator/(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) /= rhs; }
friend bool operator==(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val != rhs.val; }
dynamic_modint pow(ull n) const {
dynamic_modint res = 1, a = *this;
while (n) {
if (n & 1) res *= a;
a *= a;
n >>= 1;
}
return res;
}
dynamic_modint inv() const {
ll a = val, b = mod(), u = 1, v = 0;
while (b) {
ll t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
return dynamic_modint(u);
}
friend std::ostream& operator<<(std::ostream& os, const dynamic_modint& m) {
return os << m.val;
}
friend std::istream& operator>>(std::istream& is, dynamic_modint& m) {
ll x;
is >> x;
m = dynamic_modint(x);
return is;
}
};
#endif /* MODINT__HPP_ */
#ifndef INV_TABLE_HPP_
#define INV_TABLE_HPP_
#include <cassert>
using uint = unsigned;
using ull = unsigned long long;
using ll = long long;
using namespace std;
/// @brief mod M上での階乗,逆元テーブルを保持する構造体
/// @tparam M
template<uint M>
class mod_table{
using mint = constant_modint<M>;
list<mint> invmodlist;
list<mint> factorialmodlist;
list<mint> factorialmodinvlist;
uint N_MAX;
public:
constexpr mod_table(const uint N_MAX__){
N_MAX = max(1u, N_MAX__);
invmodlist = list<mint>(N_MAX+1);
factorialmodlist = list<mint>(N_MAX+1);
factorialmodinvlist = list<mint>(N_MAX+1);
invmodlist[1] = 1;
for (uint i = 2; i <= N_MAX; i++){
invmodlist[i] = (M-M/i)*invmodlist[M%i];
}
factorialmodinvlist[0] = 1;
factorialmodlist[0] = 1;
for (ull i = 1; i <= N_MAX; i++){
factorialmodinvlist[i] = invmodlist[i]*factorialmodinvlist[i-1];
factorialmodlist[i] = factorialmodlist[i-1]*i;
}
}
inline constexpr mint inverse_mod(ll x) const {
assert(0 <= x and x <= N_MAX);
return invmodlist[x];
}
inline constexpr mint factorialmod(ll x) const {
assert(0 <= x and x <= N_MAX);
return factorialmodlist[x];
}
inline constexpr mint factorialmodinv(ll x) const {
assert(0 <= x and x <= N_MAX);
return factorialmodinvlist[x];
}
inline constexpr mint inverse_mod(mint x) const {return inverse_mod(x.val);}
inline constexpr mint factorialmod(mint x) const {return factorialmod(x.val);}
inline constexpr mint factorialmodinv(mint x) const {return factorialmodinv(x.val);}
};
#endif /* INV_TABLE_HPP_ */
using namespace std;
using ll = long long;
using uint = unsigned;
#ifndef LEN_DEFINED
#define LEN_DEFINED
template<typename T> inline constexpr ll len(const T& A){return ssize(A);}
#endif
template<uint T>
struct FormalPowerSeries{
using mint = constant_modint<T>;
list<mint> f;
ll sz;
constexpr FormalPowerSeries(){
sz = 0;
}
/// @brief 指定したサイズのfpsを作成する。
constexpr FormalPowerSeries(uint _init_sz) : f(_init_sz), sz(_init_sz) {}
constexpr FormalPowerSeries(uint _init_sz, mint _init_val) : f(_init_sz,_init_val), sz(_init_sz) {}
template<typename U> constexpr FormalPowerSeries(const list<U>& _init){
f.assign(_init.begin(), _init.end());
sz = f.size();
}
template<typename U> inline constexpr FormalPowerSeries(initializer_list<U> _init) : FormalPowerSeries(list<U>(_init)){}
inline constexpr uint size() const {return sz;}
inline constexpr bool empty() const {return size() == 0;}
inline constexpr mint* data() {return empty() ? nullptr : &f[0];}
inline constexpr const mint* data() const {return empty() ? nullptr : &f[0];}
constexpr void resize(uint _new_size){
if (sz == _new_size){return;}
f.resize(_new_size, 0);
sz = _new_size;
}
inline constexpr mint& operator[](uint deg) {return f[deg];}
inline constexpr const mint& operator[](uint deg) const {return f[deg];}
template<typename U> constexpr operator list<U>() const {
list<U> res(sz);
for (int i = 0; i < sz; i++){
res[i] = f[i].val;
}
return res;
}
constexpr void differential(){
if (sz == 0){return;}
for (int i = 0; i < sz-1; i++){
f[i] = mint(i+1)*f[i+1];
}
f.pop();
sz--;
}
constexpr void integral(const mod_table<T>& mtable){
if (sz == 0){return;}
sz++;
f.append(0);
for (int i = sz-1; i > 0; i--){
f[i] = mtable.inverse_mod(i)*f[i-1];
}
f[0] = 0;
}
constexpr mint evaluate(const mint& x){
mint res = 0;
for (int i = 0; i < sz; i++){
res *= x;
res += f[sz-1-i];
}
return res;
}
};
#endif /* FPS__HPP_ */
#ifndef PRIME_AND_DIVISORS_HPP_
#define PRIME_AND_DIVISORS_HPP_
#include <array>
#include <algorithm>
#include <cassert>
#include <numeric>
#include <iostream>
using namespace std;
using ll = long long;
using ulll = __uint128_t;
using pll = array<ll,2>;
using pull = array<ull,2>;
using pii = array<int,2>;
/// @brief 試し割り法で正の整数Nを素因数分解する
/// @return list<array<ll,2>>{{素因数1,個数}, {素因数2,個数}, {素因数3,個数}...}
template<typename T> constexpr list<array<T,2>> trial_division(T N){
assert(N > 0);
if (N == 1){
return list<array<T,2>> {{1,0}};
}
list<array<T,2>> R;//戻り値用リスト
const T M = isqrt(static_cast<ull>(N));
for (T i = 2; i <= M; i++){
if (N % i == 0){
T divide_count = 0;
while (N % i == 0){
divide_count++;
N /= i;
}
R.append({i,divide_count});
}
}
if (N != 1){
R.append({N,1});
}
return R;
}
/// @brief 素因数分解リストを受け取って約数関数の値を求める
/// @return 約数のK乗和
template<typename T> constexpr T divisor_function(const list<array<T,2>>& vv, T K, const T MOD = -1){
if (vv[0][0] == 1){
return 1;
}
T R = 1;
if (K == 0){
for (auto x : vv){
R *= x[1]+1;
if (MOD > 0){
R %= MOD;
}
}
}
else{
for (auto x : vv){
ll r = powll(x[0],K);
R *= (powll(r,x[1]+1) - 1)/(r - 1);
if (MOD > 0){
R %= MOD;
}
}
}
return R;
}
/// @brief 素因数分解の結果pを受け取って、約数リストを生成する。
template<typename T> constexpr list<T> enumerate_divisor(const list<array<T,2>>& p){
list<T> d{1};
if (p[0][0] == 1){
return d;
}
for (auto &v : p){
T t = d.size();
T temp = 1;
for (T w = 0; w < v[1]; w++){
temp *= v[0];
for (T i = 0; i < t; i++){
d.append(d[i]*temp);
}
}
}
sort(vall(d));
return d;
}
/// @brief 線形篩
/// @attention コンストラクタに整数Nを渡すことでN以下の整数を扱うことができる。
struct LinearSieve{
list<int> p_list;
list<int> lpf;
//Nを渡すことで1以上N以下の整数を扱うことができる
LinearSieve(int N): lpf(N+1,-1){
lpf[1] = 1;
int p_list_size = 0;
for (int i = 2; i <= N; i++){
if (lpf[i] < 0){
p_list.append(i);
p_list_size++;
lpf[i] = i;
}
for (int j = 0; j < p_list_size && p_list[j] <= lpf[i] && p_list[j]*i <= N; j++){
lpf[p_list[j]*i] = p_list[j];
}
}
}
list<pii> factorize(int x){
if (x == 1){return {{1,0}};}
list<pii> r;
do{
if (r.empty() || lpf[x] != r.back()[0]){
r.append({lpf[x], 1});
}
else{
r.back()[1]++;
}
x /= lpf[x];
}while(x > 1);
return r;
}
list<list<pii>> factorize_all(int N){
list<list<pii>> r(N+1);
r[1].append({1,0});
for (int i = 2; i <= N; i++){
r[i] = factorize(i);
}
return r;
}
/// @brief N以下の整数に対するメビウス関数の値を列挙する。
list<int> enumerate_mobius(int N){
list<int> ret(N+1);
ret[1] = 1;
for (int i = 2; i <= N; i++){
if (lpf[i] == i){ret[i] = -1; continue;}
int temp = i/lpf[i];
if (lpf[i] == lpf[temp]){
ret[i] = 0;
}
else{
ret[i] = ret[lpf[i]]*ret[temp];
}
}
return ret;
}
};
constexpr ull __MillerRabin_small[3] = {2,7,61};
constexpr ull __MillerRabin_large[7] = {2,325,9375,28178,450775,9780504,1795265022};
struct Montgomery64 {
ull n, ni, r2;
constexpr Montgomery64(ull n) : n(n), ni(n), r2(-ulll(n) % n){
for (int i = 0; i < 5; ++i) ni *= 2ull - n * ni;
ni *= -1;
}
constexpr ull reduce(ulll x) const {
ull m = (ull)x * ni;
ull res = (x + ulll(m) * n) >> 64;
return res >= n ? res - n : res;
}
constexpr ull mul(ull x, ull y) const {
return reduce(ulll(x) * y);
}
constexpr ull pow_not_reduced(ull a, ull b)const{
ull res = reduce(r2);
a = reduce(ulll(a) * r2);
while (b) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}
};
///@brief ミラーラビン素数判定法
template<typename T> constexpr bool MillerRabin(T N){
if (N <= 1) return false;
if (N == 2 || N == 3) return true;
if (N % 2 == 0) return false;
T d = N - 1;
T s = 0;
while (d % 2 == 0){
d >>= 1;
s++;
}
auto modpow_T = [&](T a,T b,T m){
T t = a%m;
T ans = (m == 1 ? 0 : 1);
while (b > 0){
if (b&1){
ans = (ans*t)%m;
}
b >>= 1;
t = (t*t)%m;
}
return ans;
};
auto check = [&](T a){
if (a >= N) a %= N;
if (a == 0) return true;
T x = modpow_T(a, d, N);
if (x == 1 || x == N-1) return true;
for (int r = 1; r < s; r++){
x = x*x%N;
if (x == N - 1) return true;
}
return false;
};
if (N < 4759123141){
for (auto a : __MillerRabin_small){
if (!check(static_cast<T>(a))) return false;
}
}
else if (N < 1000000000000000000){
for (auto a : __MillerRabin_large){
if (!check(static_cast<T>(a))) return false;
}
}
else{
for (auto a : __MillerRabin_large){
if (!check(static_cast<T>(a))) return false;
}
//50個の整数でテスト
for (T a = 638245987265732141; a < 638245987265732191; a++){
if (!check(a)) return false;
}
}
return true;
}
/// @brief ミラーラビン素数判定法のull範囲限定版
template<> constexpr bool MillerRabin(ull N){
if (N <= 1) return false;
if (N == 2 || N == 3) return true;
if (N % 2 == 0) return false;
ull d = N - 1;
int s = 0;
while (d % 2 == 0){
d >>= 1;
s++;
}
Montgomery64 mg(N);
auto check = [&](ull a){
if (a >= N) a %= N;
if (a == 0) return true;
ull x = mg.pow_not_reduced(a, d);
ull y = mg.reduce(x);
if (y == 1 || y == N-1) return true;
for (int r = 1; r < s; r++){
x = mg.mul(x, x);
if (mg.reduce(x) == N - 1) return true;
}
return false;
};
if (N < 4759123141ull){
for (auto a : __MillerRabin_small){
if (!check(a)) return false;
}
}
else{
for (auto a : __MillerRabin_large){
if (!check(a)) return false;
}
}
return true;
}
/// @brief ポラード・ロー法による素因数分解
template<typename T> constexpr list<array<T, 2>> Pollard_rho(T N){
assert(N > 0);
if (N == 1){return {{1,0}};}
list<T> res;
list<T> dq;
dq.append(N);
while (!dq.empty()){
T n = dq.back();
dq.pop();
bool next_loop = false;
while (!MillerRabin<T>(n)){
if ((n&1) == 0){
while ((n&1) == 0){res.append(2); n>>=1;}
if (n == 1){next_loop = true; break;}
continue;
}
T blocksize = sqrt(sqrt(sqrt(n)));
for (T c = 1; c <= n; c++){
T y = (0b1010110101101011110011ull^c^n)%n;
T r = 1;
T k = 0;
T k_mod_blocksize = 0;
bool next_loop_inner = false;
for (int _ = 0;;){
T q = 1;
T q_old = q;
T y_old = y;
T x = y;
while (k < r){
k++;
k_mod_blocksize++;
y = y*y%n+c; if (y >= n){ y -= n;}
q = q*(y >= x ? (y-x) : (x-y))%n;
if (k_mod_blocksize == blocksize){
k_mod_blocksize -= blocksize;
T g = gcd(q,n);
if (g > 1){
if (g != n){
dq.append(g);
dq.append(n/g);
next_loop = true; break;
}
else{
y = y_old;
q = q_old;
for (int i = 0; i < blocksize; i++){
y = y*y%n+c; if (y >= n){ y -= n;}
q = q*(y >= x ? (y-x) : (x-y))%n;
g = gcd(q,n);
if (g > 1){
if (g != n){
dq.append(g);
dq.append(n/g);
next_loop = true; break;
}
next_loop_inner = true; break;
}
}
if (next_loop || next_loop_inner) break;
}
}
q_old = q;
y_old = y;
}
}
if (next_loop || next_loop_inner) break;
q_old = gcd(q,n);
if (q_old > 1){
if (q_old != n){
dq.append(q_old);
dq.append(n/q_old);
next_loop = true; break;
}
next_loop_inner = true; break;
}
r <<= 1;
}
if (next_loop) break;
}
if (next_loop) break;
}
if (!next_loop) {
res.append(n);
}
}
sort(res.begin(), res.end());
list<array<T,2>> res2;
for (ull l = 0, r = 0, sz = res.size(); l < sz;){
while (r < sz && res[l] == res[r]){
r++;
}
res2.append({res[l], r-l});
l = r;
}
return res2;
}
/// @brief ポラード・ロー法による素因数分解のull範囲限定版
template<> constexpr list<pull> Pollard_rho(ull N){
assert(N > 0);
if (N == 1){return {{1,0}};}
list<ull> res;
list<ull> dq;
dq.append(N);
while (!dq.empty()){
ull n = dq.back();
dq.pop();
bool next_loop = false;
while (!MillerRabin<ull>(n)){
if ((n&1) == 0){
while ((n&1) == 0){res.append(2); n>>=1;}
if (n == 1){next_loop = true; break;}
continue;
}
Montgomery64 mg(n);
int blocksize = max<int>(1, pow(n, 0.125));
for (ull c = 1; c <= n; c++){
ull reduced_c = mg.reduce(c*(ulll)mg.r2);
ull y = (0b1010110101101011110011ull^c^n)%n; y = mg.reduce(y*(ulll)mg.r2);
ull r = 1;
ull k = 0;
int k_mod_blocksize = 0;
bool next_loop_inner = false;
for (int _ = 0; _ < 63; _++){
ull q = mg.reduce(mg.r2);
ull q_old = q;
ull y_old = y;
ull x = y;
while (k < r){
k++;
k_mod_blocksize++;
y = mg.reduce(y*(ulll)y)+reduced_c; if (y >= n){ y -= n;}
q = mg.reduce(q*(y >= x ? (ulll)(y-x) : (ulll)(x-y)));
if (k_mod_blocksize == blocksize){
k_mod_blocksize -= blocksize;
ull g = gcd(mg.reduce(q),n);
if (g > 1){
if (g != n){
dq.append(g);
dq.append(n/g);
next_loop = true; break;
}
else{
y = y_old;
q = q_old;
for (int i = 0; i < blocksize; i++){
y = mg.reduce(y*(ulll)y)+reduced_c; if (y >= n){ y -= n;}
q = mg.reduce(q*(y >= x ? (ulll)(y-x) : (ulll)(x-y)));
g = gcd(mg.reduce(q),n);
if (g > 1){
if (g != n){
dq.append(g);
dq.append(n/g);
next_loop = true; break;
}
next_loop_inner = true; break;
}
}
if (next_loop || next_loop_inner) break;
}
}
q_old = q;
y_old = y;
}
}
if (next_loop || next_loop_inner) break;
q_old = gcd(mg.reduce(q),n);
if (q_old > 1){
if (q_old != n){
dq.append(q_old);
dq.append(n/q_old);
next_loop = true; break;
}
next_loop_inner = true; break;
}
r <<= 1;
}
if (next_loop) break;
}
if (next_loop) break;
}
if (!next_loop) {
res.append(n);
}
}
sort(res.begin(), res.end());
list<pull> res2;
for (int l = 0, r = 0, sz = res.size(); l < sz;){
while (r < sz && res[l] == res[r]){
r++;
}
res2.append({res[l], static_cast<ull>(r-l)});
l = r;
}
return res2;
}
/// @brief 素因数分解を行う。`N`の大きさによって、試し割り法とポラード・ロー法を使い分けてくれる。
template<typename T> constexpr list<array<T, 2>> factorize(T N){
if (N < 640000){
return trial_division<T>(N);
}
if (N < 1000000000000000000){
return Pollard_rho<T>(N);
}
return Pollard_rho<T>(N);
}
template<> constexpr list<pull> factorize(ull N){
if (N < 640000ull){
return trial_division<ull>(N);
}
return Pollard_rho<ull>(N);
}
template<> constexpr list<pll> factorize(ll N){
auto res = factorize<ull>(N);
list<pll> res2(res.size());
for (int i = 0, sz = res.size(); i < sz; i++){
res2[i][0] = res[i][0];
res2[i][1] = res[i][1];
}
return res2;
}
template<> constexpr list<pii> factorize(int N){
auto res = factorize<ull>(N);
list<pii> res2(res.size());
for (int i = 0, sz = res.size(); i < sz; i++){
res2[i][0] = res[i][0];
res2[i][1] = res[i][1];
}
return res2;
}
struct Xorshift64{
ull state;
constexpr Xorshift64(ull seed = 881726454633252252ull) : state(seed) {}
constexpr ull next() {
state ^= state << 13;
state ^= state >> 7;
state ^= state << 17;
return state;
}
};
/// @brief mod Pにおける原始根を1つ探す。Pはull範囲内の素数でなければならない。
constexpr ull primitive_root(ull P){
assert(MillerRabin<ull>(P));
if (P == 2){
return 1;
}
Montgomery64 mg(P);
Xorshift64 rng(0x1234567E89ABCDEF);
auto factorized = factorize<ull>(P-1);
while (true) {
ull r = rng.next()%P;
if (r == 0){continue;}
bool ok = true;
for (auto& q : factorized){
if (mg.reduce(mg.pow_not_reduced(r, (P-1)/q[0])) == 1){
ok = false;
break;
}
}
if (ok) return r;
}
}
#endif /* PRIME_AND_DIVISORS_HPP_ */
using namespace std;
using ll = long long;
using uint = unsigned;
using ull = unsigned long long;
#define vall(A) A.begin(), A.end()
template<typename T> inline bool btest_for_fps(T K, int i){return K&(1ull<<i);}
/// @brief mod M上での形式的冪級数の計算を行う構造体
template<uint M>
struct fps_operator{
template<uint T> using fps = FormalPowerSeries<T>;
template<uint T> using mint = constant_modint<T>;
uint sum_e[30];
uint sum_ie[30];
uint log_max_length;
uint last_powroot;
// --- モンゴメリ乗算用の定数と関数 ---
static constexpr uint get_r() {
uint res = M;
for (int i = 0; i < 4; ++i) res *= 2 - M * res;
return res;
}
static constexpr uint R = get_r();
static constexpr uint R2 = -ull(M) % M;
static constexpr uint NEG_INV = 0 - R;
static inline constexpr uint reduce(ull x) {
uint res = (x + ull(uint(x) * NEG_INV) * M) >> 32;
return res >= M ? res - M : res;
}
static inline constexpr uint to_montgomery(uint x) {
return reduce(ull(x) * R2);
}
static inline constexpr uint from_montgomery(uint x) {
return reduce(x);
}
static inline constexpr uint montgomery_mul(uint x, uint y) {
return reduce(ull(x) * y);
}
static inline constexpr uint mod_add(uint x, uint y) {
return x + y >= M ? x + y - M : x + y;
}
static inline constexpr uint mod_sub(uint x, uint y) {
return x < y ? x + M - y : x - y;
}
// --- AVX-512 / AVX2 SIMD Intrinsics ---
#if defined(__AVX512F__) && defined(__AVX512DQ__)
static inline __m512i montgomery_mul_simd(__m512i a, __m512i b, __m512i mod, __m512i neg_inv) {
__m512i mul0 = _mm512_mul_epu32(a, b);
__m512i mul1 = _mm512_mul_epu32(_mm512_srli_epi64(a, 32), _mm512_srli_epi64(b, 32));
__m512i q0 = _mm512_mul_epu32(_mm512_mullo_epi32(mul0, neg_inv), mod);
__m512i q1 = _mm512_mul_epu32(_mm512_mullo_epi32(mul1, neg_inv), mod);
__m512i res0 = _mm512_srli_epi64(_mm512_add_epi64(mul0, q0), 32);
__m512i res1 = _mm512_srli_epi64(_mm512_add_epi64(mul1, q1), 32);
__m512i res = _mm512_or_si512(res0, _mm512_slli_epi64(res1, 32));
__mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
return _mm512_mask_sub_epi32(res, cmp, res, mod);
}
static inline __m512i mod_add_simd(__m512i a, __m512i b, __m512i mod) {
__m512i res = _mm512_add_epi32(a, b);
__mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
return _mm512_mask_sub_epi32(res, cmp, res, mod);
}
static inline __m512i mod_sub_simd(__m512i a, __m512i b, __m512i mod) {
__mmask16 cmp = _mm512_cmplt_epu32_mask(a, b);
__m512i res = _mm512_sub_epi32(a, b);
return _mm512_mask_add_epi32(res, cmp, res, mod);
}
#elif defined(__AVX2__)
// AVX2用 モンゴメリ乗算・剰余加減算
static inline __m256i montgomery_mul_simd(__m256i a, __m256i b, __m256i mod, __m256i neg_inv) {
__m256i mul0 = _mm256_mul_epu32(a, b);
__m256i mul1 = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
__m256i q0 = _mm256_mul_epu32(_mm256_mullo_epi32(mul0, neg_inv), mod);
__m256i q1 = _mm256_mul_epu32(_mm256_mullo_epi32(mul1, neg_inv), mod);
__m256i res0 = _mm256_srli_epi64(_mm256_add_epi64(mul0, q0), 32);
__m256i res1 = _mm256_srli_epi64(_mm256_add_epi64(mul1, q1), 32);
__m256i res = _mm256_blend_epi32(res0, _mm256_slli_epi64(res1, 32), 0xAA);
__m256i diff = _mm256_sub_epi32(res, mod);
// diff < 0 なら -1 (全ビット1)、そうでないなら 0 となるマスク
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
static inline __m256i mod_add_simd(__m256i a, __m256i b, __m256i mod) {
__m256i diff = _mm256_sub_epi32(_mm256_add_epi32(a, b), mod);
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
static inline __m256i mod_sub_simd(__m256i a, __m256i b, __m256i mod) {
__m256i diff = _mm256_sub_epi32(a, b);
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
#endif
constexpr fps_operator(){
list<ll> powroot{1};
list<ll> powrootinv;
while (powroot.back() >= 0){
powroot.append(cipolla<M>(powroot.back()));
if (powroot.back() == powroot[powroot.size()-2]){
powroot.back() = M-powroot.back();
}
}
powroot.pop();
log_max_length = powroot.size()-1;
last_powroot = powroot.back();
for (auto v : powroot){
powrootinv.append(inverse_mod<ll,ll>(v,M));
}
int cnt2 = powroot.size()-1;
uint now = reduce(R2);
for (int i = 0; i <= cnt2-2; i++){
sum_e[i] = montgomery_mul(to_montgomery(powroot[i+2]), now);
now = montgomery_mul(now, to_montgomery(powrootinv[i+2]));
}
uint inow = reduce(R2);
for (int i = 0; i <= cnt2-2; i++){
sum_ie[i] = montgomery_mul(to_montgomery(powrootinv[i+2]), inow);
inow = montgomery_mul(inow, to_montgomery(powroot[i+2]));
}
for (int i = cnt2-1; i < 30; i++){
sum_e[i] = 0;
sum_ie[i] = 0;
}
}
// Public transforms return canonical residues in the original bit-reversed
// order. IDFT remains unnormalized. Transposes share the same kernels.
private:
template<bool InverseButterfly, bool InverseRoots>
void transform_scalar(mint<M>* a, int n) const {
const uint* rate = InverseRoots ? sum_ie : sum_e;
const int h = __builtin_ctz(unsigned(n));
for (int stage = 0; stage < h; ++stage) {
const int p = InverseButterfly ? 1 << stage : n >> (stage + 1);
uint w = to_montgomery(1);
for (int offset = 0, s = 0; offset < n; offset += 2*p, ++s) {
for (int i = 0; i < p; ++i) {
const uint l = a[offset+i].val, r = a[offset+i+p].val;
if constexpr (InverseButterfly) {
a[offset+i].val = mod_add(l, r);
a[offset+i+p].val = montgomery_mul(mod_sub(l, r), w);
} else {
const uint t = montgomery_mul(r, w);
a[offset+i].val = mod_add(l, t);
a[offset+i+p].val = mod_sub(l, t);
}
}
w = montgomery_mul(w, rate[__builtin_ctz(~unsigned(s))]);
}
}
}
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
// Internal SIMD residues are in [0, 2M). M < 2^30 guarantees that
// additions fit uint32_t and Montgomery products reduce to [0, 2M).
__attribute__((target("avx2")))
static inline __m256i mul_lazy256(__m256i a, __m256i b) {
const __m256i mod = _mm256_set1_epi32(M);
const __m256i ni = _mm256_set1_epi32(NEG_INV);
const __m256i x = _mm256_mul_epu32(a, b);
const __m256i y = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
// Only the low 32 bits of each 64-bit lane are consumed below; a
// widening multiply avoids computing the unused odd 32-bit products.
const __m256i qx = _mm256_mul_epu32(_mm256_mul_epu32(x, ni), mod);
const __m256i qy = _mm256_mul_epu32(_mm256_mul_epu32(y, ni), mod);
const __m256i rx = _mm256_srli_epi64(_mm256_add_epi64(x, qx), 32);
const __m256i ry = _mm256_add_epi64(y, qy);
const __m256i r = _mm256_blend_epi32(rx, ry, 0xaa);
return r;
}
__attribute__((target("avx2")))
static inline __m256i mul256(__m256i a, __m256i b) {
const __m256i r = mul_lazy256(a,b);
return _mm256_min_epu32(r,_mm256_sub_epi32(r,_mm256_set1_epi32(M)));
}
__attribute__((target("avx2")))
static inline __m256i add_lazy256(__m256i a, __m256i b) {
const __m256i r = _mm256_add_epi32(a, b);
return _mm256_min_epu32(r, _mm256_sub_epi32(r, _mm256_set1_epi32(2*M)));
}
__attribute__((target("avx2")))
static inline __m256i sub_lazy256(__m256i a, __m256i b) {
const __m256i r = _mm256_sub_epi32(a, b);
return _mm256_min_epu32(r, _mm256_add_epi32(r, _mm256_set1_epi32(2*M)));
}
__attribute__((target("avx2")))
static inline __m256i shoup256(__m256i a, __m256i b, __m256i bq) {
// b is an ordinary residue; bq = floor(b*2^32/M). The result is
// congruent to a*b and lies in [0, 2M), without a final correction.
const __m256i q0 = _mm256_srli_epi64(_mm256_mul_epu32(a,bq),32);
const __m256i q1 = _mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(bq,32));
const __m256i q = _mm256_blend_epi32(q0,q1,0xaa);
const __m256i mod = _mm256_set1_epi32(M);
const __m256i r = _mm256_sub_epi32(_mm256_mullo_epi32(a,b),_mm256_mullo_epi32(q,mod));
return r;
}
template<bool InverseButterfly, bool InverseRoots, bool Canonical = true, bool ZeroUpperHalf = false>
__attribute__((target("avx2")))
void transform_avx2(mint<M>* a, int n) const {
const uint* rate = InverseRoots ? sum_ie : sum_e;
const int h = __builtin_ctz(unsigned(n));
uint rate4[29];
for (int i = 0; i+1 < h; ++i) rate4[i] = montgomery_mul(rate[0],rate[i+1]);
const __m256i imag = _mm256_set1_epi32(rate[0]);
const uint im = from_montgomery(rate[0]);
const __m256i imv = _mm256_set1_epi32(im);
const __m256i imq = _mm256_set1_epi32((ull(im)<<32)/M);
int first_stage = 0;
if constexpr (ZeroUpperHalf && !InverseButterfly) {
// The first radix-4 block has no twiddle and its upper two inputs
// are zero. Do not load padding or execute the redundant butterflies.
const int p = n/4;
for (int i = 0; i < p; i += 8) {
const __m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
const __m256i y = _mm256_loadu_si256((const __m256i*)(a+i+p));
const __m256i t = shoup256(y,imv,imq);
_mm256_storeu_si256((__m256i*)(a+i),add_lazy256(x,y));
_mm256_storeu_si256((__m256i*)(a+i+p),sub_lazy256(x,y));
_mm256_storeu_si256((__m256i*)(a+i+2*p),add_lazy256(x,t));
_mm256_storeu_si256((__m256i*)(a+i+3*p),sub_lazy256(x,t));
}
first_stage = 2;
}
for (int stage = first_stage; stage < h; ++stage) {
const int p4 = InverseButterfly ? 1 << stage : n >> (stage + 2);
if (stage + 1 < h && p4 >= 8) {
uint w = to_montgomery(1);
for (int offset = 0, s = 0; offset < n; offset += 4*p4, ++s) {
const uint w2 = montgomery_mul(w,w), w3 = montgomery_mul(w2,w);
const uint b1 = from_montgomery(w), b2 = from_montgomery(w2), b3 = from_montgomery(w3);
const __m256i v1 = _mm256_set1_epi32(b1), q1 = _mm256_set1_epi32((ull(b1)<<32)/M);
const __m256i v2 = _mm256_set1_epi32(b2), q2 = _mm256_set1_epi32((ull(b2)<<32)/M);
const __m256i v3 = _mm256_set1_epi32(b3), q3 = _mm256_set1_epi32((ull(b3)<<32)/M);
for (int i = 0; i < p4; i += 8) {
__m256i a0 = _mm256_loadu_si256((const __m256i*)(a+offset+i));
__m256i a1 = _mm256_loadu_si256((const __m256i*)(a+offset+i+p4));
__m256i a2 = _mm256_loadu_si256((const __m256i*)(a+offset+i+2*p4));
__m256i a3 = _mm256_loadu_si256((const __m256i*)(a+offset+i+3*p4));
__m256i x,y,z,t;
if constexpr (InverseButterfly) {
x = add_lazy256(a0,a1); y = sub_lazy256(a0,a1);
z = add_lazy256(a2,a3); t = shoup256(sub_lazy256(a2,a3),imv,imq);
a0 = add_lazy256(x,z); a1 = add_lazy256(y,t);
a2 = sub_lazy256(x,z); a3 = sub_lazy256(y,t);
if (s != 0) {
a1 = shoup256(a1,v1,q1); a2 = shoup256(a2,v2,q2); a3 = shoup256(a3,v3,q3);
}
} else {
if (s != 0) {
a1 = shoup256(a1,v1,q1); a2 = shoup256(a2,v2,q2); a3 = shoup256(a3,v3,q3);
}
x = add_lazy256(a0,a2); y = add_lazy256(a1,a3);
z = sub_lazy256(a0,a2); t = shoup256(sub_lazy256(a1,a3),imv,imq);
a0 = add_lazy256(x,y); a1 = sub_lazy256(x,y);
a2 = add_lazy256(z,t); a3 = sub_lazy256(z,t);
}
_mm256_storeu_si256((__m256i*)(a+offset+i),a0);
_mm256_storeu_si256((__m256i*)(a+offset+i+p4),a1);
_mm256_storeu_si256((__m256i*)(a+offset+i+2*p4),a2);
_mm256_storeu_si256((__m256i*)(a+offset+i+3*p4),a3);
}
w = montgomery_mul(w,rate4[__builtin_ctz(~unsigned(s))]);
}
++stage;
continue;
}
if (stage + 1 < h && n >= 32) {
const int p = p4, groups = 8/p;
alignas(32) uint lane1[8], lane2[8], lane3[8];
uint w = to_montgomery(1), last = w;
for (int s = 0; s < groups; ++s) {
const uint w2 = montgomery_mul(w,w), w3 = montgomery_mul(w2,w);
for (int i = 0; i < p; ++i) {
const int j = p == 1 ? (s/2 + (s%2)*4) : s*p+i;
lane1[j] = w; lane2[j] = w2; lane3[j] = w3;
}
last = w;
w = montgomery_mul(w,rate4[__builtin_ctz(~unsigned(s))]);
}
__m256i v1 = _mm256_load_si256((const __m256i*)lane1);
__m256i v2 = _mm256_load_si256((const __m256i*)lane2);
__m256i v3 = _mm256_load_si256((const __m256i*)lane3);
const int shift = __builtin_ctz(unsigned(groups));
uint advance1[30], advance2[30], advance3[30];
for (int i = 0; i+shift+1 < h; ++i) {
const uint d = montgomery_mul(last,rate4[i+shift]);
advance1[i] = d;
advance2[i] = montgomery_mul(d,d);
advance3[i] = montgomery_mul(advance2[i],d);
}
for (int offset = 0, s = 0; offset < n; offset += 32, ++s) {
__m256i x0 = _mm256_loadu_si256((const __m256i*)(a+offset));
__m256i x1 = _mm256_loadu_si256((const __m256i*)(a+offset+8));
__m256i x2 = _mm256_loadu_si256((const __m256i*)(a+offset+16));
__m256i x3 = _mm256_loadu_si256((const __m256i*)(a+offset+24));
__m256i a0,a1,a2,a3;
if (p == 1) {
const __m256i t0 = _mm256_unpacklo_epi32(x0,x1), t1 = _mm256_unpackhi_epi32(x0,x1);
const __m256i t2 = _mm256_unpacklo_epi32(x2,x3), t3 = _mm256_unpackhi_epi32(x2,x3);
a0 = _mm256_unpacklo_epi64(t0,t2); a1 = _mm256_unpackhi_epi64(t0,t2);
a2 = _mm256_unpacklo_epi64(t1,t3); a3 = _mm256_unpackhi_epi64(t1,t3);
} else if (p == 2) {
const __m256i t0 = _mm256_unpacklo_epi64(x0,x1), t1 = _mm256_unpackhi_epi64(x0,x1);
const __m256i t2 = _mm256_unpacklo_epi64(x2,x3), t3 = _mm256_unpackhi_epi64(x2,x3);
a0 = _mm256_permute2x128_si256(t0,t2,0x20); a1 = _mm256_permute2x128_si256(t1,t3,0x20);
a2 = _mm256_permute2x128_si256(t0,t2,0x31); a3 = _mm256_permute2x128_si256(t1,t3,0x31);
} else {
a0 = _mm256_permute2x128_si256(x0,x2,0x20); a1 = _mm256_permute2x128_si256(x0,x2,0x31);
a2 = _mm256_permute2x128_si256(x1,x3,0x20); a3 = _mm256_permute2x128_si256(x1,x3,0x31);
}
__m256i x,y,z,t;
if constexpr (InverseButterfly) {
x = add_lazy256(a0,a1); y = sub_lazy256(a0,a1);
z = add_lazy256(a2,a3); t = mul_lazy256(sub_lazy256(a2,a3),imag);
a0 = add_lazy256(x,z); a1 = mul_lazy256(add_lazy256(y,t),v1);
a2 = mul_lazy256(sub_lazy256(x,z),v2); a3 = mul_lazy256(sub_lazy256(y,t),v3);
} else {
a1 = mul_lazy256(a1,v1); a2 = mul_lazy256(a2,v2); a3 = mul_lazy256(a3,v3);
x = add_lazy256(a0,a2); y = add_lazy256(a1,a3);
z = sub_lazy256(a0,a2); t = mul_lazy256(sub_lazy256(a1,a3),imag);
a0 = add_lazy256(x,y); a1 = sub_lazy256(x,y);
a2 = add_lazy256(z,t); a3 = sub_lazy256(z,t);
}
if (p == 1) {
const __m256i t0 = _mm256_unpacklo_epi32(a0,a1), t1 = _mm256_unpackhi_epi32(a0,a1);
const __m256i t2 = _mm256_unpacklo_epi32(a2,a3), t3 = _mm256_unpackhi_epi32(a2,a3);
x0 = _mm256_unpacklo_epi64(t0,t2); x1 = _mm256_unpackhi_epi64(t0,t2);
x2 = _mm256_unpacklo_epi64(t1,t3); x3 = _mm256_unpackhi_epi64(t1,t3);
} else if (p == 2) {
const __m256i t0 = _mm256_permute2x128_si256(a0,a2,0x20), t1 = _mm256_permute2x128_si256(a1,a3,0x20);
const __m256i t2 = _mm256_permute2x128_si256(a0,a2,0x31), t3 = _mm256_permute2x128_si256(a1,a3,0x31);
x0 = _mm256_unpacklo_epi64(t0,t1); x1 = _mm256_unpackhi_epi64(t0,t1);
x2 = _mm256_unpacklo_epi64(t2,t3); x3 = _mm256_unpackhi_epi64(t2,t3);
} else {
x0 = _mm256_permute2x128_si256(a0,a1,0x20); x1 = _mm256_permute2x128_si256(a2,a3,0x20);
x2 = _mm256_permute2x128_si256(a0,a1,0x31); x3 = _mm256_permute2x128_si256(a2,a3,0x31);
}
_mm256_storeu_si256((__m256i*)(a+offset),x0);
_mm256_storeu_si256((__m256i*)(a+offset+8),x1);
_mm256_storeu_si256((__m256i*)(a+offset+16),x2);
_mm256_storeu_si256((__m256i*)(a+offset+24),x3);
const int k = __builtin_ctz(~unsigned(s));
v1 = mul256(v1,_mm256_set1_epi32(advance1[k]));
v2 = mul256(v2,_mm256_set1_epi32(advance2[k]));
v3 = mul256(v3,_mm256_set1_epi32(advance3[k]));
}
++stage;
continue;
}
const int p = InverseButterfly ? 1 << stage : n >> (stage + 1);
uint w = to_montgomery(1);
if (p >= 8) {
for (int offset = 0, s = 0; offset < n; offset += 2*p, ++s) {
const __m256i vw = _mm256_set1_epi32(w);
for (int i = 0; i < p; i += 8) {
__m256i l = _mm256_loadu_si256((const __m256i*)(a+offset+i));
__m256i r = _mm256_loadu_si256((const __m256i*)(a+offset+i+p));
__m256i x, y;
if constexpr (InverseButterfly) {
x = add_lazy256(l, r); y = mul_lazy256(sub_lazy256(l, r), vw);
} else {
r = mul_lazy256(r, vw); x = add_lazy256(l, r); y = sub_lazy256(l, r);
}
_mm256_storeu_si256((__m256i*)(a+offset+i), x);
_mm256_storeu_si256((__m256i*)(a+offset+i+p), y);
}
w = montgomery_mul(w, rate[__builtin_ctz(~unsigned(s))]);
}
} else {
// Pack 8 butterflies even when each block has only 1, 2 or 4 pairs.
alignas(32) uint lane[8]{};
uint t = to_montgomery(1), last = t;
const int groups = 8/p;
for (int s = 0; s < groups; ++s) {
for (int i = 0; i < p; ++i) lane[s*p+i] = t;
last = t;
t = montgomery_mul(t, rate[__builtin_ctz(~unsigned(s))]);
}
const __m256i base = _mm256_load_si256((const __m256i*)lane);
const __m256i unpack = p == 1 ? _mm256_setr_epi32(0,2,4,6,1,3,5,7)
: _mm256_setr_epi32(0,1,4,5,2,3,6,7);
const __m256i pack = p == 1 ? _mm256_setr_epi32(0,4,1,5,2,6,3,7) : unpack;
const int shift = __builtin_ctz(unsigned(groups));
uint advance[30];
for (int i = 0; i + shift < h; ++i) advance[i] = montgomery_mul(last, rate[i+shift]);
for (int offset = 0, s = 0; offset < n; offset += 16, ++s) {
__m256i x = _mm256_loadu_si256((const __m256i*)(a+offset));
__m256i y = _mm256_loadu_si256((const __m256i*)(a+offset+8));
if (p != 4) {
x = _mm256_permutevar8x32_epi32(x, unpack);
y = _mm256_permutevar8x32_epi32(y, unpack);
}
__m256i l = _mm256_permute2x128_si256(x,y,0x20);
__m256i r = _mm256_permute2x128_si256(x,y,0x31);
const __m256i vw = mul256(base, _mm256_set1_epi32(w));
if constexpr (InverseButterfly) {
x = add_lazy256(l,r); y = mul_lazy256(sub_lazy256(l,r),vw);
} else {
r = mul_lazy256(r,vw); x = add_lazy256(l,r); y = sub_lazy256(l,r);
}
l = _mm256_permute2x128_si256(x,y,0x20);
r = _mm256_permute2x128_si256(x,y,0x31);
if (p != 4) {
l = _mm256_permutevar8x32_epi32(l,pack);
r = _mm256_permutevar8x32_epi32(r,pack);
}
_mm256_storeu_si256((__m256i*)(a+offset),l);
_mm256_storeu_si256((__m256i*)(a+offset+8),r);
w = montgomery_mul(w,advance[__builtin_ctz(~unsigned(s))]);
}
}
}
const __m256i mod = _mm256_set1_epi32(M);
if constexpr (Canonical) for (int i = 0; i < n; i += 8) {
__m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
x = _mm256_min_epu32(x,_mm256_sub_epi32(x,mod));
_mm256_storeu_si256((__m256i*)(a+i),x);
}
}
template<bool Canonical = true>
__attribute__((target("avx2")))
static void multiply_avx2(mint<M>* dst, const mint<M>* a, const mint<M>* b, int n, uint scale) {
const __m256i vs = _mm256_set1_epi32(scale);
for (int i = 0; i < n; i += 8) {
__m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
const __m256i y = _mm256_loadu_si256((const __m256i*)(b+i));
x = mul_lazy256(mul_lazy256(x,y),vs);
if constexpr (Canonical)
x = _mm256_min_epu32(x,_mm256_sub_epi32(x,_mm256_set1_epi32(M)));
_mm256_storeu_si256((__m256i*)(dst+i),x);
}
}
#endif
static bool use_avx2() {
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
// This check also makes ordinary -O2/-O3 builds use SIMD on supported CPUs.
#ifdef FPS_DISABLE_SIMD
return false;
#else
return M < (1u << 30) && __builtin_cpu_supports("avx2");
#endif
#else
return false;
#endif
}
// Pointwise product with optional normalization, all operands in ordinary
// representation. dst may alias either input.
void pointwise_product(fps<M>& dst, const fps<M>& lhs, const fps<M>& rhs,
int n, mint<M> factor = 1) const {
if (n == 0) return;
auto* out = &dst[0];
const auto* a = &lhs[0];
const auto* b = &rhs[0];
int i = 0;
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
if (use_avx2()) {
i = n & ~7;
multiply_avx2(out,a,b,i,ull(R2)*factor.val%M);
}
#endif
for (; i < n; ++i) out[i] = a[i]*b[i]*factor;
}
template<bool InverseButterfly, bool InverseRoots>
void transform(fps<M>& F) const {
if (F.empty()) return;
F.f.reserve(upperpow2(F.sz));
F.resize(upperpow2(F.sz));
const int n = F.size();
assert(unsigned(n) <= (1u << log_max_length));
auto* a = &F[0];
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
if (n >= 16 && use_avx2()) {
transform_avx2<InverseButterfly, InverseRoots>(a,n);
return;
}
#endif
transform_scalar<InverseButterfly, InverseRoots>(a,n);
}
public:
void inplaceDFT(fps<M>& F) const { transform<false,false>(F); }
void inplaceIDFT(fps<M>& F) const { transform<true,true>(F); }
void inplaceDFT_T(fps<M>& F) const { transform<true,false>(F); }
void inplaceIDFT_T(fps<M>& F) const { transform<false,true>(F); }
private:
static fps<M> convolution_copy(const fps<M>& src, int capacity) {
assert(capacity > 0 && src.size() <= unsigned(capacity));
fps<M> result(capacity);
result.resize(src.size());
if (!src.empty()) std::copy_n(src.data(),src.size(),result.data());
return result;
}
static int convolution_capacity(const fps<M>& a, const fps<M>& b) {
if (a.empty() || b.empty() || std::min(a.size(),b.size()) <= 60) return 0;
if (std::max(a.size(),b.size())/std::min(a.size(),b.size()) >= 16) return 0;
return upperpow2(a.size()+b.size()-1);
}
fps<M> convolution_owned(fps<M> F1, fps<M> F2) const {
int n = F1.size(), m = F2.size();
if (n == 0 || m == 0) return fps<M>(0);
if (std::min(n,m) <= 60) {
if (n < m) { std::swap(n,m); std::swap(F1,F2); }
fps<M> ans(n+m-1);
auto* c = &ans[0];
const auto* a = &F1[0];
const auto* b = &F2[0];
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) c[i+j] += a[i]*b[j];
return ans;
}
if (std::max(n,m)/std::min(n,m) >= 16) {
if (n < m) { std::swap(n,m); std::swap(F1,F2); }
// Overlap-add: reuse the short operand's transform for each block
// instead of padding it to the length of the entire long operand.
const int length = upperpow2(2*m-1), step = length-m+1;
F2.f.reserve(length); F2.resize(length);
inplaceDFT(F2);
fps<M> block(length), ans(n+m-1);
const mint<M> scale = mint<M>(length).inv();
for (int start = 0; start < n; start += step) {
const int count = std::min(step,n-start);
std::copy_n(F1.data()+start,count,block.data());
std::fill(block.data()+count,block.data()+length,mint<M>(0));
inplaceDFT(block);
pointwise_product(block,block,F2,length,scale);
inplaceIDFT(block);
for (int i = 0; i < count+m-1; ++i) ans.data()[start+i] += block.data()[i];
}
return ans;
}
const int size = upperpow2(n+m-1);
assert(unsigned(size) <= (1u << log_max_length));
F1.f.reserve(size); F2.f.reserve(size);
F1.resize(size); F2.resize(size);
auto* a = &F1[0];
const auto* b = &F2[0];
// Two Montgomery reductions cancel R^2. Include 1/size here, so the
// inverse transform directly returns normalized, ordinary residues.
const uint scale = ull(R2)*modpow<M>(size,M-2)%M;
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
if (use_avx2()) {
if (n <= size/2) transform_avx2<false,false,false,true>(a,size);
else transform_avx2<false,false,false>(a,size);
if (m <= size/2) transform_avx2<false,false,false,true>(&F2[0],size);
else transform_avx2<false,false,false>(&F2[0],size);
multiply_avx2<false>(a,a,b,size,scale);
transform_avx2<true,true>(a,size);
} else
#endif
{
inplaceDFT(F1); inplaceDFT(F2);
for (int i = 0; i < size; ++i)
a[i].val = montgomery_mul(montgomery_mul(a[i].val,b[i].val),scale);
inplaceIDFT(F1);
}
F1.resize(n+m-1);
return F1;
}
public:
// Copy lvalues directly into the final NTT allocation. Rvalues retain their
// storage, so std::move still avoids copying when the caller is finished.
fps<M> convolution(const fps<M>& a, const fps<M>& b) const {
const int capacity = convolution_capacity(a,b);
if (capacity == 0) return convolution_owned(a,b);
return convolution_owned(convolution_copy(a,capacity),convolution_copy(b,capacity));
}
fps<M> convolution(fps<M>&& a, fps<M>&& b) const {
// Moving the same series twice must retain the original square semantics.
if (&a == &b) return convolution_owned(a,b);
return convolution_owned(std::move(a),std::move(b));
}
fps<M> convolution(fps<M>&& a, const fps<M>& b) const {
const int capacity = convolution_capacity(a,b);
auto copy = capacity ? convolution_copy(b,capacity) : b;
return convolution_owned(std::move(a),std::move(copy));
}
fps<M> convolution(const fps<M>& a, fps<M>&& b) const {
const int capacity = convolution_capacity(a,b);
auto copy = capacity ? convolution_copy(a,capacity) : a;
return convolution_owned(std::move(copy),std::move(b));
}
fps<M> add(const fps<M>& F1, const fps<M>& F2) const {
fps<M> ret(max(F1.size(),F2.size()));
for (uint i = 0; i < F1.size(); ++i) ret[i] = F1[i];
for (uint i = 0; i < F2.size(); ++i) ret[i] += F2[i];
return ret;
}
fps<M> subtract(const fps<M>& F1, const fps<M>& F2) const {
fps<M> ret(max(F1.size(),F2.size()));
for (uint i = 0; i < F1.size(); ++i) ret[i] = F1[i];
for (uint i = 0; i < F2.size(); ++i) ret[i] -= F2[i];
return ret;
}
/// @brief F*G == 1 mod x^n となるGを求める F[x^0] != 0 が必要
fps<M> inv(const fps<M>& F, const int n) const {
assert(n >= 0);
if (n == 0) return fps<M>(0);
assert(!F.empty() && F[0] != 0);
fps<M> res{F[0].inv()};
for (int d = 1; d < n; d <<= 1) {
fps<M> f(2*d), g(2*d);
for (int j = 0, limit = min<int>(F.size(),2*d); j < limit; ++j) f[j] = F[j];
for (int j = 0; j < d; ++j) g[j] = res[j];
inplaceDFT(f); inplaceDFT(g);
const mint<M> iz = mint<M>(2*d).inv();
pointwise_product(f,f,g,2*d,iz);
inplaceIDFT(f);
for (int j = 0; j < d; ++j) { f[j] = 0; f[j+d] = -f[j+d]; }
inplaceDFT(f);
pointwise_product(f,f,g,2*d,iz);
inplaceIDFT(f);
for (int j = 0; j < d; ++j) f[j] = res[j];
res = std::move(f);
}
res.resize(n);
return res;
}
/// @brief F÷G の商と余りを求める。{商, 余り}として返される。
array<fps<M>,2> div(fps<M> F, fps<M>G) const {
while (!F.f.empty() && F.f.back() == 0){
F.f.pop();
F.sz--;
}
while (!G.f.empty() && G.f.back() == 0){
G.f.pop();
G.sz--;
}
int n = F.sz, m = G.sz;
if (n < m){return {fps<M>{}, F};}
if (min(m, n-m+1) <= 100){
fps<M> res(n-m+1);
mint<M> t = G.f.back().inv();
for (int i = n-1; i >= m-1; i--){
auto u = t*F[i];
res[i-m+1] = u;
for (int j = 0; j < m; j++){
F[i-j] -= G[m-1-j]*u;
}
}
while (!F.f.empty() && F.f.back() == 0){
F.f.pop();
F.sz--;
}
while (!G.f.empty() && G.f.back() == 0){
G.f.pop();
G.sz--;
}
return {res, F};
}
auto F2 = F;
auto G2 = G;
reverse(vall(F.f));
reverse(vall(G.f));
while (len(F.f) > n-m+1){
F.f.pop();
F.sz--;
}
while (len(G.f) > n-m+1){
G.f.pop();
G.sz--;
}
auto res = convolution(F, inv(G, n-m+1));
while (len(res.f) > n-m+1){
res.f.pop();
res.sz--;
}
reverse(vall(res.f));
auto res2 = subtract(F2, convolution(res, G2));
while (!res2.f.empty() && res2.f.back() == 0){
res2.f.pop();
res2.sz--;
}
return {res, res2};
}
/// @brief F を G で割った余りを求める。
fps<M> remainder(fps<M> F, const fps<M>& G) const {
while (!F.f.empty() && F.f.back() == 0){
F.f.pop();
F.sz--;
}
assert(!G.f.empty() && G.f.back() != 0);
int n = F.sz, m = G.sz;
if (n < m) return F;
if (min(m, n-m+1) <= 100){
const bool monic = (G[m-1] == mint<M>(1));
const mint<M> inverse_leading_coefficient = monic ? mint<M>(1) : G[m-1].inv();
for (int i = n-1; i >= m-1; i--){
if (F[i] == 0) continue;
mint<M> q = monic ? F[i] : inverse_leading_coefficient*F[i];
for (int j = 0; j < m; j++){
F[i-j] -= G[m-1-j]*q;
}
}
F.resize(m-1);
while (!F.f.empty() && F.f.back() == 0){
F.f.pop();
F.sz--;
}
return F;
}
auto qr = div(std::move(F), G);
return std::move(qr[1]);
}
/// @brief exp(F) mod x^n を求める F[x^0] == 0 が必要
fps<M> exp(const fps<M>& F, const int n, const mod_table<M>& mtable) const {
assert(n >= 0);
if (n == 0) return fps<M>(0);
assert(F.empty() || F[0] == 0);
ull sz_f = F.sz;
fps<M> F0{1},f0{1},g0{1};
for(int d = 1; d < n; d<<=1){
auto G0 = g0;
inplaceDFT(G0);
const mint<M> iz = mint<M>(d).inv(), iz2 = mint<M>(2*d).inv();
fps<M> Delta(d);
pointwise_product(Delta,F0,G0,d,iz);
inplaceIDFT(Delta);
Delta[0] -= 1;
fps<M> delta(2*d);
for(int j = 0; j < d; j++) delta[d+j] = Delta[j];
fps<M> epsilon(2*d);
fps<M> DF0(d-1);
for (ull j = 0, ulld = d; j < ulld-1; j++) DF0[j] = f0[j+1]*mint<M>(j+1);
DF0.f.append(0);
DF0.sz++;
inplaceDFT(DF0);
pointwise_product(DF0,DF0,G0,d,iz);
inplaceIDFT(DF0);
for (ull j = 0, limit = d-1; j < limit; ++j) {
const auto df = j+1 < sz_f ? F[j+1]*mint<M>(j+1) : mint<M>(0);
epsilon[j] = df;
epsilon[j+d] = DF0[j]-df;
}
epsilon[d-1] = DF0[d-1];
Delta = std::move(delta);
inplaceDFT(Delta);
fps<M> DH0(d-1);
for (ull j = 0, limit = d-1; j < limit; j++) DH0[j] = j+1 < sz_f ? F[j+1]*mint<M>(j+1) : mint<M>(0);
DH0.resize(2*d);
inplaceDFT(DH0);
pointwise_product(Delta,Delta,DH0,2*d,iz2);
inplaceIDFT(Delta);
for(int j = 0; j < d; j++){
epsilon[j+d] -= Delta[j+d];
}
for (ull i = 2*d-1; i >= 1; i--){
epsilon[i] = epsilon[i-1]*mtable.inverse_mod(i)-(i < sz_f ? F[i] : mint<M>(0));
}
epsilon[0] = 0;
auto Epsilon = std::move(epsilon);
inplaceDFT(Epsilon);
for (int j = 0; j < d; j++){DH0[j] = f0[j]; DH0[j+d].val = 0;}
inplaceDFT(DH0);
pointwise_product(Epsilon,Epsilon,DH0,2*d,iz2);
inplaceIDFT(Epsilon);
f0.resize(2*d);
for (int j = 0; j < d; j++){
f0[j+d] -= Epsilon[j+d];
}
if(2*d >= n) break;
G0.resize(2*d);
for (int j = 0; j < d; j++) G0[j] = g0[j];
inplaceDFT(G0);
F0 = f0;
inplaceDFT(F0);
fps<M> T(2*d);
pointwise_product(T,F0,G0,2*d,iz2);
inplaceIDFT(T);
for (int j = 0; j < d; j++){
T[j].val = 0;
T[j+d] = -T[j+d];
}
inplaceDFT(T);
pointwise_product(T,T,G0,2*d,iz2);
inplaceIDFT(T);
for (int j = 0; j < d; j++) T[j] = g0[j];
g0 = std::move(T);
}
f0.resize(n);
return f0;
}
/// @brief log(F) mod x^n を求める F[x^0] == 1 が必要
fps<M> log(const fps<M>& F, const int n, const mod_table<M>& mtable) const {
assert(n >= 0);
if (n == 0) return fps<M>(0);
assert(!F.empty() && F[0] == 1);
const int degree = min<int>(F.size(),n)-1;
fps<M> derivative(degree);
for (int i = 0; i < degree; ++i) derivative[i] = F[i+1]*mint<M>(i+1);
auto res = convolution(std::move(derivative),inv(F,n));
res.resize(n);
for (int i = n-1; i >= 1; --i) res[i] = res[i-1]*mtable.inverse_mod(i);
res[0] = 0;
return res;
}
/// @brief F^k mod x^n を求める
fps<M> pow(const fps<M>& F, const ull k, const int n, const mod_table<M>& mtable) const {
assert(n >= 0);
if (n == 0) return fps<M>(0);
fps<M> ret(n);
if (k == 0) { ret[0] = 1; return ret; }
int lowest = 0;
while (lowest < min<int>(F.size(),n) && F[lowest] == 0) ++lowest;
if (lowest == min<int>(F.size(),n)) return ret;
// Check before multiplying, including k = ULLONG_MAX.
if (lowest > 0 && k > ull((n-1)/lowest)) return ret;
const int offset = ull(lowest)*k;
const int precision = n-offset;
const auto leading = F[lowest], inverse_leading = leading.inv();
fps<M> g(precision);
for (int i = 0, limit = min<int>(precision,F.size()-lowest); i < limit; ++i)
g[i] = F[lowest+i]*inverse_leading;
auto logarithm = log(g,precision,mtable);
const mint<M> exponent = k;
for (auto& x : logarithm.f) x *= exponent;
auto result = exp(logarithm,precision,mtable);
const auto scale = leading.pow(k);
for (int i = 0; i < precision; ++i) ret[offset+i] = result[i]*scale;
return ret;
}
// f(g(x)) mod x^len(f) を求める
fps<M> composition(fps<M> f, fps<M> g, const mod_table<M>& mtable) const {
f.resize(max(f.size(),g.size()));
int N = f.size();
if (N == 0) return fps<M>(0);
int n = 1;
while (n < (int)f.size()) n *= 2;
f.resize(n), g.resize(n);
fps<M> W(2*n);
{
// bit reverse order
list<int> btr(2*n);
int log = 31-__builtin_clz(2*n);
for (int i = 0; i < 2*n; i++){btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (log - 1));}
int t = log_max_length;
uint r = last_powroot;
ull dw = modpow<M>(inverse_mod<ll,ll>(r, M), (1<<t)/(4*n));
uint w = 1;
for (auto i: btr) { W[i].val = w, w = w*dw%M; }
}
auto rec = [&](auto &rec, int n, int k, fps<M> &Q) -> fps<M> {
if (n == 1) {
reverse(f.f.begin(), f.f.end());
fps<M> p(2*k);
for (int i = 0; i < k; i++) p[2 * i] = f[i];
return p;
}
Q.resize(4*n*k);
Q[2*n*k].val = 1;
inplaceDFT(Q);
fps<M> nxt_Q(2*n*k);
for (int i = 0; i < 2*n*k; i++) nxt_Q[i] = Q[2*i]*Q[2*i+1];
inplaceIDFT(nxt_Q);
ull iz = inverse_mod<ll,ll>(nxt_Q.size(), M);
for (uint i = 0; i < nxt_Q.size(); i++) nxt_Q[i].val = nxt_Q[i].val*iz%M;
for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) nxt_Q[n*j+i].val = 0;
nxt_Q[0].val = 0;
fps<M> p = rec(rec, n/2, 2*k, nxt_Q);
for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) p[n * j + i].val = 0;
iz = mtable.inverse_mod(p.size()).val;
for (uint i = 0; i < p.size(); i++) p[i].val = p[i].val*iz%M;
inplaceIDFT_T(p);
p.resize(4*n*k);
for (int i = 2*n*k-1; i >= 0; i--){
const auto t = mtable.inverse_mod(2)*W[i]*p[i];
p[2*i+1] = -t*Q[2*i];
p[2*i] = t*Q[2*i+1];
}
inplaceDFT_T(p);
p.resize(2*n*k);
return p;
};
fps<M> Q(2*n);
for(int i = 0; i < n; i++) Q[i].val = g[i].val == 0 ? 0 : M-g[i].val;
fps<M> p = rec(rec, n, 1, Q);
p.resize(n);
reverse(p.f.begin(), p.f.end());
p.resize(N);
return p;
}
fps<M> power_projection(fps<M> f, fps<M> wt, int m, const mod_table<M>& mtable) const {
assert(len(f) == len(wt));
if (f.empty()) { return fps<M>(m + 1, 0); }
if (f[0] != 0) {
auto c = f[0];
f[0] = 0;
fps<M> A = power_projection(f, wt, m, mtable);
for (int p = 0; p < m + 1; p++) A[p] *= mtable.factorialmodinv(p);
fps<M> B(m + 1);
mint<M> pow = 1;
for (int q = 0; q < m + 1; q++) B[q] = pow * mtable.factorialmodinv(q), pow *= c;
A = convolution(A, B);
A.resize(m + 1);
for (int i = 0; i < m + 1; i++) A[i] *= mtable.factorialmod(i);
return A;
}
int n = 1;
while (n < len(f)) n *= 2;
for (auto& x: f.f) x = -x;
f.resize(n), wt.resize(n);
reverse(vall(wt.f));
fps<M>&P = wt, &Q = f;
P.resize(4 * n), Q.resize(4 * n);
fps<M> W(n);
{
// bit reverse order
list<int> btr(n);
int log = 31-__builtin_clz(n);
for (int i = 1; i < n; i++) { btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (log - 1)); }
int t = log_max_length;
mint<M> r = last_powroot;
mint<M> dw = r.inv().pow((1 << t) / (2 * n));
mint<M> w = 1;
for (auto i: btr) { W[i] = w, w *= dw; }
}
int k = 1;
while (n > 1) {
auto normalized_IDFT = [&](fps<M>& A) {
inplaceIDFT(A);
mint<M> iz = mint<M>(A.size()).inv();
for (auto& x: A.f) x *= iz;
};
auto doubling_y = [&](fps<M>& A, int l, int r) {
mint<M> z = W[k / 2].inv();
fps<M> t(k);
for (int i = l; i < r; i++) {
for (int j = 0; j < k; j++) t[j] = A[2 * n * j + i];
normalized_IDFT(t);
mint<M> x = 1;
for (int j = 0; j < k; j++) t[j] *= x, x *= z;
inplaceDFT(t);
for (int j = 0; j < k; j++) A[2 * n * (k + j) + i] = t[j];
}
};
auto FFT_x = [&](fps<M>& A, int l, int r) {
fps<M> t(2 * n);
for (int j = l; j < r; j++) {
copy(A.f.begin() + 2 * n * j, A.f.begin() + 2 * n * (j + 1), t.f.begin());
inplaceDFT(t);
copy(vall(t.f), A.f.begin() + 2 * n * j);
}
};
if (n <= k) {
doubling_y(P, 0, n), doubling_y(Q, 1, n);
FFT_x(P, 0, 2 * k), FFT_x(Q, 0, 2 * k);
} else {
FFT_x(P, 0, k), FFT_x(Q, 0, k);
doubling_y(P, 0, 2 * n), doubling_y(Q, 0, 2 * n);
}
for (int i = 0; i < 2 * n * k; i++) Q[i] += 1;
for (int i = 2 * n * k; i < 4 * n * k; i++) Q[i] -= 1;
fps<M> F(2 * n), G(2 * n), p(n), q(n);
for (int j = 0; j < 2 * k; j++) {
copy(P.f.begin() + 2 * n * j, P.f.begin() + 2 * n * (j + 1), F.f.begin());
copy(Q.f.begin() + 2 * n * j, Q.f.begin() + 2 * n * (j + 1), G.f.begin());
for (int i = 0; i < n; i++) {
p[i] = W[i] * (F[2 * i] * G[2 * i + 1] - F[2 * i + 1] * G[2 * i]);
q[i] = G[2 * i] * G[2 * i + 1];
}
normalized_IDFT(p), normalized_IDFT(q);
fill(p.f.begin() + n / 2, p.f.end(), 0);
fill(q.f.begin() + n / 2, q.f.end(), 0);
copy(vall(p.f), P.f.begin() + n * j);
copy(vall(q.f), Q.f.begin() + n * j);
}
fill(P.f.begin() + 2 * n * k, P.f.end(), 0);
fill(Q.f.begin() + 2 * n * k, Q.f.end(), 0);
for (int j = 0; j < 4 * k; j++) Q[n * j] = 0;
n /= 2, k *= 2;
}
for (int i = 0; i < k; i++) P[i] = P[2 * i];
P.resize(k);
mint<M> ik = mint<M>(k).inv();
for (auto& x: P.f) x *= ik;
inplaceIDFT(P);
for (auto& x: P.f) x *= ik;
reverse(vall(P.f));
P.resize(m + 1);
return P;
}
fps<M> compositional_inverse(fps<M> f, mod_table<M>& mtable) const {
/*
仮定：[x^1]f = c は 0 ではない.
f を改めて cf と書く. f の逆関数を g とすれば
g(cf(x)/c)=x なので g(x/c) が求めるもの．よって [x^1]f=1 に帰着.
f, g は n 次とする. Lagrange inversion formula
n[x^n]f(x)^i = [x^{n-i}] i(g(x)/x)^{-n}
左辺を i=1,...,n で求めると (g(x)/x)^{-n} が n-1 次まで求まる.
[x^1]g=1 に帰着していたので -1/n 乗根がとれて, g(x)/x が n-1 次まで求まる.
*/
const int n = len(f) - 1;
if (n == -1) return {};
assert(f[0] == 0);
if (n == 0) return f;
assert(f[1] != 0);
mint<M> c = f[1];
mint<M> ic = c.inv();
for (auto& x: f.f) x *= ic;
fps<M> wt(n + 1);
wt[n] = 1;
fps<M> A = power_projection(f, wt, n, mtable);
fps<M> g(n);
for (int i = 1; i < n + 1; i++) g[n - i] = mint<M>(n) * A[i] * mint<M>(i).inv();
g = pow(g, M - mint<M>(n).inv().val, len(g), mtable);
g.f.insert(g.f.begin(), 0);
g.sz++;
mint<M> pow = 1;
for (int i = 0; i < len(g); i++) g[i] *= pow, pow *= ic;
return g;
}
// f(x+a)を求める
fps<M> polynominal_taylor_shift(const fps<M>& F, ll a, const mod_table<M>& mtable){
const int n = F.size();
const mint<M> shift = a;
if (n == 0 || shift == 0) return F;
fps<M> g(n), h(n);
mint<M> power = 1;
for (int i = 0; i < n; ++i) {
g[n-1-i] = F[i]*mtable.factorialmod(i);
h[i] = power*mtable.factorialmodinv(i);
power *= shift;
}
g = convolution(std::move(g),std::move(h));
fps<M> result(n);
for (int i = 0; i < n; ++i) result[i] = g[n-1-i]*mtable.factorialmodinv(i);
return result;
}
// fの微分を求める。
fps<M> differential(const fps<M>& F) const {
fps<M> G = F;
G.differential();
return G;
}
// fの積分を求める。積分定数は0
fps<M> integral(const fps<M>& F, const mod_table<M>& mtable) const {
fps<M> G = F;
G.integral(mtable);
return G;
}
};
/// @brief `[x^N].val(P(x)/Q(x))` をmod Mで求める
template<uint M>
constant_modint<M> Bostan_Mori(const ll N, FormalPowerSeries<M> P, FormalPowerSeries<M> Q, const fps_operator<M>& op){
assert(N >= 0);
if (N == 0){
return P[0].val*inverse_mod<ll,ll>(Q[0].val, M)%M;
}
FormalPowerSeries<M> Q_minus;
const int maxloop = (N == 1 ? 1 : 65-__builtin_clzll(N-1));
for (int _i_ = 0; _i_ < maxloop; _i_++){
Q_minus.resize(Q.size());
for (size_t i = 0; i < Q_minus.size(); i += 2){
Q_minus[i] = Q[i];
}
for (size_t i = 1; i < Q_minus.size(); i += 2){
Q_minus[i].val = (M-1)*Q[i].val%M;
}
auto A = op.convolution(P,Q_minus);
auto B = op.convolution(Q,Q_minus);
Q.resize((B.size()+1)/2);
for (size_t i = 0; i < B.size(); i += 2){
Q[i/2] = B[i];
}
P.resize((A.size()+!btest_for_fps(N,_i_))/2);
for (size_t i = btest_for_fps(N,_i_); i < A.size(); i += 2){
P[i/2] = A[i];
}
}
return P[0].val*inverse_mod<ll,ll>(Q[0].val, M)%M;
}
template<uint M>
list<FormalPowerSeries<M>> multipoint_product_tree(const list<constant_modint<M>>& samples, const fps_operator<M>& op){
const int n = len(samples);
if (n == 0) return {};
list<FormalPowerSeries<M>> tree(2*n-1);
int next_id = 0;
auto build = [&](auto&& self, int l, int r) -> int {
const int node = next_id++;
if (r-l == 1){
tree[node] = FormalPowerSeries<M>{-samples[l], constant_modint<M>(1)};
return node;
}
const int mid = (l+r)/2;
const int left = self(self, l, mid);
const int right = self(self, mid, r);
if (r-l == 2){
const auto x0 = samples[l], x1 = samples[l+1];
tree[node] = FormalPowerSeries<M>{x0*x1, -(x0+x1), constant_modint<M>(1)};
}
else if (r-l == 3){
const auto x0 = samples[l], x1 = samples[l+1], x2 = samples[l+2];
tree[node] = FormalPowerSeries<M>{
-x0*x1*x2,
x0*x1+x0*x2+x1*x2,
-(x0+x1+x2),
constant_modint<M>(1)
};
}
else{
tree[node] = op.convolution(tree[left], tree[right]);
}
return node;
};
build(build, 0, n);
return tree;
}
template<uint M>
list<constant_modint<M>> multipoint_evaluation_with_product_tree(const FormalPowerSeries<M>& F, const list<constant_modint<M>>& samples, const list<FormalPowerSeries<M>>& tree, const fps_operator<M>& op){
const int n = len(samples);
if (n == 0) return {};
list<constant_modint<M>> ans(n);
list<FormalPowerSeries<M>> remainders;
remainders.reserve(n);
list<array<int, 4>> stack;
stack.reserve(64);
remainders.append(F);
stack.append({0,n,0,0});
while (!stack.empty()){
auto [l,r,node,remainder_index] = stack.pop();
if (l >= r){continue;}
if (r-l == 1){
constant_modint<M> value = 0;
for (int i = int(remainders[remainder_index].sz)-1; i >= 0; i--){
value = value*samples[l]+remainders[remainder_index][i];
}
ans[l] = value;
continue;
}
if (remainders[remainder_index].size() >= tree[node].size()){
remainders.append(op.remainder(remainders[remainder_index], tree[node]));
remainder_index = int(remainders.size())-1;
}
const int mid = (l+r)/2;
const int left = node+1;
const int right = node+2*(mid-l);
stack.append({mid,r,right,remainder_index});
stack.append({l,mid,left,remainder_index});
}
return ans;
}
template <uint M>
list<constant_modint<M>> multipoint_evaluation(const FormalPowerSeries<M>& F, const list<constant_modint<M>>& samples, const fps_operator<M>& op){
auto tree = multipoint_product_tree(samples, op);
return multipoint_evaluation_with_product_tree(F, samples, tree, op);
}
template<uint M>
FormalPowerSeries<M> product_of_polynomial_sequence(const list<FormalPowerSeries<M>>& A, const fps_operator<M>& op){
delist<FormalPowerSeries<M>> dq;
dq.append(FormalPowerSeries<998244353>{1});
for (auto& f : A){
dq.append(std::move(f));
}
while (dq.size() > 1){
dq.append(op.convolution(std::move(*dq.begin()), std::move(*++dq.begin())));
dq.popleft();
dq.popleft();
}
return dq.front();
}
template<uint M>
FormalPowerSeries<M> polynomial_interpolation(const list<constant_modint<M>>& X, const list<constant_modint<M>>& Y, const fps_operator<M>& op){
assert(len(X) == len(Y));
const int N = len(X);
if (N == 0) return {};
const auto tree = multipoint_product_tree(X, op);
auto product_derivative = op.differential(tree[0]);
list<constant_modint<M>> Z =
multipoint_evaluation_with_product_tree(product_derivative, X, tree, op);
list<constant_modint<M>> prefix_product(N+1);
prefix_product[0] = 1;
for (int i = 0; i < N; i++){
prefix_product[i+1] = prefix_product[i]*Z[i];
}
constant_modint<M> inverse_suffix_product = prefix_product[N].inv();
for (int i = N-1; i >= 0; i--){
const auto denominator = Z[i];
Z[i] = Y[i]*prefix_product[i]*inverse_suffix_product;
inverse_suffix_product *= denominator;
}
auto interpolate = [&](auto&& self, int l, int r, int node) -> FormalPowerSeries<M> {
if (r-l == 1){
return FormalPowerSeries<M>{Z[l]};
}
if (r-l == 2){
const auto x0 = X[l], x1 = X[l+1];
const auto z0 = Z[l], z1 = Z[l+1];
return FormalPowerSeries<M>{-(z0*x1+z1*x0), z0+z1};
}
if (r-l == 3){
const auto x0 = X[l], x1 = X[l+1], x2 = X[l+2];
const auto z0 = Z[l], z1 = Z[l+1], z2 = Z[l+2];
return FormalPowerSeries<M>{
z0*x1*x2+z1*x0*x2+z2*x0*x1,
-(z0*(x1+x2)+z1*(x0+x2)+z2*(x0+x1)),
z0+z1+z2
};
}
const int mid = (l+r)/2;
const int left = node+1;
const int right = node+2*(mid-l);
auto left_interpolation = self(self, l, mid, left);
auto right_interpolation = self(self, mid, r, right);
return op.add(
op.convolution(std::move(left_interpolation), tree[right]),
op.convolution(tree[left], std::move(right_interpolation))
);
};
return interpolate(interpolate, 0, N, 0);
}
#undef vall
#endif /* FPS_OPERATION_HPP_ */
using mint = constant_modint<998244353>;
fps_operator<998244353> op;
mod_table<998244353> mtable(1000000);
fastio::FastInput in;
fastio::FastOutput out;
void solve(){
ll N,M,K;
in.read(N,M,K);
mint nCr = 1;
FormalPowerSeries<998244353> F(K+1);
for (ll i = 0; i <= min(M, K); i++){
F[i] = nCr*mint(M-i).pow(N);
nCr *= mint(M-i)*mtable.inverse_mod(i+1);
}
// out.print(F);
FormalPowerSeries<998244353> G(K+1);
G[0] = 0;
for (ll i = 1; i <= K; i++){
G[i] = mtable.factorialmodinv(i);
if (i%2){G[i] *= -1;}
}
// out.print(G);
F = op.composition(F, G, mtable);
for (ll i = 0; i <= K; i++){
G[i] = mint(M).pow(i)*mtable.factorialmodinv(i);
}
F = op.convolution(F,G);
for (ll i = 1; i <= K; i++){
out.print(F[i]*mtable.factorialmod(i));
}
}
int main(){
ios::sync_with_stdio(false);
std::cin.tie(nullptr);
ll T = 1;
// cin >> T;
while (T--){
solve();
}
}