【XSY2751】Mythological IV 线性插值

题目描述

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

　　给你\(f(0),f(1),\ldots,f(k)\)，求
\[ \sum_{i=1}^nf(i)q^i \]
　　\(k\leq 500000,n\leq {10}^{18},q\neq 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\neq 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\)时：
\[ \begin{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)\)之后线性插值插出来。
\[ \begin{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}\binom{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\n",ans);
    return 0;
}

转载于:https://www.cnblogs.com/ywwyww/p/8599489.html