bzoj3157&&3516&&4126 国王奇遇记【拉格朗日差值+高阶差分】

最新推荐文章于 2019-11-23 21:13:32 发布

原创最新推荐文章于 2019-11-23 21:13:32 发布 · 708 阅读

1 ·

CC 4.0 BY-SA版权

--------数学问题-------- 同时被 2 个专栏收录

57 篇文章

订阅专栏

拉格朗日插值法

3 篇文章

订阅专栏

介绍了一种求解特定形式多项式求和问题的高效算法，适用于k次多项式的求和，通过对多项式进行转化和利用差分性质，实现了O(k)的时间复杂度。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

这里写图片描述

Input

共一行包括两个正整数N和M。

Output

共一行为所求表达式的值对10^9+7取模的值。

Sample Input

5 3

Sample Output

36363

HINT

1<=N<=10^9,1<=M<=500000

解题思路：

说一个不仅能求原题，还能求其拓展形式 $\sum\limits_{i=0}^nf(i)q^i\ (f(i)是一个k次多项式)$ 的 $O(k)$ 算法。
(最好 $f(i)$ 可以快速算或以点值形式给出，当然dalao们可以直接上多项式多点求值)

首先 $q=1$ 时还是要特判。
当 $q\ne 1$ 时，我们设 $S(n)=\sum\limits_{i=0}^{n-1}f(i)q^i$ ，所以我们要求 $S(n+1)$ 。
有一个结论：存在一个 $\le k$ 次多项式 $g(x)$ 使得 $S(n)=q^ng(n)-g(0)$
考虑用归纳法证明：
当 $n=0$ 时，显然成立。
当 $n=1$ 时， $S(n)=f(0)=qg(1)-g(0)$ ， $g(1)=(f(0)+g(0))/q$ ，也成立。
假设 $n\le m$ 时都成立，那么当 $n=m+1$ 时
$qS(n)=\sum\limits_{i=0}^{n-1}f(i)q^{i+1}=q^{n+1}g(n)-qg(0)$
$S(n+1)=\sum\limits_{i=0}^nf(i)q^i=q^{n+1}g(n+1)-g(0)$
两式相减可得 $\sum\limits_{i=1}^n(f(i)-f(i-1))q^i+f(0) =q^{n+1}(g(n+1)-g(n))+(q-1)g(0)$
而左边 $\sum\limits_{i=1}^{n-1}(f(i)-f(i-1))q^i+f(0)+(f(n)-f(n-1))q^n\\=q^n(g(n)-g(n-1))+(q-1)g(0)+(f(n)-f(n-1))q^n$

代入原式，两边同除以 $q^n$ 可得： $g(n)-g(n-1)+f(n)-f(n-1)=q(g(n+1)-g(n))$

所以 $g(n)=((q+1)g(n)-g(n-1)+f(n)-f(n-1))/q$ 也是个 $\le k$ 次多项式。得证。

现在问题就是如何求 $g(x)$ 。

$S(n)-S(n-1)=q^ng(n)-q^{n-1}g(n-1) =q^{n-1}f(n-1)$
即是 $g(n)=(g(n-1)+f(n-1))/q$

所以如果我们知道了 $g(0)$ 的值，就能推算出 $g(1),g(2),...g(k+1)$

不妨设 $g(0)=x$ ，由递推式可知 $g(1),g(2),g(k+1)$ 都可表示为一次函数的形式。

而 $k$ 次多项式的 $k+1$ 次差分为0，所以可以列出 $k+1$ 阶差分表达式: $\sum\limits_{i=0}^{k+1}(-1)^i\binom{k+1}{i}g(k+1-i)=0$ 从而解出 $g(0)$ 。

然后直接朗格朗日差值求出 $g(n+1)$ 就可以得到 $S(n+1)$ 了。

代码如下，其实很简洁：

#include<bits/stdc++.h>
#define ll long long
using namespace std; 

ll getll()
{
    ll i=0,f=1;char c;
    for(c=getchar();(c!='-')&&(c<'0'||c>'9');c=getchar());
    if(c=='-')f=-1,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}

const int N=500005,mod=1e9+7;
ll n,m;
ll f[N],a[N][2],g[N],fac[N],inv[N],pre[N],suf[N];

int Pow(ll x,ll y)
{
    ll res=1;
    for(;y;y>>=1,x=x*x%mod)
        if(y&1)res=res*x%mod;
    return res;
}

inline int C(int x,int y)
{
    return fac[x]*inv[y]%mod*inv[x-y]%mod;
}

int G(ll x)
{
    if(x<m)return g[x];
    ll res=0,t=1;suf[m+1]=1;x%=mod;
    for(int i=0;i<=m;i++)pre[i]=t*(x-i)%mod,t=pre[i];
    for(int i=m;i>=0;i--)suf[i]=suf[i+1]*(x-i)%mod;
    for(int i=0;i<=m;i++)
    {
        ll t=suf[i+1]*(i==0?1:pre[i-1])%mod;
        t=t*g[i]%mod*inv[i]%mod*(((m-i)&1)?mod-inv[m-i]:inv[m-i])%mod;
        res=(res+t)%mod;
    }
    return (res+mod)%mod;
}

int main()
{
    //freopen("lx.in","r",stdin); 
    scanf("%lld%lld",&n,&m);
    if(m==1)
    {
        printf("%lld\n",n*(n+1)%mod*Pow(2,mod-2)%mod);
        return 0;
    }
    for(int i=0;i<=m;i++)f[i]=Pow(i,m);
    fac[0]=1;
    for(int i=1;i<=m+1;i++)fac[i]=fac[i-1]*i%mod;
    inv[m+1]=Pow(fac[m+1],mod-2);
    for(int i=m;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod;
    ll t=Pow(m,mod-2);a[0][1]=1,a[0][0]=0;
    for(int i=1;i<=m+1;i++)
        a[i][1]=a[i-1][1]*t%mod,a[i][0]=(a[i-1][0]+f[i-1])*t%mod;
    ll c=0,d=0;
    for(int i=0;i<=m+1;i++)
    {
        c=(c+a[m+1-i][1]*C(m+1,i)*((i&1)?1:-1))%mod;
        d=(d+a[m+1-i][0]*C(m+1,i)*((i&1)?1:-1))%mod;
    }
    g[0]=(-d)*Pow(c,mod-2)%mod;
    for(int i=1;i<=m+1;i++)g[i]=(g[i-1]+f[i-1])*t%mod;
    ll ans=G(n+1); 
    ans=(ans*Pow(m,n+1)%mod-g[0]+mod)%mod;
    printf("%lld\n",ans);
    return 0;
}