51nod 1822 序列求和 V5

题面

组合数没了,怎么说嘛

题解

(r = 0)(r = 1) 时是 trivial 的,下面假设 (r eq 0, 1)

(S_k(n) = sum_{i=1}^n i^k r^i)

那么有

[egin{aligned} S_k(n) &= sum_{i=1}^n i^k r^i\ &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{i=1}^n r^i((i - 1)^k - i^k) ight)\ &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{i=1}^n r^i sum_{j=0}^{k - 1} (-1)^{k-j} i^j inom kj ight)\ &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj sum_{i=1}^n r^i i^j ight)\ &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj S_j(n) ight) end{aligned} ]

不过好像只能 (mathcal O(k^2)) 算。

但是如果将 (S_k(n)) 的封闭形式展开,也许就可以发现下面这个规律:

引理 1 存在一个 (k) 次多项式 (F_k(x)) 满足 (S_k(n) = r^{n} F_k(n) - F_k(0))

(k = 0) 时显然成立。

(k > 0) 时:

[egin{aligned} S_k(n) &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj S_j(n) ight)\ &= frac 1{r - 1} left( n^kr^{n + 1} + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj (r^n F_j(n) - F_j(0)) ight)\ &= r^n left(frac{r n^k + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj F_j(n)}{r-1} ight) - left(frac{ sum_{j=0}^{k - 1} (-1)^{k-j} inom kj F_j(0)}{r-1} ight)\ end{aligned} ]

也就是说,只要 (F_k(n) = frac{r n^k + sum_{j=0}^{k - 1} (-1)^{k-j} inom kj F_j(n)}{r-1}) 就是合法的,这显然是一个 (k) 次多项式。

考虑如何计算 (F_k(n)),首先由定义立即可以得出:

[F_k(n+1)=frac{F_k(n)}r + (n + 1)^k ]

我们只需要知道 (F_k(0)) 就可以插值出答案了。

引理 2 对任意一个 (n) 次多项式 (f(x)),满足如下恒等式:

[sum_{i=0}^{n+1} (-1)^i inom {n + 1} i f(i) = 0 ]

定义平移算子 (mathrm E f(x) = f(x + 1)),恒等算子 (1 f(x) = f(x)),差分算子 (Delta f(x) = f(x + 1) - f(x))

那么有 (Delta^n = (mathrm E - 1)^n = sum_{i=0}^n (-1)^i inom ni mathrm E^i 1^{n-i})

所以 (Delta^{n+1} f(0) = sum_{i=0}^{n+1} (-1)^i inom {n+1}i f(i)),由于 (f)(n) 次多项式,所以 (Delta^{n+1} f(0) = 0),得证。

由于所有的 (F_k(n)) 都可以用 (F_k(0)) 表示,所以就可以得出一个和 (F_k(0)) 有关的方程,解出来即可。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)

template <typename T> inline T Read()
{
	T data = 0, w = 1; char ch = getchar();
	while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar();
	if (ch == '-') w = -1, ch = getchar();
	while (ch >= '0' && ch <= '9') data = data * 10 + (ch ^ 48), ch = getchar();
	return data * w;
}

#define read Read<int>
#define readll Read<ll>

using ll = long long;
const int N(2e5 + 10), Mod(985661441);
inline int upd(const int &x) { return x + (x >> 31 & Mod); }
int fastpow(int x, int y)
{
	int ans = 1;
	for (; y; y >>= 1, x = 1ll * x * x % Mod)
		if (y & 1) ans = 1ll * ans * x % Mod;
	return ans;
}

int T, K, R, prime[N], not_prime[N], cnt, Pow[N], f[N], fac[N], inv[N];
long long n;

inline int C(int n, int m) { return 1ll * fac[n] * inv[m] % Mod * inv[n - m] % Mod; }
int Lagrange(long long x, int m)
{
	if (x <= m) return f[x];
	static int g[N], h[N], y, ans; g[0] = h[m + 1] = 1, y = x % Mod, ans = 0;
	for (int i = 1; i <= m; i++) g[i] = 1ll * g[i - 1] * upd(y - i) % Mod;
	for (int i = m; i; i--) h[i] = 1ll * h[i + 1] * upd(y - i) % Mod;
	for (int i = 1, op = (m & 1) ? 1 : Mod - 1; i <= m; op = Mod - op, i++)
		ans = (ans + 1ll * op * f[i] % Mod * inv[i - 1] % Mod * inv[m - i] % Mod * g[i - 1] % Mod * h[i + 1]) % Mod;
	return ans;
}

int main()
{
#ifndef ONLINE_JUDGE
	file(cpp);
#endif
	K = 2e5 + 9, fac[0] = 1;
	for (int i = 1; i <= K; i++) fac[i] = 1ll * fac[i - 1] * i % Mod;
	inv[K] = fastpow(fac[K], Mod - 2);
	for (int i = K; i; i--) inv[i - 1] = 1ll * inv[i] * i % Mod;
	for (T = read(); T--; )
	{
		K = read(), n = readll(), R = readll() % Mod;
		Pow[1] = not_prime[1] = 1, cnt = f[0] = 0;
		for (int i = 2; i <= K + 2; i++)
		{
			if (!not_prime[i]) prime[++cnt] = i, Pow[i] = fastpow(i, K);
			for (int j = 1; j <= cnt && i * prime[j] <= K + 2; j++)
			{
				not_prime[i * prime[j]] = 1, Pow[i * prime[j]] = 1ll * Pow[i] * Pow[prime[j]] % Mod;
				if (i % prime[j] == 0) break;
			}
		}
		if (R == 0) { puts("0"); continue; }
		if (R == 1)
		{
			for (int i = 1; i <= K + 2; i++) f[i] = upd(f[i - 1] + Pow[i] - Mod);
			printf("%d
", Lagrange(n, K + 2)); continue;
		}
		int k = 1, b = 0, sk = 1, sb = 0, invR = fastpow(R, Mod - 2);
		for (int i = 1; i <= K + 1; i++)
		{
			k = 1ll * k * invR % Mod, b = (1ll * b * invR + Pow[i]) % Mod;
			if (i & 1) sk = (sk - 1ll * C(K + 1, i) * k) % Mod, sb = (sb - 1ll * C(K + 1, i) * b) % Mod;
			else sk = (sk + 1ll * C(K + 1, i) * k) % Mod, sb = (sb + 1ll * C(K + 1, i) * b) % Mod;
		}
		f[0] = 1ll * sb * fastpow(sk, Mod - 2) % Mod, f[0] = upd(-f[0]);
		for (int i = 1; i <= K + 1; i++) f[i] = (1ll * f[i - 1] * invR + Pow[i]) % Mod;
		printf("%d
", upd(1ll * Lagrange(n, K + 1) * fastpow(R, n % (Mod - 1)) % Mod - f[0]));
	}
	return 0;
}
原文地址:https://www.cnblogs.com/cj-xxz/p/14430202.html