一、基本内容
考虑一个 (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 eq t)，则存在一个 (j) 使得 (j=t)。
对于乘积项的分子 (prodlimits_{j eq i}k-x_j)，(j=t) 时 (k-x_j=0)。因此 (i eq t) 的项无贡献。
若 (i=t)，乘积项变为 (prodlimits_{j eq 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 ",(ans+mod)%mod); return 0; }
二、随时加点查询
Problem：LOJ#165. 拉格朗日插值。
再写一遍式子：
有 (mathcal O(n)) 次询问，肯定不能每次询问都按拉格朗日插值的式子 (mathcal O(n^2)) 暴力算一遍。
考虑在询问时还是枚举 (i)。预先维护好求积式中分母、分子的值。
- 分母：对于每个 (i)，维护一个 (g_i=prod_{j eq 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_icdot frac{1}{g_i}cdot frac{s}{k-x_i}=ssum_{i=1}^{n+1}y_icdot 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 ",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 ",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 eq i}(k-x_j)=pre_{i-1}cdot suf_{i+1})。
对于分母，因为 (x) 取值连续，所以 (prod_{j eq i}(x_i-x_j)=prod_{j eq i}(i-j))。分 (j<i) 和 (j>i) 讨论，得 (prod_{j eq 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) 取模。
(1leq nleq 10^9,0leq kleq 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 ",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 eq i}x_i-x_j}) 作为常量可以直接求出来。我们要把 (prodlimits_{j eq i}(k-x_j)) “展开”成一个关于 (k) 的求和的形式。考虑先将 (prodlimits_{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 o H)，除以 ((k-c)) 相当于 (H o 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 ",(ans+mod)%mod); } return 0; }
2. LOJ 6024 XLkxc
给定 (k,a,n,d)，求 (sum_{i=0}^nsum_{j=1}^{a+icdot d}sum_{l=1}^jl^kmod p)。
(kleq 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) 次多项式。
[egin{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) ight)+a_k(x^k-(x+c)^k) \&=left(sum_{i=0}^{k-1}a_i(x^i-(x+c)^i) ight)+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 ",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 ",f[m][K]); return 0; }
4. BZOJ 2655 calc
BZOJ#2655. calc。可参考 这篇博客。
朴素 DP：设 (f_{i,j}) 表示在 ([1,j]) 中选出 (i) 个互不相同的数组成集合的值的和（无序）。
答案为 (f_{n,A} imes n!)。猜想 (f_{i,j}) 是关于 (j) 的多项式，即 (f_{i,j}=F_i(j))，(F_i(x)=F_i(x-1)+F_{i-1}(x-1) imes x)。
设 (F_i(x)) 的最高次数为 (N(i))，那么 (F_{i-1}(x-1) imes x) 的次数显然为 (N(i-1)+1)。而 (F_i(x)) 相当于是 (F_{i-1}(x-1) imes x) 的前缀和，所以 (N(i)=N(i-1)+2)。
边界 (F_0(x)) 是常数 (1)，次数为 (0)，所以 (N(i)=2i)，即 (F_i(x)) 是关于 (x) 的 (2i) 次多项式。
所以我们只需要取 (f_{n,1dots 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 ",calc(A)*fac[n]%mod); return 0; }
5. BZOJ 4126 国王奇遇记
给定 (n,m)，求 (sum_{i=1}^n i^mm^i)。
(1leq nleq 10^9,1leq mleq 5 imes 10^5)。
拉格朗日插值 + 高阶差分。鸽了。