【BZOJ】3453: tyvj 1858 XLkxc 拉格朗日插值（自然数幂和）

【题意】给定k<=123，a,n,d<=10^9，求：

$$f(n)=\sum_{i=0}^{n}\sum_{j=1}^{a+id}\sum_{x=1}^{j}x^k$$

【算法】拉格朗日插值

【题解】参考：拉格朗日插值法及应用 by DZYO

虽然式子很复杂，但一点一点化简有条理的化简后就可以做了。

首先最后是一个自然数幂和：

$$\sum_{x=1}^{j}x^k$$

这是一个k+1次多项式，可以理解为k+一个Σ（一般一个Σ增加一次项）。

然后会发现最后部分和第二部分之间不需要插值，因为第二部分的前若干小项的计算只需要第一部分的前若干小项，那么：

$$g(m)=\sum_{j=1}^{m}\sum_{x=1}^{j}x^k$$

这是一个k+2次多项式，对g(m)可以O(k)计算前k+2项，因为后面的自然数幂和与j无关，所以可以处理成前缀和（不过对最终复杂度没有意义）。

剩余的部分：

$$f(n)=\sum_{i=0}^{n}g(a+id)$$

这是一个k+3次多项式，强制插值。对于前k+3项，每一项需要O(k)的插值运算，复杂度O(k^2)。

从这道题可以看出相邻插值嵌套的计算是O(k^2)，而且适用于多层嵌套复杂度不变。

这道题有一点比较坑，模数*2会爆int，要使用unsigned int。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define int unsigned int
using namespace std;
const int maxn=150,MOD=1234567891;
int f[maxn],g[maxn],a,b,k,n,mx,fac[maxn],fav[maxn];
int M(int x){return x>=MOD?x-MOD:x;}
int power(int x,int k){int ans=1;while(k){if(k&1)ans=1ll*ans*x%MOD;x=1ll*x*x%MOD;k>>=1;}return ans;}
int calc(int *g,int u){
    if(u<=mx)return g[u];//?
    int v=1;
    for(int i=1;i<=mx;i++)v=1ll*v*(u-i+MOD)%MOD;
    int ans=0;
    for(int i=1;i<=mx;i++){
        int t=1ll*fav[i-1]*fav[mx-i]%MOD;
        if((mx-i)&1)t=M(MOD-t);
        ans=M(ans+1ll*g[i]*v%MOD*power(M(u-i+MOD),MOD-2)%MOD*t%MOD);
    }
    return ans;
}
#undef int
int main(){
#define int unsigned int
    int T;scanf("%d",&T);
    fac[0]=1;for(int i=1;i<=130;i++)fac[i]=1ll*fac[i-1]*i%MOD;
    fav[130]=power(fac[130],MOD-2);for(int i=130;i>=1;i--)fav[i-1]=1ll*fav[i]*i%MOD;
    while(T--){
        scanf("%d%d%d%d",&k,&a,&n,&b);mx=k+5;
        g[0]=0;
        for(int i=1;i<=mx;i++){
            g[i]=g[i-1];
            for(int j=1;j<=i;j++){
                g[i]=M(g[i]+power(j,k));
            }
        }
        f[0]=calc(g,a);//!!!
        for(int i=1;i<=mx;i++){
            f[i]=M(f[i-1]+calc(g,M(a+1ll*i*b%MOD)));
        }
        printf("%d\n",calc(f,n));
    }
    return 0;
}