拉格朗日插值小记

拉格朗日插值

基本式子

$$f_k=\sum _{i=0}^n{y}_i\prod _{i\ne j}\frac{k-x_j}{x_i-x_j}$$

当取点连续为时

$$f_k=\sum _{i=0}^n{f}_i\prod _{i\ne j}\frac{k-j}{i-j}$$

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

记 $pr{e}_i=pr{e}_{i-1}+k-i,su{f}_i=su{f}_{i+1}+k-i,fa{c}_i=fa{c}_{i-1}\times i$

则可得：

$$f_k=\sum _{i=0}^n{f}_i\frac{pr{e}_{i-1}\times su{f}_{i+1}}{fa{c}_{i-1}\times fa{c}_{n-i}}$$

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

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

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

一些题目

CF622F The Sum of the k-th Powers

结论：$\sum _{i=1}^n{i}^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 _{i=1}^{n-p}{i}^k$ 的贡献去掉空位的贡献。

容斥处理，先计算所有 $\sum _{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 _{j=1}^i{f}_{v,j}$，则 $f_{u,i}=\sum _v{g}_{v,i}$。

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

由于 $g_{u,i}-g_{u,i-1}=f_{u,i}=\sum _v\sum _{j=1}^i{f}_{v,j}$，所以不难发现 $g_{u,i}$ 是关于 $i$ 的 $si{z}_u$ 次多项式。

对于 $g_{1,i}$，次数为 $si{z}_1$，要插值则维护 $si{z}_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 (\genfrac{}{}{0ex}{}{k}{j})\times (\genfrac{}{}{0ex}{}{n-k-1}{r_i-1-(k-j)})\times g_i$$

$$g_i=\sum _{j=1}^{u_i}{j}^{n-r_i}\times (u_i-j{)}^{r_i-1}$$

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

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

时间复杂度 $O(n^{2}m)$。

代码：

#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$ 个数，非零数严格递增，则方案数为 $(\genfrac{}{}{0ex}{}{m+n}{n})$。

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

那么转移方程为：

$$f_{i,j}=\sum _{p=0}^{i-1}\sum _{k=1}^{j-1}(\genfrac{}{}{0ex}{}{L+m-1}{m})f_{p,k}$$

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

$$g_{i,j}=\sum _{k=1}^{j-1}{f}_{i,k}$$

$$f_{i,j}=\sum _{p=0}^{i-1}(\genfrac{}{}{0ex}{}{L+m-1}{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;
}

未完待续...