QOJ1285 Stirling Number

QOJ 传送门

因为 $x^{\overline p} \equiv x^p - x \pmod p$，所以设 $n = pq + r$，其中 $r \in [0, p - 1]$，则有：

$$\begin{aligned} x^{\overline n} & = (\prod\limits_{i = 0}^{q - 1} (x + ip)^{\overline p}) (x + pq)^{\overline r} \\ & = (x^p - x)^q (x + pq)^{\overline r} \\ & = (x^p - x)^q x^{\overline r} \\ & = x^q (x^{p - 1} - 1)^q x^{\overline r} \end{aligned}$$

设 $k - q = a(p - 1) + b$，其中 $b \in [0, p - 2]$，那么：

$$\begin{aligned} \begin{bmatrix} n \\ k \end{bmatrix} & = [x^k] x^{\overline n} \\ & = [x^k] x^q (x^{p - 1} - 1)^q x^{\overline r} \\ & = [x^{k - q}] (x^{p - 1} - 1)^q x^{\overline r} \\ & = (-1)^{q - a} \binom{q}{a} \begin{bmatrix} r \\ b \end{bmatrix} \end{aligned}$$

设 $m - q = k(p - 1) + t$，其中 $t \in [0, p - 2]$，那么：

$$\sum\limits_{i = 0}^m \begin{bmatrix} n \\ i \end{bmatrix} = \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} \sum\limits_{j = 0}^k (-1)^{q - j} \binom{q}{j} + \sum\limits_{i = t + 1}^r \begin{bmatrix} r \\ i \end{bmatrix} \sum\limits_{j = 0}^{k - 1} (-1)^{q - j} \binom{q}{j}$$

又因为：

$$\sum\limits_{i = 0}^m (-1)^i \binom{n}{i} = (-1)^m \binom{n - 1}{m}$$

这个式子由 $\binom{n}{m} = \binom{n - 1}{m} + \binom{n - 1}{m - 1}$ 易证。

那么：

$$\begin{aligned} \sum\limits_{i = 0}^m \begin{bmatrix} n \\ i \end{bmatrix} & = \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} (-1)^{q - k} \binom{q - 1}{k} + \sum\limits_{i = t + 1}^r \begin{bmatrix} r \\ i \end{bmatrix} (-1)^{q - k + 1} \binom{q - 1}{k - 1} \\ & = \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} (-1)^{q - k} \binom{q - 1}{k} + \sum\limits_{i = 0}^r \begin{bmatrix} r \\ i \end{bmatrix} (-1)^{q - k + 1} \binom{q - 1}{k - 1} + \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} (-1)^{q - k} \binom{q - 1}{k - 1} \\ & = (-1)^{q - k} (\binom{q}{k} \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} - \binom{q - 1}{k - 1} \sum\limits_{i = 0}^r \begin{bmatrix} r \\ i \end{bmatrix}) \\ & = (-1)^{q - k} (\binom{q}{k} \sum\limits_{i = 0}^t \begin{bmatrix} r \\ i \end{bmatrix} - \binom{q - 1}{k - 1} r!) \end{aligned}$$

那么剩下的问题是给定 $r, t$，求 $x^{\overline r}$ 的前 $t$ 次项系数和。

考虑一个多项式 $f(x)$，设模数的原根为 $g$，若我们求出 $f(g^0), f(g^1), \ldots, f(g^{p - 2})$，那么我们可以线性地求出 $f(x)$ 的前 $t$ 次项系数和。

具体地，考虑 $\sum\limits_{i = 0}^{p - 2} g^{ki}$，只有当 $k = 0$ 时这个式子等于 $p - 1 \equiv -1 \pmod p$，当 $k \ne 0$ 时由等比数列求和公式可知它等于 $0$。

那么考虑若要求 $f(x)$ 的第 $j$ 次项，$-\sum\limits_{i = 0}^{p - 2} g^{-ij} f(g^i)$ 即为答案。

所以：

$$\begin{aligned} \sum\limits_{i = 0}^t [x^i] x^{\overline r} & = \sum\limits_{i = 0}^t \sum\limits_{j = 0}^{p - 2} g^{-ij} f(g^j) \\ & = \sum\limits_{j = 0}^{p - 2} f(g^j) (\sum\limits_{i = 0}^t g^{-ij}) \end{aligned}$$

求 $f(g^i)$ 可以预处理阶乘及其逆元然后 $O(1)$ 单次计算，后面那个 $\sum\limits_{i = 0}^t g^{-ij}$ 可以用等比数列求和公式然后预处理原根的次方做到 $O(1)$ 单次计算。

总时间复杂度 $O(p)$。

code

#include <bits/stdc++.h>
#define pb emplace_back
#define fst first
#define scd second
#define mkp make_pair
#define mems(a, x) memset((a), (x), sizeof(a))
 
using namespace std;
typedef long long ll;
typedef double db;
typedef unsigned long long ull;
typedef long double ldb;
typedef pair<ll, ll> pii;
 
const int maxn = 1000100;
 
ll n, L, R, P, fac[maxn], ifac[maxn], N, G, q, r, inv[maxn], pw[maxn];
 
inline ll qpow(ll b, ll p) {
	ll res = 1;
	while (p) {
		if (p & 1) {
			res = res * b % P;
		}
		b = b * b % P;
		p >>= 1;
	}
	return res;
}
 
inline ll C(ll n, ll m) {
	if (n < m || n < 0 || m < 0) {
		return 0;
	}
	if (n >= P) {
		return C(n / P, m / P) * C(n % P, m % P) % P;
	} else {
		return fac[n] * ifac[m] % P * ifac[n - m] % P;
	}
}
 
inline bool check(ll x) {
	ll t = P - 1;
	vector<ll> vc;
	for (ll i = 2; i * i <= t; ++i) {
		if (t % i == 0) {
			vc.pb(i);
			while (t % i == 0) {
				t /= i;
			}
		}
	}
	if (t > 1) {
		vc.pb(t);
	}
	for (ll y : vc) {
		if (qpow(x, (P - 1) / y) == 1) {
			return 0;
		}
	}
	return 1;
}
 
inline ll g(ll n, ll m) {
	ll p = 1, ip = 1, iG = qpow(G, P - 2), res = 0;
	for (int i = 0; i <= P - 2; ++i, p = p * G % P, ip = ip * iG % P) {
		if (p + n - 1 >= P) {
			continue;
		}
		if (ip == 1) {
			res = (res + (m + 1) * fac[p + n - 1] % P * ifac[p - 1]) % P;
		} else {
			res = (res + (pw[((-(m + 1) * i) % (P - 1) + P - 1) % (P - 1)] + P - 1) % P * inv[(ip + P - 1) % P] % P * fac[p + n - 1] % P * ifac[p - 1]) % P;
		}
	}
	return (P - res) % P;
}
 
inline ll h(ll m) {
	if (m < q) {
		return 0;
	}
	ll k = (m - q) / (P - 1);
	ll t = (m - q) % (P - 1);
	ll ans = C(q, k) * g(r, t) % P;
	ans = (ans - C(q - 1, k - 1) * fac[r] % P + P) % P;
	if ((q - k) & 1) {
		ans = (P - ans) % P;
	}
	return ans;
}
 
void solve() {
	scanf("%lld%lld%lld%lld", &n, &L, &R, &P);
	N = P - 1;
	fac[0] = 1;
	for (int i = 1; i <= N; ++i) {
		fac[i] = fac[i - 1] * i % P;
	}
	ifac[N] = qpow(fac[N], P - 2);
	for (int i = N - 1; ~i; --i) {
		ifac[i] = ifac[i + 1] * (i + 1) % P;
	}
	inv[1] = 1;
	for (int i = 2; i <= N; ++i) {
		inv[i] = (P - P / i) * inv[P % i] % P;
	}
	for (G = 1;; ++G) {
		if (check(G)) {
			break;
		}
	}
	pw[0] = 1;
	for (int i = 1; i <= P - 2; ++i) {
		pw[i] = pw[i - 1] * G % P;
	}
	q = n / P;
	r = n - P * q;
	printf("%lld\n", (h(R) - h(L - 1) + P) % P);
}
 
int main() {
	int T = 1;
	// scanf("%d", &T);
	while (T--) {
		solve();
	}
	return 0;
}