「算法笔记」拉格朗日插值
拉格朗日插值入门。
一、基本内容
考虑一个 \(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)\) 的值。
拉格朗日插值的结论是：
可以感性理解其正确性：
将 \(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. 拉格朗日插值。
再写一遍式子：
有 \(\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)\)。
首先，对于分子，预处理出 \(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)!\)。此时有：
预处理阶乘，按这个式子直接算，时间复杂度 \(\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\)。
上式中只有 \(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] 教科书般的亵渎
注意：题目中的 \(k\) 是全局需要的总张数，不是当前还需要几张。并且每张卡的效果可能被多次施放（直到没有怪死为止），但每个怪在每张卡的回合内只会对答案产生一次贡献，而不是每次施放都有贡献。
发现若血量的区间 \([1,x]\) 连续，则只需要一张亵渎就可以杀死区间 \([1,x]\) 内所有怪物，所以 \(k = m+1\)。答案为：
\(\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
给定 \(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]成绩比较
设 \(f_{i,j}\) 表示考虑了前 \(i\) 门课，被 B 神全面碾压的同学数为 \(j\) 的方案数。下式中，\(d_i\) 表示在第 \(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 神的分数：
但是 \(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_{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 国王奇遇记
给定 \(n,m\)，求 \(\sum_{i=1}^n i^mm^i\)。
\(1\leq n\leq 10^9,1\leq m\leq 5\times 10^5\)。
拉格朗日插值 + 高阶差分。鸽了。
