【51nod1822】序列求和 V5（拉格朗日插值）（组合数学推导）

最新推荐文章于 2022-05-04 10:41:24 发布

zxyoi_dreamer

最新推荐文章于 2022-05-04 10:41:24 发布

阅读量970

点赞数

分类专栏：拉格朗日插值组合数学

本文链接：https://blog.csdn.net/zxyoi_dreamer/article/details/89048235

版权

组合数学同时被 2 个专栏收录

66 篇文章 2 订阅

订阅专栏

拉格朗日插值

10 篇文章 1 订阅

订阅专栏

传送门

解析：

如果 $r\equiv 0$ 显然答案就是0，如果 $r\equiv 1$ ，则我们要求的就是自然数的幂和，可以看一下我的博客：自然数求幂和的几种方法。

现在考虑 $r\not \equiv 1$ 的情况。

设我们要求的是 $s_k(n)=\sum\limits_{i=1}^{n}i^kr^i$

利用序列求和 V2的推导方法我们得到：

$s_k(n)=\frac{r^{n+1}n^k+\sum_{j=0}^{k-1}(-1)^{k-j}{k\choose j}s_j(n)}{r-1}$

对这个式子是用来 $O(k^2)$ 递推的，但是这道题显然要求 $O (k)$ 解决每一个询问。

我们现在考虑暴力展开一下：

$\begin{aligned} &s_0(n)=\frac{r^{n+1}-r}{r-1}\\ &s_1(n)=\frac{r^{n+1}n}{r-1}-\frac{s_0(n)}{r-1}\\ &s_2(n)=\frac{r^{n+1}n^2}{r-1}-\frac{2s_1(n)}{r-1}+\frac{s_0(n)}{r-1} \end{aligned}$

看不出来规律？展开到底：

$\begin{aligned} &s_0(n)=r^{n+1}\cdot \frac{1}{r-1}-r\cdot (\frac{1}{r-1})\\ &s_1(n)=r^{n+1}\cdot(\frac{n}{r-1}-\frac{1}{(r-1)^2})-r\cdot(-\frac{1}{(r-1)^2})\\ &s_2(n)=r^{n+1}\cdot(\frac{n^2}{r-1}-\frac{2n}{(r-1)^2}+\frac{2}{(r-1)^3}+\frac{1}{(r-1)^2})-r\cdot(\frac{2}{(r-1)^3}+\frac{1}{(r-1)^2}) \end{aligned}$

我们猜想对于 $\forall k$ ，存在 $k$ 次多项式 $S_k(n)$ ，使得： $s_k(n)=r^{n+1}S_k(n)-rS_k(0)$

现在证明 $S_k(n)$ 总是存在：

1.k=0，由上面的展开形式显然存在

2.k>0，假设结论对 $0, 1, 2, \dots \dots k - 1$ 全部成立。
$\begin{aligned} s_k(n)&=&&\frac{r^{n+1}n^k+\sum_{j=0}^{k-1}(-1)^{k-j}{k\choose j}s_j(n)}{r-1}\\ &=&&\frac{r^{n+1}n^k+\sum_{j=0}^{k-1}(-1)^{k-j}{k\choose j}[r^{n+1}S_j(n)-rS_j(0)]}{r-1}\\ &=&&r^{n+1}\cdot (\frac{n^k+\sum_{j=0}^{k-1}(-1)^{k-j}{k\choose j}S_j(n)}{r-1})+r\cdot(\frac{\sum_{j=0}^{k-1}(-1)^{k-j}{k\choose j}S_j(0)}{r-1}) \end{aligned}$

所以合法的 $S_k(n)=\frac{n^k+\sum\limits_{j=0}^{k-1}(-1)^{k-j}{k\choose j}S_j(n)}{r-1}$

现在我们将 ${n\choose k}$ 看做一个 $k$ 确定的关于 $n$ 的 $k$ 次多项式，既 ${n\choose k}=\frac{n^{\underline{k}}}{k!}$ ，其中 $n^{\underline{k}}$ 即 $k$ 次下降幂，定义为 $n^{\underline{k}}=n(n-1)(n-2)…(n-k+1)$

显然，所有的 ${n\choose k}$ 线性无关。

现在给一个引理， $\forall 0\leq m\leq k,\sum\limits_{i=0}^k(-1)^{i}{k\choose i}{i\choose m}=0$ ，证明如下：

$\begin{aligned} \sum\limits_{i=0}^k(-1)^{i}{k\choose i}{i\choose m}&=&&\sum_{i=0}^k(-1)^{i}{k\choose m}{k-m\choose i-m}\\ &=&&{k\choose m}\sum_{i=0}^{k-m}(-1)^{i+m}{k-m\choose i}\\ &=&&{k\choose m}(-1)^m(1-1)^{k-m}\\ &=&&0 \end{aligned}$

通过将 ${n\choose k}$ 看作一个多项式我们可以将 $S_k(n)$ 看作是 ${n\choose 0}{n\choose 1}{n\choose 2}…{n\choose k}$ 的线性组合。

设 $S_k(n)=\sum\limits_{i=0}^k{n\choose i}a_i$ 。

则我们有：

$\begin{aligned} \sum_{i=0}^{k+1}(-1)^i{k+1\choose i}S_k(i)&=&&\sum_{i=0}^{k+1}(-1)^i{k+1\choose i}\sum_{j=0}^{k}{i\choose j}a_j\\ &=&&\sum_{j=0}^ka_j\sum_{i=0}^{k+1}(-1)^i{k+1\choose i}{i\choose j}\\ &=&&0 \end{aligned}$

根据 $S_k(n)$ 的递推式，我们有 $S_k(n+1)=\frac{S_k(n)+(n+1)^k}{r}$

手动解方程一路回代可以得到 $S_k(0)$ ，按照最后的式子一路向上递推可以得到 $S_k(i)$ 。

然后直接拉格朗日插值就行了。

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define gc get_char
#define cs const

namespace IO{
	inline char get_char(){
		static cs int Rlen=1<<20|1;
		static char buf[Rlen],*p1,*p2;
		return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;
	}
	
	inline ll getint(){
		re char c;
		while(!isdigit(c=gc()));re ll num=c^48;
		while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
		return num;
	}
}
using namespace IO;

using std::cout;
using std::cerr;

cs int mod=985661441;
inline int add(int a,int b){return (a+b>=mod)?a+b-mod:a+b;}
inline void Inc(int &a,int b){(a+=b)>=mod?a-=mod:a;}
inline int dec(int a,int b){return a<b?a-b+mod:a-b;}
inline int mul(int a,int b){return (ll)a*b%mod;}
inline int quickpow(int a,ll b,int res=1){
	b%=mod-1;
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
template<class ...Argc>
inline int add(int a,Argc ...args){return add(a,add(args...));}
template<class ...Argc>
inline int mul(int a,Argc ...args){return mul(a,mul(args...));}

cs int N=2e5+5;

int k,r;
ll n;

int fac[N],ifac[N],inv[N];
inline int C(int n,int m){return mul(fac[n],ifac[n-m],ifac[m]);}

std::vector<int> prime;
int minpri[N];

inline void linear_sieves(int len=N-5){
	fac[0]=fac[1]=ifac[0]=ifac[1]=inv[0]=inv[1]=1;
	for(int re i=2;i<N;++i){
		fac[i]=mul(fac[i-1],i);
		inv[i]=mul(inv[mod%i],mod-mod/i);
		ifac[i]=mul(ifac[i-1],inv[i]);
		if(!minpri[i])prime.push_back(i),minpri[i]=i;
		for(re int j:prime){
			if(i*j>len)break;
			minpri[i*j]=j;
		}
	}
}

int x[N],f[N];

inline void init(int lim){
	int now=r;
	x[1]=1,f[1]=r,f[0]=0;
	for(int re i=2;i<=lim;++i){
		now=mul(now,r);
		if(minpri[i]==i)x[i]=quickpow(i,k);
		else x[i]=mul(x[minpri[i]],x[i/minpri[i]]);
		f[i]=add(f[i-1],mul(x[i],now));
	}
}

int pre[N],suf[N];
inline int calc1(int *f,int k){
	int ans=0;int n=::n%mod;
	pre[0]=suf[k+1]=1;
	for(int re i=1;i<=k+1;++i)
	pre[i]=mul(pre[i-1],dec(n+1,i));
	for(int re i=k;~i;--i)
	suf[i]=mul(suf[i+1],dec(n,i+1));
	for(int re i=0;i<=k+1;++i){
		int val=mul(f[i],ifac[i],ifac[k-i+1],pre[i],suf[i]);
		if((k^i)&1)Inc(ans,val);
		else ans=dec(ans,val);
	}
	return ans;
}

inline void calc2(){
	int inv=quickpow(r,mod-2);
	int a=1,b=0;
	pre[0]=1,suf[0]=0;
	for(int re i=1;i<=k+1;++i){
		pre[i]=mul(pre[i-1],inv);
		suf[i]=mul(add(suf[i-1],x[i]),inv);
		int coef=C(k+1,i);
		if(i&1)coef=dec(0,coef);
		Inc(a,mul(coef,pre[i]));
		Inc(b,mul(coef,suf[i]));
	}
	int x=mul(quickpow(a,mod-2),dec(0,b));
	for(int re i=0;i<=k;++i)f[i]=add(mul(pre[i],x),suf[i]);
	int ans=calc1(f,k-1);
	ans=dec(mul(ans,quickpow(r,n)),f[0]);
	cout<<mul(ans,r)<<"\n";
}

signed main(){
	linear_sieves();
	int T=getint();
	while(T--){
		k=getint(),n=getint(),r=getint()%mod;
		if(r==0){
			cout<<0<<"\n";
			continue;
		}
		int lim=std::min(n,k+1ll);
		init(lim);
		if(n==lim){
			cout<<f[n]<<"\n";
			continue;
		}
		if(r==1)cout<<calc1(f,k)<<'\n';
		else calc2();
	}
	return 0;
}

zxyoi_dreamer

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
【51nod1822】序列求和 V5（拉格朗日插值）（组合数学推导）

传送门解析：如果r≡0r\equiv 0r≡0显然答案就是0，如果r≡1r\equiv 1r≡1，则我们要求的就是自然数的幂和，可以看一下我的博客：自然数求幂和的几种方法。现在考虑r̸≡1r\not \equiv 1r̸≡1的情况。设我们要求的是sk(n)=∑i=1nikris_k(n)=\sum\limits_{i=1}^{n}i^kr^isk(n)=i=1∑nikri利用序列...
复制链接

扫一扫

专栏目录