stirling_small_prime.hpp¶
在较小素数模下查询下标可达 \(10^{18}\) 的第一类或第二类 Stirling 数。
\[
\displaystyle x^n=\sum_{k=0}^n S(n,k)x^{\underline{k}}
\]
Complexity: Time: O(p^2) preprocessing and O(log_p n) per query. Space: O(p^2) for each requested Stirling kind.
AC 记录：stirling_number_of_the_first_kind_small_p_large_n, stirling_number_of_the_second_kind_small_p_large_n。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(p^2) preprocessing and O(log_p n) per query.
/// Space: O(p^2) for each requested Stirling kind.
#include <cassert>
#include <cstdint>
#include <vector>
namespace noya {
/// @brief Query signed first-kind and second-kind Stirling numbers modulo a
/// small fixed prime for indices as large as 10^18. Over F_p, the identities
/// (x)_p=x^p-x and x^p=x split a large falling factorial or power into p-sized
/// blocks. Repeating the split leaves one small Stirling number with indices
/// below p and one binomial coefficient. Lucas' theorem evaluates the latter
/// digit by digit. The constructor can omit either Stirling table when only one
/// kind is needed.
class stirling_small_prime_table {
public:
explicit stirling_small_prime_table(std::uint32_t p, bool b1 = true,
bool b2 = true)
: p_(p) {
assert(p >= 2);
build_binomial();
if (b1) {
build_first();
}
if (b2) {
build_second();
}
}
std::uint32_t first_kind(std::uint64_t n, std::uint64_t k) const {
assert(!a_.empty());
if (k > n) {
return 0;
}
std::uint64_t bs = n / p_;
int rem = int(n % p_);
if (bs > k) {
return 0;
}
std::uint64_t quo = (k - bs) / (p_ - 1);
int sk = int((k - bs) % (p_ - 1));
if (sk == 0 && rem > 0) {
if (quo == 0) {
return 0;
}
sk = int(p_ - 1);
quo--;
}
if (bs < quo || sk > rem) {
return 0;
}
std::uint32_t res = multiply(lucas_choose(bs, quo), a_[rem][sk]);
if (((bs + quo) & 1) != 0 && res != 0) {
res = p_ - res;
}
return res;
}
std::uint32_t second_kind(std::uint64_t n, std::uint64_t k) const {
assert(!b_.empty());
if (k > n) {
return 0;
}
if (n == 0) {
return k == 0 ? 1 : 0;
}
std::uint64_t bs = k / p_;
int sk = int(k % p_);
if (n < bs) {
return 0;
}
std::uint64_t quo = (n - bs) / (p_ - 1);
int rem = int((n - bs) % (p_ - 1));
if (rem == 0) {
if (quo == 0) {
return 0;
}
rem = int(p_ - 1);
quo--;
}
if (sk > rem) {
return 0;
}
if (rem < int(p_ - 1)) {
return multiply(lucas_choose(quo, bs), b_[rem][sk]);
}
if (sk == 0) {
return bs == 0 ? 0 : lucas_choose(quo, bs - 1);
}
return multiply(lucas_choose(quo, bs), b_[p_ - 1][sk]);
}
private:
std::uint32_t p_;
std::vector<std::vector<std::uint32_t>> bin;
std::vector<std::vector<std::uint32_t>> a_;
std::vector<std::vector<std::uint32_t>> b_;
std::uint32_t multiply(std::uint64_t a, std::uint64_t b) const {
return std::uint32_t(a * b % p_);
}
std::uint32_t lucas_choose(std::uint64_t n, std::uint64_t k) const {
if (k > n) {
return 0;
}
std::uint32_t res = 1;
while (n > 0 || k > 0) {
int hi = int(n % p_);
int lo = int(k % p_);
if (lo > hi) {
return 0;
}
res = multiply(res, bin[hi][lo]);
n /= p_;
k /= p_;
}
return res;
}
void build_binomial() {
bin.resize(p_);
bin[0] = {1};
for (int n = 1; n < int(p_); n++) {
bin[n] = bin[n - 1];
bin[n].push_back(0);
for (int k = 1; k <= n; k++) {
bin[n][k] += bin[n - 1][k - 1];
if (bin[n][k] >= p_) {
bin[n][k] -= p_;
}
}
}
}
void build_first() {
a_.resize(p_);
a_[0] = {1};
for (int n = 1; n < int(p_); n++) {
a_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
a_[n][k] += a_[n - 1][k - 1];
}
if (k < n) {
a_[n][k] += multiply(a_[n - 1][k], p_ - n + 1);
}
if (a_[n][k] >= p_) {
a_[n][k] -= p_;
}
}
}
}
void build_second() {
b_.resize(p_);
b_[0] = {1};
for (int n = 1; n < int(p_); n++) {
b_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
b_[n][k] += b_[n - 1][k - 1];
}
if (k < n) {
b_[n][k] += multiply(b_[n - 1][k], k);
}
if (b_[n][k] >= p_) {
b_[n][k] -= p_;
}
}
}
}
};
} // namespace noya
#ifndef NOYA_STIRLING_SMALL_PRIME_HPP
#define NOYA_STIRLING_SMALL_PRIME_HPP 1
/// @complexity Time: O(p^2) preprocessing and O(log_p n) per query.
/// Space: O(p^2) for each requested Stirling kind.
#include <cassert>
#include <cstdint>
#include <vector>
namespace noya {
/// @brief Query signed first-kind and second-kind Stirling numbers modulo a
/// small fixed prime for indices as large as 10^18. Over F_p, the identities
/// (x)_p=x^p-x and x^p=x split a large falling factorial or power into p-sized
/// blocks. Repeating the split leaves one small Stirling number with indices
/// below p and one binomial coefficient. Lucas' theorem evaluates the latter
/// digit by digit. The constructor can omit either Stirling table when only one
/// kind is needed.
class stirling_small_prime_table {
public:
explicit stirling_small_prime_table(std::uint32_t p, bool b1 = true,
bool b2 = true)
: p_(p) {
assert(p >= 2);
build_binomial();
if (b1) {
build_first();
}
if (b2) {
build_second();
}
}
std::uint32_t first_kind(std::uint64_t n, std::uint64_t k) const {
assert(!a_.empty());
if (k > n) {
return 0;
}
std::uint64_t bs = n / p_;
int rem = int(n % p_);
if (bs > k) {
return 0;
}
std::uint64_t quo = (k - bs) / (p_ - 1);
int sk = int((k - bs) % (p_ - 1));
if (sk == 0 && rem > 0) {
if (quo == 0) {
return 0;
}
sk = int(p_ - 1);
quo--;
}
if (bs < quo || sk > rem) {
return 0;
}
std::uint32_t res = multiply(lucas_choose(bs, quo), a_[rem][sk]);
if (((bs + quo) & 1) != 0 && res != 0) {
res = p_ - res;
}
return res;
}
std::uint32_t second_kind(std::uint64_t n, std::uint64_t k) const {
assert(!b_.empty());
if (k > n) {
return 0;
}
if (n == 0) {
return k == 0 ? 1 : 0;
}
std::uint64_t bs = k / p_;
int sk = int(k % p_);
if (n < bs) {
return 0;
}
std::uint64_t quo = (n - bs) / (p_ - 1);
int rem = int((n - bs) % (p_ - 1));
if (rem == 0) {
if (quo == 0) {
return 0;
}
rem = int(p_ - 1);
quo--;
}
if (sk > rem) {
return 0;
}
if (rem < int(p_ - 1)) {
return multiply(lucas_choose(quo, bs), b_[rem][sk]);
}
if (sk == 0) {
return bs == 0 ? 0 : lucas_choose(quo, bs - 1);
}
return multiply(lucas_choose(quo, bs), b_[p_ - 1][sk]);
}
private:
std::uint32_t p_;
std::vector<std::vector<std::uint32_t>> bin;
std::vector<std::vector<std::uint32_t>> a_;
std::vector<std::vector<std::uint32_t>> b_;
std::uint32_t multiply(std::uint64_t a, std::uint64_t b) const {
return std::uint32_t(a * b % p_);
}
std::uint32_t lucas_choose(std::uint64_t n, std::uint64_t k) const {
if (k > n) {
return 0;
}
std::uint32_t res = 1;
while (n > 0 || k > 0) {
int hi = int(n % p_);
int lo = int(k % p_);
if (lo > hi) {
return 0;
}
res = multiply(res, bin[hi][lo]);
n /= p_;
k /= p_;
}
return res;
}
void build_binomial() {
bin.resize(p_);
bin[0] = {1};
for (int n = 1; n < int(p_); n++) {
bin[n] = bin[n - 1];
bin[n].push_back(0);
for (int k = 1; k <= n; k++) {
bin[n][k] += bin[n - 1][k - 1];
if (bin[n][k] >= p_) {
bin[n][k] -= p_;
}
}
}
}
void build_first() {
a_.resize(p_);
a_[0] = {1};
for (int n = 1; n < int(p_); n++) {
a_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
a_[n][k] += a_[n - 1][k - 1];
}
if (k < n) {
a_[n][k] += multiply(a_[n - 1][k], p_ - n + 1);
}
if (a_[n][k] >= p_) {
a_[n][k] -= p_;
}
}
}
}
void build_second() {
b_.resize(p_);
b_[0] = {1};
for (int n = 1; n < int(p_); n++) {
b_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
b_[n][k] += b_[n - 1][k - 1];
}
if (k < n) {
b_[n][k] += multiply(b_[n - 1][k], k);
}
if (b_[n][k] >= p_) {
b_[n][k] -= p_;
}
}
}
}
};
} // namespace noya
#endif // NOYA_STIRLING_SMALL_PRIME_HPP
#include <cassert>
#include <cstdint>
#include <vector>
/// @complexity Time: O(p^2) preprocessing and O(log_p n) per query.
/// Space: O(p^2) for each requested Stirling kind.
namespace noya {
/// @brief Query signed first-kind and second-kind Stirling numbers modulo a
/// small fixed prime for indices as large as 10^18. Over F_p, the identities
/// (x)_p=x^p-x and x^p=x split a large falling factorial or power into p-sized
/// blocks. Repeating the split leaves one small Stirling number with indices
/// below p and one binomial coefficient. Lucas' theorem evaluates the latter
/// digit by digit. The constructor can omit either Stirling table when only one
/// kind is needed.
class stirling_small_prime_table {
public:
explicit stirling_small_prime_table(std::uint32_t p, bool b1 = true,
bool b2 = true)
: p_(p) {
assert(p >= 2);
build_binomial();
if (b1) {
build_first();
}
if (b2) {
build_second();
}
}
std::uint32_t first_kind(std::uint64_t n, std::uint64_t k) const {
assert(!a_.empty());
if (k > n) {
return 0;
}
std::uint64_t bs = n / p_;
int rem = int(n % p_);
if (bs > k) {
return 0;
}
std::uint64_t quo = (k - bs) / (p_ - 1);
int sk = int((k - bs) % (p_ - 1));
if (sk == 0 && rem > 0) {
if (quo == 0) {
return 0;
}
sk = int(p_ - 1);
quo--;
}
if (bs < quo || sk > rem) {
return 0;
}
std::uint32_t res = multiply(lucas_choose(bs, quo), a_[rem][sk]);
if (((bs + quo) & 1) != 0 && res != 0) {
res = p_ - res;
}
return res;
}
std::uint32_t second_kind(std::uint64_t n, std::uint64_t k) const {
assert(!b_.empty());
if (k > n) {
return 0;
}
if (n == 0) {
return k == 0 ? 1 : 0;
}
std::uint64_t bs = k / p_;
int sk = int(k % p_);
if (n < bs) {
return 0;
}
std::uint64_t quo = (n - bs) / (p_ - 1);
int rem = int((n - bs) % (p_ - 1));
if (rem == 0) {
if (quo == 0) {
return 0;
}
rem = int(p_ - 1);
quo--;
}
if (sk > rem) {
return 0;
}
if (rem < int(p_ - 1)) {
return multiply(lucas_choose(quo, bs), b_[rem][sk]);
}
if (sk == 0) {
return bs == 0 ? 0 : lucas_choose(quo, bs - 1);
}
return multiply(lucas_choose(quo, bs), b_[p_ - 1][sk]);
}
private:
std::uint32_t p_;
std::vector<std::vector<std::uint32_t>> bin;
std::vector<std::vector<std::uint32_t>> a_;
std::vector<std::vector<std::uint32_t>> b_;
std::uint32_t multiply(std::uint64_t a, std::uint64_t b) const {
return std::uint32_t(a * b % p_);
}
std::uint32_t lucas_choose(std::uint64_t n, std::uint64_t k) const {
if (k > n) {
return 0;
}
std::uint32_t res = 1;
while (n > 0 || k > 0) {
int hi = int(n % p_);
int lo = int(k % p_);
if (lo > hi) {
return 0;
}
res = multiply(res, bin[hi][lo]);
n /= p_;
k /= p_;
}
return res;
}
void build_binomial() {
bin.resize(p_);
bin[0] = {1};
for (int n = 1; n < int(p_); n++) {
bin[n] = bin[n - 1];
bin[n].push_back(0);
for (int k = 1; k <= n; k++) {
bin[n][k] += bin[n - 1][k - 1];
if (bin[n][k] >= p_) {
bin[n][k] -= p_;
}
}
}
}
void build_first() {
a_.resize(p_);
a_[0] = {1};
for (int n = 1; n < int(p_); n++) {
a_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
a_[n][k] += a_[n - 1][k - 1];
}
if (k < n) {
a_[n][k] += multiply(a_[n - 1][k], p_ - n + 1);
}
if (a_[n][k] >= p_) {
a_[n][k] -= p_;
}
}
}
}
void build_second() {
b_.resize(p_);
b_[0] = {1};
for (int n = 1; n < int(p_); n++) {
b_[n].assign(n + 1, 0);
for (int k = 0; k <= n; k++) {
if (k > 0) {
b_[n][k] += b_[n - 1][k - 1];
}
if (k < n) {
b_[n][k] += multiply(b_[n - 1][k], k);
}
if (b_[n][k] >= p_) {
b_[n][k] -= p_;
}
}
}
}
};
} // namespace noya