显然答案为:

$$ \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;
}
分类: 文章

Qiuly

QAQ

0 条评论

发表回复

Avatar placeholder

您的电子邮箱地址不会被公开。 必填项已用 * 标注