【XSY2751】Mythological IV 线性插值

题目描述

　　已知(f(x))为(k)次多项式。

　　给你(f(0),f(1),ldots,f(k))，求

[sum_{i=1}^nf(i)q^i ]

　　(kleq 500000,nleq {10}^{18},q eq 1)

题解

　　当(q=0)时答案为(f(0))

　　当(q=1)时：记(S(n)=sum_{i=0}^nf(i))，易证(S(n))是一个(k+1)次多项式。直接求出(S(0)ldots S(k+1))然后线性插值即可。

　　当(q eq 1)时：记(S(n)=sum_{i=0}^{n-1}f(i)q^i=q^nG(n)-G(0))，其中(G(n))是一个(k)次多项式。

　　证明：

　　当(k=0)时显然成立。

　　假设当(k=d-1)时成立。

　　当(k=d)时：

[egin{align} S(n)&=sum_{i=0}^{n-1}f(i)q^i\ qS(n)&=sum_{i=0}^{n-1}f(i)q^{i+1}=sum_{i=1}^nf(i-1)q^i\ (q-1)S(n)&=f(n)q^n+sum_{i=0}^{n-1}(f(i)-f(i-1))q^i+f(-1) end{align} ]

　　因为(f(n)-f(n-1))是一个(d-1)次多项式，所以(sum_{i=0}^{n-1}(f(i)-f(i-1))q^i)可以被表示成(q^nP(n)-P(0))

　　所以(S(n))一定能被表示为(q^nG(n)-c)，其中(G(n)=frac{f(n-1)+P(n)}{q-1})，(c)为一个常数。

　　考虑当(n=0)时(S(n)=0)，所以(c=f(0))

　　因为(f(n-1))是一个(d)次多项式，(P(n))是一个(d)次多项式，所以(G(n))也是一个(d)次多项式。

　　

　　现在要算(G(n))，可以算出(G(0)ldots G(k))之后线性插值插出来。

[egin{align} S(n)&=sum_{i=0}^{n-1}f(i)q^i\ S(n+1)-S(n-1)&=q^{n+1}G(n+1)-q^nG(n)=f(n)q^n\ qG(n+1)&=G(n)+f(n)\ G(n+1)&=frac{G(n)+f(n)}{q} end{align} ]

　　所以每个(G(n))都可以被表示为(a_iG(0)+b_i)

　　由于(G(n))是一个(k)次多项式，那么就满足(k+1)次差分之后的值为(0)

[sum_{i=0}^{k+1}{(-1)}^{k+1-i}inom{k+1}{i}G(i)=0 ]

　　这是一个关于(G(0))的一元一次方程，可以解出(G(0))的值。

　　然后递推得到(G(1)ldots G(k))，线性插值插出(G(n+1))

　　时间复杂度：(O(k+log n))

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll p=1000000007;
ll n,q;
int k;
ll fp(ll a,ll b)
{
	ll s=1;
	for(;b;b>>=1,a=a*a%p)
		if(b&1)
			s=s*a%p;
	return s;
}
ll inv[500010];
ll fac[500010];
ll ifac[500010];
ll f1[500010];
ll f2[500010];
ll g[500010];
ll f[500010];
ll getc(int x,int y)
{
	return fac[x]*ifac[y]%p*ifac[x-y]%p;
}
ll pre[500010];
ll suf[500010];
int main()
{
#ifndef ONLINE_JUDGE
	freopen("b.in","r",stdin);
	freopen("b.out","w",stdout);
#endif
	scanf("%lld%d%lld",&n,&k,&q);
	n++;
	for(int i=0;i<=k;i++)
		scanf("%lld",&f[i]);
	inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
	for(int i=2;i<=k+1;i++)
	{
		inv[i]=-p/i*inv[p%i]%p;
		ifac[i]=ifac[i-1]*inv[i]%p;
	}
	for(int i=2;i<=k+1;i++)
		fac[i]=fac[i-1]*i%p;
	ll invq=fp(q,p-2);
	f1[0]=1;
	f2[0]=0;
	for(int i=1;i<=k+1;i++)
	{
		f1[i]=f1[i-1]*invq%p;
		f2[i]=(f2[i-1]+f[i-1])*invq%p;
	}
	ll v1=0,v2=0;
	for(int i=0;i<=k+1;i++)
	{
		v1=(v1+getc(k+1,i)*f1[i]%p*((k+1-i)&1?-1:1))%p;
		v2=(v2+getc(k+1,i)*f2[i]%p*((k+1-i)&1?-1:1))%p;
	}
	g[0]=-v2*fp(v1,p-2)%p;
	for(int i=1;i<=k+1;i++)
		g[i]=(f1[i]*g[0]+f2[i])%p;
	ll ans=0;
	ll n2=n%p;
	for(int i=0;i<=k;i++)
	{
		pre[i]=n2-i;
		if(i)
			pre[i]=pre[i-1]*pre[i]%p;
	}
	for(int i=k;i>=0;i--)
	{
		suf[i]=n2-i;
		if(i!=k)
			suf[i]=suf[i+1]*suf[i]%p;
	}
	for(int i=0;i<=k;i++)
	{
		ll v=g[i];
		if(i)
			v=v*pre[i-1]%p;
		if(i!=k)
			v=v*suf[i+1]%p;
		v=v*ifac[i]%p;
		v=v*ifac[k-i]*((k-i)&1?-1:1)%p;
		ans=(ans+v)%p;
	}
	ans=ans*fp(q,n)%p;
	ans=(ans-g[0])%p;
	ans=(ans+p)%p;
	printf("%lld
",ans);
	return 0;
}