[TJOI2018]教科书般的亵渎

最新推荐文章于 2022-04-29 00:17:33 发布

__zyn__

最新推荐文章于 2022-04-29 00:17:33 发布

阅读量217

点赞数

文章标签：高精度计算拉格朗日插值多项式复杂度分析算法实现

本文链接：https://blog.csdn.net/zyn_615/article/details/113834234

版权

Statement

statement

Solution

经过推导发现答案为:
$\sum_{i=1}^{m+1}\left[\sum_{j=1}^{n-a_{i-1}}j^{m+1}-\sum_{j=i}^m(a_j-a_{i-1})^{m+1}\right]$
复杂度瓶颈在 $\sum_{i=1}^ni^k$

对以上和式用扰动法

记 $S(k,n)=\sum_{i=1}^ni^k$

$S(k+1,n)=\sum_{i=1}^ni^{k+1}$
$S(k+1,n)+(n+1)^{k+1}=\sum_{i=1}^{n}(i+1)^{k+1}+1$
$S(k+1,n)+(n+1)^{k+1}=\sum_{i=1}^n\sum_{j=0}^{k+1}\dbinom{k+1}{j}i^j+1$
$S(k+1,n)+(n+1)^{k+1}=\sum_{i=0}^{k+1}\dbinom{k+1}{i}S(i,n)+1$
$S(k+1,n)+(n+1)^{k+1}=S(k+1,n)+1+\dbinom{k+1}{k}S(k,n)+\sum_{i=0}^{k-1}\dbinom{k+1}{i}S(i,n)$
$S(k,n)=\dfrac{(n+1)^{k+1}-1-\sum\limits_{i=0}^{k-1}\dbinom{k+1}{i}S(i,n)}{\dbinom{k+1}{k}}$
发现其为一个 $k + 1$ 次多项式，拉格朗日插值即可，时间复杂度为 $\mathcal {O}(Tm^3)$ 。

$f(x)=\sum_{i=1}^{m+3}y_i\prod_{j\not ={i}}\dfrac{x-x_j}{x_i-x_j}$

Code

#include<bits/stdc++.h>
using namespace std;

const int MAXN = 55;
const int P = 1E9 + 7;
int A[MAXN];

inline int qPow(int a, int b) {
	int res = 1;
	for (; b; b >>= 1, a = 1ll * a * a % P)
		(b & 1) && (res = 1ll * res * a % P);
	return res;
}

inline int calc(long long n, int k) {
	if (n <= k + 2)
		return A[n];
	
	int res = 0;
	for (int i = 1;i <= k + 2; ++i) {
		int u = 1, v = 1;
		for (int j = 1;j <= k + 2; ++j)
			if (i != j)
				u = 1ll * u * ((n - j) % P + P) % P,
				v = 1ll * v * (i - j + P) % P;
		u = 1ll * u * qPow(v, P - 2) % P * A[i] % P,
		res = (res + u) % P;
	}
	return res;
	
}

int main(void) {
	ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
	
	int testCase = 0;
	cin >> testCase;
	
	while (testCase--) {
		long long n;
		int m = 0, ans = 0;
		cin >> n >> m;
		vector <long long> v(m);
		for (long long& a : v)
			cin >> a;
		sort(v.begin(), v.end());
		
		for (int i = 1;i <= m + 3; ++i)
			A[i] = (A[i - 1] + qPow(i, m + 1)) % P;
		
		for (int i = 0, o = 0;i <= m; ++i) {
			ans = (ans + calc(n - (i ? v[i - 1] : 0), m + 1)) % P, o = 0;
			for (int j = i;j < m; ++j)
				o = (o + qPow(v[j] - (i ? v[i - 1] : 0), m + 1)) % P;
			ans = (ans - o + P) % P;
		}
		
		cout << ans << '\n';
	}
	
	return 0;
}