「loj - 6055」「from CommonAnts」一道数学题 加强版


description

已知函数 \(f(k, x)(k, x\in \mathbb{N}_+)\) 满足:

\[f(k, x) = \begin{cases} 1 & x = 1\\ \sum_{i=1}^{x-1}f(k, i) + x^k & x > 1 \end{cases} \]

现在给定 \(n, k\),求 \(f(k, n)\)

由于答案很大,你只需计算答案对 \(10^9 + 7\) 取模后的值。

problem link。

solution

考虑每个 \(x^k\) 的贡献,可以得到 \(f(k, n) = n^k + \sum_{i=0}^{n-1}2^{n-i-1}i^k=n^k + 2^{n-1}\sum_{i=0}^{n-1}\frac{i^k}{2^i}\)

\(a = \frac{1}{2}\),尝试求 \(F_k(n) = \sum_{i=0}^{n-1}a^ii^k\)

通过转下降幂构造错位相减。不妨记 \(S_k(n) =\sum_{i=0}^{n-1}a^{i}i^{\underline k}\),则 \(F_k(n) = \sum_{j=0}^{k}{k\brace j}S_{j}(n)\)

\[\begin{aligned} S_k(n) &=\sum_{i=0}^{n-1}a^{i}i^{\underline k}\\ aS_k(n) &= \sum_{i=0}^{n-1}a^{i+1}i^{\underline k} = \sum_{i=1}^{n}a^i(i-1)^{\underline k} \\ (1-a)S_k(n) &= \sum_{i=1}^{n}a^i(i^{\underline k} - (i-1)^{\underline k})-a^nn^{\underline k} + a^{0}0^{\underline k}\\ &= ak\sum_{i=0}^{n-1}a^ii^{\underline {k-1}}-a^nn^{\underline k}\\ S_k(n) &= \frac{ak}{1-a}S_{k-1}(n)-\frac{a^nn^{\underline k}}{1-a}\\ \end{aligned} \]

\(k = 0\) 时,有 \(S_0(n) = \sum_{i=0}^{n-1} a^i = \frac{1-a^n}{1-a}\)


如果硬上任意模数 fft 求斯特林数是过不了的,考虑进一步优化。

考虑 \(S_{k}(n)\) 的形式,发现存在最高次数为 \(k\) 的多项式 \(G_k(x)=\sum g_ix^i\) 满足 \(S_k(n)=G_k(n)\times a^n-G_k(0)\)

进一步地,存在最高次数为 \(k\) 的多项式 \(G'_k(x)\) 满足 \(F_k(n)=G'_k(n)\times a^n-G'_k(0)\)

可以处理出 \(G'_k(i) = \alpha_i\times G_{k}'(0)+\beta_i\) 中的 \(\alpha_i, \beta_i\)。再利用如下的 \(k+1\) 阶差分公式:

\[\Delta^{k+1}G'(0)=\sum_{i=0}^{k+1}\binom{k+1}{i}(-1)^{k+1-i}G'(i)=0 \]

解出 \(G'_k(0)\),再拉格朗日插值就可以求 \(G'_k(n)\)

code

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int MAXN = 1000000;
const int MOD = int(1E9) + 7;
const int INV2 = (MOD + 1) >> 1;

inline int add(int x, int y) {x += y; return x >= MOD ? x - MOD : x;}
inline int sub(int x, int y) {x -= y; return x < 0 ? x + MOD : x;}
inline int mul(int x, int y) {return (int) (1LL * x * y % MOD);}
int mpow(int b, int p) {
	int ret;
	for(ret=1;p;p>>=1,b=mul(b,b))
		if( p & 1 ) ret = mul(ret, b);
	return ret;
}

char s[MAXN + 5]; int len, n, p, k;
void read() {
	scanf("%s%d", s, &k), len = strlen(s);
	for(int i=0;i<len;i++) {
		n = (10LL*n + s[i] - '0') % (MOD);
		p = (10LL*p + s[i] - '0') % (MOD - 1);
	}
}

int prm[MAXN + 5], pcnt; bool nprm[MAXN + 5];
int fct[MAXN + 5], ifct[MAXN + 5], pwk[MAXN + 5];
void init() {
	fct[0] = 1; for(int i=1;i<=k+3;i++) fct[i] = mul(fct[i - 1], i);
	ifct[k+3] = mpow(fct[k+3], MOD - 2);
	for(int i=k+2;i>=0;i--) ifct[i] = mul(ifct[i + 1], i + 1);
	
	pwk[1] = 1;
	for(int i=2;i<=k+3;i++) {
		if( !nprm[i] ) prm[++pcnt] = i, pwk[i] = mpow(i, k);
		for(int j=1;i*prm[j]<=k+3;j++) {
			nprm[i*prm[j]] = false;
			pwk[i*prm[j]] = mul(pwk[i], pwk[prm[j]]);
			if( i % prm[j] == 0 ) break;
		}
	}
}

int lf[MAXN + 5], rf[MAXN + 5];
int gety(int x, int k, int *y) {
	lf[0] = 1; for(int i=1;i<=k;i++) lf[i] = mul(lf[i - 1], sub(x, i));
	rf[k + 1] = 1; for(int i=k;i>=1;i--) rf[i] = mul(rf[i + 1], sub(x, i));
	
	int ret = 0;
	for(int i=1;i<=k;i++) {
		int coef = mul(mul(lf[i - 1], rf[i + 1]), mul(ifct[i - 1], ifct[k - i]));
		if( (k - i) & 1 ) ret = sub(ret, mul(coef, y[i]));
		else ret = add(ret, mul(coef, y[i]));
	}
	return ret;
}

int a[MAXN + 5], b[MAXN + 5], g[MAXN + 5];
int main() {
	read(), init();
	
	a[0] = 1;
	for(int i=1;i<=k+1;i++)
		a[i] = mul(2, a[i - 1]), b[i] = mul(2, add(b[i - 1], pwk[i - 1]));
	
	int A = 0, B = 0;
	for(int i=0;i<=k+1;i++) {
		int del = mul(ifct[i], ifct[k + 1 - i]);
		if( (k + 1 - i) & 1 )
			A = sub(A, mul(del, a[i])), B = sub(B, mul(del, b[i]));
		else A = add(A, mul(del, a[i])), B = add(B, mul(del, b[i]));
	}
		
	g[0] = mul(mpow(A, MOD - 2), MOD - B);
	for(int i=1;i<=k+1;i++) g[i] = add(mul(a[i], g[0]), b[i]);
	
	int ans = sub(mul(mpow(INV2, p), gety(n, k + 1, g)), g[0]);
	printf("%d\n", add(mul(mpow(2, p+(MOD-1)-1), ans), mpow(n, k)));
}

details

lca给的solution 是根据递推式直接构造了一个多项式,然后转成了常系数线性递推(虽然最后算法的形式差不多)。

不过感觉 \(\sum_{i=0}^{n-1}a_ii^k\) 这个形式更普遍一些(不清楚我乱说的)?

posted @ 2020-07-02 21:44  Tiw_Air_OAO  阅读(134)  评论(0编辑  收藏
把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end