PROBLEM LINK:
Practice
Div-4 Contest
Div-3 Contest
Div-2 Contest
Div-1 Contest
Author: Aryan Choudhary
Tester: Aryan Choudhary, Istvan Nagy
Editorialist: Daanish Mahajan
DIFFICULTY:
Medium
PREREQUISITES:
Exponential Generating Functions, NTT
PROBLEM:
Given 3 integers N, P, and M.
- Let f(A) be the MEX of the array A.
- Let g(A) = (f(A))^P for the array A .
There are (M+1)^N distinct arrays A = [A_1, A_2, \ldots, A_N] such that (0 \le A_i \le M) \forall (1\le i \le N).
Calculate the sum of g(A) over all such (M+1)^N arrays. Since the answer can be huge, print it modulo 754974721.
EXPLANATION:
The problem is an application of contribution technique.
answer = \sum_{i = 0}^{M + 1} i^p \cdot number of ways i can be the mex of the array.
Number of ways
For i to be mex, following conditions must be satisfied:
- All numbers less than i must appear in the array.
- i must not be present in the array.
- The numbers greater than i may / maynot be present in the array.
Suppose x numbers less than i appear, which means N - x numbers greater than i appear in the array.
Assume the frequency of elements less than i appearing is F_0, F_1, \ldots, F_{i - 1} such that \sum_{j = 0}^{i - 1} F_j = x and F_j \gt 0.
So number of ways
= \sum_{x = i}^{N} \binom{N}{x} \cdot \frac{x !}{\prod_{j = 0}^{i - 1}F_j !} \cdot (M - i)^{N - x}
= \sum_{x = i}^{N} \frac{N !}{x ! \cdot (N - x) !} \cdot \frac{x !}{\prod_{j = 0}^{i - 1}F_j !} \cdot (M - i)^{N - x}
= N ! \cdot \sum_{x = i}^{N} \frac{1}{\prod_{j = 0}^{i - 1}F_j !} \cdot \frac{(M - i)^{N - x}}{(N - x) !}
= N ! \cdot \sum_{x = i}^{N} ([y ^ x](e ^ y - 1) ^ i) \cdot ([y ^ {N - x}] e ^ {y \cdot (M - i)}) ((e^y - 1) because F_j \gt 0)
= N ! \cdot \sum_{x = 0}^{N} ([y ^ x](e ^ y - 1) ^ i) \cdot ([y ^ {N - x}] e ^ {y \cdot (M - i)}) (limits changed from x = i to x = 0 since the coefficients from y^0 to y^{i - 1} will be zero in the first expression)
= N ! \cdot ([y ^ N]((e ^ y - 1) ^ i \cdot e ^ {y \cdot (M - i)}))
= N ! \cdot ((\sum_{i = 0}^{M} i ^ p \cdot ([y ^ N]((e ^ y - 1) ^ i \cdot e ^ {y \cdot (M - i)}))) + (M + 1) ^ p \cdot ([y ^ N]((e ^ y - 1) ^ {M + 1})) (M + 1 case is to be handled separately, you can easily check why?)
= N ! \cdot ((\sum_{i = 0}^{M} i ^ p \cdot ([y ^ N](\sum_{j = 0}^{i} \binom{i}{j} \cdot (-1)^j \cdot e ^ {y \cdot (i - j)} \cdot e ^ {y \cdot (M - i)}))) + (M + 1) ^ p \cdot ([y ^ N] (\sum_{j = 0}^{M + 1} \binom{M + 1}{j} \cdot (-1) ^ j \cdot e^{y \cdot (M + 1 - j)})))
= N ! \cdot ((\sum_{i = 0}^{M} i ^ p \cdot (\sum_{j = 0}^{i} \binom{i}{j} \cdot (-1)^j \cdot \frac{(M - j) ^ N}{N !})) + (M + 1) ^ p \cdot (\sum_{j = 0}^{M + 1} \binom{M + 1}{j} \cdot (-1)^j \cdot \frac{(M + 1 - j) ^ N}{N !}))
= (\sum_{i = 0}^{M} i ^ p \cdot (\sum_{j = 0}^{i} \binom{i}{j} \cdot (-1)^j \cdot (M - j) ^ N)) + (M + 1) ^ p \cdot (\sum_{j = 0}^{M + 1} \binom{M + 1}{j} \cdot (-1)^j \cdot (M + 1 - j) ^ N)
C_i
\sum_{i = 0}^{M} i ^ p \cdot (\sum_{j = 0}^{i} \binom{i}{j} \cdot (-1)^j \cdot (M - j) ^ N)
= \sum_{i = 0}^{M} i ^ p \cdot (\sum_{j = 0}^{i} \frac{i !} {(i - j) ! \cdot j !} \cdot (-1)^j \cdot (M - j) ^ N)
= \sum_{i = 0}^{M} i ^ p \cdot i ! (\sum_{j = 0}^{i} \frac{1}{(i - j) !} \cdot \frac{(-1)^j \cdot (M - j) ^ N}{j !})
Now, C_i = [x^i](A(x) \cdot B(x))
where:
A(x) = \sum_{j = 0}^{M} \frac{1}{j !} \cdot x ^ j
B(x) = \sum_{j = 0}^{M} \frac{(-1)^j \cdot (M - j) ^ N}{j !} \cdot x ^ j
= (\sum_{i = 0}^{M} i ^ p \cdot i ! \cdot C_i) + (M + 1) ^ p \cdot (\sum_{j = 0}^{M + 1} \binom{M + 1}{j} \cdot (-1)^j \cdot (M + 1 - j) ^ N)
where C(x) can be calculated using NTT in \mathcal{O}(M\log{}M) and the second term in summation can be calculated in \mathcal{O}(M).
COMPLEXITY ANALYSIS:
\mathcal{O}(M\log{}M)
SOLUTIONS:
Setter's Solution
#include <cassert>
#include <numeric>
#include <type_traits>
#ifdef _MSC_VER
#include <intrin.h>
#endif
#include <utility>
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder {
namespace internal {
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 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; // |m1 * u| <= |m1| * s <= b
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);
unsigned long long floor_sum_unsigned(unsigned long long n,
unsigned long long m,
unsigned long long a,
unsigned long long b) {
unsigned long long ans = 0;
while (true) {
if (a >= m) {
ans += n * (n - 1) / 2 * (a / m);
a %= m;
}
if (b >= m) {
ans += n * (b / m);
b %= m;
}
unsigned long long y_max = a * n + b;
if (y_max < m) break;
n = (unsigned long long)(y_max / m);
b = (unsigned long long)(y_max % m);
std::swap(m, a);
}
return ans;
}
} // namespace internal
} // namespace atcoder
#include <cassert>
#include <numeric>
#include <type_traits>
namespace atcoder {
namespace internal {
#ifndef _MSC_VER
template <class T>
using is_signed_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value ||
std::is_same<T, __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int128 =
typename std::conditional<std::is_same<T, __uint128_t>::value ||
std::is_same<T, unsigned __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using make_unsigned_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value,
__uint128_t,
unsigned __int128>;
template <class T>
using is_integral = typename std::conditional<std::is_integral<T>::value ||
is_signed_int128<T>::value ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_signed_int = typename std::conditional<(is_integral<T>::value &&
std::is_signed<T>::value) ||
is_signed_int128<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) ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<
is_signed_int128<T>::value,
make_unsigned_int128<T>,
typename std::conditional<std::is_signed<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type>::type;
#else
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;
#endif
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;
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
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>;
} // namespace internal
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());
}
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());
}
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>;
namespace internal {
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>;
} // namespace internal
} // namespace atcoder
using mint2 = atcoder::static_modint<754974721>;
#ifdef ARYANC403
#include <header.h>
#else
#pragma GCC optimize ("Ofast")
#pragma GCC target ("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx")
//#pragma GCC optimize ("-ffloat-store")
#include<bits/stdc++.h>
#define dbg(args...) 42;
#endif
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
using namespace std;
using ll = long long int;
mt19937_64 rng(chrono::high_resolution_clock::now().time_since_epoch().count());
using namespace std;
#define fo(i,n) for(i=0;i<(n);++i)
#define repA(i,j,n) for(i=(j);i<=(n);++i)
#define repD(i,j,n) for(i=(j);i>=(n);--i)
#define all(x) begin(x), end(x)
#define sz(x) ((lli)(x).size())
#define pb push_back
#define mp make_pair
#define X first
#define Y second
#define endl "\n"
typedef long long int lli;
typedef long double mytype;
typedef pair<lli,lli> ii;
typedef vector<ii> vii;
typedef vector<lli> vi;
// https://judge.yosupo.jp/submission/69895
struct IOPre {
static constexpr int TEN = 10, SZ = TEN * TEN * TEN * TEN;
std::array<char, 4 * SZ> num;
constexpr IOPre() : num{} {
for (int i = 0; i < SZ; i++) {
int n = i;
for (int j = 3; j >= 0; j--) {
num[i * 4 + j] = static_cast<char>(n % TEN + '0');
n /= TEN;
}
}
}
};
struct IO {
#if !HAVE_DECL_FREAD_UNLOCKED
#define fread_unlocked fread
#endif
#if !HAVE_DECL_FWRITE_UNLOCKED
#define fwrite_unlocked fwrite
#endif
static constexpr int SZ = 1 << 17, LEN = 32, TEN = 10, HUNDRED = TEN * TEN,
THOUSAND = HUNDRED * TEN, TENTHOUSAND = THOUSAND * TEN,
MAGIC_MULTIPLY = 205, MAGIC_SHIFT = 11, MASK = 15,
TWELVE = 12, SIXTEEN = 16;
static constexpr IOPre io_pre = {};
std::array<char, SZ> input_buffer, output_buffer;
int input_ptr_left, input_ptr_right, output_ptr_right;
IO()
: input_buffer{},
output_buffer{},
input_ptr_left{},
input_ptr_right{},
output_ptr_right{} {}
IO(const IO&) = delete;
IO(IO&&) = delete;
IO& operator=(const IO&) = delete;
IO& operator=(IO&&) = delete;
~IO() { flush(); }
template <class T>
struct is_char {
static constexpr bool value = std::is_same_v<T, char>;
};
template <class T>
struct is_bool {
static constexpr bool value = std::is_same_v<T, bool>;
};
template <class T>
struct is_string {
static constexpr bool value =
std::is_same_v<T, std::string> || std::is_same_v<T, const char*> ||
std::is_same_v<T, char*> || std::is_same_v<std::decay_t<T>, char*>;
;
};
template <class T, class D = void>
struct is_custom {
static constexpr bool value = false;
};
template <class T>
struct is_custom<T, std::void_t<typename T::internal_value_type>> {
static constexpr bool value = true;
};
template <class T>
struct is_default {
static constexpr bool value = is_char<T>::value || is_bool<T>::value ||
is_string<T>::value ||
std::is_integral_v<T>;
};
template <class T, class D = void>
struct is_iterable {
static constexpr bool value = false;
};
template <class T>
struct is_iterable<
T, typename std::void_t<decltype(std::begin(std::declval<T>()))>> {
static constexpr bool value = true;
};
template <class T, class D = void, class E = void>
struct is_applyable {
static constexpr bool value = false;
};
template <class T>
struct is_applyable<T, std::void_t<typename std::tuple_size<T>::type>,
std::void_t<decltype(std::get<0>(std::declval<T>()))>> {
static constexpr bool value = true;
};
template <class T>
static constexpr bool needs_newline = (is_iterable<T>::value ||
is_applyable<T>::value) &&
(!is_default<T>::value);
template <typename T, typename U>
struct any_needs_newline {
static constexpr bool value = false;
};
template <typename T>
struct any_needs_newline<T, std::index_sequence<>> {
static constexpr bool value = false;
};
template <typename T, std::size_t I, std::size_t... Is>
struct any_needs_newline<T, std::index_sequence<I, Is...>> {
static constexpr bool value =
needs_newline<decltype(std::get<I>(std::declval<T>()))> ||
any_needs_newline<T, std::index_sequence<Is...>>::value;
};
inline void load() {
memmove(std::begin(input_buffer),
std::begin(input_buffer) + input_ptr_left,
input_ptr_right - input_ptr_left);
input_ptr_right =
input_ptr_right - input_ptr_left +
static_cast<int>(fread_unlocked(
std::begin(input_buffer) + input_ptr_right - input_ptr_left, 1,
SZ - input_ptr_right + input_ptr_left, stdin));
input_ptr_left = 0;
}
inline void read_char(char& c) {
if (input_ptr_left + LEN > input_ptr_right) load();
c = input_buffer[input_ptr_left++];
}
inline void read_string(std::string& x) {
char c;
while (read_char(c), c < '!') continue;
x = c;
while (read_char(c), c >= '!') x += c;
}
template <class T>
inline std::enable_if_t<std::is_integral_v<T>, void> read_int(T& x) {
if (input_ptr_left + LEN > input_ptr_right) load();
char c = 0;
do c = input_buffer[input_ptr_left++];
while (c < '-');
[[maybe_unused]] bool minus = false;
if constexpr (std::is_signed<T>::value == true)
if (c == '-') minus = true, c = input_buffer[input_ptr_left++];
x = 0;
while (c >= '0')
x = x * TEN + (c & MASK), c = input_buffer[input_ptr_left++];
if constexpr (std::is_signed<T>::value == true)
if (minus) x = -x;
}
inline void skip_space() {
if (input_ptr_left + LEN > input_ptr_right) load();
while (input_buffer[input_ptr_left] <= ' ') input_ptr_left++;
}
inline void flush() {
fwrite_unlocked(std::begin(output_buffer), 1, output_ptr_right, stdout);
output_ptr_right = 0;
}
inline void write_char(char c) {
if (output_ptr_right > SZ - LEN) flush();
output_buffer[output_ptr_right++] = c;
}
inline void write_bool(bool b) {
if (output_ptr_right > SZ - LEN) flush();
output_buffer[output_ptr_right++] = b ? '1' : '0';
}
inline void write_string(const std::string& s) {
for (auto x : s) write_char(x);
}
inline void write_string(const char* s) {
while (*s) write_char(*s++);
}
inline void write_string(char* s) {
while (*s) write_char(*s++);
}
template <typename T>
inline std::enable_if_t<std::is_integral_v<T>, void> write_int(T x) {
if (output_ptr_right > SZ - LEN) flush();
if (!x) {
output_buffer[output_ptr_right++] = '0';
return;
}
if constexpr (std::is_signed<T>::value == true)
if (x < 0) output_buffer[output_ptr_right++] = '-', x = -x;
int i = TWELVE;
std::array<char, SIXTEEN> buf{};
while (x >= TENTHOUSAND) {
memcpy(std::begin(buf) + i,
std::begin(io_pre.num) + (x % TENTHOUSAND) * 4, 4);
x /= TENTHOUSAND;
i -= 4;
}
if (x < HUNDRED) {
if (x < TEN) {
output_buffer[output_ptr_right++] = static_cast<char>('0' + x);
} else {
std::uint32_t q =
(static_cast<std::uint32_t>(x) * MAGIC_MULTIPLY) >>
MAGIC_SHIFT;
std::uint32_t r = static_cast<std::uint32_t>(x) - q * TEN;
output_buffer[output_ptr_right] = static_cast<char>('0' + q);
output_buffer[output_ptr_right + 1] =
static_cast<char>('0' + r);
output_ptr_right += 2;
}
} else {
if (x < THOUSAND) {
memcpy(std::begin(output_buffer) + output_ptr_right,
std::begin(io_pre.num) + (x << 2) + 1, 3),
output_ptr_right += 3;
} else {
memcpy(std::begin(output_buffer) + output_ptr_right,
std::begin(io_pre.num) + (x << 2), 4),
output_ptr_right += 4;
}
}
memcpy(std::begin(output_buffer) + output_ptr_right,
std::begin(buf) + i + 4, TWELVE - i);
output_ptr_right += TWELVE - i;
}
template <typename T_>
IO& operator<<(T_&& x) {
using T = typename std::remove_cv<
typename std::remove_reference<T_>::type>::type;
static_assert(is_custom<T>::value or is_default<T>::value or
is_iterable<T>::value or is_applyable<T>::value);
if constexpr (is_custom<T>::value) {
write_int(x.get());
} else if constexpr (is_default<T>::value) {
if constexpr (is_bool<T>::value) {
write_bool(x);
} else if constexpr (is_string<T>::value) {
write_string(x);
} else if constexpr (is_char<T>::value) {
write_char(x);
} else if constexpr (std::is_integral_v<T>) {
write_int(x);
}
} else if constexpr (is_iterable<T>::value) {
// strings are immune
using E = decltype(*std::begin(x));
constexpr char sep = needs_newline<E> ? '\n' : ' ';
int i = 0;
for (const auto& y : x) {
if (i++) write_char(sep);
operator<<(y);
}
} else if constexpr (is_applyable<T>::value) {
// strings are immune
constexpr char sep =
(any_needs_newline<
T, std::make_index_sequence<std::tuple_size_v<T>>>::value)
? '\n'
: ' ';
int i = 0;
std::apply(
[this, &sep, &i](auto const&... y) {
(((i++ ? write_char(sep) : void()), this->operator<<(y)),
...);
},
x);
}
return *this;
}
template <typename T>
IO& operator>>(T& x) {
static_assert(is_custom<T>::value or is_default<T>::value or
is_iterable<T>::value or is_applyable<T>::value);
static_assert(!is_bool<T>::value);
if constexpr (is_custom<T>::value) {
typename T::internal_value_type y;
read_int(y);
x = y;
} else if constexpr (is_default<T>::value) {
if constexpr (is_string<T>::value) {
read_string(x);
} else if constexpr (is_char<T>::value) {
read_char(x);
} else if constexpr (std::is_integral_v<T>) {
read_int(x);
}
} else if constexpr (is_iterable<T>::value) {
for (auto& y : x) operator>>(y);
} else if constexpr (is_applyable<T>::value) {
std::apply([this](auto&... y) { ((this->operator>>(y)), ...); }, x);
}
return *this;
}
IO* tie(std::nullptr_t) { return this; }
void sync_with_stdio(bool) {}
};
IO io;
#define cin io
#define cout io
namespace ntt {
template <class T, class F = multiplies<T>>
T power(T a, long long n, F op = multiplies<T>(), T e = {1}) {
// assert(n >= 0);
T res = e;
while (n) {
if (n & 1) res = op(res, a);
if (n >>= 1) a = op(a, a);
}
return res;
}
constexpr int mod = int(1e9) + 7;
constexpr int nttmod = 754'974'721;
template <std::uint32_t P>
struct ModInt32 {
public:
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using m32 = ModInt32;
using internal_value_type = u32;
private:
u32 v;
static constexpr u32 get_r() {
u32 iv = P;
for (u32 i = 0; i != 4; ++i) iv *= 2U - P * iv;
return -iv;
}
static constexpr u32 r = get_r(), r2 = -u64(P) % P;
static_assert((P & 1) == 1);
static_assert(-r * P == 1);
static_assert(P < (1 << 30));
static constexpr u32 pow_mod(u32 x, u64 y) {
u32 res = 1;
for (; y != 0; y >>= 1, x = u64(x) * x % P)
if (y & 1) res = u64(res) * x % P;
return res;
}
static constexpr u32 reduce(u64 x) {
return (x + u64(u32(x) * r) * P) >> 32;
}
static constexpr u32 norm(u32 x) { return x - (P & -(x >= P)); }
public:
static constexpr u32 get_pr() {
u32 tmp[32] = {}, cnt = 0;
const u64 phi = P - 1;
u64 m = phi;
for (u64 i = 2; i * i <= m; ++i)
if (m % i == 0) {
tmp[cnt++] = i;
while (m % i == 0) m /= i;
}
if (m != 1) tmp[cnt++] = m;
for (u64 res = 2; res != P; ++res) {
bool flag = true;
for (u32 i = 0; i != cnt && flag; ++i)
flag &= pow_mod(res, phi / tmp[i]) != 1;
if (flag) return res;
}
return 0;
}
constexpr ModInt32() : v(0){};
~ModInt32() = default;
constexpr ModInt32(u32 _v) : v(reduce(u64(_v) * r2)) {}
constexpr ModInt32(i32 _v) : v(reduce(u64(_v % P + P) * r2)) {}
constexpr ModInt32(u64 _v) : v(reduce((_v % P) * r2)) {}
constexpr ModInt32(i64 _v) : v(reduce(u64(_v % P + P) * r2)) {}
constexpr ModInt32(const m32& rhs) : v(rhs.v) {}
constexpr u32 get() const { return norm(reduce(v)); }
explicit constexpr operator u32() const { return get(); }
explicit constexpr operator i32() const { return i32(get()); }
constexpr m32& operator=(const m32& rhs) { return v = rhs.v, *this; }
constexpr m32 operator-() const {
m32 res;
return res.v = (P << 1 & -(v != 0)) - v, res;
}
constexpr m32 inv() const { return pow(P - 2); }
constexpr m32& operator+=(const m32& rhs) {
return v += rhs.v - (P << 1), v += P << 1 & -(v >> 31), *this;
}
constexpr m32& operator-=(const m32& rhs) {
return v -= rhs.v, v += P << 1 & -(v >> 31), *this;
}
constexpr m32& operator*=(const m32& rhs) {
return v = reduce(u64(v) * rhs.v), *this;
}
constexpr m32& operator/=(const m32& rhs) {
return this->operator*=(rhs.inv());
}
friend m32 operator+(const m32& lhs, const m32& rhs) {
return m32(lhs) += rhs;
}
friend m32 operator-(const m32& lhs, const m32& rhs) {
return m32(lhs) -= rhs;
}
friend m32 operator*(const m32& lhs, const m32& rhs) {
return m32(lhs) *= rhs;
}
friend m32 operator/(const m32& lhs, const m32& rhs) {
return m32(lhs) /= rhs;
}
friend bool operator==(const m32& lhs, const m32& rhs) {
return norm(lhs.v) == norm(rhs.v);
}
friend bool operator!=(const m32& lhs, const m32& rhs) {
return norm(lhs.v) != norm(rhs.v);
}
friend std::istream& operator>>(std::istream& is, m32& rhs) {
return is >> rhs.v, rhs.v = reduce(u64(rhs.v) * r2), is;
}
friend std::ostream& operator<<(std::ostream& os, const m32& rhs) {
return os << rhs.get();
}
constexpr m32 pow(i64 y) const {
// assumes P is a prime
i64 rem = y % (P - 1);
if (y > 0 && rem == 0)
y = P - 1;
else
y = rem;
m32 res(1), x(*this);
for (; y != 0; y >>= 1, x *= x)
if (y & 1) res *= x;
return res;
}
};
using mint = ModInt32<nttmod>;
void ntt(vector<mint>& a, bool inverse) {
static array<mint, 30> dw{}, idw{};
if (dw[0] == 0) {
mint root = 2;
while (power(root, (nttmod - 1) / 2) == 1) root += 1;
for (int i = 0; i < 30; ++i)
dw[i] = -power(root, (nttmod - 1) >> (i + 2)),
idw[i] = 1 / dw[i];
}
int n = (int)a.size();
assert((n & (n - 1)) == 0);
if (not inverse) {
for (int m = n; m >>= 1;) {
mint w = 1;
for (int s = 0, k = 0; s < n; s += 2 * m) {
for (int i = s, j = s + m; i < s + m; ++i, ++j) {
auto x = a[i], y = a[j] * w;
a[i] = x + y;
a[j] = x - y;
}
w *= dw[__builtin_ctz(++k)];
}
}
} else {
for (int m = 1; m < n; m *= 2) {
mint w = 1;
for (int s = 0, k = 0; s < n; s += 2 * m) {
for (int i = s, j = s + m; i < s + m; ++i, ++j) {
auto x = a[i], y = a[j];
a[i] = x + y;
a[j] = (x - y) * w;
}
w *= idw[__builtin_ctz(++k)];
}
}
auto inv = 1 / mint(n);
for (auto&& e : a) e *= inv;
}
}
vector<mint> operator*(vector<mint> l, vector<mint> r) {
if (l.empty() or r.empty()) return {};
int n = (int)l.size(), m = (int)r.size(),
sz = 1 << __lg(2 * (n + m - 1) - 1);
if (min(n, m) < 30) {
vector<mint> res(n + m - 1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[i + j] += (l[i] * r[j]);
return {begin(res), end(res)};
}
bool eq = l == r;
l.resize(sz), ntt(l, false);
if (eq)
r = l;
else
r.resize(sz), ntt(r, false);
for (int i = 0; i < sz; ++i) l[i] *= r[i];
ntt(l, true), l.resize(n + m - 1);
return l;
}
vector<mint>& operator*=(vector<mint>& l, vector<mint> r) {
if (l.empty() or r.empty()) {
l.clear();
return l;
}
int n = (int)l.size(), m = (int)r.size(),
sz = 1 << __lg(2 * (n + m - 1) - 1);
if (min(n, m) < 30) {
vector<mint> res(n + m - 1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[i + j] += (l[i] * r[j]);
l = {begin(res), end(res)};
return l;
}
bool eq = l == r;
l.resize(sz), ntt(l, false);
if (eq)
r = l;
else
r.resize(sz), ntt(r, false);
for (int i = 0; i < sz; ++i) l[i] *= r[i];
ntt(l, true), l.resize(n + m - 1);
return l;
}
} // namespace ntt
using namespace ntt;
using vm = vector<mint2>;
// std::ostream& operator << (std::ostream& out, const mint2& rhs) {
// return out<<rhs.val();
// }
// Ref - https://www.codechef.com/viewsolution/41909444 Line 827 - 850.
const int maxnCr=2e6+5;
array<mint2,maxnCr+1> fac,inv;
mint2 nCr(lli n,lli r)
{
if(n<0||r<0||r>n)
return 0;
return fac[n]*inv[r]*inv[n-r];
}
void pre(lli n)
{
fac[0]=1;
for(int i=1;i<=n;++i)
fac[i]=i*fac[i-1];
inv[n]=1/fac[n];
for(int i=n;i>0;--i)
inv[i-1]=i*inv[i];
assert(inv[0]==mint2(1));
}
vm multiply(const vm &a,const vm &b){
vector<mint> aa,bb,cc;
for(auto x:a)
aa.pb(x.val());
for(auto x:b)
bb.pb(x.val());
cc=aa*bb;
vm c;
for(auto x:cc)
c.pb(x.get());
return c;
};
int main(void) {
ios_base::sync_with_stdio(false);cin.tie(NULL);
// freopen("txt.in", "r", stdin);
// freopen("txt.out", "w", stdout);
// cout<<std::fixed<<std::setprecision(35);
lli T;
cin>>T;
pre(maxnCr);
while(T--)
{
lli n,p,m;
cin>>n>>p>>m;
mint2 ans=0;
// mex <= m
vm a(m+1),b(m+1);
for(lli i=0;i<=m;i++)
a[i]=inv[i];
for(lli i=0;i<=m;i++){
b[i]=mint2(m-i).pow(n)*inv[i];
if(i&1)
b[i]*=-1;
}
const auto c = multiply(a,b);
for(lli t=0;t<=m;t++)
ans+=c[t]*fac[t]*(mint2(t).pow(p));
// dbg(mint(m+1).pow(p));
// mex is m+1
for(lli i=0;i<=m+1;i++){
mint2 cur=mint2(m+1).pow(p);
if(i&1)
cur*=-1;
cur*=nCr(m+1,i);
cur*=mint2(m+1-i).pow(n);
ans+=cur;
}
cout<<ans.val()<<endl;
}
return 0;
}
Editorialist's Solution
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <deque>
#include <forward_list>
#include <fstream>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;
using lint = long long;
using pint = pair<int, int>;
using plint = pair<lint, lint>;
struct fast_ios { fast_ios(){ cin.tie(nullptr), ios::sync_with_stdio(false), cout << fixed << setprecision(20); }; } fast_ios_;
#define ALL(x) (x).begin(), (x).end()
#define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i<i##_end_;i++)
#define IFOR(i, begin, end) for(int i=(end)-1,i##_begin_=(begin);i>=i##_begin_;i--)
#define REP(i, n) FOR(i,0,n)
#define IREP(i, n) IFOR(i,0,n)
template <typename T, typename V>
void ndarray(vector<T>& vec, const V& val, int len) { vec.assign(len, val); }
template <typename T, typename V, typename... Args> void ndarray(vector<T>& vec, const V& val, int len, Args... args) { vec.resize(len), for_each(begin(vec), end(vec), [&](T& v) { ndarray(v, val, args...); }); }
template <typename T> bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; }
template <typename T> bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; }
int floor_lg(long long x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); }
template <typename T1, typename T2> pair<T1, T2> operator+(const pair<T1, T2> &l, const pair<T1, T2> &r) { return make_pair(l.first + r.first, l.second + r.second); }
template <typename T1, typename T2> pair<T1, T2> operator-(const pair<T1, T2> &l, const pair<T1, T2> &r) { return make_pair(l.first - r.first, l.second - r.second); }
template <typename T> vector<T> sort_unique(vector<T> vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; }
template <typename T> int arglb(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); }
template <typename T> int argub(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); }
template <typename T> istream &operator>>(istream &is, vector<T> &vec) { for (auto &v : vec) is >> v; return is; }
template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <typename T, size_t sz> ostream &operator<<(ostream &os, const array<T, sz> &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; }
#if __cplusplus >= 201703L
template <typename... T> istream &operator>>(istream &is, tuple<T...> &tpl) { std::apply([&is](auto &&... args) { ((is >> args), ...);}, tpl); return is; }
template <typename... T> ostream &operator<<(ostream &os, const tuple<T...> &tpl) { os << '('; std::apply([&os](auto &&... args) { ((os << args << ','), ...);}, tpl); return os << ')'; }
#endif
template <typename T> ostream &operator<<(ostream &os, const deque<T> &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <typename T> ostream &operator<<(ostream &os, const set<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T, typename TH> ostream &operator<<(ostream &os, const unordered_set<T, TH> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T> ostream &operator<<(ostream &os, const multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T> ostream &operator<<(ostream &os, const unordered_multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T1, typename T2> ostream &operator<<(ostream &os, const pair<T1, T2> &pa) { os << '(' << pa.first << ',' << pa.second << ')'; return os; }
template <typename TK, typename TV> ostream &operator<<(ostream &os, const map<TK, TV> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
template <typename TK, typename TV, typename TH> ostream &operator<<(ostream &os, const unordered_map<TK, TV, TH> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
#ifdef HITONANODE_LOCAL
const string COLOR_RESET = "\033[0m", BRIGHT_GREEN = "\033[1;32m", BRIGHT_RED = "\033[1;31m", BRIGHT_CYAN = "\033[1;36m", NORMAL_CROSSED = "\033[0;9;37m", RED_BACKGROUND = "\033[1;41m", NORMAL_FAINT = "\033[0;2m";
#define dbg(x) cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << endl
#define dbgif(cond, x) ((cond) ? cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << endl : cerr)
#else
#define dbg(x) (x)
#define dbgif(cond, x) 0
#endif
template <int mod> struct ModInt {
#if __cplusplus >= 201402L
#define MDCONST constexpr
#else
#define MDCONST
#endif
using lint = long long;
MDCONST static int get_mod() { return mod; }
static int get_primitive_root() {
static int primitive_root = 0;
if (!primitive_root) {
primitive_root = [&]() {
std::set<int> fac;
int v = mod - 1;
for (lint i = 2; i * i <= v; i++)
while (v % i == 0) fac.insert(i), v /= i;
if (v > 1) fac.insert(v);
for (int g = 1; g < mod; g++) {
bool ok = true;
for (auto i : fac)
if (ModInt(g).pow((mod - 1) / i) == 1) {
ok = false;
break;
}
if (ok) return g;
}
return -1;
}();
}
return primitive_root;
}
int val;
MDCONST ModInt() : val(0) {}
MDCONST ModInt &_setval(lint v) { return val = (v >= mod ? v - mod : v), *this; }
MDCONST ModInt(lint v) { _setval(v % mod + mod); }
MDCONST explicit operator bool() const { return val != 0; }
MDCONST ModInt operator+(const ModInt &x) const { return ModInt()._setval((lint)val + x.val); }
MDCONST ModInt operator-(const ModInt &x) const { return ModInt()._setval((lint)val - x.val + mod); }
MDCONST ModInt operator*(const ModInt &x) const { return ModInt()._setval((lint)val * x.val % mod); }
MDCONST ModInt operator/(const ModInt &x) const { return ModInt()._setval((lint)val * x.inv() % mod); }
MDCONST ModInt operator-() const { return ModInt()._setval(mod - val); }
MDCONST ModInt &operator+=(const ModInt &x) { return *this = *this + x; }
MDCONST ModInt &operator-=(const ModInt &x) { return *this = *this - x; }
MDCONST ModInt &operator*=(const ModInt &x) { return *this = *this * x; }
MDCONST ModInt &operator/=(const ModInt &x) { return *this = *this / x; }
friend MDCONST ModInt operator+(lint a, const ModInt &x) { return ModInt()._setval(a % mod + x.val); }
friend MDCONST ModInt operator-(lint a, const ModInt &x) { return ModInt()._setval(a % mod - x.val + mod); }
friend MDCONST ModInt operator*(lint a, const ModInt &x) { return ModInt()._setval(a % mod * x.val % mod); }
friend MDCONST ModInt operator/(lint a, const ModInt &x) { return ModInt()._setval(a % mod * x.inv() % mod); }
MDCONST bool operator==(const ModInt &x) const { return val == x.val; }
MDCONST bool operator!=(const ModInt &x) const { return val != x.val; }
MDCONST bool operator<(const ModInt &x) const { return val < x.val; } // To use std::map<ModInt, T>
friend std::istream &operator>>(std::istream &is, ModInt &x) {
lint t;
return is >> t, x = ModInt(t), is;
}
MDCONST friend std::ostream &operator<<(std::ostream &os, const ModInt &x) { return os << x.val; }
MDCONST ModInt pow(lint n) const {
lint ans = 1, tmp = this->val;
while (n) {
if (n & 1) ans = ans * tmp % mod;
tmp = tmp * tmp % mod, n /= 2;
}
return ans;
}
static std::vector<long long> facs, facinv, invs;
MDCONST static void _precalculation(int N) {
int l0 = facs.size();
if (N <= l0) return;
facs.resize(N), facinv.resize(N); invs.resize(N);
for (int i = l0; i < N; i++) facs[i] = facs[i - 1] * i % mod;
facinv[N - 1] = ModInt(facs.back()).pow(mod - 2).val;
for (int i = N - 1; i >= l0; i--) {
invs[i] = facinv[i] * facs[i - 1] % mod;
if(i != 0)facinv[i - 1] = facinv[i] * i % mod;
}
}
MDCONST lint inv() const {
if (this->val < std::min(mod >> 1, 1 << 21)) {
while (this->val >= int(facs.size())) _precalculation(facs.size() * 2);
return invs[this->val];
} else {
return this->pow(mod - 2).val;
}
}
MDCONST ModInt fac() const {
while (this->val >= int(facs.size())) _precalculation(facs.size() * 2);
return facs[this->val];
}
MDCONST ModInt doublefac() const {
lint k = (this->val + 1) / 2;
return (this->val & 1) ? ModInt(k * 2).fac() / (ModInt(2).pow(k) * ModInt(k).fac())
: ModInt(k).fac() * ModInt(2).pow(k);
}
MDCONST ModInt nCr(const ModInt &r) const {
return (this->val < r.val) ? 0 : this->fac() / ((*this - r).fac() * r.fac());
}
ModInt sqrt() const {
if (val == 0) return 0;
if (mod == 2) return val;
if (pow((mod - 1) / 2) != 1) return 0;
ModInt b = 1;
while (b.pow((mod - 1) / 2) == 1) b += 1;
int e = 0, m = mod - 1;
while (m % 2 == 0) m >>= 1, e++;
ModInt x = pow((m - 1) / 2), y = (*this) * x * x;
x *= (*this);
ModInt z = b.pow(m);
while (y != 1) {
int j = 0;
ModInt t = y;
while (t != 1) j++, t *= t;
z = z.pow(1LL << (e - j - 1));
x *= z, z *= z, y *= z;
e = j;
}
return ModInt(std::min(x.val, mod - x.val));
}
};
template <int mod> std::vector<long long> ModInt<mod>::facs = {1};
template <int mod> std::vector<long long> ModInt<mod>::facinv = {0};
template <int mod> std::vector<long long> ModInt<mod>::invs = {0};
using mint = ModInt<754974721>;
template <typename MODINT> std::vector<MODINT> nttconv(std::vector<MODINT> a, std::vector<MODINT> b);
// Integer FFT (Fast Fourier Transform) for ModInt class
// (Also known as Number Theoretic Transform, NTT)
// is_inverse: inverse transform
// ** Input size must be 2^n **
template <typename MODINT> void ntt(std::vector<MODINT> &a, bool is_inverse = false) {
int n = a.size();
if (n == 1) return;
static const int mod = MODINT::get_mod();
static const MODINT root = MODINT::get_primitive_root();
assert(__builtin_popcount(n) == 1 and (mod - 1) % n == 0);
static std::vector<MODINT> w{1}, iw{1};
for (int m = w.size(); m < n / 2; m *= 2) {
MODINT dw = root.pow((mod - 1) / (4 * m)), dwinv = 1 / dw;
w.resize(m * 2), iw.resize(m * 2);
for (int i = 0; i < m; i++) w[m + i] = w[i] * dw, iw[m + i] = iw[i] * dwinv;
}
if (!is_inverse) {
for (int m = n; m >>= 1;) {
for (int s = 0, k = 0; s < n; s += 2 * m, k++) {
for (int i = s; i < s + m; i++) {
MODINT x = a[i], y = a[i + m] * w[k];
a[i] = x + y, a[i + m] = x - y;
}
}
}
} else {
for (int m = 1; m < n; m *= 2) {
for (int s = 0, k = 0; s < n; s += 2 * m, k++) {
for (int i = s; i < s + m; i++) {
MODINT x = a[i], y = a[i + m];
a[i] = x + y, a[i + m] = (x - y) * iw[k];
}
}
}
int n_inv = MODINT(n).inv();
for (auto &v : a) v *= n_inv;
}
}
template <typename MODINT> std::vector<MODINT> nttconv(std::vector<MODINT> a, std::vector<MODINT> b) {
int sz = 1, n = a.size(), m = b.size();
while (sz < n + m) sz <<= 1;
if (sz <= 16) {
std::vector<MODINT> ret(n + m - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) ret[i + j] += a[i] * b[j];
}
return ret;
}
a.resize(sz), b.resize(sz);
if (a == b) {
ntt(a, false);
b = a;
} else
ntt(a, false), ntt(b, false);
for (int i = 0; i < sz; i++) a[i] *= b[i];
ntt(a, true);
a.resize(n + m - 1);
return a;
}
int main()
{
ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);
mint::_precalculation(1000'10);
int t;
cin >> t;
int n, p, m;
while (t--)
{
cin >> n >> p >> m;
vector<mint> a(m + 1), b(m + 1);
mint m1 = m, mul = 1;
for(int i = 0; i <= m; i++){
a[i] = mint::facinv[i];
b[i] = a[i] * mul * m1.pow(n);
mul *= -1;
m1 = m1 - 1;
}
vector<mint> c = nttconv(a, b);
mint ans = 0;
for(int i = 0; i <= m; i++){
ans += mint(i).pow(p) * mint::facs[i] * c[i];
}
m1 = m + 1;
mint eans = 0, m2 = m + 1;
mul = 1;
for(int i = 0; i <= m + 1; i++){
eans += m2.nCr(i) * mul * m1.pow(n);
mul *= -1;
m1 = m1 - 1;
}
ans += eans * m2.pow(p);
cout << ans << endl;
}
return 0;
}