結果
| 問題 | No.502 階乗を計算するだけ |
| コンテスト | |
| ユーザー |
 tonegawa
|
| 提出日時 | 2023-10-11 22:03:12 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.89.0) |
| 結果 |
AC
|
| 実行時間 | 785 ms / 1,000 ms |
| コード長 | 65,254 bytes |
| 記録 | |
| コンパイル時間 | 3,451 ms |
| コンパイル使用メモリ | 196,196 KB |
| 最終ジャッジ日時 | 2025-02-17 06:45:18 |
|
ジャッジサーバーID （参考情報） |
judge2 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 52 |
ソースコード
#line 1 ".lib/template.hpp"
#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <tuple>
#include <stack>
#include <queue>
#include <deque>
#include <algorithm>
#include <set>
#include <map>
#include <unordered_set>
#include <unordered_map>
#include <bitset>
#include <cmath>
#include <functional>
#include <cassert>
#include <climits>
#include <iomanip>
#include <numeric>
#include <memory>
#include <random>
#include <thread>
#include <chrono>
#define allof(obj) (obj).begin(), (obj).end()
#define range(i, l, r) for(int i=l;i<r;i++)
#define unique_elem(obj) obj.erase(std::unique(allof(obj)), obj.end())
#define bit_subset(i, S) for(int i=S, zero_cnt=0;(zero_cnt+=i==S)<2;i=(i-1)&S)
#define bit_kpop(i, n, k) for(int i=(1<<k)-1,x_bit,y_bit;i<(1<<n);x_bit=(i&-i),y_bit=i+x_bit,i=(!i?(1<<n):((i&~y_bit)/x_bit>>1)|y_bit))
#define bit_kth(i, k) ((i >> k)&1)
#define bit_highest(i) (i?63-__builtin_clzll(i):-1)
#define bit_lowest(i) (i?__builtin_ctzll(i):-1)
#define sleepms(t) std::this_thread::sleep_for(std::chrono::milliseconds(t))
using ll = long long;
using ld = long double;
using ul = uint64_t;
using pi = std::pair<int, int>;
using pl = std::pair<ll, ll>;
using namespace std;
template<typename F, typename S>
std::ostream &operator<<(std::ostream &dest, const std::pair<F, S> &p){
dest << p.first << ' ' << p.second;
return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::vector<std::vector<T>> &v){
int sz = v.size();
if(sz==0) return dest;
for(int i=0;i<sz;i++){
int m = v[i].size();
for(int j=0;j<m;j++) dest << v[i][j] << (i!=sz-1&&j==m-1?'\n':' ');
}
return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::vector<T> &v){
int sz = v.size();
if(sz==0) return dest;
for(int i=0;i<sz-1;i++) dest << v[i] << ' ';
dest << v[sz-1];
return dest;
}
template<typename T, size_t sz>
std::ostream &operator<<(std::ostream &dest, const std::array<T, sz> &v){
if(sz==0) return dest;
for(int i=0;i<sz-1;i++) dest << v[i] << ' ';
dest << v[sz-1];
return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::set<T> &v){
for(auto itr=v.begin();itr!=v.end();){
dest << *itr;
itr++;
if(itr!=v.end()) dest << ' ';
}
return dest;
}
template<typename T, typename E>
std::ostream &operator<<(std::ostream &dest, const std::map<T, E> &v){
for(auto itr=v.begin();itr!=v.end();){
dest << '(' << itr->first << ", " << itr->second << ')';
itr++;
if(itr!=v.end()) dest << '\n';
}
return dest;
}
std::ostream &operator<<(std::ostream &dest, __int128_t value) {
std::ostream::sentry s(dest);
if (s) {
__uint128_t tmp = value < 0 ? -value : value;
char buffer[128];
char *d = std::end(buffer);
do {
--d;
*d = "0123456789"[tmp % 10];
tmp /= 10;
} while (tmp != 0);
if (value < 0) {
--d;
*d = '-';
}
int len = std::end(buffer) - d;
if (dest.rdbuf()->sputn(d, len) != len) {
dest.setstate(std::ios_base::badbit);
}
}
return dest;
}
template<typename T>
vector<T> make_vec(size_t sz, T val){return std::vector<T>(sz, val);}
template<typename T, typename... Tail>
auto make_vec(size_t sz, Tail ...tail){
return std::vector<decltype(make_vec<T>(tail...))>(sz, make_vec<T>(tail...));
}
template<typename T>
vector<T> read_vec(size_t sz){
std::vector<T> v(sz);
for(int i=0;i<(int)sz;i++) std::cin >> v[i];
return v;
}
template<typename T, typename... Tail>
auto read_vec(size_t sz, Tail ...tail){
auto v = std::vector<decltype(read_vec<T>(tail...))>(sz);
for(int i=0;i<(int)sz;i++) v[i] = read_vec<T>(tail...);
return v;
}
void io_init(){
std::cin.tie(nullptr);
std::ios::sync_with_stdio(false);
}
#line 1 ".lib/math/combinatorics.hpp"
#line 1 ".lib/math/mod.hpp"
#line 6 ".lib/math/mod.hpp"
#include <type_traits>
#line 8 ".lib/math/mod.hpp"
#include <ostream>
#line 1 ".lib/math/minior/mod_base.hpp"
#line 4 ".lib/math/minior/mod_base.hpp"
// @param m `1 <= m`
constexpr long long safe_mod(long long x, long long m){
x %= m;
if (x < 0) x += m;
return x;
}
struct barrett{
unsigned int _m;
unsigned long long im;
explicit barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1){}
unsigned int umod()const{return _m;}
unsigned int mul(unsigned int a, unsigned int b)const{
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x = (unsigned long long)(((unsigned __int128)(z) * im) >> 64);
#endif
unsigned long long y = x * _m;
return (unsigned int)(z - y + (z < y ? _m : 0));
}
};
// @param n `0 <= n`
// @param m `1 <= m`
constexpr long long pow_mod_constexpr(long long x, long long n, int m){
if(m == 1) return 0;
unsigned int _m = (unsigned int)(m);
unsigned long long r = 1;
unsigned long long y = safe_mod(x, m);
while(n){
if (n & 1) r = (r * y) % _m;
y = (y * y) % _m;
n >>= 1;
}
return r;
}
constexpr bool is_prime_constexpr(int n) {
if (n <= 1) return false;
if (n == 2 || n == 7 || n == 61) return true;
if (n % 2 == 0) return false;
long long d = n - 1;
while (d % 2 == 0) d /= 2;
constexpr long long bases[3] = {2, 7, 61};
for(long long a : bases){
long long t = d;
long long y = pow_mod_constexpr(a, t, n);
while(t != n - 1 && y != 1 && y != n - 1){
y = y * y % n;
t <<= 1;
}
if(y != n - 1 && t % 2 == 0){
return false;
}
}
return true;
}
template<int n>
constexpr bool is_prime = is_prime_constexpr(n);
constexpr int primitive_root_constexpr(int m){
if(m == 2) return 1;
if(m == 167772161) return 3;
if(m == 469762049) return 3;
if(m == 754974721) return 11;
if(m == 998244353) return 3;
int divs[20] = {};
divs[0] = 2;
int cnt = 1;
int x = (m - 1) / 2;
while (x % 2 == 0) x /= 2;
for(int i = 3; (long long)(i)*i <= x; i += 2){
if(x % i == 0){
divs[cnt++] = i;
while(x % i == 0){
x /= i;
}
}
}
if(x > 1) divs[cnt++] = x;
for(int g = 2;; g++){
bool ok = true;
for(int i = 0; i < cnt; i++){
if(pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1){
ok = false;
break;
}
}
if(ok)return g;
}
}
template <int m>
constexpr int primitive_root = primitive_root_constexpr(m);
int ceil_pow2(int n){
int x = 0;
while ((1U << x) < (unsigned int)(n)) x++;
return x;
}
int bsf(unsigned int n){
return __builtin_ctz(n);
}
// @param b `1 <= b`
// @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
constexpr std::pair<long long, long long> inv_gcd(long long a, long long b){
a = safe_mod(a, b);
if(a == 0) return {b, 0};
long long s = b, t = a;
long long m0 = 0, m1 = 1;
while (t){
long long u = s / t;
s -= t * u;
m0 -= m1 * u;
auto tmp = s;
s = t;
t = tmp;
tmp = m0;
m0 = m1;
m1 = tmp;
}
if(m0 < 0) m0 += b / s;
return {s, m0};
}
#line 13 ".lib/math/mod.hpp"
template<int m>
long long modpow(long long a, long long b){
assert(0 <= b);
assert(0 < m);
a = safe_mod(a, m);
long long ret = 1;
while(b){
if(b & 1) ret = (ret * a) % m;
a = (a * a) % m;
b >>= 1;
}
return ret;
}
// @param 0 <= b, 0 < m
long long modpow(long long a, long long b, int m){
assert(0 <= b);
assert(0 < m);
a = safe_mod(a, m);
long long ret = 1;
while(b){
if(b & 1) ret = (ret * a) % m;
a = (a * a) % m;
b >>= 1;
}
return ret;
}
struct modint_base {};
struct static_modint_base : modint_base {};
template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : static_modint_base{
using mint = static_modint;
public:
static constexpr int mod(){return m;}
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
static_modint(): _v(0){}
template <class T>
static_modint(T v){
long long x = v % (long long)umod();
if (x < 0) x += umod();
_v = x;
}
unsigned int val()const{return _v;}
mint& operator++(){
_v++;
if (_v == umod()) _v = 0;
return *this;
}
mint& operator--(){
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int){
mint result = *this;
++*this;
return result;
}
mint operator--(int){
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs){
_v += rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs){
_v -= rhs._v;
if (_v >= umod()) _v += umod();
return *this;
}
mint& operator*=(const mint& rhs){
unsigned long long z = _v;
z *= rhs._v;
_v = (unsigned int)(z % umod());
return *this;
}
mint& operator/=(const mint& rhs){return *this = *this * rhs.inv();}
mint operator+()const{return *this;}
mint operator-()const{return mint() - *this;}
mint pow(long long n)const{
assert(0 <= n);
mint x = *this, r = 1;
while(n){
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv()const{
if(prime){
assert(_v);
return pow(umod() - 2);
}else{
auto eg = inv_gcd(_v, m);
assert(eg.first == 1);
return eg.second;
}
}
friend mint operator+(const mint& lhs, const mint& rhs){return mint(lhs) += rhs;}
friend mint operator-(const mint& lhs, const mint& rhs){return mint(lhs) -= rhs;}
friend mint operator*(const mint& lhs, const mint& rhs){return mint(lhs) *= rhs;}
friend mint operator/(const mint& lhs, const mint& rhs){return mint(lhs) /= rhs;}
friend bool operator==(const mint& lhs, const mint& rhs){return lhs._v == rhs._v;}
friend bool operator!=(const mint& lhs, const mint& rhs){return lhs._v != rhs._v;}
private:
unsigned int _v;
static constexpr unsigned int umod(){return m;}
static constexpr bool prime = is_prime<m>;
};
template<int id>
struct dynamic_modint : modint_base{
using mint = dynamic_modint;
public:
static int mod(){return (int)(bt.umod());}
static void set_mod(int m){
assert(1 <= m);
bt = barrett(m);
}
static mint raw(int v){
mint x;
x._v = v;
return x;
}
dynamic_modint(): _v(0){}
template <class T>
dynamic_modint(T v){
long long x = v % (long long)(mod());
if (x < 0) x += mod();
_v = x;
}
unsigned int val()const{return _v;}
mint& operator++(){
_v++;
if(_v == umod()) _v = 0;
return *this;
}
mint& operator--(){
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int){
mint result = *this;
++*this;
return result;
}
mint operator--(int){
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs){
_v += rhs._v;
if(_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs){
_v += mod() - rhs._v;
if(_v >= umod()) _v -= umod();
return *this;
}
mint& operator*=(const mint& rhs){
_v = bt.mul(_v, rhs._v);
return *this;
}
mint& operator/=(const mint& rhs){return *this = *this * rhs.inv();}
mint operator+()const{return *this;}
mint operator-()const{return mint() - *this;}
mint pow(long long n)const{
assert(0 <= n);
mint x = *this, r = 1;
while(n){
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv()const{
auto eg = inv_gcd(_v, mod());
assert(eg.first == 1);
return eg.second;
}
friend mint operator+(const mint& lhs, const mint& rhs){return mint(lhs) += rhs;}
friend mint operator-(const mint& lhs, const mint& rhs){return mint(lhs) -= rhs;}
friend mint operator*(const mint& lhs, const mint& rhs){return mint(lhs) *= rhs;}
friend mint operator/(const mint& lhs, const mint& rhs){return mint(lhs) /= rhs;}
friend bool operator==(const mint& lhs, const mint& rhs){return lhs._v == rhs._v;}
friend bool operator!=(const mint& lhs, const mint& rhs){return lhs._v != rhs._v;}
private:
unsigned int _v;
static barrett bt;
static unsigned int umod(){return bt.umod();}
};
template <int id>
barrett dynamic_modint<id>::bt(998244353);
using modint = dynamic_modint<-1>;
using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
template <class T>
using is_modint = std::is_base_of<modint_base, T>;
template <class T>
using is_modint_t = std::enable_if_t<is_modint<T>::value>;
template <class T>
using is_static_modint = std::is_base_of<static_modint_base, T>;
template <class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template <class> struct is_dynamic_modint : public std::false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template <class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
template<int m>
std::ostream &operator<<(std::ostream &dest, const static_modint<m> &a){
dest << a.val();
return dest;
}
template<int id>
std::ostream &operator<<(std::ostream &dest, const dynamic_modint<id> &a){
dest << a.val();
return dest;
}
// 0 <= n < m <= int_max
// 前処理 O(n + log(m))
// 各種計算 O(1)
// 変数 <= n
template<typename mint, is_modint<mint>* = nullptr>
struct modcomb{
private:
int n;
std::vector<mint> f, i, fi;
void init(int _n){
assert(0 <= _n && _n < mint::mod());
if(_n < f.size()) return;
n = _n;
f.resize(n + 1), i.resize(n + 1), fi.resize(n + 1);
f[0] = fi[0] = mint(1);
if(n) f[1] = fi[1] = i[1] = mint(1);
for(int j = 2; j <= n; j++) f[j] = f[j - 1] * j;
fi[n] = f[n].inv();
for(int j = n; j >= 2; j--){
fi[j - 1] = fi[j] * j;
i[j] = f[j - 1] * fi[j];
}
}
public:
modcomb(): n(-1){}
modcomb(int _n){
init(_n);
}
void recalc(int _n){
init(std::min(mint::mod() - 1, 1 << ceil_pow2(_n)));
}
mint comb(int a, int b){
if((a < 0) || (b < 0) || (a < b)) return 0;
return f[a] * fi[a - b] * fi[b];
}
mint perm(int a, int b){
if((a < 0) || (b < 0) || (a < b)) return 0;
return f[a] * fi[a - b];
}
mint fac(int x){
assert(0 <= x && x <= n);
return f[x];
}
mint inv(int x){
assert(0 < x && x <= n);
return i[x];
}
mint finv(int x){
assert(0 <= x && x <= n);
return fi[x];
}
};
template<typename mint, is_modint<mint>* = nullptr>
struct modpow_table{
std::vector<mint> v;
// x^maxkまで計算できる
modpow_table(){}
void init(int x, int maxk){
v.resize(maxk + 1);
v[0] = 1;
for(int i = 1; i <= maxk; i++) v[i] = v[i - 1] * x;
}
mint pow(int k){
assert(0 <= k && k < v.size());
return v[k];
}
};
#line 1 ".lib/math/fps_extra.hpp"
#line 1 ".lib/math/fps.hpp"
#line 1 ".lib/math/convolution.hpp"
#line 5 ".lib/math/convolution.hpp"
template<typename mint>
void butterfly(std::vector<mint>& a){
static constexpr int g = primitive_root<mint::mod()>;
int n = int(a.size());
int h = ceil_pow2(n);
static bool first = true;
static mint sum_e[30]; // sum_e[i] = ies[0] * ... * ies[i - 1] * es[i]
if(first){
first = false;
mint es[30], ies[30]; // es[i]^(2^(2+i)) == 1
int cnt2 = bsf(mint::mod() - 1);
mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
for(int i = cnt2; i >= 2; i--){
// e^(2^i) == 1
es[i - 2] = e;
ies[i - 2] = ie;
e *= e;
ie *= ie;
}
mint now = 1;
for(int i = 0; i <= cnt2 - 2; i++){
sum_e[i] = es[i] * now;
now *= ies[i];
}
}
for(int ph = 1; ph <= h; ph++){
int w = 1 << (ph - 1), p = 1 << (h - ph);
mint now = 1;
for(int s = 0; s < w; s++){
int offset = s << (h - ph + 1);
for(int i = 0; i < p; i++){
auto l = a[i + offset];
auto r = a[i + offset + p] * now;
a[i + offset] = l + r;
a[i + offset + p] = l - r;
}
now *= sum_e[bsf(~(unsigned int)(s))];
}
}
}
template<typename mint>
void butterfly_inv(std::vector<mint>& a){
static constexpr int g = primitive_root<mint::mod()>;
int n = int(a.size());
int h = ceil_pow2(n);
static bool first = true;
static mint sum_ie[30]; // sum_ie[i] = es[0] * ... * es[i - 1] * ies[i]
if(first){
first = false;
mint es[30], ies[30]; // es[i]^(2^(2+i)) == 1
int cnt2 = bsf(mint::mod() - 1);
mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
for(int i = cnt2; i >= 2; i--){
// e^(2^i) == 1
es[i - 2] = e;
ies[i - 2] = ie;
e *= e;
ie *= ie;
}
mint now = 1;
for(int i = 0; i <= cnt2 - 2; i++){
sum_ie[i] = ies[i] * now;
now *= es[i];
}
}
for(int ph = h; ph >= 1; ph--){
int w = 1 << (ph - 1), p = 1 << (h - ph);
mint inow = 1;
for(int s = 0; s < w; s++){
int offset = s << (h - ph + 1);
for(int i = 0; i < p; i++){
auto l = a[i + offset];
auto r = a[i + offset + p];
a[i + offset] = l + r;
a[i + offset + p] = (unsigned long long)(mint::mod() + l.val() - r.val()) * inow.val();
}
inow *= sum_ie[bsf(~(unsigned int)(s))];
}
}
}
template<typename mint>
std::vector<mint> convolution_mod(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){
if(n < m){
std::swap(n, m);
std::swap(a, b);
}
std::vector<mint> ans(n + m - 1);
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
ans[i + j] += a[i] * b[j];
}
}
return ans;
}
int z = 1 << ceil_pow2(n + m - 1);
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(int i = 0; i < n + m - 1; i++) a[i] *= iz;
return a;
}
template <unsigned int mod = 998244353, typename T>
std::vector<T> convolution_mod(const std::vector<T>& a, const std::vector<T>& b) {
int n = int(a.size()), m = int(b.size());
if(!n || !m) return {};
using mint = static_modint<mod>;
std::vector<mint> a2(n), b2(m);
for(int i = 0; i < n; i++) a2[i] = mint(a[i]);
for(int i = 0; i < m; i++) b2[i] = mint(b[i]);
auto c2 = convolution_mod<mint>(move(a2), move(b2));
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++) c[i] = c2[i].val();
return c;
}
template<typename mint>
std::vector<mint> square_mod(std::vector<mint> a){
int n = int(a.size());
if(!n) return {};
if(n <= 60){
std::vector<mint> ans(2 * n - 1);
for(int i = 0; i < n; i++){
for(int j = i + 1; j < n; j++){
ans[i + j] += 2 * a[i] * a[j];
}
ans[2 * i] += a[i] * a[i];
}
return ans;
}
int z = 1 << ceil_pow2(2 * n - 1);
a.resize(z);
butterfly(a);
for(int i = 0; i < z; i++) a[i] *= a[i];
butterfly_inv(a);
a.resize(2 * n - 1);
mint iz = mint(z).inv();
for (int i = 0; i < 2 * n - 1; i++) a[i] *= iz;
return a;
}
template <unsigned int mod = 998244353, class T>
std::vector<T> square_mod(const std::vector<T> &a){
int n = int(a.size());
if(!n) return {};
using mint = static_modint<mod>;
std::vector<mint> a2(n);
for(int i = 0; i < n; i++) a2[i] = mint(a[i]);
auto c2 = square_mod<mint>(move(a2));
std::vector<T> c(2 * n - 1);
for(int i = 0; i < 2 * n - 1; i++) c[i] = c2[i].val();
return c;
}
// ntt-friendlyでないmodで畳み込みを行う
template<typename mint>
std::vector<mint> convolution_int_mod(const std::vector<mint>& a, const std::vector<mint>& b){
if(mint::mod() == 998244353) return convolution_mod<mint>(a, b);
int n = int(a.size()), m = int(b.size());
if(!n || !m) return {};
if(std::min(n, m) <= 60){
std::vector<mint> ans(n + m - 1, 0);
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
ans[i + j] += a[i] * b[j];
}
}
return ans;
}
static constexpr long long MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
static constexpr long long M1M2 = MOD1 * MOD2, ix = inv_gcd(MOD1, MOD2).second, i3 = inv_gcd(MOD1 * MOD2, MOD3).second;
std::vector<long long> a2(n), b2(m);
for(int i = 0; i < n; i++) a2[i] = a[i].val();
for(int i = 0; i < m; i++) b2[i] = b[i].val();
auto c1 = convolution_mod<MOD1>(a2, b2);
auto c2 = convolution_mod<MOD2>(a2, b2);
auto c3 = convolution_mod<MOD3>(a2, b2);
std::vector<mint> c(n + m - 1);
int M1M2m = M1M2 % mint::mod();
for(int i = 0; i < n + m - 1; i++){
long long v = (((long long)c2[i] - c1[i]) * ix) % MOD2;
if(v < 0) v += MOD2;
long long xxv = c1[i] + MOD1 * v;
v = ((c3[i] - (xxv % MOD3)) * i3) % MOD3;
if(v < 0) v += MOD3;
c[i] = mint(xxv + M1M2m * v);
}
return c;
}
// ntt-friendlyでないmodでa * aを計算する
template<typename mint>
std::vector<mint> square_int_mod(const std::vector<mint>& a){
if(mint::mod() == 998244353) return square_mod<mint>(a);
int n = int(a.size());
if(!n) return {};
if(n <= 60){
std::vector<mint> ans(2 * n - 1);
for(int i = 0; i < n; i++){
for(int j = i + 1; j < n; j++){
ans[i + j] += 2 * a[i] * a[j];
}
ans[2 * i] += a[i] * a[i];
}
return ans;
}
static constexpr long long MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
static constexpr long long M1M2 = MOD1 * MOD2, ix = inv_gcd(MOD1, MOD2).second, i3 = inv_gcd(MOD1 * MOD2, MOD3).second;
std::vector<long long> a2(n);
for(int i = 0; i < n; i++) a2[i] = a[i].val();
auto c1 = square_mod<MOD1>(a2);
auto c2 = square_mod<MOD2>(a2);
auto c3 = square_mod<MOD3>(a2);
std::vector<mint> c(2 * n - 1);
int M1M2m = M1M2 % mint::mod();
for(int i = 0; i < 2 * n - 1; i++){
long long v = (((long long)c2[i] - c1[i]) * ix) % MOD2;
if(v < 0) v += MOD2;
long long xxv = c1[i] + MOD1 * v;
v = ((c3[i] - (xxv % MOD3)) * i3) % MOD3;
if(v < 0) v += MOD3;
c[i] = mint(xxv + M1M2m * v);
}
return c;
}
// @param 答えがlong longに収まる
std::vector<long long> convolution_ll(const std::vector<long long>& a, const std::vector<long long>& b){
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
if(std::min(n, m) <= 60){
std::vector<long long> ans(n + m - 1, 0);
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
ans[i + j] += a[i] * b[j];
}
}
return ans;
}
// 2^24, 2^25, 2^26
static constexpr unsigned long long MOD1 = 754974721, MOD2 = 167772161, MOD3 = 469762049;
static constexpr unsigned long long M2M3 = MOD2 * MOD3, M1M3 = MOD1 * MOD3, M1M2 = MOD1 * MOD2, M1M2M3 = MOD1 * MOD2 * MOD3;
static constexpr unsigned long long i1 = inv_gcd(MOD2 * MOD3, MOD1).second, i2 = inv_gcd(MOD1 * MOD3, MOD2).second, i3 = inv_gcd(MOD1 * MOD2, MOD3).second;
auto c1 = convolution_mod<MOD1>(a, b);
auto c2 = convolution_mod<MOD2>(a, b);
auto c3 = convolution_mod<MOD3>(a, b);
std::vector<long long> c(n + m - 1);
for(int i = 0; i < n + m - 1; i++){
unsigned long long x = 0;
x += (c1[i] * i1) % MOD1 * M2M3;
x += (c2[i] * i2) % MOD2 * M1M3;
x += (c3[i] * i3) % MOD3 * M1M2;
long long diff = c1[i] - safe_mod((long long)(x), (long long)(MOD1));
if (diff < 0) diff += MOD1;
static constexpr unsigned long long offset[5] = {0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3};
x -= offset[diff % 5];
c[i] = x;
}
return c;
}
// @param 答えがlong longに収まる
std::vector<long long> square_ll(const std::vector<long long>& a){
int n = int(a.size());
if (!n) return {};
if(n <= 60){
std::vector<long long> ans(2 * n - 1, 0);
for(int i = 0; i < n; i++){
for(int j = i + 1; j < n; j++){
ans[i + j] += 2 * a[i] * a[j];
}
ans[2 * i] += a[i] * a[i];
}
return ans;
}
// 2^24, 2^25, 2^26
static constexpr unsigned long long MOD1 = 754974721, MOD2 = 167772161, MOD3 = 469762049;
static constexpr unsigned long long M2M3 = MOD2 * MOD3, M1M3 = MOD1 * MOD3, M1M2 = MOD1 * MOD2, M1M2M3 = MOD1 * MOD2 * MOD3;
static constexpr unsigned long long i1 = inv_gcd(MOD2 * MOD3, MOD1).second, i2 = inv_gcd(MOD1 * MOD3, MOD2).second, i3 = inv_gcd(MOD1 * MOD2, MOD3).second;
auto c1 = square_mod<MOD1>(a);
auto c2 = square_mod<MOD2>(a);
auto c3 = square_mod<MOD3>(a);
std::vector<long long> c(2 * n - 1);
for(int i = 0; i < 2 * n - 1; i++){
unsigned long long x = 0;
x += (c1[i] * i1) % MOD1 * M2M3;
x += (c2[i] * i2) % MOD2 * M1M3;
x += (c3[i] * i3) % MOD3 * M1M2;
long long diff = c1[i] - safe_mod((long long)(x), (long long)(MOD1));
if (diff < 0) diff += MOD1;
static constexpr unsigned long long offset[5] = {0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3};
x -= offset[diff % 5];
c[i] = x;
}
return c;
}
#line 7 ".lib/math/fps.hpp"
template<typename mint>
struct formal_power_series;
// 全ての要素の次数が同じだと仮定して総積を求める
// 全ての次数が1の場合 O(Nlog^2N)
template<typename mint>
formal_power_series<mint> multiply_constant_degree(std::vector<formal_power_series<mint>> v, int max_size = -1){
if(max_size == -1) max_size = std::numeric_limits<int>::max();
int n = v.size();
for(int d = 1; d < n; d <<= 1){
for(int j = 0; j < n; j += 2 * d){
if(j + d < n){
v[j] *= v[j + d];
if(v[j].size() > max_size) v[j].resize(max_size);
}
}
}
return v[0];
}
// 次数が低い順に計算
// 次数の総和をNとしてO(Nlog^2N)
template<typename mint>
formal_power_series<mint> multiply_lower_degree(std::vector<formal_power_series<mint>> v, int max_size = std::numeric_limits<int>::max()){
using p = std::pair<int, int>;
std::priority_queue<p, std::vector<p>, std::greater<p>> que;
int n = v.size();
if(n == 0) return {mint(1)};
for(int i = 0; i < n; i++){
if(v[i].size() > max_size) v[i].resize(max_size);
que.push({v[i].size(), i});
}
while(que.size() > 1){
int a = que.top().second;
que.pop();
int b = que.top().second;
que.pop();
v[a] *= v[b];
if(v[a].size() > max_size) v[a].resize(max_size);
que.push({v[a].size(), a});
}
return v[que.top().second];
}
template<typename mint>
struct formal_power_series : std::vector<mint>{
using std::vector<mint>::vector;
using fps = formal_power_series<mint>;
formal_power_series(const std::vector<mint> &v){
this->resize(v.size());
std::copy(v.begin(), v.end(), this->begin());
}
fps operator *= (const fps &vr){
return *this = convolution_int_mod<mint>(*this, vr);
}
fps operator /= (fps &vr){
return (*this) *= vr.inv();
}
fps operator += (const fps &vr){
int n = this->size(), m = vr.size();
if(n < m) this->resize(m);
for(int i = 0; i < m; i++) (*this)[i] += vr[i];
return *this;
}
fps operator -= (const fps &vr){
int n = this->size(), m = vr.size();
if(n < m) this->resize(m);
for(int i = 0; i < m; i++) (*this)[i] -= vr[i];
return *this;
}
fps operator %= (const fps &vr){
int n = this->size(), m = vr.size();
if(n < m) return *this;
n = n - m + 1;
fps r = ((rev().prefix(n) * vr.rev().inv(n)).prefix(n).rev(n)) * vr;
return (*this - r).prefix(m - 1);
}
// 係数が0でない最大の次数, deg(0) = -1とする
int deg()const{
int n = this->size();
for(int i = n - 1; i >= 0; i--) if((*this)[i].val()) return i;
return -1;
}
// 係数が0でない最大の次数, deg(0) = -1とする, 余分な0を消す
int deg_fix(){
int n = this->size();
for(int i = n - 1; i >= 0; i--){
if(this->back().val() != 0) return i;
else this->pop_back();
}
return -1;
}
// f(x) = q(x) * g(x) + r(x)となるような商q(x)と余りr(x)を求める
// f(x)がn - 1次, g(x)がm - 1次のとき, q(x)はn - m次, r(x)はm - 2次以下になる
std::pair<fps, fps> polynomial_div(fps g){
int n = this->size(), m = g.size();
if(deg_fix() < g.deg_fix()) return {fps{0}, *this};
n = n - m + 1;
fps q = (rev().prefix(n) * g.rev().inv(n)).prefix(n).rev(n);
return {q, ((*this) - q * g).prefix(m - 1)};
}
// どちらかの次数が-1の場合空のfpsを返す, O(nm)
static fps gcd_naive(const fps &f, const fps &g){
if(g.deg() == -1) return g;
if(f.deg() == -1) return f;
fps st[2] = {g, f};
bool flip = false;
while(st[!flip].deg_fix() != -1){
int n = st[flip].deg_fix(), m = st[!flip].deg_fix(), deg_diff = n - m;
if(n >= m){
mint c = st[flip][n] / st[!flip][m];
for(int i = deg_diff; i <= n; i++) st[flip][i] -= c * st[!flip][i - deg_diff];
n = st[flip].deg_fix(), m = st[!flip].deg_fix();
}
if(n < m) flip ^= 1;
}
return st[flip];
}
// f(x) * h(x)がmod g(x)でgcd(f(x), g(x))になるような{h(x), gcd(f(x), g(x))}を求める
// O(nm)
static std::pair<fps, fps> polynomial_inv_naive(const fps &f, const fps &g){
if(f.deg() == -1) return std::make_pair(g, fps{0}); // fが0
if(g.deg() == -1) return std::make_pair(fps{}, fps{}); // gで割れない
fps st[2] = {g, f};
fps ms[2] = {fps{0}, fps{1}};
bool flip = false;
while(st[!flip].deg_fix() != -1){
int n = st[flip].deg_fix(), m = st[!flip].deg_fix(), deg_diff = n - m;
if(n >= m){
mint c = st[flip][n] / st[!flip][m];
for(int i = deg_diff; i <= n; i++) st[flip][i] -= c * st[!flip][i - deg_diff];
int o = ms[!flip].deg_fix() + 1;
if(o + deg_diff > ms[flip].size()) ms[flip].resize(o + deg_diff, 0);
for(int i = 0; i < o; i++) ms[flip][i + deg_diff] -= ms[!flip][i] * c;
n = st[flip].deg_fix(), m = st[!flip].deg_fix();
}
if(n < m) flip ^= 1;
}
return {st[flip], ms[flip]};
}
fps operator += (const mint &vr){
int n = this->size();
if(n == 0) this->resize(1, 0);
(*this)[0] += vr;
return *this;
}
fps operator -= (const mint &vr){
int n = this->size();
if(n == 0) this->resize(1, 0);
(*this)[0] -= vr;
return *this;
}
fps operator *= (const mint &vr){
int n = this->size();
for(int i = 0; i < n; i++) (*this)[i] *= vr;
return *this;
}
fps operator /= (const mint &vr){
mint ir = vr.inv();
int n = this->size();
for(int i = 0; i < n; i++) (*this)[i] *= ir;
return *this;
}
fps operator + (const fps& vr)const{return fps(*this) += vr;}
fps operator - (const fps& vr)const{return fps(*this) -= vr;}
fps operator * (const fps& vr)const{return fps(*this) *= vr;}
fps operator / (const fps& vr)const{return fps(*this) /= vr;}
fps operator % (const fps& vr)const{return fps(*this) %= vr;}
fps operator + (const mint& vr)const{return fps(*this) += vr;}
fps operator - (const mint& vr)const{return fps(*this) -= vr;}
fps operator * (const mint& vr)const{return fps(*this) *= vr;}
fps operator / (const mint& vr)const{return fps(*this) /= vr;}
void print(int printsize = -1)const{
if(printsize == -1) printsize = (int)this->size();
printsize = std::min(printsize, (int)this->size());
if(printsize == 0){
std::cout << '\n';
return;
}
for(int i = 0; i < printsize; i++){
if(i == printsize - 1) std::cout << (*this)[i].val() << '\n';
else std::cout << (*this)[i].val() << ' ';
}
}
fps rev(int deg = -1)const{
fps res(*this);
if(deg != -1) res.resize(deg, 0);
std::reverse(res.begin(), res.end());
return res;
}
fps prefix(int deg){
int n = std::min((int)this->size(), deg);
return fps(this->begin(), this->begin() + n);
}
fps rshift(int s){
fps res(*this);
if(res.size() <= s) return {};
res.erase(res.begin(), res.begin() + s);
return res;
}
fps lshift(int s){
fps res(*this);
res.insert(res.begin(), s, mint(0));
return res;
}
fps inv_any_mod(int deg = -1){
assert((*this)[0].val());
int n = this->size();
if(deg == -1) deg = n;
fps res{(*this)[0].inv()};
for(int i = 1; i < deg; i <<= 1){
res = (res + res - (res * res * prefix(i << 1))).prefix(i << 1);
}
return res.prefix(deg);
}
fps inv(int deg = -1){
assert((*this)[0].val());
if(mint::mod() != 998244353) return inv_any_mod(deg);
int n = this->size();
if(deg == -1) deg = n;
fps res{(*this)[0].inv()};
for(int i = 1, z = i << 2; i < deg; i <<= 1, z <<= 1){
std::vector<mint> res_old = res, f_prefix = prefix(i << 1);
res_old.resize(z);
butterfly(res_old);
f_prefix.resize(z);
butterfly(f_prefix);
for(int j = 0; j < z; j++) res_old[j] *= res_old[j] * f_prefix[j];
butterfly_inv(res_old);
mint iz = mint(z).inv();
res_old.resize(z >> 1);
for(int j = 0; j < (z >> 1); j++) res_old[j] *= iz;
res = res + res - (fps)res_old;
}
return res.prefix(deg);
}
fps diff(){
int n = (int) this->size();
fps res(std::max(0, n - 1));
for(int i = 1; i < n; i++) res[i - 1] = (*this)[i] * i;
return res;
}
fps integral(){
static modcomb<mint> mcb;
int n = (int) this->size();
fps res(n + 1);
mcb.recalc(n);
res[0] = 0;
for(int i = 0; i < n; i++) res[i + 1] = (*this)[i] * mcb.inv(i + 1);
return res;
}
fps log(int deg = -1){
assert((*this)[0].val() == 1);
int n = (int)this->size();
if(deg == -1) deg = n;
return (this->diff() * this->inv(deg)).prefix(deg - 1).integral();
}
fps exp_any_mod(int deg = -1){
assert((*this)[0].val() == 0);
int n = this->size();
if(deg == -1) deg = n;
fps res{1};
for(int i = 1; i < deg; i <<= 1){
res = (res * (prefix(i << 1) + mint(1) - res.log(i << 1))).prefix(i << 1);
}
return res.prefix(deg);
}
fps exp(int deg = -1){
if(mint::mod() != 998244353) return exp_any_mod(deg);
assert((*this)[0].val() == 0);
int n = this->size();
if(deg == -1) deg = n;
fps res{1}, diff_res, diff_res_prev, inv_res, inv_res_old, inv_res_prev, log_res, log_res_prev;
for(int i = 1; i < deg; i <<= 1){
diff_res_prev = diff_res;
diff_res.resize(i - 1);
for(int j = std::max(1, i >> 1); j < i; j++) diff_res[j - 1] = res[j] * j;
if(i == 1){
inv_res = inv_res_old = {1};
}else{
fps i1 = inv_res_old;
i1.resize(i << 1);
butterfly(i1);
fps i2 = res.prefix(i);
i2.resize(i << 1);
butterfly(i2);
for(int j = 0; j < (i << 1); j++) i1[j] *= i1[j] * i2[j];
butterfly_inv(i1);
i1.resize(i);
mint iz = mint(i << 1).inv();
for(int j = 0; j < i; j++) i1[j] *= iz;
inv_res_old = inv_res_old + inv_res_old - i1;
i1 = inv_res_old;
i1.resize(i << 2);
butterfly(i1);
i2 = res.prefix(i << 1);
i2.resize(i << 2);
butterfly(i2);
for(int j = 0; j < (i << 2); j++) i1[j] *= i1[j] * i2[j];
butterfly_inv(i1);
i1.resize(i << 1);
iz = mint(i << 2).inv();
for(int j = 0; j < (i << 1); j++) i1[j] *= iz;
inv_res = inv_res_old + inv_res_old - i1;
}
log_res = (diff_res * inv_res).prefix((i << 1) - 1).integral();
res = (res * (prefix(i << 1) + mint(1) - log_res)).prefix(i << 1);
}
return res.prefix(deg);
}
fps pow(long long k, int deg = -1){
int n = (int) this->size();
if(deg == -1) deg = n;
if(!k){
fps res(deg, 0);
res[0] = 1;
return res;
}
for(long long i = 0; i < n; i++){
if((*this)[i] != 0) {
fps C = (*this) * (*this)[i].inv();
fps D(deg);
for(int j = i; j < n; j++) D[j - i] = C[j];
D = (D.log(deg) * mint(k)).exp(deg) * (*this)[i].pow(k);
fps E(deg);
if(i * k > deg) return E;
for(long long j = 0; j + __int128_t(i) * k < deg && j < D.size(); j++) E[j + i * k] = D[j];
return E;
}
}
return *this;
}
//初めて0以外が現れるのが奇数次数 or 初めに現れる数が平方根でない -> 解無し(空のfpsを返す)
fps sqrt(int deg = -1){
int n = (int)this->size();
if(deg == -1) deg = n;
if((*this)[0].val()==0){
for(int i = 1; i < n;i++){
if((*this)[i].val() != 0){
if(i & 1) return {};
if(deg - i / 2 <= 0) return fps(deg, 0);
fps res = rshift(i).sqrt(deg - i / 2);
if(res.empty()) return {};
res = res.lshift(i / 2);
if(res.size() < deg) res.resize(deg, 0);
return res;
}
}
return fps(deg, 0);
}
int x = modsqrt((*this)[0].val(), mint::mod());
if(x == -1) return {};
fps res({mint(x)});
mint i2 = mint(2).inv();
for(int i = 1; i < deg; i <<= 1) res = (res + prefix(i << 1) * res.inv(i << 1)) * i2;
return res;
}
// reference: http://www.eecs.harvard.edu/~htk/publication/1978-jacm-brent-kung.pdf
// N次多項式 f(x)にM次多項式 g(x)を合成した結果を求める
// deg(f) = deg(g) = deg(ans) = Nとして、
// f(x)をk := √N ブロックに平方分割すると f(x) = f_0(x) + f_1(x) x^k ... と約k項になる
// f((g(x))) = f_0(g(x)) + f_1(g(x))g(x)^k + ...
// 1. f_i(g(x))はk項とM項の合成なので、g(x)^i (0<=i<=k)を前処理しておくとO(Nk)
// 2. g(x)^kiを求める、かけるのは共に一回あたりO(NlogN)
// (1, 2)をkブロック分行うのでO(k × (Nk + NlogN)) = O(N^2 + N^1.5 * logN)
static fps composition(const fps &f, const fps &g, int deg = -1){
int n = f.size();
//if(deg == -1) deg = (n - 1) * (m - 1) + 1;
if(deg == -1) deg = n;
int k = (int)std::sqrt(n) + 1; // ブロック数
int d = (n + k - 1) / k; // 1ブロックあたりの次数
std::vector<fps> gpow(k + 1);
gpow[0] = {mint(1)};
for(int i = 1; i <= k; i++){
gpow[i] = gpow[i - 1] * g;
if(gpow[i].size() > deg) gpow[i].resize(deg);
}
std::vector<fps> fi(k, fps(deg, 0));
for(int i = 0; i < k; i++){
for(int j = 0; j < d; j++){
int idx = i * d + j;
if(idx >= n) break;
for(int t = 0; t < gpow[j].size(); t++) fi[i][t] += gpow[j][t] * f[idx];
}
}
fps res(deg, 0), gd = {1};
for(int i = 0; i < k; i++){
fi[i] *= gd;
int sz = std::min(deg, (int)fi[i].size());
for(int j = 0; j < sz; j++) res[j] += fi[i][j];
gd *= gpow[d];
if(gd.size() > deg) gd.resize(deg);
}
return res;
}
};
#line 5 ".lib/math/fps_extra.hpp"
template<typename mint>
std::vector<mint> multipoint_evaluation(const formal_power_series<mint> &f, const std::vector<mint> &v){
using fps = formal_power_series<mint>;
int m = v.size();
int N = 1;
while(N < m) N *= 2;
std::vector<fps> t(2 * N - 1, fps{1});
for(int i = 0; i < m; i++) t[N - 1 + i] = fps{-v[i], 1};
for(int i = N - 2; i >= 0; i--) t[i] = t[i * 2 + 1] * t[i * 2 + 2];
t[0] = f % t[0];
for(int i = 1; i < 2 * N - 1; i++){
t[i] = t[(i - 1) / 2] % t[i];
}
std::vector<mint> res(m);
for(int i = 0; i < m; i++){
res[i] = t[N - 1 + i][0];
}
return res;
}
template<typename mint, typename fps>
fps interpolation_divide_and_conquer(const std::vector<mint> &fx, const fps &f, const std::vector<fps> &tree, int k, int l, int r){
if(r - l == 1) return fps{fx[l] * f[0].inv()};
int mid = (l + r) / 2;
if(tree[k * 2 + 2].empty()) return interpolation_divide_and_conquer<mint, fps>(fx, f, tree, k * 2 + 1, l, mid);
fps left = interpolation_divide_and_conquer<mint, fps>(fx, f % tree[k * 2 + 1], tree, k * 2 + 1, l, mid);
fps right = interpolation_divide_and_conquer<mint, fps>(fx, f % tree[k * 2 + 2], tree, k * 2 + 2, mid, r);
return left * tree[k * 2 + 2] + right * tree[k * 2 + 1];
}
template<typename mint>
formal_power_series<mint> polynomial_interpolation(const std::vector<mint> &xi, const std::vector<mint> &fx){
using fps = formal_power_series<mint>;
int n = xi.size();
int N = 1;
while(N < n) N *= 2;
std::vector<fps> tree(2 * N - 1, fps{});
for(int i = 0; i < n; i++) tree[N - 1 + i] = fps{-xi[i], 1};
for(int i = N - 2; i >= 0; i--){
if(tree[i * 2 + 2].empty()) tree[i] = tree[i * 2 + 1];
else tree[i] = tree[i * 2 + 1] * tree[i * 2 + 2];
}
for(int i = 0; i < n; i++) tree[0][i] = tree[0][i + 1] * (i + 1);
tree[0].pop_back();
return interpolation_divide_and_conquer<mint, fps>(fx, tree[0], tree, 0, 0, N).prefix(n);
}
// a / bx + c を大量に足す
template<typename mint>
formal_power_series<mint> fractions_sum_binomial(const std::vector<std::tuple<mint, mint, mint>> &v, int max_size = -1){
using fps = formal_power_series<mint>;
if(max_size == -1) max_size = std::numeric_limits<int>::max();
int n = v.size();
std::vector<std::pair<fps, fps>> res(n);
for(int i = 0; i < n; i++){
auto [a, b, c] = v[i];
assert(b.val() || c.val());
res[i].first = fps{a};
res[i].second = fps{c, b};
}
for(int d = 1; d < n; d <<= 1){
for(int j = 0; j < n; j += 2 * d){
if(j + d < n){
v[j].first *= v[j + d].second;
v[j + d].first *= v[j].second;
v[j].second *= v[j + d].second;
v[j].first += v[j + d].first;
if(v[j].first.size() > max_size) v[j].first.resize(max_size);
if(v[j].second.size() > max_size) v[j].second.resize(max_size);
}
}
}
return res[0].first * res[0].second.inv();
}
/*
template<typename mint>
std::vector<mint> shiftof_sampling_points(mint x0, mint d, const std::vector<mint> &y, const std::vector<mint> &p){
using fps = formal_power_series<mint>;
static modcomb<mint> mcb;
int n = y.size();
mcb.recalc(n);
mint pi_xi_x0(1);
// xi == 0 or xi == xjで壊れる
// xi == 0は調べる
// 素数modだと仮定して, n >= mod or dがmodの倍数で壊れる
assert(n < mint::mod() && d.val()) ;
for(int i = 0; i < n; i++){
mint xi = x0 + d * i;
assert(xi.val());
pi_xi_x0 *= d * i;
}
std::vector<fps> tmp;
std::vector<std::tuple<mint, mint, mint>> tmp_frac;
for(int i = 0; i < n; i++){
mint xi = x0 + d * i;
tmp.push_back({-xi, 1});
}
fps x = multiply_constant_degree<mint>(tmp);
x /= pi_xi_x0;
for(int i = 0; i < n; i++){
mint iQ = mcb.finv(i) * mcb.finv(n - 1 - i);
mint xi = x0 + d * i;
if((n - i - 1) & 1) iQ = -iQ;
mint iC = y[i] * iQ;
tmp_frac.push_back({iC, 1, -xi});
}
x *= fractions_sum_binomial(tmp_frac);
x.print();
exit(0);
}
*/
/*
//f_i(x) := 1 ± a_i x^b_i を大量に掛ける
//O(N^2/th + thMlogM)
// N = M = 1e5 なら50が早かった
template<ll MOD>
StaticModFPS<MOD> multiply_binomial(vector<vector<ll>> v, ll m, ll th){
using fps = StaticModFPS<MOD>;
fps ret_big(m, 0);
vvl big;
vector<vector<fps>> small(th+1, vector<fps>());
vector<ll> inv((m/th)+2);
inv[0] = inv[1] = 1;
for(ll i=2;i<=(m/th)+1;i++){
ll u = MOD - (MOD/i);
inv[i] = (u*inv[MOD%i])%MOD;
}
for(auto e:v){
e[0] %= MOD;
if(e[0] < 0) e[0] += MOD;
if(e[1]>th) big.push_back(e);
else small[e[1]].push_back({1, e[0]});
}
for(int i=0;i<big.size();i++){
ll j = 1;
ll mul = big[i][0];
while(j * big[i][1] < m){
ll y;
if(j%2==0) y = ((MOD - abs(mul)) * inv[j])%MOD;
else y = (abs(mul) * inv[j])%MOD;
ret_big[j * big[i][1]] = (ret_big[j * big[i][1]] + y)%MOD;
mul = (mul * big[i][0])%MOD;
j++;
}
}
ret_big = ret_big.exp(m);
for(ll i=1;i<=th;i++){
if(small[i].size()==0) continue;
fps tmp = multiply_lower_degree<MOD>(small[i], m/i+1);
bool f = true;
for(int j=m-1;j>=1;j--){
if(f&&j%i==0&&tmp.size()>(j/i)&&tmp[j/i]){
tmp.resize(j+1);
break;
}
}
for(int j=tmp.size()-1;j>=1;j--){
if(j%i==0) tmp[j] = tmp[j/i];
else tmp[j] = 0;
}
ret_big *= tmp;
if(ret_big.size()>m) ret_big.resize(m);
}
return ret_big;
}
*/
/*
//1 ± x^a_i を大量に掛ける
template<ll MOD>
StaticModFPS<MOD> multiply_binomial_one(const vector<vector<ll>> &v, ll m){
using fps = StaticModFPS<MOD>;
vector<ll> pl(m, 0), mi(m, 0);
for(int i=0;i<v.size();i++){
if(v[i][0]==-1) mi[v[i][1]]++;
else pl[v[i][1]]++;
}
vector<ll> inv(m+1);
inv[0] = inv[1] = 1;
for(ll i=2;i<=m;i++){
ll u = MOD - (MOD/i);
inv[i] = (u*inv[MOD%i])%MOD;
}
fps ret(m, 0);
for(ll i=1;i<m;i++){
ll j = 1;
while(j*i<m){
ll mlt = (j%2==0?(-pl[i]-mi[i]):pl[i]-mi[i]);
mlt = (mlt * inv[j])%MOD;
mlt = (mlt + MOD)%MOD;
ret[i*j] += mlt;
if(ret[i*j]>=MOD) ret[i*j] -= MOD;
j++;
}
}
return ret.exp(m);
}
*/
template<typename mint>
mint _kth_coefficient(long long k, formal_power_series<mint> &p, formal_power_series<mint> &q){
using fps = formal_power_series<mint>;
if(k <= 100000) return (p.prefix(k + 1) * q.inv(k + 1))[k];
int n = p.deg_fix() + 1, m = q.deg_fix() + 1;
fps _q(m, 0);
for(int i = 0; i < m; i++) _q[i] = i & 1 ? -q[i] : q[i];
// p *= _q
// q *= _q
// pとqがほとんど同じサイズである場合は早い
int N = std::max(n, m), M = _q.size();
int z = 1 << ceil_pow2(N + M - 1);
p.resize(z), q.resize(z), _q.resize(z);
butterfly(p), butterfly(q), butterfly(_q);
mint iz = mint(z).inv();
for(int i = 0; i < z; i++){
p[i] *= _q[i];
q[i] *= _q[i];
}
butterfly_inv(p), butterfly_inv(q);
p.resize(n + M - 1), q.resize(m + M - 1);
//
n = p.size(), m = q.size();
for(int i = k & 1; i < n; i += 2) p[i / 2] = p[i] * iz;
for(int i = 0; i < m; i += 2) q[i / 2] = q[i] * iz;
p.resize((n + 1) / 2);
q.resize((m + 1) / 2);
return _kth_coefficient<mint>(k / 2, p, q);
}
template<typename mint>
mint kth_coefficient(long long k, formal_power_series<mint> p, formal_power_series<mint> q){
return _kth_coefficient<mint>(k, p, q);
}
// _a := 初期値, a[0, d)
// c := (a[i] = ∑{0 <= j < d} a[i - 1 - j] * c[j]を満たす)
// 0-indexed
template<typename mint>
mint kth_term_of_linear_reccurrence(long long k, const std::vector<mint> &_a, const std::vector<mint> &c){
using fps = formal_power_series<mint>;
int d = c.size();
assert(_a.size() == d);
if(k < d) return _a[k];
fps Q(d + 1);
Q[0] = 1;
for(int i = 0; i < d; i++) Q[i + 1] = -c[i];
fps P = Q * fps(_a);
P.resize(d);
return kth_coefficient<mint>(k, P, Q);
}
#line 1 ".lib/data_structure/range_query/mo_algorithm.hpp"
#line 1 ".lib/misc/random_number.hpp"
#line 5 ".lib/misc/random_number.hpp"
unsigned long long random_once(){
static std::random_device seed_gen;
static std::mt19937_64 engine(seed_gen());
static unsigned long long ret = engine();
return ret;
}
unsigned long long random_number(){
static std::random_device seed_gen;
static std::mt19937_64 engine(seed_gen());
return engine();
}
// [low, high]
unsigned long long random_number(unsigned long long low, unsigned long long high){
static std::random_device seed_gen;
static std::mt19937_64 engine(seed_gen());
assert(high >= low);
unsigned long long diff = high - low + 1;
return engine() % diff + low;
}
#line 5 ".lib/data_structure/range_query/mo_algorithm.hpp"
#include <math.h>
template<typename mo_struct>
struct mo_algorithm{
private:
struct query{
int i, b, l, r, label;
query(int i, int b, int l, int r, int label) : i(i), b(b), l(l), r(r), label(label){}
};
std::vector<query> q;
const int block_size, random_shift;
int pq, pl, pr;
mo_struct &_st;
void build(){
std::sort(q.begin(), q.end(), [](query &a, query &b){
if(a.b != b.b) return a.b < b.b;
return a.b & 1 ? a.r > b.r : a.r < b.r;
});
}
public:
mo_algorithm(int n, mo_struct &_st) : block_size(round(sqrt(n))), random_shift(0 % block_size), pq(0), pl(0), pr(0), _st(_st){}
inline int insert(int l, int r, int label = -1){
int qsz = q.size();
q.emplace_back(query(qsz, (l + random_shift) / block_size, l, r, label));
return qsz;
}
// すでに全てのクエリを処理し終わっている場合は{-1, -1}
// 今見ているクエリが追加された順番, ラベルを返す
inline std::pair<int, int> process(){
if(pq == 0) build();
if(pq == q.size()) return {-1, -1};
query &qi = q[pq];
while(pl > qi.l) _st.add_left(--pl);
while(pr < qi.r) _st.add_right(pr++);
while(pl < qi.l) _st.del_left(pl++);
while(pr > qi.r) _st.del_right(--pr);
pq++;
return {qi.i, qi.label};
}
// [l, r)が欲しい場合
inline std::pair<int, int> process2(){
if(pq == 0) build();
if(pq == q.size()) return {-1, -1};
query &qi = q[pq];
while(pl > qi.l) _st.add_left(--pl, pr);
while(pr < qi.r) _st.add_right(pl, pr++);
while(pl < qi.l) _st.del_left(pl++, pr);
while(pr > qi.r) _st.del_right(pl, --pr);
pq++;
return {qi.i, qi.label};
}
};
/*
struct range_inversion{
binary_indexed_tree<int> b;
std::vector<int> a;
int n;
long long sum = 0;
range_inversion(int n): n(n){}
void add_left(int i){
sum += b.query(0, a[i]);
b.update(a[i], 1);
}
void del_left(int i){
sum -= b.query(0, a[i]);
b.update(a[i], -1);
}
void add_right(int i){
sum += b.query(a[i] + 1, n);
b.update(a[i], 1);
}
void del_right(int i){
sum -= b.query(a[i] + 1, n);
b.update(a[i], -1);
}
};
*/
#line 7 ".lib/math/combinatorics.hpp"
// f(i, j) := comb(i, j) * comb(n - i, m - j), {0 <= i <= n, 0 <= j <= m}
// h(i, j) := ∑{0 <= k <= j} f(i, j) <-jについての和
// f(i, j)は(n - m) * mグリッドにおける点(i - j, j)を通る最短経路{左上(0, 0)と右下(n, m)の最短経路}の数
// f(i, j)とf(i, j + 1)は最短経路の数え上げにおいて独立{(i - j, j)と(i - j - 1, j + 1)を両方通ることが不可能}
// h(i, j)は(i - k, k) {0 <= k <= i}を通る場合の数(ナナメ)である
// 遷移
// h(i, j + 1) = h(i, j) + f(i, j + 1)
// h(i, j - 1) = h(i, j) - f(i, j)
// h(i, j) -> h(i + 1, j) : 点(i - j, j)および点(i - j, j + 1)を通る最短経路の数を引く -comb(i, j) * comb(n - i, m - j)
// h(i, j) -> h(i - 1, j) : 点(i - j - 1, j)および点(i - j - 1, j + 1)を通る最短経路の数を足す +comb(i - 1, j) * comb(n - i + 1, m - j)
/* #を通る最短経路の数
h(3, 1)
s..........
...........
.#.........
#.........t
h(3, 1)の言い換え
s..........
...........
##.........
..........t
s..........
...........
.#.........
.#........t
*/
template<typename mint>
struct diagonal_path_sum_struct{
int n, m;
modcomb<mint> mcb;
// f(0, 0) = comb(n, m)
mint x;
diagonal_path_sum_struct(int n, int m): n(n), m(m), mcb(n + 1), x(mcb.comb(n, m)){}
void add_left(int l, int r){
x -= mcb.comb(r - 1, l + 1) * mcb.comb(n - r + 1, m - l - 1);
}
void add_right(int l, int r){
r--;
x -= mcb.comb(r, l) * mcb.comb(n - r - 1, m - l - 1);
}
void del_left(int l, int r){
x += mcb.comb(r - 1, l + 1) * mcb.comb(n - r + 1, m - l - 1);
}
void del_right(int l, int r){
r--;
x += mcb.comb(r, l) * mcb.comb(n - r - 1, m - l - 1);
}
};
template<typename mint>
std::vector<mint> diagonal_path_sum(int n, int m, std::vector<std::pair<int, int>> Q){
diagonal_path_sum_struct<mint> dp(n, m);
mo_algorithm mo(n + 1, dp);
int q = Q.size();
std::vector<mint> ans(q);
for(int i = 0; i < q; i++) mo.insert(Q[i].first, Q[i].second + 1);
for(int i = 0; i < q; i++) ans[mo.process2().first] = dp.x;
return ans;
}
/*
// f(i, j) = ∑{0 <= k <= j} nCi * (n - i)Ck
// n個のものを{i, j, n - i - j}の3つに分ける場合の数のjに関する和
// f(i, j) = nCi * (n - i + 1)C(j + 1)
template<typename mint>
mint split_3_row(int n, int i, int j){
static modcomb<mint> mcb;
mcb.recalc(n + 1);
return mcb.comb(n, i) * mcb.comb(n - i + 1, j + 1);
}
*/
// f(i, j) = iCj + iC(j - 2) + iC(j - 4) ....
template<typename mint>
struct stripe_sum_struct{
modcomb<mint> mcb;
mint i2;
// f(0, 0) = 1
mint x, y, tmp;
stripe_sum_struct(): i2(mint(2).inv()), x(1), y(0){}
void add_left(int l, int r){
tmp = x - mcb.comb(r - 1, l + 1);
x = y;
y = tmp;
}
void del_left(int l, int r){
tmp = y + mcb.comb(r - 1, l + 1);
y = x;
x = tmp;
}
void add_right(int l, int r){
tmp = x + y; // ∑{0 <= k <= j}iCk
x = tmp;
y = tmp - mcb.comb(r, l);
}
void del_right(int l, int r){
r--;
tmp = x + y; // ∑{0 <= k <= j}iCk
tmp = (tmp + mcb.comb(r, l)) * i2; // ∑{0 <= k <= j}(i - 1)Ck
tmp = (tmp + mcb.comb(r - 1, l)) * i2; // ∑{0 <= k <= j}(i - 2)Ck
x = tmp;
y = tmp - mcb.comb(r, l);
}
};
template<typename mint>
std::vector<mint> stripe_sum(std::vector<std::pair<int, int>> Q){
int max_n = 0;
for(auto [i, j] : Q) max_n = std::max(max_n, i); // nCrの最大のn
stripe_sum_struct<mint> st;
st.mcb.recalc(max_n + 1);
mo_algorithm mo(max_n + 1, st);
int q = Q.size();
std::vector<mint> ans(q);
for(int i = 0; i < q; i++){
if(Q[i].first < Q[i].second) Q[i].second = Q[i].first;
mo.insert(Q[i].second, Q[i].first + 1);
}
for(int i = 0; i < q; i++){
auto [idx, id] = mo.process2();
ans[idx] = st.x;
}
return ans;
}
// f(i, j) := ∑{0 <= k <= j} comb(i, k)
// パスカルの三角形におけるi行目[0, j]列目の和
// 遷移
// f(i, j + 1) = f(i, j) + comb(i, j + 1)を足す
// f(i, j - 1) = f(i, j) - comb(i, j)を引く
// f(i + 1, j) = f(i, j) + f(i, j - 1)
// = 2 * f(i, j) - comb(i, j) = f(i + 1, j)
// f(i, j) -> f(i - 1, j)も同様
template<typename mint>
struct pascal_sum_row_struct{
modcomb<mint> mcb;
mint i2 = mint(2).inv();
mint x = i2;
pascal_sum_row_struct(){}
void add_left(int l, int r){
x -= mcb.comb(r - 1, l + 1);
}
void del_left(int l, int r){
x += mcb.comb(r - 1, l + 1);
}
void add_right(int l, int r){
r--;
x = 2 * x - mcb.comb(r, l);
}
void del_right(int l, int r){
r--;
x = (x + mcb.comb(r, l)) * i2;
}
};
/*
i < jの場合i = jにする
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
*/
template<typename mint>
std::vector<mint> pascal_sum_row(std::vector<std::pair<int, int>> Q){
int max_n = 0;
for(auto [i, j] : Q) max_n = std::max(max_n, i); // nCrの最大のn
pascal_sum_row_struct<mint> pas;
pas.mcb.recalc(max_n + 1);
mo_algorithm mo(max_n + 1, pas);
int q = Q.size();
std::vector<mint> ans(q);
for(int i = 0; i < q; i++){
if(Q[i].first < Q[i].second) Q[i].second = Q[i].first;
mo.insert(Q[i].second, Q[i].first + 1);
}
for(int i = 0; i < q; i++){
auto [idx, id] = mo.process2();
ans[idx] = pas.x;
}
return ans;
}
// 等差数列の重み付き
// g(i, j) := ∑{0 <= k <= j} (ak + b) * comb(i, k)
// 遷移
// g(i, j) -> g(i, j + 1): (a(j + 1) + b) * comb(i, j + 1)を足す
// g(i, j) -> g(i, j - 1): (aj + b) * comb(i, j)を引く
// g(i, j) -> g(i + 1, j):
// 各項(ak + b) * comb(i + 1, k)
// = (ak + b) * (comb(i, k) + comb(i, k - 1))
// = (ak + b) * comb(i, k) + (a(k - 1) + b) * comb(i, k - 1) + a * comb(i, k - 1)
// より
// g(i + 1, j) = g(i, j) + g(i, j - 1) + a * f(i, j - 1)
// = 2 * g(i, j) - (aj + b) * comb(i, j) + a * f(i, j) - a * comb(i, j)
// = 2 * g(i, j) + a * f(i, j) - (aj + a + b) * comb(i, j)
// fはpascal_sum_row, fとgを同時に計算する
// g(i - 1, j) = (g(i, j) - a * f(i - 1, j) + (aj + a + b) * comb(i - 1, j)) / 2
template<typename mint>
struct pascal_sum_row_arithmetic_struct{
modcomb<mint> mcb;
mint i2 = mint(2).inv();
// [0, 1)が合うように初期値を設定
// f(0, 0) = 1
// g(0, 0) = b
mint a, b, g, f = i2;
pascal_sum_row_arithmetic_struct(mint _a, mint _b): a(_a), b(_b), g(b * i2 - a * i2 * i2){}
void add_left(int l, int r){
l++;
g -= (a * l + b) * mcb.comb(r - 1, l);
f -= mcb.comb(r - 1, l);
}
void del_left(int l, int r){
l++;
g += (a * l + b) * mcb.comb(r - 1, l);
f += mcb.comb(r - 1, l);
}
void add_right(int l, int r){
r--;
g = 2 * g + a * f - (a * l + a + b) * mcb.comb(r, l);
f = 2 * f - mcb.comb(r, l);
}
void del_right(int l, int r){
r--;
f = (f + mcb.comb(r, l)) * i2;
g = (g - a * f + (a * l + a + b) * mcb.comb(r - 1, l)) * i2;
}
};
template<typename mint>
std::vector<mint> pascal_sum_row_arithmetic_weighted(mint a, mint b, std::vector<std::pair<int, int>> Q){
int max_n = 0;
for(auto [i, j] : Q) max_n = std::max(max_n, i); // nCrの最大のn
pascal_sum_row_arithmetic_struct<mint> pas(a, b);
pas.mcb.recalc(max_n + 1);
mo_algorithm mo(max_n + 1, pas);
int q = Q.size();
std::vector<mint> ans(q);
for(int i = 0; i < q; i++){
if(Q[i].first < Q[i].second) Q[i].second = Q[i].first;
mo.insert(Q[i].second, Q[i].first + 1);
}
for(int i = 0; i < q; i++){
auto [idx, id] = mo.process2();
ans[idx] = pas.g;
}
return ans;
}
// 等比級数の重み付き
// f(i, j) := ∑{0 <= k <= j} ab^k * comb(i, k)
// 遷移
// f(i, j) -> f(i, j + 1): ab^(j + 1) * comb(i, j + 1) を足す
// f(i, j) -> f(i, j - 1): ab^j * comb(i, j)を引く
// f(i, j) -> f(i + 1, j)
// f(i + 1, j) = f(i, j) + b * f(i, j - 1)
// = (1 + b) * f(i, j) - b * ab^j * comb(i, j)
//
// f(i - 1, j) = (f(i, j) + b * ab^j * comb(i - 1, j)) / (1 + b)
template<typename mint>
struct pascal_sum_row_geometric_struct{
modcomb<mint> mcb;
// [0, 1)が合うように初期値を設定
// f(0, 0) = a
const mint a, b, binv, b_plus_oneinv;
mint a_b_pow_l, f;
pascal_sum_row_geometric_struct(mint _a, mint _b): a(_a), b(_b), binv(b.inv()), b_plus_oneinv((b + 1).inv()), a_b_pow_l(a), f(a * b_plus_oneinv){}
void add_left(int l, int r){
l++;
f -= a_b_pow_l * mcb.comb(r - 1, l);
a_b_pow_l *= binv;
}
void del_left(int l, int r){
l++;
a_b_pow_l *= b;
f += a_b_pow_l * mcb.comb(r - 1, l);
}
void add_right(int l, int r){
r--;
f = (1 + b) * f - b * a_b_pow_l * mcb.comb(r, l);
}
void del_right(int l, int r){
r--;
f = (f + b * a_b_pow_l * mcb.comb(r - 1, l)) * b_plus_oneinv;
}
};
template<typename mint>
std::vector<mint> pascal_sum_row_geometric_weighted(mint a, mint b, std::vector<std::pair<int, int>> Q){
int max_n = 0;
for(auto [i, j] : Q) max_n = std::max(max_n, i); // nCrの最大のn
pascal_sum_row_geometric_struct<mint> pas(a, b);
pas.mcb.recalc(max_n + 1);
mo_algorithm mo(max_n + 1, pas);
int q = Q.size();
std::vector<mint> ans(q);
for(int i = 0; i < q; i++){
if(Q[i].first < Q[i].second) Q[i].second = Q[i].first;
mo.insert(Q[i].second, Q[i].first + 1);
}
for(int i = 0; i < q; i++){
auto [idx, id] = mo.process2();
ans[idx] = pas.f;
}
return ans;
}
// f(i, j) := ∑{0 <= k <= i} comb(k, j)
// パスカルの三角形における[0, i]行目j列目の和
// f(i, j) = comb(i + 1, j + 1)
// comb(i + 1, j + 1) = comb(i, j) + comb(i, j + 1)
// = comb(i, j) + comb(i - 1, j) + comb(i - 1, j + 1)
// ...
// = f(i, j)
template<typename mint>
mint pascal_sum_column(int i, int j){
static modcomb<mint> mcb;
mcb.recalc(i + 1);
return mcb.comb(i + 1, j + 1);
}
// f(i, j) := ∑{0 <= k <= i} (ak + b) * comb(k, j)
// ガリガリ計算するとこの式になる
template<typename mint>
mint pascal_sum_column_arithmetic_weighted(mint a, mint b, int i, int j){
static modcomb<mint> mcb;
mcb.recalc(i + 1);
return (a * i + b) * mcb.comb(i + 1, j + 1) - a * mcb.comb(i + 1, j + 2);
}
// ガリガリ計算すると
// f(i, j + 1) = (b * f(i, j) - a * b^(i + 1) * comb(i + 1, j + 1)) / (1 - b)
// f(i, j - 1) = ((1 - b) * f(i, j) + a * b ^ (i + 1) * comb(i + 1, j)) / b
// になる. moじゃなくてもできそうな気がする
/*
template<typename mint>
mint pascal_sum_column_geometric_weighted(mint a, mint b, int i, int j){
}
*/
// f(i, j) = ∑{0 <= k <= i}∑{0 <= l <= j} comb(k, l)
// 遷移
// f(i, j) = ∑{1 <= k <= j + 1} comb(i + 1, k)
// O(N√Q)
template<typename mint>
std::vector<mint> pascal_sum(std::vector<std::pair<int, int>> Q){
int q = Q.size();
for(int i = 0; i < q; i++){
Q[i].first++;
Q[i].second++;
}
auto res = pascal_sum_row<mint>(Q);
for(int i = 0; i < q; i++) res[i] -= 1;
return res;
}
template<typename mint>
mint factorial_mod_large(long long n){
if(n >= mint::mod()) return 0;
if(!n) return 1;
using fps = formal_power_series<mint>;
static constexpr int block_size = round(sqrt(mint::mod())), block_n = mint::mod() / block_size;
static std::vector<mint> x, fx;
// init O(k * log^2 k), k = max(block_size, mod / block_size)
// O(√mod * log^2 mod) if block_size = block_n = √mod
if(fx.empty()){
x.resize(block_n);
for(int i = 0; i < block_n; i++) x[i] = (i + 1) * block_size;
std::vector<fps> tmp(block_size);
for(int i = 0; i < block_size; i++) tmp[i] = {-i, 1};
fps f = multiply_constant_degree<mint>(tmp); // x(x - 1)...(x - block_size + 1)
fx = multipoint_evaluation<mint>(f, x);
for(int i = 1; i < block_n; i++) fx[i] *= fx[i - 1];
}
// O(block_size)
auto restore_fact = [](long long a)->mint{
int b = a / block_size;
mint res = (!b ? 1 : fx[b - 1]);
for(long long i = (long long)b * block_size + 1; i <= a; i++) res *= i;
return res;
};
return restore_fact(n);
}
template<typename mint>
mint combination_mod_large(long long n, long long k){
if(n < 0 || k < 0 || n < k) return 0;
if(n == k || k == 0) return 1;
// [1]
// comb(n, k) = (n * n - 1 .... n - k + 1) / k !
// n - k + 1 > k となるように変形(nCk = nC(n - k))したとき, 分子がmod上で0を跨いでいるなら0
k = std::min(k, n - k);
long long numerator_min = n - k + 1;
if(numerator_min % mint::mod() == 0 || (n / mint::mod() != numerator_min / mint::mod())) return 0;
// [1]より1 <= k < mod かつ nPk ≡ (n%mod)Pk
assert(1 <= k && k < mint::mod());
n %= mint::mod();
return factorial_mod_large<mint>(n) * (factorial_mod_large<mint>(k) * factorial_mod_large<mint>(n - k)).inv();
}
template<typename mint>
mint permutation_mod_large(long long n, long long k){
if(n < 0 || k < 0 || n < k) return 0;
if(k == 0) return 1;
// [1]
// perm(n, k) = n * n - 1 .... n - k + 1
// mod上で0を跨いでいるなら0
long long numerator_min = n - k + 1;
if(numerator_min % mint::mod() == 0 || (n / mint::mod() != numerator_min / mint::mod())) return 0;
// [1]より1 <= k < mod かつ nPk ≡ (n%mod)Pk
assert(1 <= k && k < mint::mod());
n %= mint::mod();
return factorial_mod_large<mint>(n) / factorial_mod_large<mint>(n - k);
}
template<typename mint>
formal_power_series<mint> stirling_number_first_kind(int n){
if(!n) return {1};
std::vector<formal_power_series<mint>> f(n);
for(int i = 0; i < n; i++) f[i] = {-i, 1};
return multiply_constant_degree<mint>(f);
}
struct four_square_theorem{
int n;
std::vector<int> A;
std::vector<int> sz;
// 0 <= i < Nについて以下のテーブルを作る
// sum(A) = iを満たしlen(A)が最小(常に4つ以下になる)
// O(N√N)
// N = 10^6で1sec程度
void init(int _n){
n = _n;
A.resize(n, -1);
sz.resize(n, 5);
sz[0] = 0;
for(int i = 1; i < n; i++){
for(int j = 1; j * j <= i; j++){
if(sz[i - j * j] + 1 < sz[i]){
A[i] = j;
sz[i] = sz[i - j * j] + 1;
}
}
}
}
four_square_theorem(int _n){
init(_n);
}
// x = ∑A[i]^2 を満たす要素数最小のA(高々4つ)
std::vector<int> get_square(int x){
assert(n);
assert(0 <= x && x < n);
std::vector<int> ret;
while(x){
int y = A[x];
ret.push_back(y);
x -= y * y;
}
return ret;
}
// x = ∑A[i]^2 を満たすA(要素数最小とは限らない, N = 10^5, x < 10^18なら高々6要素になる)
// @param x <= 10^18
std::vector<unsigned long long> get_square_large(unsigned long long x){
std::vector<unsigned long long> ret;
while(n <= x){
unsigned long long y = sqrtl(x);
ret.push_back(y);
x -= y * y;
}
while(x){
int y = A[x];
ret.push_back(y);
x -= y * y;
}
return ret;
}
// x = ∑ A[i] * (A[i] + 1) / 2 を満たすA(3つ以下)
std::vector<int> get_trianglar(int x){
auto v = get_square(8 * x + 3);
assert(v.size() == 3);
std::vector<int> ret;
for(int i = 0; i < 3; i++) if(v[i] != 1) ret.push_back(v[i] / 2);
return ret;
}
// x = ∑ A[i] * (A[i] + 1) / 2 を満たすA(要素数最小とは限らない, N = 10^5, x < 10^18なら高々5要素になる)
// @param x <= 10^18
std::vector<unsigned long long> get_trianglar_large(unsigned long long x){
std::vector<unsigned long long> ret;
while(n <= 8 * x + 3){
unsigned long long y = sqrtl(2 * x);
while(y * (y + 1) > 2 * x) y--;
ret.push_back(y);
x -= (y & 1 ? y * ((y >> 1) + 1) : (y >> 1) * (y + 1));
}
auto v = get_square(8 * x + 3);
for(int i = 0; i < 3; i++) if(v[i] != 1) ret.push_back(v[i] / 2);
return ret;
}
};
// ピタゴラス数について
// a^2 + b^2 = c^3を満たす自然数の組(a, b, c)をピタゴラス数と呼ぶ
// bが偶数, a, cが奇数でなければならない
// 自然数x, y(x > y)を用いて(a, b, c) = (x^2 - y^2, 2xy, x^2 + y^2)を表すことができる
#line 3 "b.cpp"
using mint = modint1000000007;
int main(){
io_init();
ll n;
std::cin >> n;
if(n >= mint::mod()){
std::cout << 0 << '\n';
return 0;
}
mint ans = factorial_mod_large<mint>(n);
std::cout << ans << '\n';
}