拉格朗日插值小记

拉格朗日插值

基本式子

\[f_k=\sum\limits_{i=0}^n y_i\prod\limits_{i\not=j}\frac{k-x_j}{x_i-x_j} \]

当取点连续为时

\[f_k=\sum\limits_{i=0}^nf_i\prod\limits_{i\not=j}\frac{k-j}{i-j} \]

观察发现，\(k-j,i-j\) 是连续一段，可以预处理为维护。

记 \(pre_i=pre_{i-1}+k-i,suf_i=suf_{i+1}+k-i,fac_i=fac_{i-1}\times i\)

则可得：

\[f_k=\sum\limits_{i=0}^nf_i\frac{pre_{i-1}\times suf_{i+1}}{fac_{i-1}\times fac_{n-i}} \]

注意当 \(n-i\) 为奇数时，分母应该取负。

时间复杂度 \(O(n)\)。

当然还有重心拉格朗日插值，但用的少；至于多项式快速插值，太毒瘤了（

一些题目

CF622F The Sum of the k-th Powers

结论：\(\sum\limits_{i=1}^ni^k\) 是关于 \(n\) 的 \(k+1\) 次多项式。

直接插值维护。

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
const int N=1e6+10;
int n,k,y,ans;
int l[1000005],r[1000005],jc[1000005];
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;y>>=1;
    }return res;
}
signed main() 
{
    
    scanf("%d%d", &n, &k);l[0]=r[k+3]=jc[0]=1;
    for(int i=1;i<=k+2;++i)l[i]=l[i-1]*(n-i)%mod;
    for(int i=k+2;i>=1;--i)r[i]=r[i+1]*(n-i)%mod;
    for(int i=1;i<=k+2;i++)jc[i]=jc[i-1]*i%mod;
    for(int i=1;i<=k+2;i++)
    {
        y=(y+qpow(i,k))%mod;
        int a=l[i-1]*r[i+1]%mod;
        int b=jc[i-1]*((k-i)&1?-1:1)*jc[k+2-i]%mod;
        ans=(ans+1ll*y*a%mod*qpow(b,mod-2)%mod)%mod;
    }write((ans+mod)%mod),pc('\n');
    return 0;
}

P4593 [TJOI2018]教科书般的亵渎

每有一个空位 \(p\)，则会产生 \(\sum\limits_{i=1}^{n-p}i^k\) 的贡献去掉空位的贡献。

容斥处理，先计算所有 \(\sum\limits_{i=1}^{n-p}i^k\)，最后空位贡献统一处理。

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
int T,n,m,k,ans;
int a[55],f[55];
int pre[55],suf[55],jc[55];
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;y>>=1;
    }return res;
}
int inv(int x){return qpow(x,mod-2);}
map<int,bool>mp;
int calc(int x)
{
    if(x<=k+2)return f[x];
    int res=0;pre[0]=suf[k+3]=1;
    for(int i=1;i<=k+2;++i)pre[i]=pre[i-1]*(x-i)%mod;
    for(int i=k+2;i>=1;--i)suf[i]=suf[i+1]*(x-i)%mod;
    for(int i=1;i<=k+2;++i)
    {
        int res1=pre[i-1]*suf[i+1]%mod;
        int res2=jc[i-1]*jc[k+2-i]*(((k+2-i)&1)?-1:1)%mod;
        res=(res+f[i]*res1%mod*inv(res2)%mod)%mod;
    }return (res+mod)%mod;
}
signed main()
{
    T=read();jc[0]=1;
    for(int i=1;i<=52;++i)jc[i]=jc[i-1]*i%mod;
    while(T-->0)
    {
        n=read(),m=read();k=m+1;mp.clear();
        for(int i=1;i<=m;++i)a[i]=read(),mp[a[i]]=true;
        sort(a+1,a+m+1);while(mp[n])--n,--k,--m;
        for(int i=1;i<=k+2;++i)f[i]=(f[i-1]+qpow(i,k))%mod;
        ans=calc(n);for(int i=1;i<=m;++i)ans=(ans-qpow(a[i],k))%mod;
        for(int i=1;i<=m;++i)ans=(ans+calc(n-a[i]))%mod;
        for(int i=1;i<=m;++i)
            for(int j=i-1;j>=1;--j)
                ans=(ans-qpow(a[i]-a[j],k))%mod;
        write((ans+mod)%mod),pc('\n');
    }return 0;
}

CF995F Cowmpany Cowmpensation

首先容易想到 \(O(nD)\) 的 DP，记 \(f_{u,i}\) 为 \(u\) 节点工资为 \(i\) 时的子树方案数量。

\[f_{u,i}=\sum_v\sum_{j=1}^{i}f_{v,j} \]

后面可以优化，记 \(g_{v,i}=\sum\limits_{j=1}^if_{v,j}\)，则 \(f_{u,i}=\sum\limits_vg_{v,i}\)。

叶子节点的 \(f_{u,i}=1\)，所以可以看成 \(0\) 次多项式。

由于 \(g_{u,i}-g_{u,i-1}=f_{u,i}=\sum\limits_v\sum\limits_{j=1}^if_{v,j}\)，所以不难发现 \(g_{u,i}\) 是关于 \(i\) 的 \(siz_u\) 次多项式。

对于 \(g_{1,i}\)，次数为 \(siz_1\)，要插值则维护 \(siz_1+1\) 个点值。答案即为插值出来的 \(g_{1,D}\)。

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
int n,d,ans,siz[3005],f[3005][3005],g[3005][3005];
int pre[3005],suf[3005],fac[3005];
vector<int>e[3005];
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;y>>=1;
    }return res;
}
void dfs(int x)
{
    siz[x]=1;
    for(int y:e[x]){dfs(y);siz[x]+=siz[y];}
    for(int i=1;i<=n+1;++i)
    {
        int res=1;g[x][i]=1;
        for(int y:e[x])res=res*g[y][i]%mod;
        g[x][i]=(g[x][i-1]+res)%mod;
    }
}
signed main()
{
    n=read(),d=read();
    for(int i=2;i<=n;++i)e[read()].push_back(i);
    dfs(1);if(d<=n+1){write(g[1][d]),pc('\n');return 0;}
    pre[0]=suf[n+2]=fac[0]=1;
    for(int i=1;i<=n+1;++i)
    {
        pre[i]=pre[i-1]*(d-i+mod)%mod;
        fac[i]=fac[i-1]*i%mod;
    }
    for(int i=n+1;i>=1;--i)suf[i]=suf[i+1]*(d-i+mod)%mod;
    for(int i=1;i<=n+1;++i)
    {
        int res=fac[i-1]*fac[n+1-i]%mod;res=(n-i)&1?res:mod-res;
        res=qpow(res,mod-2);res=res*pre[i-1]%mod*suf[i+1]%mod;
        ans=(ans+g[1][i]*res)%mod;
    }write(ans),pc('\n');
    return 0;
}

P3270 [JLOI2016]成绩比较

记 \(f_{i,j}\) 是 B神 在前 \(i\) 门课比 \(j\) 个人强。

\[f_{i,j}=\sum_{k=j}^{n}f_{i-1,k}\times {k \choose j}\times {n-k-1 \choose r_i-1-(k-j)}\times g_i \]

\[g_i=\sum\limits_{j=1}^{u_i}j^{n-r_i}\times (u_i-j)^{r_i-1} \]

可得：\(G(x)=\sum\limits_{j=1}^{x}j^{n-r_i}\times (u_i-j)^{r_i-1}\)

\(G(x)\) 是关于 \(x\) 的不超过 \(n\) 次的多项式，维护 \(n+1\) 个点值后暴力插值求出 \(g_i\)。

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

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
int n,m,K ,u[105],r[105];
int pre[105],suf[105],fac[105];
int f[105][105],g[105],G[105],C[105][105];
int c(int x,int y){return (y<0||y>x)?0:C[x][y];}
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;y>>=1;
    }return res;
}
int calc(int u,int r)
{
    int res=0;
    for(int i=1;i<=n+1;++i)G[i]=0;
    for(int i=1;i<=n+1;++i)
            G[i]=(G[i-1]+qpow(i,n-r)*qpow(u-i+mod,r-1))%mod;
    pre[0]=suf[n+2]=1;
    for(int i=1;i<=n+1;++i)pre[i]=pre[i-1]*(u-i+mod)%mod;
    for(int i=n+1;i>=1;--i)suf[i]=suf[i+1]*(u-i+mod)%mod;
    for(int i=1;i<=n+1;++i)
    {
        int tmp=fac[i-1]*fac[n+1-i]%mod;tmp=(n-i)&1?tmp:mod-tmp;
        tmp=qpow(tmp,mod-2);tmp=tmp*pre[i-1]%mod*suf[i+1]%mod;
        res=(res+G[i]*tmp)%mod;
    }return res;
}
signed main()
{
    n=read(),m=read(),K=read();fac[0]=f[0][n-1]=C[0][0]=1;
    for(int i=1;i<=m;++i)u[i]=read();
    for(int i=1;i<=m;++i)r[i]=read();
    for(int i=1;i<=n;C[i][0]=1,++i)
        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<=n+1;++i)fac[i]=fac[i-1]*i%mod;
    for(int i=1;i<=m;++i)g[i]=calc(u[i],r[i]);
    for(int i=1;i<=m;++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,j)%mod*c(n-k-1,r[i]-1-(k-j))%mod*g[i])%mod;
    write(f[m][K]),pc('\n');
    return 0;
}

P4463 [集训队互测 2012] calc

不妨设 \(a_i<a_{i+1}\)，那么答案即为 \(f_{n,k}\times n!\)。

设 \(f_{i,j}\) 为前 \(i\) 个数在 \([1,j]\) 范围内的方案乘积和。

可得：\(f_{i,j}=f_{i,j-1}+f_{i-1,j-1}\times j\)

假设 \(f_{n,j}\) 是关于 \(j\) 的 \(g(n)\) 次多项式。

因为 \(f_{n,j}-f_{n,j-1}=f_{n-1,j-1}\times j\)，所以 \(g(n)-1=g(n-1)+1\)。

又 \(g(0)=0\)，所以 \(g(n)=2n\)。

对 \(f_{n,j}\) 维护 \(2n+1\) 个点值，插值求出 \(f_{n,k}\)。

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
int n,k,p,ans;
int f[505][1005],pre[1005],suf[1005],fac[1005];
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%p;
        x=x*x%p;y>>=1;
    }return res;
}
signed main()
{
    k=read(),n=read(),p=read();
    for(int i=0;i<=2*n+1;++i)f[0][i]=1;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=2*n+1;++j)
            f[i][j]=(f[i][j-1]+f[i-1][j-1]*j)%p;
    fac[0]=pre[0]=suf[2*n+2]=1;
    for(int i=1;i<=2*n+1;++i)fac[i]=fac[i-1]*i%p;
    for(int i=1;i<=2*n+1;++i)pre[i]=pre[i-1]*(k-i+p)%p;
    for(int i=2*n+1;i>=1;--i)suf[i]=suf[i+1]*(k-i+p)%p;
    for(int i=1;i<=2*n+1;++i)
    {
        int res=fac[i-1]*fac[2*n+1-i]%p;res=(2*n-i)&1?res:p-res;
        res=qpow(res,p-2)*pre[i-1]%p*suf[i+1]%p;
        ans=(ans+res*f[n][i])%p;
    }write(ans*fac[n]%p),pc('\n');
    return 0;
}

P3643 [APIO2016] 划艇

组合法

值域 \(10^9\)，直接 DP 必爆，先离散化所有区间。

记 \(f_{i,j}\) 为前 \(i\) 个学校，第 \(i\) 个学校派出的划艇数在区间 \(j\) 中。

那么这 \(i\) 个学校分为 \(2\) 种：划艇数在区间 \(j\)，和不在区间 \(j\)。

不在区间 \(j\) 可以直接转移，在区间 \(j\) 需要特殊处理。

注意到，从区间 \([0,m]\) 中取 \(n\) 个数，非零数严格递增，则方案数为 \(m+n\choose n\)。

而取到第 \(i\) 个 \(0\) 其实就代表着第 \(i\) 个学校不参赛。

那么转移方程为：

\[f_{i,j}= \sum\limits_{p=0}^{i-1}\sum\limits_{k=1}^{j-1}{L+m-1\choose m}f_{p,k} \]

注意到 \(L,m\) 不随 \(k\) 变化而变化，所以可以前缀和优化下。

\[g_{i,j}=\sum\limits_{k=1}^{j-1}f_{i,k} \]

\[ f_{i,j}=\sum\limits_{p=0}^{i-1}{L+m-1\choose m}g_{p,j-1} \]

然后优化空间，发现第一层按照 \(j\) 枚举，第二层按照 \(i\) 倒序枚举，可以不用开 \(f\) 数组，直接用前缀和数组 \(g\) 就行。第三层倒序枚举 \(p\) ，组合数的 \(m\) 是递增的；在第一层预处理出当前要用的所有组合数就行。

时间复杂度 \(O(n^3)\)，空间复杂度 \(O(n)\)。

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
int n,m,ans,g[505];
struct node
{
    int l,r;
}a[505];
int b[1005],tot;
int inv[505],C[505];
signed main()
{
    n=read();inv[1]=1;
    for(int i=2;i<=n;++i)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for(int i=1;i<=n;++i)
    {
        a[i].l=read(),a[i].r=read()+1;
        b[++tot]=a[i].l;b[++tot]=a[i].r;
    }sort(b+1,b+tot+1);tot=unique(b+1,b+tot+1)-(b+1);
    for(int i=1;i<=n;++i)
    {
        a[i].l=lower_bound(b+1,b+tot+1,a[i].l)-b;
        a[i].r=lower_bound(b+1,b+tot+1,a[i].r)-b;
    }C[0]=g[0]=1;
    for(int j=1;j<=tot;++j)//区间
    {
        int L=b[j+1]-b[j];
        for(int i=1;i<=n;++i)C[i]=C[i-1]*(L+i-1)%mod*inv[i]%mod;
        for(int i=n;i>=1;--i)//学校
            if(a[i].l<=j&&j+1<=a[i].r)
            {
                int f=0,m=1;
                for(int p=i-1;p>=0;--p)
                {
                    f=f+C[m]*g[p]%mod;
                    if(a[p].l<=j&&j+1<=a[p].r)++m;
                }g[i]=(g[i]+f)%mod;
            }
    }for(int i=1;i<=n;++i)ans=(ans+g[i])%mod;
    write(ans),pc('\n');
    return 0;
}

插值法

记 \(f_{i,j}\) 为前 \(i\) 个学校（可选可不选），选择的最后一个学校划艇数为 \(j\) 的方案数。

那么转移方程为：

\[f_{i,j}=f_{i-1,j}+\sum_{k=1}^{j-1}f_{i-1,k} \]

若把第二维看作自变量的话，可以发现其实是分段多项式。

每次前缀和，多项式次数 \(+1\)，所以每一段的次数最多为 \(n\)。

我们可以将每一段离散后进行处理，每一段维护 \(n+1\) 个点值，插值求出每一段的答案。

时间复杂度 \(O(n^3)\)，空间复杂度 \(O(n^2)\)。

超级大常数，比组合法慢将近20倍（

代码：

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
const int mod=1e9+7;
int n,l[505],r[505],a[1005],tot;
int fac[505],pre[505],suf[505];
int sum[1005],g[1005][505];
int val[1005];
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;y>>=1;
    }return res;
}
int LGRG(int k,int *y)
{
    pre[0]=suf[n+2]=1;int res=0;
    for(int i=1;i<=n+1;++i)pre[i]=pre[i-1]*(k-i+mod)%mod;
    for(int i=n+1;i>=1;--i)suf[i]=suf[i+1]*(k-i+mod)%mod;
    for(int i=1;i<=n+1;++i)
    {
        int tmp=fac[i-1]*fac[n+1-i]%mod;tmp=(n-i)&1?tmp:mod-tmp;
        tmp=tmp*pre[i-1]%mod*suf[i+1]%mod;
        res=(res+tmp*y[i])%mod;
    }return res;
}
signed main()
{
    n=read();
    for(int i=1;i<=n;++i)
    {
        l[i]=read(),r[i]=read()+1;
        a[++tot]=l[i],a[++tot]=r[i];
    }sort(a+1,a+tot+1);tot=unique(a+1,a+1+tot)-(a+1);
    for(int i=1;i<=n;++i) 
        l[i]=lower_bound(a+1,a+1+tot,l[i])-a,
        r[i]=lower_bound(a+1,a+1+tot,r[i])-a;
    fac[0]=1;for(int i=1;i<=n+1;++i)fac[i]=fac[i-1]*i%mod;
    for(int i=2;i<=n+1;++i)fac[i]=qpow(fac[i],mod-2);
    for(int i=0;i<=tot;++i)sum[i]=1;
    for(int i=1;i<=n;++i)
    {
        for(int j=l[i];j<r[i];++j)
        {
            for(int k=1;k<=n+1;++k)g[j][k]=(g[j][k]+sum[j-1])%mod;
            for(int k=2;k<=n+1;++k)g[j][k]=(g[j][k]+g[j][k-1])%mod;
            val[j]=LGRG(a[j+1]-a[j],g[j]);
        }for(int j=1;j<=tot;++j)sum[j]=(sum[j-1]+val[j])%mod;
    }write(sum[tot]-1),pc('\n');
    return 0;
}

未完待续...