結果
| 問題 | No.665 Bernoulli Bernoulli |
| ユーザー |
 tonegawa
|
| 提出日時 | 2020-10-07 06:39:31 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 378 ms / 2,000 ms |
| コード長 | 28,195 bytes |
| コンパイル時間 | 4,109 ms |
| コンパイル使用メモリ | 167,068 KB |
| 実行使用メモリ | 7,056 KB |
| 最終ジャッジ日時 | 2023-09-27 09:07:20 |
| 合計ジャッジ時間 | 11,000 ms |
|
ジャッジサーバーID （参考情報） |
judge15 / judge12 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 2 ms
4,380 KB |
| testcase_01 | AC | 2 ms
4,380 KB |
| testcase_02 | AC | 378 ms
7,056 KB |
| testcase_03 | AC | 367 ms
6,860 KB |
| testcase_04 | AC | 364 ms
6,788 KB |
| testcase_05 | AC | 349 ms
6,800 KB |
| testcase_06 | AC | 353 ms
6,640 KB |
| testcase_07 | AC | 344 ms
6,760 KB |
| testcase_08 | AC | 354 ms
6,760 KB |
| testcase_09 | AC | 358 ms
6,716 KB |
| testcase_10 | AC | 348 ms
6,916 KB |
| testcase_11 | AC | 369 ms
6,872 KB |
| testcase_12 | AC | 361 ms
6,824 KB |
| testcase_13 | AC | 375 ms
6,904 KB |
| testcase_14 | AC | 376 ms
6,844 KB |
| testcase_15 | AC | 355 ms
6,712 KB |
| testcase_16 | AC | 354 ms
6,984 KB |
| testcase_17 | AC | 358 ms
6,720 KB |
| testcase_18 | AC | 347 ms
6,748 KB |
ソースコード
#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <queue>
#include <deque>
#include <algorithm>
#include <set>
#include <map>
#include <bitset>
#include <cmath>
#include <functional>
#include <cassert>
#include <iomanip>
#define vll vector<ll>
#define vvvl vector<vvl>
#define vvl vector<vector<ll>>
#define VV(a, b, c, d) vector<vector<d>>(a, vector<d>(b, c))
#define VVV(a, b, c, d) vector<vvl>(a, vvl(b, vll (c, d)));
#define re(c, b) for(ll c=0;c<b;c++)
#define all(obj) (obj).begin(), (obj).end()
typedef long long int ll;
typedef long double ld;
using namespace std;
#include <numeric>
#include <type_traits>
namespace internal{
template <class T> using is_integral = typename std::is_integral<T>;
template <class T>
using is_signed_int =
typename std::conditional<is_integral<T>::value && std::is_signed<T>::value, std::true_type, std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<is_integral<T>::value && std::is_unsigned<T>::value,
std::true_type, std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<is_signed_int<T>::value,
std::make_unsigned<T>, std::common_type<T>>::type;
template <class T>
using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
template <class T>
using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
//<internal_math>
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;
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;
unsigned long long x = (unsigned long long)(((unsigned __int128)(z)*im) >> 64);
unsigned int v = (unsigned int)(z - x * _m);
if (_m <= v) v += _m;
return v;
}
};
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 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};
}
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);
//<internal_bit>
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);
}
//<modint>
struct modint_base {};
struct static_modint_base : modint_base {};
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 <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : internal::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, internal::is_signed_int_t<T>* = nullptr>
static_modint(T v) {
long long x = (long long)(v % (long long)(umod()));
if (x < 0) x += umod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
static_modint(T v) {
_v = (unsigned int)(v % umod());
}
static_modint(bool v) { _v = ((unsigned int)(v) % umod()); }
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 = internal::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 = internal::is_prime<m>;
};
template <int id> struct dynamic_modint : internal::modint_base {
using mint = dynamic_modint;
public:
static int mod() { return (int)(bt.umod()); }
static void set_mod(int m) {
assert(1 <= m);
bt = internal::barrett(m);
}
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
dynamic_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
dynamic_modint(T v) {
long long x = (long long)(v % (long long)(mod()));
if (x < 0) x += mod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
dynamic_modint(T v) {
_v = (unsigned int)(v % mod());
}
dynamic_modint(bool v) { _v = ((unsigned int)(v) % mod()); }
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 = internal::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 internal::barrett bt;
static unsigned int umod() { return bt.umod(); }
};
template <int id> internal::barrett dynamic_modint<id>::bt = 998244353;
//using modint998244353 = static_modint<998244353>;
//using modint1000000007 = static_modint<1000000007>;
//using modint = dynamic_modint<-1>;
template <class T>
using is_static_modint = std::is_base_of<internal::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>;
//<conbovution>
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly(std::vector<mint>& a) {
static constexpr int g = internal::primitive_root<mint::mod()>;
int n = int(a.size());
int h = internal::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[internal::bsf(~(unsigned int)(s))];
}
}
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly_inv(std::vector<mint>& a) {
static constexpr int g = internal::primitive_root<mint::mod()>;
int n = int(a.size());
int h = internal::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[internal::bsf(~(unsigned int)(s))];
}
}
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> _convolution(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 << internal::ceil_pow2(n + m - 1);
a.resize(z);
internal::butterfly(a);
b.resize(z);
internal::butterfly(b);
for (int i = 0; i < z; i++) {
a[i] *= b[i];
}
internal::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, class T,
std::enable_if_t<internal::is_integral<T>::value>* = nullptr>
std::vector<T> _convolution(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 = internal::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(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;
}
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 {};
static constexpr unsigned long long MOD1 = 754974721; // 2^24
static constexpr unsigned long long MOD2 = 167772161; // 2^25
static constexpr unsigned long long MOD3 = 469762049; // 2^26
static constexpr unsigned long long M2M3 = MOD2 * MOD3;
static constexpr unsigned long long M1M3 = MOD1 * MOD3;
static constexpr unsigned long long M1M2 = MOD1 * MOD2;
static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3;
static constexpr unsigned long long i1 =
internal::inv_gcd(MOD2 * MOD3, MOD1).second;
static constexpr unsigned long long i2 =
internal::inv_gcd(MOD1 * MOD3, MOD2).second;
static constexpr unsigned long long i3 =
internal::inv_gcd(MOD1 * MOD2, MOD3).second;
auto c1 = internal::_convolution<MOD1>(a, b);
auto c2 = internal::_convolution<MOD2>(a, b);
auto c3 = internal::_convolution<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] - internal::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;
}
}
ll mpow(ll a, ll b, ll MOD = -1){
ll ret = 1, num = a;
while(b>0){
if(b%2) ret = (ret*num)%MOD;
num = (num*num)%MOD;
b /= 2;
}
return ret;
}
vector<ll> int32mod_conv(vector<ll> a, vector<ll> b, ll MOD, int MAXSIZE=-1){
if(MAXSIZE!=-1){
if(a.size()>MAXSIZE) a.resize(MAXSIZE);
if(b.size()>MAXSIZE) b.resize(MAXSIZE);
}
if(MOD==998244353) return internal::_convolution<998244353, ll>(a, b);
vector<ll> A = internal::_convolution<167772161>(a, b);
vector<ll> B = internal::_convolution<469762049>(a, b);
vector<ll> C = internal::_convolution<1224736769>(a, b);
ll N = A.size();
vector<ll> ret(N);
ll x = 167772161, y = 469762049, z = 1224736769;
ll ix = mpow(x, y-2, y);
ll ixy = mpow((x*y)%z, z-2, z);
for(int i=0;i<N;i++){
ll v = ((B[i] - A[i])*ix)%y;
if(v<0) v += y;
ll xxv = A[i]+x*v;
v = ((C[i] - (xxv%z))*ixy)%z;
if(v<0) v += z;
ret[i] = ((xxv%MOD) + ((x*y)%MOD)*v)%MOD;
}
if(MAXSIZE!=-1&&(int)ret.size()>MAXSIZE) ret.resize(MAXSIZE);
return ret;
}
//using modint998244353 = internal::static_modint<998244353>;
//using modint1000000007 = internal::static_modint<1000000007>;
//using modint = internal::dynamic_modint<-1>;
template<ll MOD, int MAXSIZE=-1>
struct StaticModFPS: vector<ll>{
using vector<ll>::vector;
using fps = StaticModFPS<MOD, MAXSIZE>;
StaticModFPS(vector<ll> v){
int n = v.size();
this->resize(n);
for(int i=0;i<n;i++){
(*this)[i] = v[i] % MOD;
if((*this)[i] < 0) (*this)[i] += MOD;
}
}
fps operator *= (const fps &vr) {
*this = int32mod_conv(*this, vr, MOD);
return *this;
}
fps operator /= (fps &vr){
return (*this) *= vr.inv();
}
fps operator += (const fps &vr){
int n = this->size();
int m = vr.size();
if(n < m) this->resize(m);
for(int i=0;i<m;i++) {
(*this)[i] += vr[i];
if((*this)[i] >= MOD) (*this)[i] -= MOD;
}
return *this;
}
fps operator -= (const fps &vr){
int n = this->size();
int m = vr.size();
if(n < m) this->resize(m);
for(int i=0;i<m;i++) {
(*this)[i] -= vr[i];
if((*this)[i] < 0) (*this)[i] += MOD;
}
return *this;
}
fps operator += (const ll &vr){
int n = this->size();
ll r = vr % MOD;
if(r < 0) r += MOD;
for(int i=0;i<n;i++){
(*this)[i] += r;
if((*this)[i] >= MOD) (*this)[i] -= MOD;
}
return *this;
}
fps operator -= (const ll &vr){
int n = this->size();
ll r = vr % MOD;
if(r<0) r += MOD;
for(int i=0;i<n;i++){
(*this[i]) -= r;
if((*this[i]) < 0) (*this)[i] += MOD;
}
return *this;
}
fps operator *= (const ll &vr){
int n = this->size();
ll r = vr % MOD;
if(r<0) r += MOD;
for(int i=0;i<n;i++){
(*this)[i] = ((*this)[i] * r)%MOD;
}
return *this;
}
fps operator /= (const ll &vr){
ll r = vr % MOD;
if(r < 0) r += MOD;
assert(r!=0);
r = mpow(r, MOD-2, MOD);
int n = (int)this->size();
for(int i=0;i<n;i++) (*this)[i] = ((*this)[i] * r)%MOD;
return *this;
}
fps operator + (const fps& vr){return fps(*this) += vr;}
fps operator - (const fps& vr){return fps(*this) -= vr;}
fps operator * (const fps& vr){return fps(*this) *= vr;}
fps operator / (const fps& vr){return fps(*this) /= vr;}
fps operator + (const ll& vr){return fps(*this) += vr;}
fps operator - (const ll& vr){return fps(*this) -= vr;}
fps operator * (const ll& vr){return fps(*this) *= vr;}
fps operator / (const ll& vr){return fps(*this) /= vr;}
void debug(int printsize = 20){
int n = min(20, (int)this->size());
for(int i=0;i<n;i++){
if(i==n-1) printf("%lld\n", (*this)[i]);
else printf("%lld ", (*this)[i]);
}
}
void print(){
int n = (int)this->size();
for(int i=0;i<n;i++){
if(i==n-1) printf("%lld\n", (*this)[i]);
else printf("%lld ", (*this)[i]);
}
}
fps rev(int deg=-1) const {
fps ret(*this);
if(deg != -1) ret.resize(deg, 0);
reverse(ret.begin(), ret.end());
return ret;
}
fps prefix(int deg){
int n = min((int)this->size(), deg);
return fps(this->begin(), this->begin() + n);
}
// https://en.wikipedia.org/wiki/Formal_power_series#Multiplicative_inverse
// invertible in R[X] if and only if its constant coefficient a_{0} is invertible in R
fps inv(int deg=-1){
assert((*this)[0]);
int n = this->size();
if(deg==-1) deg = n;
fps ret({mpow((*this)[0], MOD-2, MOD)});
for(int i=1;i<deg;i<<=1){
ret = ((ret + ret) - (ret * ret * prefix(i << 1))).prefix(i << 1);
}
return ret.prefix(deg);
}
fps modulo(const fps &vr){
int n = this->size();
int m = vr.size();
if(n<m) return (*this);
n = n - m + 1;
//n-m次多項式を返す
fps r = ((rev().prefix(n) * vr.rev().inv(n)).prefix(n).rev(n)) * vr;
return (*this - r).prefix(m-1);
}
};
// n: 多項式の次数
// m: 評価する点の数
// n >> mの場合, n, m共に小さい場合は愚直に計算したほうが早い
// n < mの場合はn<-mとなるとして、
// n * m < 10^6 なら愚直
// n * m < 2*10^8 かつ n / m > 2000
template<ll MOD>
vector<ll> simpleMultipointEvaluation(const StaticModFPS<MOD> &F, const vector<ll> &v){
int n = F.size();
int m = v.size();
vector<ll> ret(m, 0);
for(int i=0;i<m;i++){
ll x = 1;
for(int j=0;j<n;j++){
ret[i] = (ret[i] + (x * F[j])%MOD)%MOD;
if(ret[i] < 0) ret[i] += MOD;
x = (x * v[i])%MOD;
}
}
return ret;
}
template<ll MOD>
vector<ll> MultipointEvaluation(StaticModFPS<MOD> F, const vector<ll> &v){
typedef StaticModFPS<MOD> fps;
ll n = F.size();
ll m = v.size();
int add_coef = -1;
if(MOD==998244353){
if(max(n, m)*m < 2000000LL || (max(n, m)*m<100000000LL&&n/m>2500LL)) {
return simpleMultipointEvaluation<MOD>(F, v);
}
}else{
if(max(n, m)*m < 15000000LL || (max(n, m)*m<100000000LL&&n/m>1000LL)) {
return simpleMultipointEvaluation<MOD>(F, v);
}
}
if(n < m) {
add_coef = m-1;//m-1次に1を足す
F.resize(m);
F[m-1] = 1;
n = m;
}
int N = 1;
while(N < m) N *= 2;
vector<fps> tree(2*N-1, fps{1});
for(int i=0;i<m;i++) tree[N-1+i] = fps{(MOD-v[i]%MOD)%MOD, 1};
for(int i=N-2;i>=0;i--){
tree[i] = tree[i*2+1] * tree[i*2+2];
}
tree[0] = F.modulo(tree[0]);
for(int i=1;i<2*N-1;i++){
int par = (i-1)/2;
tree[i] = tree[par].modulo(tree[i]);
}
vector<ll> ret(m);
for(int i=0;i<m;i++) {
ret[i] = tree[N-1+i][0];
if(add_coef!=-1) ret[i] = (ret[i] - mpow(v[i], add_coef, MOD) + MOD)%MOD;
}
return ret;
}
//N次多項式f(x)について, x0, x1....xNの時のf(x0), f(x1)...f(xN)からfを求める
//z = (x-x0)(x-x1).....(x-xN)としてQi[i] = z/(x-xi)
// 係数C[i]は f(xi) = C[i]Qi[i](xi)よりC[i] = f(xi)/Qi[i](xi) <- 標本点が連続していれば高速化可能
// f(T) = ∑C[i]*Qi(T)
template<ll MOD>
ll rag(vvl sample, ll T){
// sample: xi, f(x1)
ll n = sample.size();
ll zT = 1, ret = 0;
vll QT(n), Qi(n), C(n);
for(ll i=0;i<n;i++) zT = (zT * (T-sample[i][0]+MOD)%MOD)%MOD;
for(ll i=0;i<n;i++){
ll tmp = 1;
for(ll j=0;j<n;j++){
if(i==j) continue;
tmp = (tmp * (sample[i][0] - sample[j][0] + MOD)%MOD)%MOD;
}
Qi[i] = tmp;
}
for(int j=0;j<n;j++) std::cout << Qi[j] << (j==n-1?"\n":" ");
for(ll i=0;i<n;i++) {
QT[i] = (zT * mpow((T-sample[i][0]+MOD)%MOD, MOD-2, MOD))%MOD;
C[i] = (sample[i][1] * mpow(Qi[i], MOD-2, MOD))%MOD;
}
for(ll i=0;i<n;i++) ret = (ret + (C[i]*QT[i])%MOD)%MOD;
return ret;
}
template<ll MOD>
StaticModFPS<MOD> InterpolationbyDivideandConquer(
const vector<ll> &fx, StaticModFPS<MOD> F, const vector<StaticModFPS<MOD>> &tree, int k, int l, int r){
if(r-l==1) return StaticModFPS<MOD>{(fx[l]*mpow(F[0], MOD-2, MOD))%MOD};
int mid = (l+r)/2;
if(tree[k*2+2].size()==0) return InterpolationbyDivideandConquer<MOD>(fx, F, tree, k*2+1, l, mid);
StaticModFPS<MOD> left = InterpolationbyDivideandConquer<MOD>(fx, F.modulo(tree[k*2+1]), tree, k*2+1, l, mid);
StaticModFPS<MOD> right = InterpolationbyDivideandConquer<MOD>(fx, F.modulo(tree[k*2+2]), tree, k*2+2, mid, r);
return left * tree[k*2+2] + right * tree[k*2+1];
}
template<ll MOD>
StaticModFPS<MOD> PolynomialInterpolation(const vector<ll> &xi, const vector<ll> &fx){
typedef StaticModFPS<MOD> fps;
int n = xi.size();
int N = 1;
while(N<n) N*=2;
vector<fps> tree(2*N-1, fps{});
for(int i=0;i<n;i++) tree[N-1+i] = fps{(MOD - xi[i]%MOD)%MOD, 1};
for(int i=N-2;i>=0;i--){
if(tree[i*2+2].size()==0) tree[i] = tree[i*2+1];
else tree[i] = tree[i*2+1] * tree[i*2+2];
}
for(ll i=0;i<n;i++) tree[0][i] = (tree[0][i+1] * (i+1))%MOD;
tree[0].pop_back();
return InterpolationbyDivideandConquer<MOD>(fx, tree[0], tree, 0, 0, N).prefix(n);
}
// O(N^2logN) (deg(f)=deg(g)=deg(ans)=Nとして)
template<ll MOD>
StaticModFPS<MOD> simpleComposition(const StaticModFPS<MOD> &f, const StaticModFPS<MOD> &g, int deg){
int n = f.size();
StaticModFPS<MOD> c(deg, 0), gpower{1};
for(int i=0;i<min(deg, n);i++){
int d = min(deg, (int)gpower.size());
for(int j=0;j<d;j++){
c[j] = (c[j] + f[i] * gpower[j])%MOD;
}
gpower *= g;
if(gpower.size() > deg) gpower.resize(deg);
}
return c;
}
// reference: http://www.eecs.harvard.edu/~htk/publication/1978-jacm-brent-kung.pdf
// N-1次多項式 f(x)にM-1次多項式 g(x)を合成した結果を求める
// deg(f) = deg(g) = deg(ans) = Nとして、
// f(x)をk:=ceil(√N+1)ブロックに平方分割すると 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項とN項の合成なので、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)
template<ll MOD>
StaticModFPS<MOD> Composition(const StaticModFPS<MOD> &f, const StaticModFPS<MOD> &g, int deg){
typedef StaticModFPS<MOD> fps;
int n = f.size();
int k = (int)sqrt(n);
if(k*k<n) k++;
int d = n / k;
if(k*d<n) d++;
vector<fps> gpower(k+1, {1});
for(int i=1;i<=k;i++){
gpower[i] = gpower[i-1] * g;
if(gpower[i].size()>deg) gpower[i].resize(deg);
}
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;
int sz = gpower[j].size();
for(int t=0;t<sz;t++){
fi[i][t] = (fi[i][t] + gpower[j][t] * f[idx])%MOD;
}
}
}
fps ret(deg, 0), gd = {1};
for(int i=0;i<k;i++){
fi[i] *= gd;
int sz = min(deg, (int)fi[i].size());
for(int j=0;j<sz;j++) ret[j] = (ret[j] + fi[i][j])%MOD;
gd *= gpower[d];
if(gd.size() > deg) gd.resize(deg);
}
return ret;
}
int main(){
ll P = 1000000007;
ll n, k;scanf("%lld %lld", &n, &k);
vector<ll> x(k+2), y(k+2);
for(ll i=0;i<=k+1;i++){
x[i] = i+1;
y[i] = ((i==0?0:y[i-1]) + mpow(i+1, k, P))%P;
}
StaticModFPS<1000000007> f = PolynomialInterpolation<1000000007>(x, y);
ll ans = 0;
for(int i=0;i<f.size();i++) ans = (ans + mpow(n%P, i, P)*f[i])%P;
printf("%lld\n", ans);
}