51nod1229 序列求和 V2 数论+自然数幂和

本文介绍了一种求解形如∑i=1nik⋅ri的序列和的方法，针对不同的参数范围采用不同策略，包括利用二项式定理展开、递推公式及拉格朗日插值等数学工具。

Description

求 $∑i=1nik⋅ri\sum_{i=1}^n{i^k\cdot r^i}$
其中 $n,r≤1018,k≤2000n,r\le10^{18},k\le2000$

Solution

不会做(艹皿艹 )，怒orz题解
记 $S(k)=∑i=1nik⋅riS(k)=\sum\limits_{i=1}^{n}i^k\cdot r^i$ ，那么 $(r−1)S(k)=rn+1⋅nk−r+∑i=0k−1ri[(i−1)k−ik](r-1)S(k)=r^{n+1}\cdot n^k-r+\sum\limits_{i=0}^{k-1}{r^i\left[{(i-1)}^k-i^k\right]}$
我们把 ${(i-1)}^k$ 二项式定理展开，那么就是 $(r−1)S(k)=rn+1⋅nk−r+∑j=0k−1(kj)(−1)k−j∑i=2nri⋅ij(r-1)S(k)=r^{n+1}\cdot n^k-r+\sum\limits_{j=0}^{k-1}{\binom{k}{j}{\left(-1\right)}^{k-j}\sum\limits_{i=2}^{n}r^i\cdot i^j}$
然后发现最右边就是 $S (j) - r$ ，于是我们就可以愉快地k²递推这个东西了。。

然后当r=1的时候上面就不成立了，这个时候求自然数幂和就行。
拉格朗日插值： $F(x)=∑i=0k+1yi∏j=0,j≠ik+1x−xjxi−xjF(x)=\sum_{i=0}^{k+1}y_i\prod_{j=0,j =\not i}^{k+1}\frac{x-x_j}{x_i-x_j}$
对于自然数幂和可以发现x都是连续的，那么就变成了下面这个样子 $F(x)=∑i=0k+1yi∏j=0i−1(n−j)∏j=i+1k+1(n−j)i!(k+1−i)!(−1)k+1−iF(x)=\sum_{i=0}^{k+1}y_i\frac{\prod_{j=0}^{i-1}{(n-j)}\prod_{j=i+1}^{k+1}{(n-j)}}{i!(k+1-i)!{(-1)}^{k+1-i}}$
分子刚好就是前缀后缀积，分母刚好就是阶乘了。自然数幂和的前k项可以线性筛，阶乘和逆元可以线性求，前缀后缀积同理，那么整个做下来都是O(n)的

Code

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define drp(i,st,ed) for (int i=st;i>=ed;--i)

typedef long long LL;
const int MOD=1e9+7;
const int N=200005;

LL s[N],fac[N],inv[N],pre[N],nex[N];

LL ksm(LL x,LL dep) {
	LL res=1;
	x=(x%MOD+MOD)%MOD;
	for (;dep;dep>>=1) {
		(dep&1)?(res=res*x%MOD):0;
		x=x*x%MOD;
	}
	return res;
}

void solve1(LL n,int k) {
	LL ans=0; n%=MOD;
	rep(i,1,k+1) s[i]=(s[i-1]+ksm(i,k))%MOD;
	pre[0]=n; nex[k+2]=1;
	rep(i,1,k+1) pre[i]=pre[i-1]*(n+MOD-i)%MOD;
	drp(i,k+1,0) nex[i]=nex[i+1]*(n+MOD-i)%MOD;
	rep(i,0,k+1) {
		LL fz=((i>0)?(pre[i-1]):1)*(nex[i+1])%MOD;
		LL fm=fac[i]*fac[k+1-i]%MOD;
		LL tmp=fz*ksm(fm,MOD-2)%MOD*s[i]%MOD;
		if ((k+1-i)&1) tmp=(MOD-tmp);
		ans=(ans+tmp)%MOD;
	}
	printf("%lld\n", ans);
}

void solve2(LL n,LL k,LL r) {
	s[0]=(ksm(r,n+1)+MOD-r)*ksm(r-1,MOD-2)%MOD;
	rep(i,1,k) {
		s[i]=(ksm(r,n+1)*ksm(n,i)%MOD+MOD-r)%MOD;
		rep(j,0,i-1) {
			LL tmp=fac[i]*inv[j]%MOD*inv[i-j]%MOD*(s[j]+MOD-r)%MOD;
			if ((i-j)&1) tmp=(MOD-tmp)%MOD;
			s[i]=((s[i]+tmp)%MOD+MOD)%MOD;
		}
		s[i]=s[i]*ksm(r-1,MOD-2)%MOD;
	}
	s[k]=(s[k]%MOD+MOD)%MOD;
	printf("%lld\n", s[k]);
}

int main(void) {
	freopen("data.in","r",stdin);
	fac[0]=fac[1]=inv[0]=inv[1]=1;
	rep(i,2,N-1) {
		fac[i]=fac[i-1]*i%MOD;
		inv[i]=inv[MOD%i]*(MOD-MOD/i)%MOD;
	}
	rep(i,2,N-1) inv[i]=inv[i-1]*inv[i]%MOD;
	int T; scanf("%d",&T);
	for (;T--;) {
		LL n,r,k; scanf("%lld%lld%lld",&n,&k,&r);
		r%=MOD;
		if (r==1) solve1(n,k);
		else solve2(n,k,r);
	}
	return 0;
}