2023 USP Try-outs-A-Lagrange插值

2023 USP Try-outs-A-Lagrange插值

题目链接:https://codeforces.com/gym/104505/problem/A

简化题意:求

\[\frac{\sum_{i=1}^n i^k (n+1-i)^k}{\sum_{i=1}^n i^k} \]

其中\(1\leq n\leq 10^9,1\leq k\leq 2\times 10^5\),模数是\(10^9+7\)


题解:

自然数幂求和

分母是经典的自然数幂求和,常见方法挺多的,比如Lagrange插值,或者用生成函数的方法,可以得到就是\(\frac{e^{(k+1)x}-e^x}{e^x-1}\)的展开:

\[f_n=\sum_{i=1}^m i^n,F(x)\xlongequal{\Delta}\sum_{n=0}^\infty \frac{x^n}{n!}f_n \]

\[F(x)=\sum_{i=1}^m \sum_{n=0}^\infty \frac{(xi)^n}{n!}=\sum_{i=1}^m e^{xi}=\frac{e^{(m+1)x}-e^x}{e^x-1} \]

这样一来只需要多项式exp和求逆,对于每个固定的上界\(m\),就可以快速得到对应的自然数\(n\)次幂的和。

分母

分母:

\[=\sum_{i=1}^n i^k \sum_{j=0}^k\binom{k}{j}(n+1)^{k-j}(-1)^j i^j=\sum_{j=0}^k\binom{k}{j}(n+1)^{k-j}(-1)^j \sum_{i=1}^n i^{k+j}=\sum_{j=0}^k\binom{k}{j}(n+1)^{k-j}(-1)^j {\color{red}f_{k+j}} \]

好写到这里我们就很开心地想上去写一发了,结果发现模数是\(10^9+7\)!!!我会写任意模数ntt!!!

再看一眼这式子,\(f_{k+j}\)是关于 \(n\)\(k+j+1\)次多项式,算上后面的\((n+1)^{k-j}\),最后的结果一定能是一个\(2k+1\)次多项式,完全可以插值,当然插值的时候肯定不能算 \(k\) 次自然数幂和…

回归原来的式子,只要算

\[\sum_{i=1}^n i^k (n+1-i)^k \]

这是关于 \(n\)\(2k+1\) 次多项式,只要求出前 \(2k+2\) 项就行,赛时我盯了好久没想到怎么简单做(只想到了任意模数ntt…)

但其实我们把问题拆开看,作为上界的 \(n\) 是一个变量,而对于固定问题的\(n\) ,求和符号里的 \((n+1-i)^k\) 里的 \(n\) 其实和一个常数没什么区别,不妨直接看成是一个 \(g_n=\sum_{i=1}^n i^k (C-i)^k\),这样一来要求前 \(2k+2\) 项就可以直接计算了~,算完再插值

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
using namespace std;
typedef long long ll;
const int MOD=1e9+7;
const int N=4e5+5;
ll ksm(ll a,ll b){
	ll ret=1;a%=MOD;
	for(;b;b>>=1,a=a*a%MOD)if(b&1)ret=ret*a%MOD;
	return ret;
}
ll inv(ll x){return ksm(x,MOD-2);}
ll Lagrange(int *g,int k,int n){//g是k次多项式,需要准备1~k+1处的值,求g[n]
	vector<ll> fact,pre,suf;
	fact=pre=suf=vector<ll>(k+3,0);
	fact[0]=1;pre[0]=1;suf[k+2]=1;
	for(int i=1;i<=k+2;i++)fact[i]=fact[i-1]*i%MOD;
	for(int i=1;i<=k+2;i++)pre[i]=pre[i-1]*(n-i)%MOD;
	for(int i=k+1;i>=1;i--)suf[i]=suf[i+1]*(n-i)%MOD;
	ll ret=0;
	for(int i=1;i<=k+1;i++){
		ll t=(ll)g[i]*pre[i-1]%MOD*suf[i+1]%MOD;
		t=t*inv(fact[i-1]*fact[k-i+1]%MOD)%MOD;
		if((k-i+1)&1)t=(MOD-t)%MOD;
		ret=(ret+t)%MOD;
	}
	return ret;
}
int n,k,f[N],g[N];
int main(){
	cin>>n>>k;
	int len=2*k+1;
	rep(i,1,len+1){
		f[i]=(f[i-1]+ksm(i,k))%MOD;
		g[i]=(g[i-1]+ksm(i,k)*ksm(n+1-i,k)%MOD)%MOD;
	}
	cout<<Lagrange(g,2*k+1,n)*inv(Lagrange(f,k+1,n))%MOD;
	return 0;
}
posted @ 2023-10-04 00:26  yoshinow2001  阅读(53)  评论(0编辑  收藏  举报