显然答案为：
$$ \sum_{i=1}^{n}F(i)\left({n-i+1\choose m}^2-{n-i\choose m}^2\right) $$
因为 ${n\choose m}$ 是一个关于 $n$ 的 $m$ 次多项式，又因为 $n$ 次多项式的前缀和是 $n+1$ 次多项式。所以上式其实是一个 $3m+1$ 次多项式。
不难想到令：
$$ f(k)=\sum_{i=1}^{k}F(i)\left({n-i+1\choose m}^2-{n-i\choose m}^2\right) $$
然后求 $f(0),f(1),f(2),\cdots,f(3m+1)$ 后线性插值即可。
$f$ 是个次数不超过 $m$ 的多项式，现在要快速求出 $f(m+1)$ 至 $f(3m+1)$ 。
考虑拉格朗日插值的式子：
$$ \begin{aligned} f(k)&=\sum_{i=0}^{m}f(i)\prod_{i\not= j}\frac{x-j}{i-j}\\
&=\sum_{i=0}^{m}f(i)\frac{\prod_{i\not= j}(x-j)}{\prod_{i\not= j}(i-j)}\\ \end{aligned} $$
显然下面的 $\prod_{i\not= j}(i-j)$ 预处理阶乘后可以 $O(1)$ 算，令 $g(i)=\frac{f(i)}{\prod_{i\not= j}(i-j)}$ ，有：
$$ \begin{aligned} f(k)&=\sum_{i=0}^{m}g(i)\prod_{i\not= j}(x-j)\\ &=\sum_{i=0}^{m}g(i)\prod_{j=0}^{m}(x-j)\frac{1}{(x-i)}\\ &=\prod_{j=0}^{m}(x-j)\sum_{i=0}^{m}g(i)\frac{1}{x-i}\\ &=x^{\underline{m+1}}\sum_{i=0}^{m}g(i)\frac{1}{x-i}\\ \end{aligned} $$
（式子没列清除导致这里想了好久，最后想出来个线段树做法，感觉不行。瞅了眼题解才发现就是个简单卷积 >_<
后面的直接 $\rm{NTT}$ 卷就好了。 时间复杂度 $O(m\log m)$ 。
// {{{ Poly
const int LogN=26;
namespace Poly {
int tim,lim,inv[LogN],w[1<<24|5];
poly rev[LogN];
inline void init_r(int len) {
for(lim=1,tim=0;lim<len;lim<<=1,++tim);
if(rev[tim].size()) return ;
rev[tim].rez(lim),inv[tim]=modpow(lim,mod-2);
for(int i=0;i<lim;++i) rev[tim][i]=(rev[tim][i>>1]>>1)|((i&1)<<(tim-1));
}
inline int W(int a,int b) {
return modpow(modpow(3,(mod-1)>>(a+1)),b);
}
inline void NTT(poly &f,int typ) {
static ull g[1<<24|5];
for(int i=0;i<lim;++i) g[rev[tim][i]]=f[i];
for(int p=1,s=0,t=0;p<lim;p<<=1,++t) {
w[0]=1,w[1]=modpow(3,(mod-1)>>(t+1));
lep(i,2,p) w[i]=mul(w[i-1],w[1]);
for(int k=0;k<lim;k+=p<<1)
for(int l=k;l<k+p;++l)
s=mul(g[l+p]%mod,w[l-k]),g[l+p]=g[l]+mod-s,g[l]+=s;
}
for(int i=0;i<lim;++i) f[i]=g[i]%mod;
if(~typ) return ;
std::reverse(++f.begin(),f.end());
for(int i=0;i<lim;++i) f[i]=mul(f[i],inv[tim]);
}
inline poly operator * (poly a,poly b) {
int n=a.size(),m=b.size(),t=n+m-1;
init_r(t);
a.rez(lim),NTT(a,1);
b.rez(lim),NTT(b,1);
for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
return NTT(a,-1),a.rez(t),a;
}
}
using Poly::operator *;
// }}}
const int N=4e6+5;
int n,m,l,f[N],inv[N],fac[N];
inline void init() {
fac[0]=1;
lep(i,1,l) fac[i]=mul(fac[i-1],i);
inv[l]=modpow(fac[l],mod-2);
rep(i,l,1) inv[i-1]=mul(inv[i],i);
}
inline void solve1() { /*get f(m+1~3m+1)*/
poly p,g; g.rez(m+1),p.rez(l+1);
lep(i,0,m) g[i]=mul(f[i],mul(inv[i],((m-i)&1)?mod-inv[m-i]:inv[m-i]));
lep(i,0,l) p[i]=modpow(i,mod-2);
poly h=p*g;
int t=fac[m+1];
lep(k,m+1,l) f[k]=mul(h[k],t),t=mul(t,mul(k+1,p[k-m]));
}
inline void solve2() {
static int d[N];
int s1=1,s2=1;
lep(i,0,m-1) s1=mul(s1,n-i+1),s2=mul(s2,n-i);
#define S(x) mul(x,x)
lep(i,0,l) {
if(i) d[i]=d[i-1];
pls(d[i],mul(f[i],dec(S(mul(s1,inv[m])),S(mul(s2,inv[m])))));
s1=s2,s2=(n-i>m)?mul(s2,mul(n-m-i,modpow(n-i,mod-2))):0;
}
#undef S
if(n<=l) return printf("%d\n",d[n]),void();
int ans=0,res=1;
lep(i,0,l) res=mul(res,dec(n,i));
lep(i,0,l) d[i]=mul(d[i],mul(inv[i],((l-i)&1)?mod-inv[l-i]:inv[l-i]));
lep(i,0,l) pls(ans,mul(mul(res,modpow(dec(n,i),mod-2)),d[i]));
printf("%d\n",ans);
}
int main() {
IN(n,m),l=3*m+1,--n;
lep(i,0,m) IN(f[i]);
init();
solve1();
solve2();
return 0;
}
0 条评论