結果
| 問題 | No.1145 Sums of Powers |
| コンテスト | |
| ユーザー |
 kakao743
|
| 提出日時 | 2026-05-01 23:44:11 |
| 言語 | C++23 (gcc 15.2.0 + boost 1.89.0) |
| 結果 |
AC
|
| 実行時間 | 217 ms / 2,000 ms |
| コード長 | 65,160 bytes |
| 記録 | |
| コンパイル時間 | 14,562 ms |
| コンパイル使用メモリ | 553,864 KB |
| 実行使用メモリ | 11,284 KB |
| 最終ジャッジ日時 | 2026-05-01 23:44:29 |
| 合計ジャッジ時間 | 15,188 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / judge3_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 6 |
ソースコード
#ifdef LOCAL
#include "my_header.h"
#else
#define PRAGMA_OPTIMIZE(s) _Pragma(#s)
PRAGMA_OPTIMIZE(GCC optimize("Ofast"))
PRAGMA_OPTIMIZE(GCC optimize("unroll-loops")) // ループ
#pragma GCC optimize("fast-math", "no-stack-protector")
// #pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,fma,tune=native")//浮動小数点
#pragma GCC target("avx2,bmi,bmi2,popcnt")//四則演算
#include <bits/stdc++.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <unistd.h>
using namespace std;
#if __has_include(<atcoder/all>)
#include <atcoder/all>
using namespace atcoder;
using minta = modint998244353;
using mintb = modint1000000007;
std::ostream&operator<<(std::ostream& os,const minta& v){os << v.val();return os;}
std::istream&operator>>(std::istream&is,minta &v){long long t;is >> t;v=t;return is;}
std::ostream&operator<<(std::ostream& os,const mintb& v){os << v.val();return os;}
std::istream&operator>>(std::istream&is,mintb &v){long long t;is >> t;v=t;return is;}
istream &operator>>(istream &is, __int128_t &x) {
string S;
is >> S;
x = 0;
int flag = 0;
for (auto &c : S) {
if (c == '-') {
flag = true;
continue;
}
x *= 10;
x += c - '0';
}
if (flag) x = -x;
return is;
}
ostream &operator<<(ostream &os, __int128_t x) {
if (x == 0) return os << 0;
if (x < 0) os << '-', x = -x;
string S;
while (x) S.push_back('0' + x % 10), x /= 10;
reverse(begin(S), end(S));
return os << S;
}
template<typename T1,typename T2>ostream&operator<<(ostream&os,const pair<T1,T2>& v){os << '(' << v.first << ',' << v.second << ')';return os;}
template<typename T1,typename T2>istream&operator>>(istream&is,pair<T1,T2> &v){is >> v.first >> v.second;return is;}
#define HAS_ACL 1
#else
#define HAS_ACL 0
#endif
using namespace std;
using ll=long long;
template <typename T1, typename T2>
struct ob2 {
T1 a; T2 b;
ob2() : a(), b() {}
ob2(T1 a, T2 b) : a(a), b(b) {}
friend auto operator<=>(const ob2&, const ob2&) = default;
ob2 operator+(const ob2& o) const { return {a + o.a, b + o.b}; }
ob2 operator-(const ob2& o) const { return {a - o.a, b - o.b}; }
friend std::ostream& operator<<(std::ostream& os, const ob2& p) {
return os << "(" << p.a << ", " << p.b << ")";
}
friend std::istream& operator>>(std::istream& is, ob2& p) {
return is >> p.a >> p.b;
}
};
template <typename T1, typename T2, typename T3>
struct ob3 {
T1 a; T2 b; T3 c;
ob3() : a(), b(), c() {}
ob3(T1 a, T2 b, T3 c) : a(a), b(b), c(c) {}
friend auto operator<=>(const ob3&, const ob3&) = default;
ob3 operator+(const ob3& o) const { return {a + o.a, b + o.b, c + o.c}; }
ob3 operator-(const ob3& o) const { return {a - o.a, b - o.b, c - o.c}; }
friend std::ostream& operator<<(std::ostream& os, const ob3& p) {
return os << "(" << p.a << ", " << p.b << ", " << p.c << ")";
}
friend std::istream& operator>>(std::istream& is, ob3& p) {
return is >> p.a >> p.b >> p.c;
}
};
template <typename T1, typename T2, typename T3, typename T4>
struct ob4 {
T1 a; T2 b; T3 c; T4 d;
ob4() : a(), b(), c(), d() {}
ob4(T1 a, T2 b, T3 c, T4 d) : a(a), b(b), c(c), d(d) {}
friend auto operator<=>(const ob4&, const ob4&) = default;
ob4 operator+(const ob4& o) const { return {a + o.a, b + o.b, c + o.c, d + o.d}; }
ob4 operator-(const ob4& o) const { return {a - o.a, b - o.b, c - o.c, d - o.d}; }
friend std::ostream& operator<<(std::ostream& os, const ob4& p) {
return os << "(" << p.a << ", " << p.b << ", " << p.c << ", " << p.d << ")";
}
friend std::istream& operator>>(std::istream& is, ob4& p) {
return is >> p.a >> p.b >> p.c >> p.d;
}
};
template <typename T1, typename T2, typename T3, typename T4, typename T5>
struct ob5 {
T1 a; T2 b; T3 c; T4 d; T5 e;
ob5() : a(), b(), c(), d(), e() {}
ob5(T1 a, T2 b, T3 c, T4 d, T5 e) : a(a), b(b), c(c), d(d), e(e) {}
friend auto operator<=>(const ob5&, const ob5&) = default;
ob5 operator+(const ob5& o) const { return {a + o.a, b + o.b, c + o.c, d + o.d, e + o.e}; }
ob5 operator-(const ob5& o) const { return {a - o.a, b - o.b, c - o.c, d - o.d, e - o.e}; }
friend std::ostream& operator<<(std::ostream& os, const ob5& p) {
return os << "(" << p.a << ", " << p.b << ", " << p.c << ", " << p.d << ", " << p.e << ")";
}
friend std::istream& operator>>(std::istream& is, ob5& p) {
return is >> p.a >> p.b >> p.c >> p.d >> p.e;
}
};
template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
struct ob6 {
T1 a; T2 b; T3 c; T4 d; T5 e; T6 f;
ob6() : a(), b(), c(), d(), e(), f() {}
ob6(T1 a, T2 b, T3 c, T4 d, T5 e, T6 f) : a(a), b(b), c(c), d(d), e(e), f(f) {}
friend auto operator<=>(const ob6&, const ob6&) = default;
ob6 operator+(const ob6& o) const { return {a + o.a, b + o.b, c + o.c, d + o.d, e + o.e, f + o.f}; }
ob6 operator-(const ob6& o) const { return {a - o.a, b - o.b, c - o.c, d - o.d, e - o.e, f - o.f}; }
friend std::ostream& operator<<(std::ostream& os, const ob6& p) {
return os << "(" << p.a << ", " << p.b << ", " << p.c << ", " << p.d << ", " << p.e << ", " << p.f << ")";
}
friend std::istream& operator>>(std::istream& is, ob6& p) {
return is >> p.a >> p.b >> p.c >> p.d >> p.e >> p.f;
}
};
using ll = long long;
using ll2 = ob2<ll, ll>;
using ll3 = ob3<ll, ll, ll>;
using ll4 = ob4<ll, ll, ll, ll>;
using ll5 = ob5<ll, ll, ll, ll, ll>;
using ll6 = ob6<ll, ll, ll, ll, ll, ll>;
template<typename T> struct is_ob : std::false_type {};
template<typename... Ts> struct is_ob<ob2<Ts...>> : std::true_type {};
template<typename... Ts> struct is_ob<ob3<Ts...>> : std::true_type {};
template<typename... Ts> struct is_ob<ob4<Ts...>> : std::true_type {};
template<typename... Ts> struct is_ob<ob5<Ts...>> : std::true_type {};
template<typename... Ts> struct is_ob<ob6<Ts...>> : std::true_type {};
namespace for_debugging{
struct subscript_and_location{
int sub;
std::source_location loc;
template<class T> subscript_and_location(T sub_,std::source_location loc_=std::source_location::current()){
if(!std::is_integral<T>::value){
std::clog << loc_.file_name() << ":(" << loc_.line() << ":" << loc_.column() << "):" << loc_.function_name() << std::endl;
std::clog << "subscript is not integer: subscript = " << sub_ << std::endl;
exit(EXIT_FAILURE);
}
sub=sub_;
loc=loc_;
}
void check_out_of_range(size_t sz){
if(sub<0||(int)sz<=sub){
std::clog << loc.file_name() << ":(" << loc.line() << ":" << loc.column() << "):" << loc.function_name() << std::endl;
std::clog << "out of range: subscript = " << sub << ", vector_size = " << sz << std::endl;
exit(EXIT_FAILURE);
}
}
};
}
namespace std{
template<class T,class Allocator=std::allocator<T>> class vector_for_debugging:public std::vector<T,Allocator>{
using std::vector<T,Allocator>::vector;
public:
[[nodiscard]] constexpr std::vector<T,Allocator>::reference operator[](for_debugging::subscript_and_location n) noexcept(!std::is_same<T,bool>::value){
n.check_out_of_range(this->size());
return std::vector<T,Allocator>::operator[](n.sub);
}
[[nodiscard]] constexpr std::vector<T,Allocator>::const_reference operator[](for_debugging::subscript_and_location n) const noexcept(!std::is_same<T,bool>::value){
n.check_out_of_range(this->size());
return std::vector<T,Allocator>::operator[](n.sub);
}
};
namespace pmr{
template<class T> using vector_for_debugging=std::vector_for_debugging<T,std::pmr::polymorphic_allocator<T>>;
}
}
#define vfd vector_for_debugging
/*//多倍長整数
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;
// 任意長整数型
using Bint = mp::cpp_int;
//仮数部が10進数で1024桁の浮動小数点数型(TLEしたら小さくする)
using Real = mp::number<mp::cpp_dec_float<1024>>;
*/
#define rep(i,n) for(long long i=0;i<(long long)n;i++)
#define reps(i,n) for(long long i=1;i<=(long long)n;i++)
#define repi(i,n) for(int i=0;i<(int)n;i++)
#define loop(i,l,r) for(long long i=(long long)l;i<=(long long)r;i++)
#define loopi(i,l,r) for(int i=(int)l;i<=(int)r;i++)
#define drep(i,n) for(long long i=(long long)n-1;i>=0;i--)
#define drepi(i,n) for(int i=(int)n-1;i>=0;i--)
#define dreps(i,n) for(int i=(int)n;i>=1;i--)
#define dloop(i,l,r) for(long long i=(long long)l;i>=(long long)r;i--)
#define dloopi(i,l,r) for(int i=(int)l;i>=(int)r;i--)
#define all(v) v.begin(), v.end()
#define rall(v) v.rbegin(), v.rend()
#define yna(x) cout << (x? "Yes":"No") << endl;
#define yn(x) out(bool(x));
#define cou(x) cout << x << endl;
#define emp emplace_back
#define mp make_pair
const long long moda=998244353LL;
const long long modb=1000000007LL;
const int kaz=1000000005;
long long yab=2500000000000000000LL;
const long long aho =-yab;
const long double eps=1.0e-14L;
const long double pi=acosl(-1.0L);
using st=string;
using P=pair<ll,ll>;
using tup=tuple<ll,ll,ll>;
using vi=vector<ll>;
using vin=vector<int>;
using vc=vector<char>;
using vb=vector<bool>;
using vd=vector<double>;
using vs=vector<string>;
using vp=vector<P>;
using sp=set<P>;
using si=set<ll>;
using vvi=vector<vector<ll>>;
using vvin=vector<vin>;
using vvc=vector<vc>;
using vvb=vector<vb>;
using vvvi=vector<vvi>;
using vvvin=vector<vvin>;
const int dx[4]={0,1,0,-1};
const int dy[4]={1,0,-1,0};
const vector<int> ex = {-1, -1, -1, 0, 0, 1, 1, 1};
const vector<int> ey = {-1, 0, 1, -1, 1, -1, 0, 1};
template<typename T>istream&operator>>(istream&is,vector<T>&v){for(T&in:v)is>>in;return is;}
template<typename T>ostream&operator<<(ostream&os,const vector<T>&v){for(int i=0;i<int(v.size());i++)os<<v[i]<<(i+1!=v.size()?" ":"\n");return os;}
struct FastIO {
char *p;
char out_buf[1 << 21];
int out_pos = 0;
std::vector<char> fallback_buf;
template<typename T> struct is_pair : std::false_type {};
template<typename T, typename U> struct is_pair<std::pair<T,U>> : std::true_type {};
FastIO() {
struct stat st;
if (isatty(0)) { p = nullptr; return; }
if (fstat(0, &st) == 0 && st.st_size > 0) {
void *res = mmap(0, st.st_size, PROT_READ, MAP_PRIVATE, 0, 0);
if (res != MAP_FAILED) { p = (char*)res; return; }
}
setup_fallback();
}
void setup_fallback() {
char buf[4096];
while (true) {
ssize_t len = read(0, buf, sizeof(buf));
if (len <= 0) break;
fallback_buf.insert(fallback_buf.end(), buf, buf + len);
}
if (!fallback_buf.empty()) { fallback_buf.push_back('\0'); p = fallback_buf.data(); }
else p = nullptr;
}
~FastIO() { flush(); }
void flush() { if (out_pos > 0) { write(1, out_buf, out_pos); out_pos = 0; } }
void write_char(char c) { out_buf[out_pos++] = c; if (out_pos >= (1<<20)) flush(); }
// ---------- 入力 ----------
template <typename T>
void read_recursive(T &x) {
if (!p) { cin >> x; return; }
if constexpr (is_same_v<T, bool>) {
long long t; read_recursive(t); x = (t != 0);
} else if constexpr (
#if HAS_ACL
is_same_v<T, minta> || is_same_v<T, mintb>
#else
false
#endif
) {
long long t; read_recursive(t); x = t;
} else if constexpr (is_same_v<T, __int128_t>) {
x = 0; bool neg = false;
while (p && *p && (*p < '0' || *p > '9') && *p != '-') p++;
if (!p || !*p) return;
if (*p == '-') { neg = true; p++; }
while (*p >= '0' && *p <= '9') x = x * (__int128_t)10 + (*p++ - '0');
if (neg) x = -x;
} else if constexpr (is_same_v<T, uint64_t>) {
x = 0;
while (p && *p && (*p < '0' || *p > '9')) p++;
if (!p || !*p) return;
while (*p >= '0' && *p <= '9') x = x * (uint64_t)10 + (uint64_t)(*p++ - '0');
} else if constexpr (is_same_v<T, char>) {
while (p && *p && *p <= ' ') p++;
if (p && *p) x = *p++;
} else if constexpr (is_integral_v<T>) {
x = 0; bool neg = false;
while (p && *p && (*p < '0' || *p > '9') && *p != '-') p++;
if (!p || !*p) return;
if (*p == '-') { neg = true; p++; }
while (*p >= '0' && *p <= '9') x = x * 10LL + (*p++ - '0');
if (neg) x = -x;
} else if constexpr (is_floating_point_v<T>) {
while (p && *p && *p <= ' ') p++;
if (!p || !*p) return;
char *end; x = (T)strtod(p, &end); p = end;
} else if constexpr (is_same_v<T, string>) {
x.clear();
while (p && *p && *p <= ' ') p++;
if (!p || !*p) return;
while (p && *p && *p > ' ') x += *p++;
} else if constexpr (is_pair<T>::value) {
read_recursive(x.first); read_recursive(x.second);
} else if constexpr (is_same_v<T, vector<char>>) {
// string と同様に空白区切りで1トークン読む
x.clear();
while (p && *p && *p <= ' ') p++;
if (!p || !*p) return;
while (p && *p && *p > ' ') x.push_back(*p++);
} else if constexpr (is_same_v<T, vector<bool>>) {
for (int i = 0; i < (int)x.size(); i++) { int t; read_recursive(t); x[i] = (bool)t; }
} else if constexpr (is_ob<T>::value) {
if constexpr (requires { x.a; }) read_recursive(x.a);
if constexpr (requires { x.b; }) read_recursive(x.b);
if constexpr (requires { x.c; }) read_recursive(x.c);
if constexpr (requires { x.d; }) read_recursive(x.d);
if constexpr (requires { x.e; }) read_recursive(x.e);
if constexpr (requires { x.f; }) read_recursive(x.f);
}
else if constexpr (requires { typename T::fp_io_tag; }) {
long long t; read_recursive(t); x = t;
}else {
for (auto &it : x) read_recursive(it);
}
}
// ---------- 出力 ----------
void write_int64(long long n) {
if (n == 0) { write_char('0'); return; }
if (n < 0) { write_char('-'); n = -n; }
char tmp[20]; int t = 0;
while (n > 0) { tmp[t++] = (char)((n % 10) + '0'); n /= 10; }
for (int i = t-1; i >= 0; i--) write_char(tmp[i]);
}
void write_uint64(uint64_t n) {
if (n == 0) { write_char('0'); return; }
char tmp[20]; int t = 0;
while (n > 0) { tmp[t++] = (char)((n % 10) + '0'); n /= 10; }
for (int i = t-1; i >= 0; i--) write_char(tmp[i]);
}
void write_int128(__int128_t n) {
if (n == 0) { write_char('0'); return; }
if (n < 0) { write_char('-'); n = -n; }
char tmp[40]; int t = 0;
while (n > 0) { tmp[t++] = (char)((n % 10) + '0'); n /= 10; }
for (int i = t-1; i >= 0; i--) write_char(tmp[i]);
}
void write_double(double x) {
if (std::isnan(x)) { write_char('0'); return; }
if (std::isinf(x)) { if (x < 0) write_char('-'); write_char('i'); write_char('n'); write_char('f'); return; }
if (x < 0) { write_char('-'); x = -x; }
if (x >= 9.2e18) {
write_uint64((uint64_t)x); write_char('.');
for (int i = 0; i < 15; i++) write_char('0');
return;
}
double offset = 0.5;
for (int i = 0; i < 15; i++) offset /= 10.0;
x += offset;
long long int_part = (long long)x;
write_int64(int_part); write_char('.');
double fraction = x - (double)int_part;
for (int i = 0; i < 15; i++) {
fraction *= 10; int d = (int)fraction;
write_char('0' + d); fraction -= d;
}
}
// const char* と char[] 用（セパレータなしで文字列として出力）
void write_recursive(const char *x) { while (*x) write_char(*x++); }
template <size_t N>
void write_recursive(const char (&x)[N]) { for (size_t i = 0; i+1 < N && x[i]; i++) write_char(x[i]); }
template <typename T>
void write_recursive(const T &x) {
if constexpr (is_same_v<T, bool>) {
if (x) { write_char('Y'); write_char('e'); write_char('s'); }
else { write_char('N'); write_char('o'); }
} else if constexpr (
#if HAS_ACL
is_same_v<T, minta> || is_same_v<T, mintb>
#else
false
#endif
) {
write_int64((long long)x.val());
} else if constexpr (is_same_v<T, __int128_t>) {
write_int128(x);
} else if constexpr (is_same_v<T, uint64_t>) {
write_uint64(x);
} else if constexpr (is_same_v<T, char>) {
write_char(x);
} else if constexpr (is_floating_point_v<T>) {
write_double((double)x);
} else if constexpr (is_integral_v<T>) {
write_int64((long long)x);
} else if constexpr (is_same_v<T, string> || is_same_v<T, vector<char>>) {
// vector<char> も string と同様セパレータなしで出力
for (char c : x) write_char(c);
} else if constexpr (is_pair<T>::value) {
write_char('('); write_recursive(x.first);
write_char(','); write_char(' '); write_recursive(x.second);
write_char(')');
} else if constexpr (is_ob<T>::value) {
write_char('(');
if constexpr (requires { x.a; }) { write_recursive(x.a); }
if constexpr (requires { x.b; }) { write_char(','); write_char(' '); write_recursive(x.b); }
if constexpr (requires { x.c; }) { write_char(','); write_char(' '); write_recursive(x.c); }
if constexpr (requires { x.d; }) { write_char(','); write_char(' '); write_recursive(x.d); }
if constexpr (requires { x.e; }) { write_char(','); write_char(' '); write_recursive(x.e); }
if constexpr (requires { x.f; }) { write_char(','); write_char(' '); write_recursive(x.f); }
write_char(')');
}
else if constexpr (requires { typename T::fp_io_tag; }) {
write_uint64(x.val());
}else {
for (int i = 0; i < (int)x.size(); i++) {
write_recursive(x[i]);
if (i+1 != (int)x.size()) write_char(' ');
}
}
}
} io;
template <typename... Args> void in(Args &...args) { (io.read_recursive(args), ...); }
template <typename... Args> void out(const Args &...args) {
int n = sizeof...(Args), i = 0;
( (io.write_recursive(args), io.write_char(++i == n ? '\n' : ' ')), ... );
}
template<typename T1,typename T2>
void co(bool x,T1 y,T2 z){
if(x)cout << y << endl;
else cout << z << endl;
}
long long isqrt(long long n){
long long ok=0,ng=1000000000;//1e9
while(ng-ok>1){
long long mid=(ng+ok)/2;
if(mid*mid<=n)ok=mid;
else ng=mid;
}
return ok;
}
template<typename T>
bool chmax(T &a, T b){
if(a<b){
a=b;
return true;
}
return false;
}
template <typename T, typename U>
T ceil(T x, U y) {
return (x > 0 ? (x + y - 1) / y : x / y);
}
template <typename T, typename U>
T floor(T x, U y) {
return (x > 0 ? x / y : (x - y + 1) / y);
}
template<typename T,typename U>
pair<pair<T,T>,T> unit(T x, T y,U k){
T s=floor(x-1,k);
T t=floor(y,k);
if(s==t)return make_pair(make_pair(0,0),-1);
return make_pair(make_pair(s+1,t),1);
}
pair<pair<ll,ll>,int> jufuku(ll a,ll b,ll c,ll d){
//a<=b,c<=dが保証されているとする
if(a>c){
swap(a,c);
swap(b,d);
}
if(c>b)return make_pair(make_pair(0,0),-1);
return make_pair(make_pair(c,min(b,d)),1);
}
template<typename T>
bool chmin(T &a, T b){
if(a>b){
a=b;
return true;
}
return false;
}
template<typename T>
void her(vector<T> &a){
for(auto &g:a)g--;
}
template <typename T>
void dec(vector<T> &t,T k=1){
for(auto &i:t)i-=k;
}
template <typename T>
void inc(vector<T> &t,T k=1){
for(auto &i:t)i+=k;
}
ll mypow(ll x,ll y,ll MOD){
if(MOD==-1){
MOD=9223372036854775807LL;
}
ll ret=1;
while(y>0){
if(y&1)ret=ret*x%MOD;
x=x*x%MOD;
y>>=1;
}
return ret;
}
struct UnionFind {
vector<int> par,siz,mi,ma;
UnionFind(int n) : par(n,-1), siz(n,1) {mi.resize(n);ma.resize(n);iota(mi.begin(),mi.end(),0);iota(ma.begin(),ma.end(),0); }
int root(int x) {
if(par[x]==-1)return x;
else return par[x]=root(par[x]);
}
bool same(int x, int y) {
return root(x)==root(y);
}
bool marge(int x, int y) {
int rx = root(x), ry = root(y);
if (rx==ry) return false;
if(siz[rx]<siz[ry])swap(rx,ry);
par[ry] = rx;
siz[rx] += siz[ry];
mi[rx]=min(mi[rx],mi[ry]);
ma[rx]=max(ma[rx],ma[ry]);
return true;
}
int size(int x) {
return siz[root(x)];
}
int mini(int x){
return mi[root(x)];
}
int maxi(int x){
return ma[root(x)];
}
};
ll nckmod;
vector<long long> fac,finv,invs;
// テーブルを作る前処理
void COMinit(int MAX,ll MOD) {
nckmod=MOD;
fac.resize(MAX);
finv.resize(MAX);
invs.resize(MAX);
fac[0] = fac[1] = 1;
finv[0] = finv[1] = 1;
invs[1] = 1;
for (int i = 2; i < MAX; i++){
fac[i] = fac[i - 1] * i % MOD;
invs[i] = MOD - invs[MOD%i] * (MOD / i) % MOD;
finv[i] = finv[i - 1] * invs[i] % MOD;
}
}
// 二項係数計算
long long binom(int n, int k){
if (n < k) return 0;
if (n < 0 || k < 0) return 0;
return fac[n] * (finv[k] * finv[n - k] % nckmod) % nckmod;
}
template<typename Container>
auto assyuku(const Container&v){
using T=typename Container::value_type;
vector<ob3<int,int,T>> ans;
if(v.empty())return ans;
int sum=0;
T pos=v[0];
int n=v.size();
for(int i=0;i<n;i++){
if(v[i]==pos)sum++;
else {
ans.emplace_back(ob3(i-sum,i-1,pos));
sum=1;
pos=v[i];
}
}
ans.emplace_back(ob3(n-sum,n-1,pos));
return ans;
};
template<typename Container,typename F>
vector<ob3<int,int,bool>> assyuku2(const Container &v,F f){
vector<ob3<int,int,bool>> ans;
if(v.empty())return ans;
int sum=0;
bool pos=f(v[0]);
int n=v.size();
for(int i=0;i<n;i++){
if(f(v[i])==pos)sum++;
else {
ans.emplace_back(ob3(i-sum,i-1,pos));
sum=1;
pos=f(v[i]);
}
}
ans.emplace_back(ob3(n-sum,n-1,pos));
return ans;
};
#endif
//けんちょんさんの
// modint
//NTTのmulのあれは絶対にやろう！
// min non-negative i such that (x & (1 << i)) != 0
int bsf(int x) { return __builtin_ctz(x); }
int bsf(unsigned int x) { return __builtin_ctz(x); }
int bsf(long long x) { return __builtin_ctzll(x); }
int bsf(unsigned long long x) { return __builtin_ctzll(x); }
// max non-negative i such that (x & (1 << i)) != 0
int bsr(int x) { return 8 * (int)sizeof(int) - 1 - __builtin_clz(x); }
int bsr(unsigned int x) { return 8 * (int)sizeof(unsigned int) - 1 - __builtin_clz(x); }
int bsr(long long x) { return 8 * (int)sizeof(long long) - 1 - __builtin_clzll(x); }
int bsr(unsigned long long x) { return 8 * (int)sizeof(unsigned long long) - 1 - __builtin_clzll(x); }
// ============================================================
// ★追加: mod 2^64 演算用 FpU64 (FastIO より前に定義が必要)
// Fp<-1> として使う。+,-,* はすべて mod 2^64 (自然なwrap around)
// ============================================================
struct FpU64 {
using fp_io_tag = void;
uint64_t _v;
constexpr FpU64() : _v(0) {}
constexpr FpU64(long long v) : _v((uint64_t)v) {}
constexpr uint64_t val() const { return _v; }
constexpr uint64_t get() const { return _v; }
static constexpr int get_mod() { return -1; }
constexpr FpU64 operator+() const { return *this; }
constexpr FpU64 operator-() const { return FpU64((long long)(-_v)); }
constexpr FpU64 operator+(const FpU64 &r) const { return FpU64((long long)(_v + r._v)); }
constexpr FpU64 operator-(const FpU64 &r) const { return FpU64((long long)(_v - r._v)); }
constexpr FpU64 operator*(const FpU64 &r) const { return FpU64((long long)(_v * r._v)); }
constexpr FpU64& operator+=(const FpU64 &r) { _v += r._v; return *this; }
constexpr FpU64& operator-=(const FpU64 &r) { _v -= r._v; return *this; }
constexpr FpU64& operator*=(const FpU64 &r) { _v *= r._v; return *this; }
constexpr bool operator==(const FpU64 &r) const { return _v == r._v; }
constexpr bool operator!=(const FpU64 &r) const { return _v != r._v; }
// fft_info<FpU64> のインスタンス化を通すためのダミー (if constexpr で実際には呼ばれない)
constexpr FpU64 pow(long long) const { return *this; }
constexpr FpU64 inv() const { return *this; }
friend istream& operator>>(istream &is, FpU64 &x) { uint64_t v; is >> v; x._v = v; return is; }
friend ostream& operator<<(ostream &os, const FpU64 &x) { return os << x._v; }
};
template<int MOD>
struct FpEven {
using fp_io_tag = void;
long long _v;
constexpr FpEven() : _v(0) {}
constexpr FpEven(long long v) : _v(v % MOD) {
if (_v < 0) _v += MOD;
}
// val() : 真の値を返す (偶数MODでは内部値がそのまま真の値)
constexpr long long val() const { return _v; }
constexpr long long get() const { return val(); }
static constexpr int get_mod() { return MOD; }
constexpr FpEven operator+() const { return FpEven(*this); }
constexpr FpEven operator-() const { return FpEven(0) - FpEven(*this); }
constexpr FpEven operator+(const FpEven &r) const { return FpEven(*this) += r; }
constexpr FpEven operator-(const FpEven &r) const { return FpEven(*this) -= r; }
constexpr FpEven operator*(const FpEven &r) const { return FpEven(*this) *= r; }
constexpr FpEven operator/(const FpEven &r) const { return FpEven(*this) /= r; }
constexpr FpEven& operator+=(const FpEven &r) {
_v += r._v;
if (_v >= MOD) _v -= MOD;
return *this;
}
constexpr FpEven& operator-=(const FpEven &r) {
_v -= r._v;
if (_v < 0) _v += MOD;
return *this;
}
constexpr FpEven& operator*=(const FpEven &r) {
_v = _v * r._v % MOD;
return *this;
}
constexpr FpEven& operator/=(const FpEven &r) {
long long a = r._v, b = MOD, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
_v = _v * u % MOD;
if (_v < 0) _v += MOD;
return *this;
}
constexpr FpEven pow(long long n) const {
FpEven res(1), mul(*this);
while (n > 0) {
if (n & 1) res *= mul;
mul *= mul;
n >>= 1;
}
return res;
}
constexpr FpEven inv() const {
return FpEven(1) / (*this);
}
constexpr FpEven sqrt() {
if (val() < 2) return *this;
if (pow((MOD - 1) >> 1).val() != 1) return FpEven(0);
FpEven b = 1, one = 1;
while (b.pow((MOD - 1) >> 1) == FpEven(1)) b += one;
long long m = MOD - 1, e = 0;
while (m % 2 == 0) m >>= 1, e += 1;
FpEven x = pow((m - 1) >> 1);
FpEven y = (*this) * x * x;
x *= (*this);
FpEven z = b.pow(m);
while (y.val() != 1) {
int j = 0;
FpEven t = y;
while (t != one) j += 1, t *= t;
z = z.pow(1LL << (e - j - 1));
x *= z; z *= z; y *= z; e = j;
}
return x;
}
constexpr bool operator==(const FpEven &r) const { return _v == r._v; }
constexpr bool operator!=(const FpEven &r) const { return _v != r._v; }
constexpr FpEven& operator++() {
++_v;
if (_v >= MOD) _v -= MOD;
return *this;
}
constexpr FpEven& operator--() {
if (_v == 0) _v += MOD;
--_v;
return *this;
}
constexpr FpEven operator++(int) {
FpEven res = *this; ++(*this); return res;
}
constexpr FpEven operator--(int) {
FpEven res = *this; --(*this); return res;
}
friend istream& operator>>(istream &is, FpEven<MOD> &x) {
long long v; is >> v;
x = FpEven(v);
return is;
}
friend ostream& operator<<(ostream &os, const FpEven<MOD> &x) {
return os << x.val();
}
friend FpEven<MOD> pow(const FpEven<MOD> &r, long long n) { return r.pow(n); }
friend FpEven<MOD> inv(const FpEven<MOD> &r) { return r.inv(); }
};
// ============================================================
// 奇数MOD用 FpOdd (Montgomery乗算)
// ============================================================
template<int MOD>
struct FpOdd {
static_assert(MOD % 2 == 1, "MOD must be odd for Montgomery");
using fp_io_tag = void;
using u32 = uint32_t;
using u64 = uint64_t;
using i64 = int64_t;
// -MOD^{-1} mod 2^32 (Newton法: 4回で32bit収束)
static constexpr u32 calc_mod_inv() {
u32 x = (u32)get_mod();
for (int i = 0; i < 4; i++) x *= 2 - (u32)get_mod() * x;
return (u32)(-(u64)x);
}
// R^2 mod MOD (R = 2^32)
static constexpr u32 calc_R2() {
u64 r = ((u64)1 << 32) % (u32)get_mod();
return (u32)(r * r % (u32)get_mod());
}
static constexpr int get_mod() { return MOD; }
static constexpr u32 MOD_INV = calc_mod_inv();
static constexpr u32 R2 = calc_R2();
u32 _v; // Montgomery表現: a * R mod MOD
// Montgomery reduction: T * R^{-1} mod MOD
static constexpr u32 reduce(u64 T) {
u32 m = (u32)T * MOD_INV;
u32 t = (u32)((T + (u64)m * (u32)get_mod()) >> 32);
return t >= (u32)get_mod() ? t - (u32)get_mod() : t;
}
// --- constructors ---
constexpr FpOdd() : _v(0) {}
constexpr FpOdd(long long v) {
v %= get_mod();
if (v < 0) v += get_mod();
_v = reduce((u64)(u32)v * R2); // aR mod MOD
}
static constexpr FpOdd raw(u32 v) { FpOdd x; x._v = v; return x; }
constexpr long long val() const { return reduce(_v); }
constexpr long long get() const { return val(); }
constexpr FpOdd operator+() const { return *this; }
constexpr FpOdd operator-() const { return FpOdd(0) - (*this); }
constexpr FpOdd operator+(const FpOdd &r) const { return FpOdd(*this) += r; }
constexpr FpOdd operator-(const FpOdd &r) const { return FpOdd(*this) -= r; }
constexpr FpOdd operator*(const FpOdd &r) const { return FpOdd(*this) *= r; }
constexpr FpOdd operator/(const FpOdd &r) const { return FpOdd(*this) /= r; }
constexpr FpOdd& operator+=(const FpOdd &r) {
_v += r._v;
if (_v >= (u32)get_mod()) _v -= (u32)get_mod();
return *this;
}
constexpr FpOdd& operator-=(const FpOdd &r) {
if (_v < r._v) _v += (u32)get_mod();
_v -= r._v;
return *this;
}
constexpr FpOdd& operator*=(const FpOdd &r) {
_v = reduce((u64)_v * r._v);
return *this;
}
constexpr FpOdd& operator/=(const FpOdd &r) {
return *this *= r.inv();
}
constexpr FpOdd pow(long long n) const {
FpOdd res(1), mul(*this);
while (n > 0) {
if (n & 1) res *= mul;
mul *= mul;
n >>= 1;
}
return res;
}
constexpr FpOdd inv() const {
long long a = val(), b = get_mod(), u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
return FpOdd(u);
}
constexpr FpOdd sqrt() {
if (val() < 2) return *this;
if (pow((get_mod() - 1) >> 1).val() != 1) return FpOdd(0);
FpOdd b(1), one(1);
while (b.pow((get_mod() - 1) >> 1) == FpOdd(1)) b += one;
long long m = get_mod() - 1, e = 0;
while (m % 2 == 0) m >>= 1, e += 1;
FpOdd x = pow((m - 1) >> 1);
FpOdd y = (*this) * x * x;
x *= (*this);
FpOdd z = b.pow(m);
while (y.val() != 1) {
int j = 0;
FpOdd t = y;
while (t != one) j += 1, t *= t;
z = z.pow(1LL << (e - j - 1));
x *= z; z *= z; y *= z; e = j;
}
return x;
}
constexpr bool operator==(const FpOdd &r) const { return _v == r._v; }
constexpr bool operator!=(const FpOdd &r) const { return _v != r._v; }
constexpr FpOdd& operator++() { return *this += FpOdd(1); }
constexpr FpOdd& operator--() { return *this -= FpOdd(1); }
constexpr FpOdd operator++(int) { FpOdd res = *this; ++(*this); return res; }
constexpr FpOdd operator--(int) { FpOdd res = *this; --(*this); return res; }
friend istream& operator>>(istream &is, FpOdd<MOD> &x) {
long long v; is >> v;
x = FpOdd(v);
return is;
}
friend ostream& operator<<(ostream &os, const FpOdd<MOD> &x) {
return os << x.val();
}
friend FpOdd<MOD> pow(const FpOdd<MOD> &r, long long n) { return r.pow(n); }
friend FpOdd<MOD> inv(const FpOdd<MOD> &r) { return r.inv(); }
};
// ============================================================
// Fp<MOD> : MODの偶奇で自動切り替え
// ★追加: Fp<-1> → FpU64 (mod 2^64)
template<int MOD, bool is_odd = (MOD % 2 == 1)>
struct FpSelector { using type = FpOdd<MOD>; };
template<int MOD>
struct FpSelector<MOD, false> { using type = FpEven<MOD>; };
// ★追加: Fp<-1> → FpU64
template<>
struct FpSelector<-1, false> { using type = FpU64; };
template<int MOD>
using Fp = typename FpSelector<MOD>::type;
namespace NTT {
// --- 基本ユーティリティ ---
long long modpow(long long a, long long n, int mod) {
long long res = 1;
a %= mod;
while (n > 0) {
if (n & 1) res = res * a % mod;
a = a * a % mod;
n >>= 1;
}
return res;
}
long long modinv(long long a, int mod) {
long long b = mod, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; std::swap(a, b);
u -= t * v; std::swap(u, v);
}
u %= mod;
if (u < 0) u += mod;
return u;
}
// 2のべき乗への切り上げ (旧 get_fft_size の代わり)
int bit_ceil(int n) {
int x = 1;
while (x < n) x <<= 1;
return x;
}
constexpr int calc_primitive_root(long long mod) {
if (mod == 1) return -1;
if (mod == 2) return 1;
if (mod == 998244353) return 3;
if (mod == 167772161) return 3;
if (mod == 469762049) return 3;
if (mod == 754974721) return 11;
if (mod == 645922817) return 3;
if (mod == 897581057) return 3;
if (mod == 2013265921) return 31;
long long divs[20] = {};
divs[0] = 2;
long long cnt = 1;
long long x = (mod - 1) / 2;
while (x % 2 == 0) x /= 2;
for (long long i = 3; i * i <= x; i += 2) {
if (x % i == 0) {
divs[cnt++] = (int)i;
while (x % i == 0) x /= i;
}
}
if (x > 1) divs[cnt++] = (int)x;
for (int g = 2;; g++) {
bool ok = true;
for (int i = 0; i < cnt; i++) {
if (modpow(g, (mod - 1) / divs[i], mod) == 1) {
ok = false; break;
}
}
if (ok) return g;
}
}
// --- NTT 内部構造 ---
template <class mint> struct fft_info {
static constexpr int rank2 = __builtin_ctz(mint::get_mod() - 1);
mint root[rank2 + 1], iroot[rank2 + 1];
mint rate2[std::max(0, rank2 - 2 + 1)], irate2[std::max(0, rank2 - 2 + 1)];
mint rate3[std::max(0, rank2 - 3 + 1)], irate3[std::max(0, rank2 - 3 + 1)];
fft_info() {
int g = calc_primitive_root(mint::get_mod());
root[rank2] = mint(g).pow((mint::get_mod() - 1) >> rank2);
iroot[rank2] = root[rank2].inv();
for (int i = rank2 - 1; i >= 0; i--) {
root[i] = root[i + 1] * root[i + 1];
iroot[i] = iroot[i + 1] * iroot[i + 1];
}
mint prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - 2; i++) {
rate2[i] = root[i + 2] * prod;
irate2[i] = iroot[i + 2] * iprod;
prod *= iroot[i + 2];
iprod *= root[i + 2];
}
prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - 3; i++) {
rate3[i] = root[i + 3] * prod;
irate3[i] = iroot[i + 3] * iprod;
prod *= iroot[i + 3];
iprod *= root[i + 3];
}
}
};
template <class mint> void butterfly(std::vector<mint>& a) {
int n = int(a.size());
int h = __builtin_ctz((unsigned int)n);
static const fft_info<mint> info;
int len = 0;
while (len < h) {
if (h - len == 1) {
int p = 1 << (h - len - 1);
mint rot = 1;
for (int s = 0; s < (1 << len); s++) {
int offset = s << (h - len);
for (int i = 0; i < p; i++) {
auto l = a[i + offset], r = a[i + offset + p] * rot;
a[i + offset] = l + r; a[i + offset + p] = l - r;
}
if (s + 1 != (1 << len)) rot *= info.rate2[__builtin_ctz(~(unsigned int)(s))];
}
len++;
} else {
int p = 1 << (h - len - 2);
mint rot = 1, imag = info.root[2];
for (int s = 0; s < (1 << len); s++) {
mint rot2 = rot * rot, rot3 = rot2 * rot;
int offset = s << (h - len);
for (int i = 0; i < p; i++) {
auto a0 = a[i + offset], a1 = a[i + offset + p] * rot;
auto a2 = a[i + offset + 2 * p] * rot2, a3 = a[i + offset + 3 * p] * rot3;
auto a1na3imag = (a1 - a3) * imag;
a[i + offset] = a0 + a2 + a1 + a3;
a[i + offset + p] = a0 + a2 - (a1 + a3);
a[i + offset + 2 * p] = a0 - a2 + a1na3imag;
a[i + offset + 3 * p] = a0 - a2 - a1na3imag;
}
if (s + 1 != (1 << len)) rot *= info.rate3[__builtin_ctz(~(unsigned int)(s))];
}
len += 2;
}
}
}
template <class mint> void butterfly_inv(std::vector<mint>& a) {
int n = int(a.size());
int h = __builtin_ctz((unsigned int)n);
static const fft_info<mint> info;
int len = h;
while (len) {
if (len == 1) {
int p = 1 << (h - len);
mint irot = 1;
for (int s = 0; s < (1 << (len - 1)); s++) {
int offset = s << (h - len + 1);
for (int i = 0; i < p; i++) {
auto l = a[i + offset], r = a[i + offset + p];
a[i + offset] = l + r; a[i + offset + p] = (l - r) * irot;
}
if (s + 1 != (1 << (len - 1))) irot *= info.irate2[__builtin_ctz(~(unsigned int)(s))];
}
len--;
} else {
int p = 1 << (h - len);
mint irot = 1, iimag = info.iroot[2];
for (int s = 0; s < (1 << (len - 2)); s++) {
mint irot2 = irot * irot, irot3 = irot2 * irot;
int offset = s << (h - len + 2);
for (int i = 0; i < p; i++) {
auto a0 = a[i + offset], a1 = a[i + offset + p];
auto a2 = a[i + offset + 2 * p], a3 = a[i + offset + 3 * p];
auto a2na3iimag = (a2 - a3) * iimag;
a[i + offset] = a0 + a1 + a2 + a3;
a[i + offset + p] = (a0 - a1 + a2na3iimag) * irot;
a[i + offset + 2 * p] = (a0 + a1 - a2 - a3) * irot2;
a[i + offset + 3 * p] = (a0 - a1 - a2na3iimag) * irot3;
}
if (s + 1 != (1 << (len - 2))) irot *= info.irate3[__builtin_ctz(~(unsigned int)(s))];
}
len -= 2;
}
}
}
// power_projection 専用: 正規化込みの変換
template<class mint>
void trans(std::vector<mint>& v) {
butterfly(v);
}
template<class mint>
void trans_inv(std::vector<mint>& v) {
butterfly_inv(v);
int n = (int)v.size();
mint iv = mint(n).inv();
for (auto& x : v) x *= iv;
}
// --- 畳み込みの実装 ---
template<class T> std::vector<T> naive_mul(const std::vector<T> &A, const std::vector<T> &B) {
if (A.empty() || B.empty()) return {};
int N = (int)A.size(), M = (int)B.size();
std::vector<T> res(N + M - 1);
for (int i = 0; i < N; ++i)
for (int j = 0; j < M; ++j) res[i + j] += A[i] * B[j];
return res;
}
// ★追加: mod 2^64 畳み込み (32bit split + 3-mod NTT)
// a[i], b[i] を上位32bit/下位32bitに分割し、各サブ畳み込みをGarnerで復元
// 各サブ係数の最大値 <= N * (2^32)^2 = N*2^64 で、N<=2^19 なら < p0*p1*p2 ≈ 2^89
// ※ mul より前に定義するため、butterfly/butterfly_inv を直接呼ぶ
std::vector<FpU64> mul_u64(std::vector<FpU64> a, std::vector<FpU64> b) {
int n = (int)a.size(), m = (int)b.size();
if (!n || !m) return {};
int sz = n + m - 1;
static constexpr int M0 = 167772161; // 2^25
static constexpr int M1 = 469762049; // 2^26
static constexpr int M2 = 2013265921; // g=31, 2^27
using mint0 = Fp<M0>; using mint1 = Fp<M1>; using mint2 = Fp<M2>;
// NTT-friendly素数での単一畳み込みヘルパ (butterfly を直接利用)
auto ntt_conv = [&]<class mint>(std::vector<mint> x, std::vector<mint> y) {
int nx = (int)x.size(), ny = (int)y.size();
int z = bit_ceil(nx + ny - 1);
x.resize(z); butterfly(x);
y.resize(z); butterfly(y);
for (int i = 0; i < z; i++) x[i] *= y[i];
butterfly_inv(x);
x.resize(nx + ny - 1);
mint iz = mint(z).inv();
for (auto& v : x) v *= iz;
return x;
};
// a, b を lo(下位32bit) / hi(上位32bit) に分割
std::vector<mint0> alo0(n),ahi0(n),blo0(m),bhi0(m);
std::vector<mint1> alo1(n),ahi1(n),blo1(m),bhi1(m);
std::vector<mint2> alo2(n),ahi2(n),blo2(m),bhi2(m);
for (int i = 0; i < n; ++i) {
uint32_t lo = (uint32_t)(a[i]._v), hi = (uint32_t)(a[i]._v >> 32);
alo0[i]=lo; ahi0[i]=hi; alo1[i]=lo; ahi1[i]=hi; alo2[i]=lo; ahi2[i]=hi;
}
for (int i = 0; i < m; ++i) {
uint32_t lo = (uint32_t)(b[i]._v), hi = (uint32_t)(b[i]._v >> 32);
blo0[i]=lo; bhi0[i]=hi; blo1[i]=lo; bhi1[i]=hi; blo2[i]=lo; bhi2[i]=hi;
}
// 3つの素数でサブ畳み込み (ll=lo*lo, lh=lo*hi, hl=hi*lo)
// hh=hi*hi は mod 2^64 で消える (係数に 2^64 がかかる)
auto cll0=ntt_conv.template operator()<mint0>(alo0,blo0);
auto clh0=ntt_conv.template operator()<mint0>(alo0,bhi0);
auto chl0=ntt_conv.template operator()<mint0>(ahi0,blo0);
auto cll1=ntt_conv.template operator()<mint1>(alo1,blo1);
auto clh1=ntt_conv.template operator()<mint1>(alo1,bhi1);
auto chl1=ntt_conv.template operator()<mint1>(ahi1,blo1);
auto cll2=ntt_conv.template operator()<mint2>(alo2,blo2);
auto clh2=ntt_conv.template operator()<mint2>(alo2,bhi2);
auto chl2=ntt_conv.template operator()<mint2>(ahi2,blo2);
// Garnerで各サブ係数の真値を復元してから mod 2^64 で合成
// 各サブ係数 <= N*(2^32-1)^2 < 2^83 < M0*M1*M2 ≈ 2^89 → Garner正確
static const mint1 iM0_1 = mint1(M0).inv();
static const mint2 iM0_2 = mint2(M0).inv();
static const mint2 iM1_2 = mint2(M1).inv();
// Garner係数 (mod 2^64 での復元用)
static const uint64_t m0 = (uint64_t)M0;
static const uint64_t m01 = (uint64_t)M0 * (uint64_t)M1; // mod 2^64 自動
auto garner3_u64 = [&](mint0 v0, mint1 v1, mint2 v2) -> uint64_t {
// x ≡ v0 (mod M0), x ≡ v1 (mod M1), x ≡ v2 (mod M2)
// CRT → x = y0 + y1*M0 + y2*M0*M1 (mod M0*M1*M2)
uint64_t y0 = (uint64_t)(int)v0.val();
uint64_t y1 = (uint64_t)(int)((v1 - mint1((long long)y0)) * iM0_1).val();
uint64_t y2 = (uint64_t)(int)((v2 - mint2((long long)y0) - mint2((long long)y1) * mint2(M0)) * iM0_2 * iM1_2).val();
return y0 + y1 * m0 + y2 * m01; // mod 2^64 自動
};
std::vector<FpU64> res(sz);
for (int i = 0; i < sz; ++i) {
uint64_t vll = garner3_u64(cll0[i], cll1[i], cll2[i]);
uint64_t vlh = garner3_u64(clh0[i], clh1[i], clh2[i]);
uint64_t vhl = garner3_u64(chl0[i], chl1[i], chl2[i]);
// C[k] = vll + (vlh + vhl) * 2^32 (hh項は 2^64 係数なので消える)
res[i]._v = vll + ((vlh + vhl) << 32);
}
return res;
}
template<class mint> std::vector<mint> mul(std::vector<mint> a, std::vector<mint> b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
if (std::min(n, m) <= 60) return naive_mul(a, b);
int MOD = mint::get_mod();
// ★追加: MOD==-1 は FpU64 → mod 2^64 畳み込み
// if constexpr で butterfly<FpU64> のインスタンス化を防ぐ
if constexpr (mint::get_mod() == -1) {
// mint == FpU64 が保証されているのでreinterpret_castは安全
auto *pa = reinterpret_cast<std::vector<FpU64>*>(&a);
auto *pb = reinterpret_cast<std::vector<FpU64>*>(&b);
auto rc = mul_u64(std::move(*pa), std::move(*pb));
return *reinterpret_cast<std::vector<mint>*>(&rc);
}
int z = bit_ceil(n + m - 1);
if ((MOD - 1) % z == 0) {//NTT-Friendlyの場合
a.resize(z); butterfly(a);
b.resize(z); butterfly(b);
for (int i = 0; i < z; i++) a[i] *= b[i];
butterfly_inv(a);
a.resize(n+m-1);
mint iz = mint(z).inv();
for (auto& x : a) x *= iz;
return a;
}
static constexpr int MOD0 = 167772161; // 2^25
static constexpr int MOD1 = 469762049; // 2^26
static constexpr int MOD2 = 2013265921; // g=31, 2^27
using mint0 = Fp<MOD0>; using mint1 = Fp<MOD1>; using mint2 = Fp<MOD2>;
std::vector<mint0> a0(n), b0(m); std::vector<mint1> a1(n), b1(m); std::vector<mint2> a2(n), b2(m);
for(int i=0; i<n; ++i) { a0[i] = a[i].val(); a1[i] = a[i].val(); a2[i] = a[i].val(); }
for(int i=0; i<m; ++i) { b0[i] = b[i].val(); b1[i] = b[i].val(); b2[i] = b[i].val(); }
auto c0 = mul(std::move(a0), std::move(b0));
auto c1 = mul(std::move(a1), std::move(b1));
auto c2 = mul(std::move(a2), std::move(b2));
static const mint1 imod0 = mint1(MOD0).inv();
static const mint2 imod1 = mint2(MOD1).inv(), imod01 = imod1 * mint2(MOD0).inv();
int res_size = n + m - 1;
std::vector<mint> res(res_size);
long long m0 = MOD0 % MOD, m01 = (long long)MOD0 * MOD1 % MOD;
for (int i = 0; i < res_size; ++i) {
int y0 = c0[i].val();
int y1 = ((c1[i] - y0) * imod0).val();
int y2 = ((c2[i] - y0) * imod01 - mint2(y1) * imod1).val();
res[i] = (long long)y0 + (long long)y1 * m0 + (long long)y2 * m01;
}
return res;
}
template<class mint> void ntt(vector<mint>&a) {
butterfly(a);
}
template<class mint> void intt(vector<mint>&a) {
butterfly_inv(a);
mint iz = mint(a.size()).inv();
for (auto& x : a) x *= iz;
}
template<class mint> void ntt_doubling(vector<mint> &a) {
int M = (int)a.size();
auto b = a;
intt(b);
mint r = 1, zeta = mint(calc_primitive_root(mint::get_mod())).pow((mint::get_mod() - 1) / (M << 1));
for (int i = 0; i < M; i++) b[i] *= r, r *= zeta;
ntt(b);
copy(begin(b), end(b), back_inserter(a));
}
};
template<class T> struct BiCoef {
vector<T> fact_, inv_, finv_;
constexpr BiCoef() {}
constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
init(n);
}
constexpr void init(int n) noexcept {
fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
int MOD = fact_[0].get_mod();
for(int i = 2; i < n; i++){
fact_[i] = fact_[i-1] * i;
inv_[i] = -inv_[MOD%i] * (MOD/i);
finv_[i] = finv_[i-1] * inv_[i];
}
}
constexpr T com(int n, int k) const noexcept {
if (n < k || n < 0 || k < 0) return 0;
return fact_[n] * finv_[k] * finv_[n-k];
}
constexpr T fact(int n) const noexcept {
if (n < 0) return 0;
return fact_[n];
}
constexpr T inv(int n) const noexcept {
if (n < 0) return 0;
return inv_[n];
}
constexpr T finv(int n) const noexcept {
if (n < 0) return 0;
return finv_[n];
}
};
// Formal Power Series
template<typename mint> struct FPS : vector<mint> {
using vector<mint>::vector;
// constructor
constexpr FPS(const vector<mint> &r) : vector<mint>(r) {}
// core operator
constexpr FPS pre(int siz) const {
return FPS(begin(*this), begin(*this) + min((int)this->size(), siz));
}
constexpr FPS rev() const {
FPS res = *this;
reverse(begin(res), end(res));
return res;
}
constexpr FPS& normalize() {
while (!this->empty() && this->back() == 0) this->pop_back();
return *this;
}
// basic operator
constexpr FPS operator - () const noexcept {
FPS res = (*this);
for (int i = 0; i < (int)res.size(); ++i) res[i] = -res[i];
return res;
}
constexpr FPS operator + (const mint &v) const { return FPS(*this) += v; }
constexpr FPS operator + (const FPS &r) const { return FPS(*this) += r; }
constexpr FPS operator - (const mint &v) const { return FPS(*this) -= v; }
constexpr FPS operator - (const FPS &r) const { return FPS(*this) -= r; }
constexpr FPS operator * (const mint &v) const { return FPS(*this) *= v; }
constexpr FPS operator * (const FPS &r) const { return FPS(*this) *= r; }
constexpr FPS operator / (const mint &v) const { return FPS(*this) /= v; }
constexpr FPS operator / (const FPS &r) const { return FPS(*this) /= r; }
constexpr FPS operator % (const FPS &r) const { return FPS(*this) %= r; }
constexpr FPS operator << (int x) const { return FPS(*this) <<= x; }
constexpr FPS operator >> (int x) const { return FPS(*this) >>= x; }
constexpr FPS& operator += (const mint &v) {
if (this->empty()) this->resize(1);
(*this)[0] += v;
return this->normalize();
}
constexpr FPS& operator += (const FPS &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); ++i) (*this)[i] += r[i];
return *this;
}
constexpr FPS& operator -= (const mint &v) {
if (this->empty()) this->resize(1);
(*this)[0] -= v;
return *this;
}
constexpr FPS& operator -= (const FPS &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); ++i) (*this)[i] -= r[i];
return this ->normalize();
}
constexpr FPS& operator *= (const mint &v) {
for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= v;
return *this;
}
constexpr FPS& operator *= (const FPS &r) {
return *this = NTT::mul((*this), r);
}
constexpr FPS& operator*=(vector<pair<int, mint>> g) {
using T=mint;
int n = (*this).size();
auto [d, c] = g.front();
if (d == 0) g.erase(g.begin());
else c = 0;
drep(i, n) {
(*this)[i] *= c;
for (auto &[j, b] : g) {
if (j > i) break;
(*this)[i] += (*this)[i-j] * b;
}
}
return *this;
}
constexpr FPS& operator/=(vector<pair<int, mint>> g) {
using T=mint;
int n = (*this).size();
auto [d, c] = g.front();
assert(d == 0 && c != T(0));
T ic = c.inv();
g.erase(g.begin());
rep(i, n) {
for (auto &[j, b] : g) {
if (j > i) break;
(*this)[i] -= (*this)[i-j] * b;
}
(*this)[i] *= ic;
}
return *this;
}
constexpr FPS& operator /= (const mint &v) {
assert(v != 0);
mint iv = v.inv();
for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= iv;
return *this;
}
// division, r must be normalized (r.back() must not be 0)
constexpr FPS& operator /= (const FPS &r) {
assert(!r.empty());
assert(r.back() != 0);
this->normalize();
if (this->size() < r.size()) {
this->clear();
return *this;
}
int need = (int)this->size() - (int)r.size() + 1;
*this = (rev().pre(need) * r.rev().inv(need)).pre(need).rev();
return *this;
}
constexpr FPS& operator %= (const FPS &r) {
assert(!r.empty());
assert(r.back() != 0);
this->normalize();
FPS q = (*this) / r;
return *this -= q * r;
}
constexpr FPS& operator <<= (int x) {
FPS res(x, 0);
res.insert(res.end(), begin(*this), end(*this));
return *this = res;
}
constexpr FPS& operator >>= (int x) {
FPS res;
res.insert(res.end(), begin(*this) + x, end(*this));
return *this = res;
}
constexpr mint eval(const mint &v) const {
mint res = 0;
for (int i = (int)this->size()-1; i >= 0; --i) {
res *= v;
res += (*this)[i];
}
return res;
}
// advanced operation
// df/dx
constexpr FPS diff() const {
int n = (int)this->size();
FPS res(n-1);
for (int i = 1; i < n; ++i) res[i-1] = (*this)[i] * i;
return res;
}
// \int f dx
constexpr FPS integral() const {
int n = (int)this->size();
FPS res(n+1, 0);
for (int i = 0; i < n; ++i) res[i+1] = (*this)[i] / (i+1);
return res;
}
// inv(f), f[0] must not be 0
constexpr FPS inv(int deg) const {
assert((*this)[0] != 0);
if (deg < 0) deg = (int)this->size();
FPS res({mint(1) / (*this)[0]});
for (int i = 1; i < deg; i <<= 1) {
res = (res + res - res * res * pre(i << 1)).pre(i << 1);
}
res.resize(deg);
return res;
}
constexpr FPS inv() const {
return inv((int)this->size());
}
constexpr FPS inv_sparse(vector<pair<int, mint>> f){
mint f0 = f[0].second ;
mint g0 = mint(1)/f0 ;
f.erase(f.begin()) ;
const int N = (*this).size() ;
FPS g(N) ;
g[0] = g0 ;
for(int n=1; n < N ;n++){
mint rhs = 0 ;
for(auto &&[i, fi]:f){
if(i <= n) rhs -= fi * g[n-i] ;
}
g[n] = rhs * g0 ;
}
return *this = g ;
}
// log(f) = \int f'/f dx, f[0] must be 1
constexpr FPS log(int deg) const {
assert((*this)[0]==1);
FPS res = (diff() * inv(deg)).integral();
res.resize(deg);
return res;
}
constexpr FPS log() const {
return log((int)this->size());
}
constexpr FPS log_sparse(){
FPS f=(*this);
assert(f[0]==1);
const int n=f.size();
const int mod=f[0].get_mod();
vector<pair<int,mint>>f_sub;
for(int i=1;i<n;i++)if(f[i]!=0)f_sub.emplace_back(i,f[i]);
FPS g(n),inv(n);
for(auto[a,c]:f_sub)if(a)g[a]=c;
if(n>1)inv[1]=1;
for(int i=1;i<n;i++){
mint sum=0;
for(auto[a,c]:f_sub)if(a&&a<=i)sum+=c*(i-a)*g[i-a];
if(i>=2){
inv[i]=-inv[mod%i]*(mod/i);
sum*=inv[i];
}
g[i]-=sum;
}
return g;
}
// exp(f), f[0] must be 0
constexpr FPS exp(int deg) const {
assert((*this)[0] == 0);
FPS res(1, 1);
for (int i = 1; i < deg; i <<= 1) {
res = res * (pre(i<<1) - res.log(i<<1) + 1).pre(i<<1);
}
res.resize(deg);
return res;
}
constexpr FPS exp() const {
return exp((int)this->size());
}
static constexpr FPS exp_sparse(const vector<pair<int,mint>>& f_sub, int n){
int mod=f_sub[0].second.get_mod();
FPS g(n),inv(n);
g[0]=1;
if(n>1)inv[1]=1;
for(int i=1;i<n;i++){
for(auto[a,c]:f_sub)if(a<=i)g[i]+=c*g[i-a]*a;
if(i>=2){
inv[i]=-inv[mod%i]*(mod/i);
g[i]*=inv[i];
}
}
return g;
}
//pow(f) = exp(e * log f)
constexpr FPS pow(long long e, int deg) const {
if (e == 0) {
FPS res(deg, 0);
res[0] = 1;
return res;
}
long long i = 0;
while (i < (int)this->size() && (*this)[i] == 0) ++i;
if (i == (int)this->size() || i > (deg - 1) / e) return FPS(deg, 0);
mint k = (*this)[i];
FPS res = ((((*this) >> i) / k).log(deg) * e).exp(deg) * mint(k).pow(e) << (e * i);
res.resize(deg);
return res;
}
constexpr FPS pow(long long e) const {
return pow(e, (int)this->size());
}
constexpr FPS pow_sparse(const vector<pair<int,mint>>&f_sub, int n, long long k){
FPS g(n);
if(f_sub.empty()){
if(k==0)g[0]=1;
return g;
}
auto[a0,c0]=f_sub[0];
if(a0>0){
if(k>(n-1)/a0)return FPS(n);
vector<pair<int,mint>>g_sub=f_sub;
for(auto&[a,c]:g_sub)a-=a0;
auto h=pow_sparse(g_sub,n-k*a0,k);
for(int i=0;i<n-k*a0;i++)g[i+k*a0]=h[i];
return g;
}
const long long mod=(*this)[0].get_mod();
FPS inv(n);
inv[1]=1;
for(int i=2;i<n;i++)inv[i]=-inv[mod%i]*(mod/i);
g[0]=c0.pow(k);
mint cef=c0.inv();
for(int i=1;i<n;i++){
for(auto[a,c]:f_sub)if(a&&i-a>=0)g[i]+=c*g[i-a]*(mint(k+1)*a-i);
g[i]*=inv[i]*cef;
}
return g;
}
// sqrt(f), f[0] must be 1
constexpr FPS sqrt_base(int deg) const {
assert((*this)[0] == 1);
mint inv2 = mint(1) / 2;
FPS res(1, 1);
for (int i = 1; i < deg; i <<= 1) {
res = (res + pre(i << 1) * res.inv(i << 1)).pre(i << 1);
for (mint &x : res) x *= inv2;
}
res.resize(deg);
return res;
}
constexpr FPS sqrt_sparse() {
FPS f = *this ;
int n = f.size();
int di = n;
for (int i = n - 1; i >= 0; --i)if (f[i] != 0) di = i;
if (di == n) return f;
if (di & 1) return {};
mint y = f[di];
mint x = y.sqrt();
if (x * x != y) return {};
mint c = mint(1) / y;
FPS g(n - di);
for (int i = 0; i < n - di; ++i) g[i] = f[di + i] * c;
g = g.sqrt_base();
for (int i = 0; i < int(g.size()); ++i) g[i] *= x; g.resize(n);
for (int i = n - 1; i >= 0; --i) {
if (i >= di / 2) g[i] = g[i - di / 2];
else g[i] = 0;
}
return *this = g;
}
constexpr FPS sqrt_base() const {
return sqrt_base((int)this->size());
}
static constexpr FPS subset_sum(const vector<int>& s,int t) {
vector<mint> a(t+1);
for(auto g:s)if(g<=t)a[g]+=1;
vector<mint> invmemo(t+1);
reps(i,t){
mint k=i;
invmemo[i]=k.inv();
}
FPS g(t+1);
reps(i,t){
for(int j=1;j*i<=t;j++){
if(j&1)g[j*i]+=invmemo[j]*a[i];
else g[j*i]-=invmemo[j]*a[i];
}
}
g=g.exp();
return g;
}
static constexpr FPS subset_sum2(const vector<int>& s,int t) {
vector<mint> a(t+1);
for(auto g:s)if(g<=t)a[g]+=1;
vector<mint> invmemo(t+1);
reps(i,t){
mint k=i;
invmemo[i]=k.inv();
}
FPS g(t+1);
reps(i,t){
for(int j=1;j*i<=t;j++){
g[j*i]+=invmemo[j]*a[i];
}
}
g=g.exp();
return g;
}
// ntt 系ラッパー (multipoint_eval の高速版などで使う)
constexpr void ntt() {
NTT::butterfly(*this);
}
constexpr void intt() {
NTT::butterfly_inv(*this);
mint iz = mint(this->size()).inv();
for (auto& x : *this) x *= iz;
}
constexpr void ntt_doubling() {
NTT::ntt_doubling(*this);
}
static constexpr mint ntt_pr() {
return mint(NTT::calc_primitive_root(mint::get_mod()));
}
// friend operators
friend constexpr FPS diff(const FPS &f) { return f.diff(); }
friend constexpr FPS integral(const FPS &f) { return f.integral(); }
friend constexpr FPS inv(const FPS &f, int deg) { return f.inv(deg); }
friend constexpr FPS inv(const FPS &f) { return f.inv((int)f.size()); }
friend constexpr FPS log(const FPS &f, int deg) { return f.log(deg); }
friend constexpr FPS log(const FPS &f) { return f.log((int)f.size()); }
friend constexpr FPS exp(const FPS &f, int deg) { return f.exp(deg); }
friend constexpr FPS exp(const FPS &f) { return f.exp((int)f.size()); }
friend constexpr FPS pow(const FPS &f, long long e, int deg) { return f.pow(e, deg); }
friend constexpr FPS pow(const FPS &f, long long e) { return f.pow(e, (int)f.size()); }
friend constexpr FPS sqrt_base(const FPS &f, int deg) { return f.sqrt_base(deg); }
friend constexpr FPS sqrt_base(const FPS &f) { return f.sqrt_base((int)f.size()); }
};
// 有理式の和を計算する。分割統治 O(Nlog^2N)。N は次数の和。
template <typename mint>
ob2<FPS<mint>, FPS<mint>> sum_of_rationals(vector<ob2<FPS<mint>, FPS<mint>>> dat) {
if (dat.size() == 0) {
FPS<mint> f = {0}, g = {1};
return {f, g};
}
using P = ob2<FPS<mint>, FPS<mint>>;
auto add = [&](P& a, P& b) -> P {
int na = a.a.size() - 1, da = a.b.size() - 1;
int nb = b.a.size() - 1, db = b.b.size() - 1;
int n = max(na + db, da + nb);
FPS<mint> num(n + 1);
{
auto f = a.a*b.b;
repi(i, f.size()) num[i] += f[i];
}
{
auto f = a.b*b.a;
repi(i, f.size()) num[i] += f[i];
}
auto den = a.b*b.b;
return {num, den};
};
while (dat.size() > 1) {
int n = dat.size();
for(int i=1;i<n;i+=2) { dat[i - 1] = add(dat[i - 1], dat[i]); }
repi(i, ceil(n, 2)) dat[i] = dat[2 * i];
dat.resize(ceil(n, 2));
}
return dat[0];
}
// sum wt[i]/(1-A[i]x)
template <typename mint>
ob2<FPS<mint>, FPS<mint>> sum_of_rationals_1(FPS<mint> A, FPS<mint> wt) {
using poly = FPS<mint>;
if ((mint::get_mod()-1)%(1<<23)!=0) {
vector<ob2<poly, poly>> rationals;
repi(i, A.size()) rationals.emplace_back(poly({wt[i]}), poly({mint(1), -A[i]}));
return sum_of_rationals(rationals);
}
int n = 1;
while (n < A.size()) n *= 2;
int k = __builtin_ctz(n);
FPS<mint> F(n), G(n);
FPS<mint> nxt_F(n), nxt_G(n);
repi(i, A.size()) F[i] = -A[i], G[i] = wt[i];
int D = 6;
repi(d, k) {
int b = 1 << d;
if (d < D) {
fill(all(nxt_F), mint(0)), fill(all(nxt_G), mint(0));
for (int L = 0; L < n; L += 2 * b) {
repi(i, b) repi(j, b) nxt_F[L + i + j] += F[L + i] * F[L + b + j];
repi(i, b) repi(j, b) nxt_G[L + i + j] += F[L + i] * G[L + b + j];
repi(i, b) repi(j, b) nxt_G[L + i + j] += F[L + b + i] * G[L + j];
repi(i, b) nxt_F[L + b + i] += F[L + i] + F[L + b + i];
repi(i, b) nxt_G[L + b + i] += G[L + i] + G[L + b + i];
}
}
else if (d == D) {
for (int L = 0; L < n; L += 2 * b) {
poly f1 = {F.begin() + L, F.begin() + L + b};
poly f2 = {F.begin() + L + b, F.begin() + L + 2 * b};
poly g1 = {G.begin() + L, G.begin() + L + b};
poly g2 = {G.begin() + L + b, G.begin() + L + 2 * b};
f1.resize(2 * b), f2.resize(2 * b), g1.resize(2 * b), g2.resize(2 * b);
NTT::ntt(f1), NTT::ntt(f2), NTT::ntt(g1), NTT::ntt(g2);
repi(i, b) f1[i] += 1, f2[i] += 1;
for(int i=b;i<2*b;i++) f1[i] -= 1, f2[i] -= 1;
repi(i, 2 * b) nxt_F[L + i] = f1[i] * f2[i] - 1;
repi(i, 2 * b) nxt_G[L + i] = g1[i] * f2[i] + g2[i] * f1[i];
}
}
else {
for (int L = 0; L < n; L += 2 * b) {
poly f1 = {F.begin() + L, F.begin() + L + b};
poly f2 = {F.begin() + L + b, F.begin() + L + 2 * b};
poly g1 = {G.begin() + L, G.begin() + L + b};
poly g2 = {G.begin() + L + b, G.begin() + L + 2 * b};
NTT::ntt_doubling(f1), NTT::ntt_doubling(f2), NTT::ntt_doubling(g1), NTT::ntt_doubling(g2);
repi(i, b) f1[i] += 1, f2[i] += 1;
for(int i=b;i<2*b;i++) f1[i] -= 1, f2[i] -= 1;
repi(i, 2 * b) nxt_F[L + i] = f1[i] * f2[i] - 1;
repi(i, 2 * b) nxt_G[L + i] = g1[i] * f2[i] + g2[i] * f1[i];
}
}
swap(F, nxt_F), swap(G, nxt_G);
}
if (k - 1 >= D) NTT::intt(F), NTT::intt(G);
F.emplace_back(1);
reverse(all(F)), reverse(all(G));
F.resize(A.size() + 1);
G.resize(A.size());
return {G, F};
}
using mint=Fp<moda>;
int main(){
cin.tie(nullptr);
int n;in(n);int m;in(m);FPS<mint> a(n);in(a);FPS<mint> b(n,mint(1));auto v=sum_of_rationals_1<mint>(a,b);
FPS<mint> f=(v.a*(v.b.inv(m+1)));f.resize(m+1);
f.erase(f.begin());
out(f);
}