令$S(n)=\sum\limits_{i=0}^{n-1}f(i)q^i$，那么存在一个次数$\leq k$的多项式使得$S(n)=q^ng(n)-g(0)$（证明来自杜教的PPT）

设$f$的次数为$d$，对$d$归纳

当$d=0$时，取$g(n)=\dfrac{f(0)}{q-1}$可使其成立

假设当$d=k-1$时，命题成立

当$d=k$时，有$S(n)=\sum\limits_{i=0}^{n-1}f(i)q^i$和$qS(n)=\sum\limits_{i=0}^{n-1}f(i)q^{i+1}=\sum\limits_{i=1}^nf(i-1)q^i$

将两式相减，得到$(q-1)S(n)=f(n-1)q^n-f(0)+\sum\limits_{i=1}^{n-1}\left[f(i-1)-f(i)\right]q^i$

因为$f(i-1)-f(i)$是一个$k-1$次多项式，所以$\sum\limits_{i=0}^{n-1}\left[f(i-1)-f(i)\right]q^i$可以被表示为$q^nh(n)-h(0)$，其中$h(x)$是一个次数$\leq k-1$的多项式

代回原式继续处理

$\begin{align*}(q-1)S(n)&=f(n-1)q^n-f(0)+q^nh(n)-h(0)\\S(n)&=q^n\dfrac{f(n-1)+h(n)}{q-1}-\dfrac{f(0)+h(0)}{q-1}\end{align*}$

令$F(n)=\dfrac{f(n-1)+h(n)}{q-1},c=\dfrac{f(0)+h(0)}{q-1}$，那么原式变为$S(n)=q^nF(n)-c$

因为$S(0)=0$，所以$F(0)-c=0$，即$c=F(0)$

而$F(n)$是一个$k$次多项式（含$f$）

所以当$d=k$时，命题成立

有了这个定理，我们可以做一番推导

一方面，由$S(n)$的定义，$S(n)-S(n-1)=q^{n-1}f(n-1)$，另一方面，由上述定理，$S(n)-S(n-1)=q^ng(n)-q^{n-1}g(n-1)$，联立起来我们可以得到$qg(n)=f(n-1)+g(n-1)$

如果我们知道了$g(0)$，那么$g(1)\cdots g(k+1)$都可求得

不妨设$g(0)$为一个未知数，那么$g(1)\cdots g(k+1)$都是关于$f(0)$的一次函数

另一方面，$k$次多项式$g(x)$的$k+1$阶差分等于$0$，即$\triangle^{k+1}g(x)=\sum\limits_{i=0}^{k+1}(-1)^i\binom{k+1}ig(x-i)=0$

代入$k+1$，我们得到$\sum\limits_{i=0}^{k+1}(-1)^i\binom{k+1}ig(k+1-i)=0$，它是关于$g(0)$的一元一次方程，直接解即可

于是我们得到了$g(0)\cdots g(k+1)$，使用线性插值的技巧，我们可以做到在$O(k)$的时间内求出$g(n)$

关于线性插值，这里不展开讲，只放一个miskcoo的博客和公式（推荐去看他的博客，推导过程挺好玩的）

设$k$次多项式为$P(x)$，已经知道$P(0)\cdots P(k)$，则$P(x)=\sum\limits_{j=0}^k(-1)^{k-j}\binom xj\binom{x-j-1}{k-j}P(j)$

第一个组合数预处理，第二个组合数从大到小枚举$j$的时候迭代计算

于是这道题就做完了，总复杂度$O(k)$，真是奇妙

#include<stdio.h>
typedef long long ll;
const int mod=1000000007;
int mul(int a,int b){return a*(ll)b%mod;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
ll n;
int k,q,invq,fac[500010],rfac[500010],inv[500010],f[500010],G[500010],cn[500010];
int C(int n,int k){return mul(fac[n],mul(rfac[n-k],rfac[k]));}
int pow(int a,ll b){int s=1;while(b){if(b&1)s=mul(s,a);a=mul(a,a);b>>=1;}return s;
}
struct single{int x,c;single(int a=0,int b=0){x=a;c=b;}
}g[500010],s;
single operator+(single a,single b){return single(ad(a.x,b.x),ad(a.c,b.c));}
single operator+(single a,int b){return single(a.x,ad(a.c,b));}
single operator*(single a,int b){return single(mul(a.x,b),mul(a.c,b));}
int main(){int i,t,c,gn;scanf("%lld%d%d",&n,&k,&q);for(i=0;i<=k;i++)scanf("%d",f+i);fac[0]=1;for(i=1;i<=k+1;i++)fac[i]=mul(fac[i-1],i);rfac[k+1]=pow(fac[k+1],mod-2);for(i=k+1;i>0;i--)rfac[i-1]=mul(rfac[i],i);for(i=1;i<=k+1;i++)inv[i]=mul(rfac[i],fac[i-1]);g[0].x=1;invq=pow(q,mod-2);for(i=1;i<=k+1;i++)g[i]=(g[i-1]+f[i-1])*invq;t=1;for(i=0;i<=k+1;i++){s=s+((g[k+1-i]*C(k+1,i))*t);t=-t;}G[0]=mul(-s.c,pow(s.x,mod-2));for(i=1;i<=k;i++)G[i]=mul(ad(G[i-1],f[i-1]),invq);cn[0]=1;n++;for(i=0;i<k;i++)cn[i+1]=mul(cn[i],mul((n-i)%mod,inv[i+1]));c=1;t=1;gn=0;for(i=k;i>=0;i--){gn=ad(gn,mul(mul(G[i],cn[i]),mul(t,c)));t=-t;c=mul(c,mul((n-i)%mod,inv[k-i+1]));}gn=de(mul(pow(q,n),gn),G[0]);if(gn<0)gn+=mod;printf("%d",gn);
}