拉格朗日插值

标签：拉格朗 power 插值 sum int MAXN binom mod

如果答案能表示为一个 \(K\) 次多项式的形式可以考虑插值求解，\(O(k^2)\)。

比如求 \(\sum_{i=1}^ni^k\)，可以表示为一个 \(k+1\) 次多项式 \(f(n)\)，事实上次数开高是没影响的（插值出来系数为0）但是不能插低。

一个人的数论

毒瘤。

\[f_k(n)=\sum_{i=1}^{n}[gcd(i,n)=1]i^k \]

\[=\sum_{i=1}^ni^k\sum_{d|gcd(i,n)}\mu(d) \]

\[=\sum_{i=1}^ni^k\sum_{d=1}^n[d|i][d|n]\mu(d) \]

\[=\sum_{d=1}^n[d|n]\mu(d)\sum_{d|i}^ni^k \]

\[=\sum_{d|n}\mu(d)\sum_{i=ad}^n(ad)^k \]

\[=\sum_{d|n}\mu(d)d^k\sum_{i=1}^{\lfloor n/d\rfloor}i^k \]

后半部分的求和就是拉插，然后这个题的 \(n\) 特别大。再看一下式子发现前半部分是个狄卷的形式，考虑把后边的式子看作一个多项式(因为本来就是)，又设

\[n=\prod p_i^{a_i},n'=\prod p_i \]

就有

\[Ans=\sum_{d|n'}\mu(d)d^k\sum_{i=0}^{k+1}a_i(\frac{n}{d})^i \]

\[=\sum_{i=0}^{k+1}a_in^i\sum_{d|n'}\mu(d)d^{k-i} \]

考虑到后半部分仍为一个积性函数的形式，且根据莫函数性质，设

\[ F(x)=\sum_{d|x}\mu(d)d^{k-i} \]

则

\[Ans=\sum_{i=0}^{k+1}a_in^i\prod_{j=1}^w F(p_j) \]

\[=\sum_{i=0}^{k+1}a_in^i\prod_{j=1}^w(1-p_j^{k-i}) \]

阿好的现在终于得到了答案式子但是不难发现 \(a_i\) 是不知道怎么求的，事实上拉插可以跑出来多项式系数的，放一下代码。

#include<bits/stdc++.h>
#define MAXN 2005
#define MAXM 105
#define int long long
using namespace std;
const int mod=1e9+7;
int K,n,N;
int p[MAXN],t[MAXN];
inline int qpow(int base,int power){
	int res=1;
	while(power){
		if(power&1)res=res*base%mod;
		base=base*base%mod;
		power>>=1;
	}
	return res%mod;
}
inline int inv(int x){
	return qpow(x,mod-2);
}
struct Lagrange{
	int x[MAXM],y[MAXM];
	int t[MAXM],pre[MAXM],f[MAXM],g[MAXM];
	inline void INIT(int K){
		x[0]=y[0]=0;
		for(int i=1;i<=K;i++){
			x[i]=i;
			y[i]=(y[i-1]+qpow(i,K-1))%mod;
		}
		for(int i=0;i<=K;i++){
			t[i]=1;
			for(int j=0;j<=K;j++){
				if(i==j)continue;
				t[i]=t[i]*(x[i]-x[j]+mod)%mod;
			}
			t[i]=y[i]*inv(t[i])%mod;
		}
		pre[0]=1;
		for(int i=0;i<=K;i++){
			for(int j=K;j>0;j--)pre[j]=(pre[j-1]-pre[j]*x[i]%mod+mod)%mod;
			pre[0]=(pre[0]*(-x[i])%mod+mod)%mod;
		}
		for(int i=0;i<=K;i++){
			int val=inv(mod-x[i]);
			g[0]=pre[0]*val%mod;
			for(int j=1;j<=K;j++)g[j]=(pre[j]-g[j-1]+mod)%mod*val%mod;
			for(int j=0;j<=K;j++)f[j]=(f[j]+t[i]*g[j]%mod)%mod;
		}
	}
}Lag;
int ans;
signed main(){
	scanf("%lld%lld",&K,&n);N=1;
	for(int i=1;i<=n;i++)scanf("%lld%lld",&p[i],&t[i]),N=N*qpow(p[i],t[i])%mod;
	Lag.INIT(K+1);
//	printf("ced\n");
	int bas=1;
	for(int i=0;i<=K+1;i++){
	//	printf("nowi=%lld\n",i);
		int res=bas*Lag.f[i]%mod;
		bas=bas*N%mod;
		for(int j=1;j<=n;j++)res=res*(1-qpow(p[j],K-i-1+mod)+mod)%mod;
		ans=(ans+res+mod)%mod;
	}
	printf("%lld\n",ans);
	return 0;
}

教科书般的亵渎

不难发现需要 \(m+1\) 次出牌。

不难发现每次造成的贡献为阶梯状有缺口。

把缺口 \(b_i\) 从小到大排序，\(b_0\)=0。

\[Ans=\sum_{i=0}^m(\sum_{j=b_i}^n(j-b_i)^k-\sum_{j=i}^m(b_j-b_i)^k) \]

\[Ans=\sum_{i=0}^m\sum_{j=0}^{n-b_i}j^k-\sum_{i=0}^m\sum_{j=i}^m(b_j-b_i)^k \]

然后就做出来了。

XLkxc

预处理 \(f\)，\(f\) 为 \(k+1\) 次多项式，\(g\) 为 \(f\) 的前缀和，为 \(k+2\) 次多项式，插值计算。
\(Ans\) 为 \(g\) 的前缀和，为 \(k+3\) 次多项式，仍然插值计算。

多项式笑传之插插倍。

成绩比较

更是毒瘤我草。

考虑容斥。

\[Ans=\sum_{i=K}^{min(N-R)}(-1)^{i-k}\binom{i}{k}*\binom{n}{i}\prod_{j=1}^M\binom{n-1-i}{n-R_j-i}\sum_{x=1}^{U_i}x^{n-R_j}(U_j-x)^{R_j-1} \]

就是说每次从 \(n-1\) 个人里挑 \(i\) 个完虐，然后分别乘法原理每一科的答案，由于每一科都有排名所以要从剩下的 \(n-1-i\) 个没被完虐的人里挑 \(n-R_j-i\) 个来低于他这一课然后美句他这一课的分数 \(x\) 低于的有 \(x\) 种取值高于的有 \(U_j-x\) 种取值。

后面这一坨

\[\sum_{x=1}^{U_i}x^{n-R_j}(U_j-x)^{R_j-1} \]

拿二项式定理展开

\[=\sum_{x=1}^{U_j}x^{n-R_j}\sum_{k=0}^{R_j-1}\binom{R_j-1}{k}(-1)^kx^kU_j^{R_j-k-1} \]

\[=\sum_{k=0}^{R_j-1}\binom{R_j-1}{k}(-1)^kU_j^{R_j-k-1}\sum_{x=1}^{U_j}x^{n-R_j+k} \]

啊终于变成我们想看到的样子了qwq然后就到了非常蛋疼的实现环节。总之需要 \(O(n)\) 插值（x连续）和巨量预处理，我的我要爆了。

\[Ans=\sum_{i=K}^{min(N-R)}(-1)^{i-K}\binom{i}{K}*\binom{n}{i}\prod_{j=1}^M\binom{n-1-i}{n-R_j-i}\sum_{k=0}^{R_j-1}\binom{R_j-1}{k}(-1)^kU_j^{R_j-k-1}\sum_{x=1}^{U_j}x^{n-R_j+k} \]

#include<bits/stdc++.h>
#define int long long
#define MAXN 205
using namespace std;
const int mod=1e9+7;
int n,m,k;
int U[MAXN],R[MAXN];
int f[MAXN],g[MAXN];
inline int qpow(int base,int power){
	int res=1;
	while(power){
		if(power&1)res=res*base%mod;
		base=base*base%mod;
		power>>=1;
	}
	return res%mod;
}
inline void INIT(){
	f[0]=g[0]=1;
	for(int i=1;i<MAXN;i++)f[i]=f[i-1]*i%mod;
	g[MAXN-1]=qpow(f[MAXN-1],mod-2);
	for(int i=MAXN-2;i>=1;i--)g[i]=g[i+1]*(i+1)%mod;
}
inline int inv(int x){
	return qpow(x,mod-2);
}
inline int getC(int x,int y){
	if(y>x)return 0;
	if(!y)return 1;
	return f[x]*g[y]%mod*g[x-y]%mod;
}
int sav[MAXN];
struct Lagrange{
	int y[MAXN],pre[MAXN],nxt[MAXN];
	inline void INIT(int a[],int K){
		for(int i=0;i<=K;i++)y[i]=a[i];
	}
	inline int calc(int K,int x0){
		pre[0]=nxt[K+1]=1;
		for(int i=1;i<=K;i++)pre[i]=pre[i-1]*(x0-i)%mod;
		for(int i=K;i>=1;i--)nxt[i]=nxt[i+1]*(x0-i)%mod;
		int res=0;
		for(int i=1;i<=K;i++){
			int val=f[i-1]*f[K-i]%mod;
			if((K-i)&1)val=-val;
		//	printf("nowcon=%lld %lld %lld %lld\n",y[i],p)
			res=(res+y[i]*pre[i-1]%mod*nxt[i+1]%mod*qpow(val,mod-2)%mod)%mod;
		}
	//	printf("res=%lld\n",res);
		return res;
	}
}Lag[MAXN];
int Uk[MAXN][MAXN];//(Ui)^j
int lsav[MAXN][MAXN];//lag:Ui j ci
int Hb[MAXN];//j
int Lim=mod,ans;
signed main(){
	INIT();
	scanf("%lld%lld%lld",&n,&m,&k);
	for(int i=1;i<=m;i++)scanf("%lld",&U[i]);
	for(int i=1;i<=m;i++)scanf("%lld",&R[i]),Lim=min(Lim,n-R[i]);
	for(int i=1;i<=n;i++){
	//	printf("nowi=%lld\n",i);
		for(int j=1;j<=n+1;j++)
			sav[j]=(sav[j-1]+qpow(j,i))%mod;
		Lag[i].INIT(sav,n+1);
	}
	
	for(int i=1;i<=m;i++){
		Uk[i][0]=1;
		for(int j=1;j<=n;j++)Uk[i][j]=Uk[i][j-1]*U[i]%mod;
	}
	for(int i=1;i<=m;i++)for(int j=1;j<=n;j++){
		lsav[i][j]=Lag[j].calc(n+1,U[i]);
	}
	for(int j=1;j<=m;j++){
		for(int d=0,bas=1;d<=R[j]-1;d++,bas=-bas)
			Hb[j]=(Hb[j]+bas*getC(R[j]-1,d)%mod*Uk[j][R[j]-d-1]%mod*lsav[j][n-R[j]+d]%mod+mod)%mod;
	}
	for(int i=k,bas=1;i<=Lim;i++,bas=-bas){
		int P1=(bas*getC(i,k)*getC(n-1,i)%mod+mod)%mod;
		int val=1;
		for(int j=1;j<=m;j++)val=val*getC(n-i-1,n-R[j]-i)%mod*Hb[j]%mod;
		ans=(ans+P1*val%mod+mod)%mod;
	}
	printf("%lld",ans);
	return 0;
}

calc

咕咕咕

标签：拉格朗,power,插值,sum,int,MAXN,binom,mod
From： https://www.cnblogs.com/Claimo0611/p/18639546