「算法笔记」拉格朗日插值

一、基本内容

考虑一个 \(n\) 次多项式 \(f(x)=\sum_{i=0}^{n}a_ix^i\)

若已知 \(n+1\) 个点值,则可构造出多项式。具体来说,拉格朗日插值要解决的问题是:

给定 \(n+1\)\((x_i,y_i)\),其中 \(y_i=f(x_i)\)。对于某个给定的 \(k\),求 \(f(k)\) 的值。

拉格朗日插值的结论是:

\[f(k)=\sum_{i=1}^{n+1}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j} \]

可以感性理解其正确性:

\(n+1\) 个点值代入即可检验。当 \(k\) 等于某个 \(x_t\) 时:

  • \(i\neq t\),则存在一个 \(j\) 使得 \(j=t\)

    对于乘积项的分子 \(\prod\limits_{j\neq i}k-x_j\)\(j=t\)\(k-x_j=0\)。因此 \(i\neq t\) 的项无贡献。

  • \(i=t\),乘积项变为 \(\prod\limits_{j\neq i}\frac{x_i-x_j}{x_i-x_j}=1\)

因此 \(f(k)=y_t\),这与题目条件是相符的。

暴力计算这个式子,时间复杂度 \(\mathcal O(n^2)\)

//Luogu P4781
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e3+5,mod=998244353;
int n,k,x[N],y[N],ans;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
signed main(){
    scanf("%lld%lld",&n,&k);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&x[i],&y[i]);
    for(int i=1;i<=n;i++){
        int s1=1,s2=1;
        for(int j=1;j<=n;j++)
            if(j!=i) s1=s1*(k-x[j])%mod,s2=s2*(x[i]-x[j])%mod;
        ans=(ans+y[i]*s1%mod*mul(s2,mod-2,mod)%mod)%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
} 

二、随时加点查询

Problem:LOJ#165. 拉格朗日插值

再写一遍式子:

\[f(k)=\sum_{i=1}^{n+1}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j} \]

\(\mathcal O(n)\) 次询问,肯定不能每次询问都按拉格朗日插值的式子 \(\mathcal O(n^2)\) 暴力算一遍。

考虑在询问时还是枚举 \(i\)。预先维护好求积式中分母、分子的值。

  • 分母:对于每个 \(i\),维护一个 \(g_i=\prod_{j\neq i}(x_i-x_j)\),在加点时顺便更新。
  • 分子:在询问时,先特判 \(k\) 等于其中某个 \(x_t\) 的情况(\(f(k)=y_t\))。否则维护 \(s=\prod_{j=1}^{n+1}(k-x_j)\),在枚举 \(i\) 时除以 \((k-x_i)\) 即可。

于是,\(f(k)=\sum_{i=1}^{n+1}y_i\cdot \frac{1}{g_i}\cdot \frac{s}{k-x_i}=s\sum_{i=1}^{n+1}y_i\cdot \frac{1}{g_i(k-x_i)}\)。这样,若不管求逆元的复杂度,单次加点、查询的复杂度都是 \(\mathcal O(n)\),总复杂度 \(\mathcal O(n^2)\)

//LOJ 165
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=3e3+5,mod=998244353;
int n,opt,x,y,k,tot,qx[N],qy[N],g[N],s,res;
bool flag;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans; 
}
signed main(){
    scanf("%lld",&n);
    while(n--){
        scanf("%lld",&opt);
        if(opt==1){
            scanf("%lld%lld",&x,&y);
            ++tot,qx[tot]=x,qy[tot]=y,g[tot]=1;
            for(int i=1;i<tot;i++){
                g[i]=g[i]*(qx[i]-qx[tot]+mod)%mod;
                g[tot]=g[tot]*(qx[tot]-qx[i]+mod)%mod;
            } 
        }
        else{
            scanf("%lld",&k),s=1,res=0,flag=0;
            for(int i=1;i<=tot;i++){
                if(k==qx[i]){printf("%lld\n",qy[i]),flag=1;break;} 
                s=s*(k-qx[i]+mod)%mod,res=(res+qy[i]*mul(g[i]*(k-qx[i]+mod)%mod,mod-2,mod)%mod)%mod;
            }
            if(!flag) printf("%lld\n",s*res%mod);
        }
    }
    return 0;
}

三、取值连续

\(x\) 取值为连续的 \(n+1\) 个自然数时,可以把拉格朗日插值算法的时间复杂度优化至 \(\mathcal O(n)\)

\[f(k)=\sum_{i=1}^{n+1}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j} \]

首先,对于分子,预处理出 \(k-x_i\) 的前缀积 \(pre_i=\prod_{j=1}^{i}(k-x_j)\)、后缀积 \(suf_i=\prod_{j=i}^{n+1}(k-x_j)\),那么 \(\prod_{j\neq i}(k-x_j)=pre_{i-1}\cdot suf_{i+1}\)

对于分母,因为 \(x\) 取值连续,所以 \(\prod_{j\neq i}(x_i-x_j)=\prod_{j\neq i}(i-j)\)。分 \(j<i\)\(j>i\) 讨论,得 \(\prod_{j\neq i}(i-j)=(i-1)!\cdot(-1)^{(n+1)-i}\cdot((n+1)-i)!\)。此时有:

\[f(k)=\sum_{i=1}^{n+1}y_i\frac{pre_{i-1}\cdot suf_{i+1}}{(i-1)!\cdot(-1)^{n+1-i}\cdot(n+1-i)!} \]

预处理阶乘,按这个式子直接算,时间复杂度 \(\mathcal O(n)\)

int calc(int n,int k,int qx[N],int qy[N]){    //f(k)
    fac[0]=pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod,fac[i]=fac[i-1]*i%mod;
    for(int i=n+1;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=n+1;i++){
        int s=fac[i-1]*((n+1-i)&1?mod-1:1)%mod*fac[n+1-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}

四、自然数幂求和

CF622F The Sum of the k-th Powers

给定 \(n,k\),求 \(S_k(n)=\sum_{i=1}^n i^k\),对 \(10^9+7\) 取模。

\(1\leq n\leq 10^9,0\leq k\leq 10^6\)

计算 \(S_k(n)\) 的方法有多种,各有各的优缺点。拉格朗日插值法,一次只能求出一个 \(S_k(n)\),但时间复杂度只有 \(\mathcal O(k)\),简单好写。

一个结论:\(S_k(n)\) 为关于 \(n\)\(k+1\) 次多项式(将 \(n\) 看作变量,\(k\) 看作定值)。例如 \(\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6}\) 是一个关于 \(n\)\(3\) 次多项式。

那么代入 \(k+2\) 个点值 \(x_1,x_2,\cdots,x_{k+2}\),算出 \(S_k(x_1),S_k(x_2),\cdots,S_k(x_{k+2})\),然后用拉格朗日插值法直接求出 \(S_k(n)\) 即可。

任意取 \(x_i\) 复杂度为 \(\mathcal O(k^2)\),但我们可以取 \(k+2\) 个连续的自然数 \(1,2,\dots,k+2\),复杂度 \(\mathcal O(k)\)

预处理 \(k+2\) 个点时,可以用线性筛求 \(i^k\)(这是完全积性函数)。在质数处暴力快速幂,而因为质数只有 \(\mathcal O(\frac{k}{\log k})\) 个,所以预处理的复杂度也是 \(\mathcal O(\frac{k}{\log k}\cdot \log k)=\mathcal O(k)\)。代码里直接暴力算了 QAQ。

//CF622F 
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e6+5,mod=1e9+7;
int n,k,pre[N],suf[N],fac[N],qx[N],qy[N];
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans; 
}
int calc(int n,int k,int qx[N],int qy[N]){
    fac[0]=pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod,fac[i]=fac[i-1]*i%mod;
    for(int i=n+1;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=n+1;i++){
        int s=fac[i-1]*((n+1-i)&1?mod-1:1)%mod*fac[n+1-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}
signed main(){
    scanf("%lld%lld",&n,&k);
    for(int i=1;i<=k+2;i++)
        qx[i]=i,qy[i]=(qy[i-1]+mul(i,k,mod))%mod;
    printf("%lld\n",calc(k+1,n,qx,qy));
    return 0;
}

五、求系数

考虑一个 \(n\) 次多项式 \(f(x)=\sum_{i=0}^{n}a_ix^i\)。给定 \(n+1\) 个点值,我们不满足于对于某个 \(k\)\(f(k)\) 的值,还想求出 \(a_0,a_1,\cdots,a_n\)

\[f(k)=\sum_{i=1}^{n+1}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j}\Rightarrow f(k)=\sum_{i=1}^{n+1}\frac{y_i}{\prod_{j\neq i}x_i-x_j}\prod_{j\neq i}(k-x_j) \]

上式中只有 \(k\) 是变量,其他值都是常量。我们要把式子写成 \(f(k)=a_0+a_1k+a_2k^2+\cdots a_nk^n\) 的形式。

\(\frac{y_i}{\prod_{j\neq i}x_i-x_j}\) 作为常量可以直接求出来。我们要把 \(\prod\limits_{j\neq i}(k-x_j)\) “展开”成一个关于 \(k\) 的求和的形式。考虑先将 \(\prod\limits_{j=1}^{n+1}(k-x_j)\) 写成一个关于 \(k\)\(n+1\) 次的多项式(预处理),然后对于每个 \(i\),在这个 \(n+1\) 次多项式的基础上除以 \((k-x_i)\)

乘以 \((k-x_i)\) 和除以 \((k-x_i)\) 的操作都能 \(\mathcal O(n)\) 实现。具体地,考虑生成函数 \(H(k)\)\(G(k)\),其中 \(H(k)=G(k)(k-c)\)(注意 \(c\) 在此处是一个常数,\(c=x_i\)。而 \(k\) 是变量)。那么,乘以 \((k-c)\) 相当于 \(G\to H\),除以 \((k-c)\) 相当于 \(H\to G\)。考虑生成函数的 \(k^i\) 项,可知 \(h_i=g_{i-1}-cg_i\)\(g_{i-1}=h_i+cg_i\),直接递推即可。

总时间复杂度 \(\mathcal O(n^2)\)

六、Problem

1. Luogu P4593 [TJOI2018] 教科书般的亵渎

P4593 [TJOI2018]教科书般的亵渎

注意:题目中的 \(k\) 是全局需要的总张数,不是当前还需要几张。并且每张卡的效果可能被多次施放(直到没有怪死为止),但每个怪在每张卡的回合内只会对答案产生一次贡献,而不是每次施放都有贡献。

发现若血量的区间 \([1,x]\) 连续,则只需要一张亵渎就可以杀死区间 \([1,x]\) 内所有怪物,所以 \(k = m+1\)。答案为:

\[\sum_{i=1}^k\left(\sum_{j=a_{i-1}+1}^{n}(j-a_{i-1})^k-\sum_{j=i}^m(a_j-a_{i-1})^k\right)=\sum_{i=1}^k\left(\sum_{j=1}^{n-a_{i-1}}j^k-\sum_{j=i}^m(a_j-a_{i-1})^k\right) \]

\(\sum_{j=1}^{n-a_{i-1}}j^k\) 可以用拉格朗日插值算(自然数幂求和),\(\sum_{j=i}^m(a_j-a_{i-1})^k\) 直接暴力算就好了。

//Luogu P4593
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e6+5,mod=1e9+7;
int t,n,m,k,a[N],pre[N],suf[N],fac[N],qx[N],qy[N],ans=0;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans; 
}
int calc(int n,int k,int qx[N],int qy[N]){
    fac[0]=pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod,fac[i]=fac[i-1]*i%mod;
    for(int i=n+1;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=n+1;i++){
        int s=fac[i-1]*((n+1-i)&1?mod-1:1)%mod*fac[n+1-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}
signed main(){
    scanf("%lld",&t);
    while(t--){ 
        scanf("%lld%lld",&n,&m),ans=0;
        for(int i=1;i<=m;i++)
            scanf("%lld",&a[i]);
        sort(a+1,a+1+m),m=unique(a+1,a+1+m)-a-1,k=m+1;
        for(int i=1;i<=k+2;i++)
            qx[i]=i,qy[i]=(qy[i-1]+mul(i,k,mod))%mod;
        for(int i=1;i<=k;i++){
            ans=(ans+calc(k+1,n-a[i-1],qx,qy))%mod;
            for(int j=i;j<=m;j++)
                ans=(ans-mul(a[j]-a[i-1],k,mod))%mod;
        }
        printf("%lld\n",(ans+mod)%mod);
    } 
    return 0;
}

2. LOJ 6024 XLkxc

LOJ#6024. XLkxc

给定 \(k,a,n,d\),求 \(\sum_{i=0}^n\sum_{j=1}^{a+i\cdot d}\sum_{l=1}^jl^k\bmod p\)

\(k\leq 123,a,n,d<p=1234567891\)

\(f(n)=\sum_{i=1}^n i^k\)。这是自然数幂求和,是一个关于 \(n\)\(k+1\) 次多项式。

\(g(n)=\sum_{i=1}^n f(i)\)。因为它差分之后是 \(f\),所以这是一个 \(k+2\) 次多项式。

同理最后我们要求的是一个 \(k+3\) 次多项式。

\[\begin{aligned} \sum_{i=0}^k a_ix^i-\sum_{i=0}^k a_i(x+c)^i&=\left(\sum_{i=0}^{k-1}a_i(x^i-(x+c)^i)\right)+a_k(x^k-(x+c)^k) \\&=\left(\sum_{i=0}^{k-1}a_i(x^i-(x+c)^i)\right)+a_k(x-(x+c))(x^{k-1}+\\&\quad \quad x^{k-2}(x+c)+x^{k-3}(x+c)^2+\cdots+(x+c)^{k-1}) \end{aligned} \]

用了 \(k\) 次方差公式。

最初的 \(\sum_{i=0}^k a_ix^i\) 是一个 \(k\) 次多项式。最后得出的式子中,左边是 \(k-1\) 次的,右边也是 \(k-1\) 次的。因此 \(\sum_{i=0}^k a_ix^i-\sum_{i=0}^k a_i(x+c)^i\) 是一个 \(k-1\) 次多项式。

所以直接用拉格朗日插值算出 \(g\),再根据 \(g\) 算出答案即可。

//LOJ 6024
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=140,mod=1234567891;
int t,k,a,n,d,b[N],c[N],pre[N],suf[N],fac[N],qx[N],qy[N],ans=0;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans; 
}
int calc(int n,int k,int qx[N],int qy[N]){
    fac[0]=pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod,fac[i]=fac[i-1]*i%mod;
    for(int i=n+1;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=n+1;i++){
        int s=fac[i-1]*((n+1-i)&1?mod-1:1)%mod*fac[n+1-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}
signed main(){
    scanf("%lld",&t);
    for(int i=1;i<=130;i++) qx[i]=i;
    while(t--){
        scanf("%lld%lld%lld%lld",&k,&a,&n,&d);
        for(int i=1;i<=k+3;i++) b[i]=(b[i-1]+mul(i,k,mod))%mod;
        for(int i=2;i<=k+3;i++) b[i]=(b[i-1]+b[i])%mod;
        for(int i=0;i<=k+4;i++)
            c[i]=((i?c[i-1]:0)+calc(k+2,(a+i*d%mod)%mod,qx,b))%mod;
        printf("%lld\n",calc(k+3,n,qx,c));
    }
    return 0;
}

3. Luogu P3270 [JLOI2016]成绩比较

P3270 [JLOI2016]成绩比较

\(f_{i,j}\) 表示考虑了前 \(i\) 门课,被 B 神全面碾压的同学数为 \(j\) 的方案数。下式中,\(d_i\) 表示在第 \(i\) 门课上分配成绩的方案数。

\[f_{i,j}=\sum_{k=j}^nf_{i-1,k}\binom{k}{k-j}\binom{n-k-1}{r_i-1-(k-j)}d_i \]

意思是,被 B 神全面碾压的人数从 \(k\) 减为了 \(j\),说明要从(原来被全面碾压的)\(k\) 个人中,选出 \((k-j)\) 个人第 \(i\) 门课成绩比 B 神高;又由于本身有 \((r_i-1)\) 个人第 \(i\) 门课比 B 神高,所以要从剩下的(没有被 B 神全面碾压的)\((n-k-1)\) 个人中选 \(r_i-1-(k-j)\) 个。

现在考虑计算 \(d_i\)。枚举 B 神的分数:

\[d_i=\sum_{j=1}^{u_i}j^{n-r_i}(u_i-j)^{r_i-1} \]

但是 \(u_i\)\(10^9\) 级别的,显然不能直接枚举。考虑将 \(d_i\) 看作关于 \(u_i\) 的多项式(将 \(u_i\) 看作变量,\(r_i\) 看作定值),那么根据自然数幂和的结论,这个多项式的次数是不超过 \(r_i\) 的(也就是不超过 \(n\) 的)。暴力算出 \((n+1)\) 个点值然后拉格朗日插值即可。

总时间复杂度 \(\mathcal O(n^2m)\)

//Luogu P3270
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=110,mod=1e9+7;
int n,m,K,u[N],r[N],c[N][N],f[N][N],pre[N],suf[N],fac[N],qx[N],qy[N];
int C(int n,int m){
    if(n<m||n<0||m<0) return 0;
    return c[n][m];
}
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
int calc(int n,int k){
    fac[0]=pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod,fac[i]=fac[i-1]*i%mod;
    for(int i=n+1;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=n+1;i++){
        int s=fac[i-1]*((n+1-i)&1?mod-1:1)%mod*fac[n+1-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}
int get(int u,int r){
    for(int i=1;i<=n+1;i++)
        qx[i]=i,qy[i]=mul(i,n-r,mod)*mul(u-i+mod,r-1,mod)%mod;
    for(int i=1;i<=n+1;i++) qy[i]=(qy[i-1]+qy[i])%mod;
    return calc(n,u);
}
signed main(){
    scanf("%lld%lld%lld",&n,&m,&K),c[0][0]=1;
    for(int i=1;i<=n;i++){
        c[i][0]=1;
        for(int j=1;j<=i;j++)
            c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
    }
    for(int i=1;i<=m;i++)
        scanf("%lld",&u[i]);
    for(int i=1;i<=m;i++)
        scanf("%lld",&r[i]);
    f[0][n-1]=1;
    for(int i=1;i<=m;i++){
        int val=get(u[i],r[i]);
        for(int j=K;j<n;j++)
            for(int k=j;k<n;k++)
                f[i][j]=(f[i][j]+f[i-1][k]*C(k,k-j)%mod*C(n-k-1,r[i]-1-(k-j))%mod*val%mod)%mod;
    }
    printf("%lld\n",f[m][K]);
    return 0;
} 

4. BZOJ 2655 calc

BZOJ#2655. calc。可参考 这篇博客

朴素 DP:设 \(f_{i,j}\) 表示在 \([1,j]\) 中选出 \(i\) 个互不相同的数组成集合的值的和(无序)。

\[f_{i,j}=f_{i,j-1}+f_{i-1,j-1}\times j \]

答案为 \(f_{n,A}\times n!\)。猜想 \(f_{i,j}\) 是关于 \(j\) 的多项式,即 \(f_{i,j}=F_i(j)\)\(F_i(x)=F_i(x-1)+F_{i-1}(x-1)\times x\)

\(F_i(x)\) 的最高次数为 \(N(i)\),那么 \(F_{i-1}(x-1)\times x\) 的次数显然为 \(N(i-1)+1\)。而 \(F_i(x)\) 相当于是 \(F_{i-1}(x-1)\times x\) 的前缀和,所以 \(N(i)=N(i-1)+2\)

边界 \(F_0(x)\) 是常数 \(1\),次数为 \(0\),所以 \(N(i)=2i\),即 \(F_i(x)\) 是关于 \(x\)\(2i\) 次多项式。

所以我们只需要取 \(f_{n,1\dots 2n+1}\)\(2n+1\) 个点值,然后用拉格朗日插值求出 \(f_{n,A}\) 即可。

//BZOJ 2655
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e3+5;
int A,n,mod,up,fac[N],pre[N],suf[N],qx[N],qy[N],f[N][N];
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans; 
}
int calc(int k){
    if(k<=up) return qy[k];
    pre[0]=suf[up+1]=1;
    for(int i=1;i<=up;i++) pre[i]=pre[i-1]*(k-qx[i]+mod)%mod;
    for(int i=up;i>=1;i--) suf[i]=suf[i+1]*(k-qx[i]+mod)%mod;
    int ans=0;
    for(int i=1;i<=up;i++){
        int s=fac[i-1]*((up-i)&1?mod-1:1)%mod*fac[up-i]%mod;
        ans=(ans+qy[i]*(pre[i-1]*suf[i+1]%mod)%mod*mul(s,mod-2,mod)%mod)%mod;
    }
    return ans;
}
signed main(){
    scanf("%lld%lld%lld",&A,&n,&mod),up=n*2+1,fac[0]=1;
    for(int i=1;i<=up;i++) fac[i]=fac[i-1]*i%mod;
    for(int i=0;i<=up;i++) f[0][i]=1;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=up;j++)
            f[i][j]=(f[i][j-1]+f[i-1][j-1]*j%mod)%mod;
    for(int i=1;i<=up;i++)
        qx[i]=i,qy[i]=f[n][i];
    printf("%lld\n",calc(A)*fac[n]%mod);
    return 0;
}

5. BZOJ 4126 国王奇遇记

BZOJ#4126. 国王奇遇记加强版之再加强版

给定 \(n,m\),求 \(\sum_{i=1}^n i^mm^i\)

\(1\leq n\leq 10^9,1\leq m\leq 5\times 10^5\)

拉格朗日插值 + 高阶差分。鸽了。

posted @ 2021-04-16 11:29  maoyiting  阅读(149)  评论(0编辑  收藏  举报