参考资料:https://attack.blog.luogu.org/solution-p4781
因为高斯消元法是n立方的,有些鬼畜问题需要n平方的拉格朗日插值法。
通过n个不同的点直接构造出多项式。
$f(x)=\sum\limits_{i=1}^{n} y_i \prod\limits_{i\ne j}\frac{x-x_j}{x_i-x_j} $
正确性:代入 \(x_i\) 只有 \(y_i\) 右边是1,其他分子都有0,消掉了。
那么可以n平方构造。
需要小心溢出以及卡常数。假如进行同一个多项式的多次求 \(f(k)\) 的值 ,\(x_i-x_j\) 的逆元需要多次使用,可以像下面的类似做法预处理他们的差,当然kd和pkd要重新处理。不需要用map存,这样和直接求qpow一个鬼样。
开个invd数组,保存每两个项之间的差的逆元。
注意到prod里面有一项是 \(k-x_j\) 的积,这里也可以选择直接记录这个积。毕竟模运算比较费时,卡出来的常数可能差几倍
重心拉格朗日插值法
因为是要求逆元的,所以预处理 \(O(n^2logn)\) 。
然后每次求新k时,预处理重心拉格朗日插值法的g, \(O(nlogn)\) 。
多次求同一多项式的值快了不少的,比原版的拉格朗日插值。就是要 \(O(n^2)\) 的额外空间。
先卡掉重复求差的逆元的,已经可以通过了。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD=998244353;
int n,x[2005],y[2005];
int qpow(int x,int n){
int ans=1;
while(n){
if(n&1)
ans=((ll)ans*x)%MOD;
x=((ll)x*x)%MOD;
n>>=1;
}
return ans;
}
int invd[2005][2005];
int kd[2005];
void init_prod(int k){
for(int i=1;i<=n;i++){
kd[i]=(k-x[i]);
if(kd[i]<MOD)
kd[i]+=MOD;
for(int j=1;j<=n;j++){
if(i==j)
continue;
else{
int d=x[i]-x[j];
if(d<MOD)
d+=MOD;
invd[i][j]=qpow(d,MOD-2);;
}
}
}
}
int prod(int i,int k) {
ll ans=1;
for(int j=1;j<=n;j++){
if(j==i)
continue;
else{
ans=ans*(kd[j])%MOD*(invd[i][j])%MOD;
}
}
return ans;
}
int lagrange(int k) {
ll ans=0;
init_prod(k);
for(int i=1;i<=n;i++){
ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
}
return ans;
}
int main() {
#ifdef Yinku
freopen("Yinku.in","r",stdin);
#endif // Yinku
int k;
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++){
scanf("%d%d",&x[i],&y[i]);
}
printf("%d\n",lagrange(k));
}
再卡掉n倍logn,效果不是很明显,毕竟是n平方的算法:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD=998244353;
int n,x[2005],y[2005];
int qpow(int x,int n){
int ans=1;
while(n){
if(n&1)
ans=((ll)ans*x)%MOD;
x=((ll)x*x)%MOD;
n>>=1;
}
return ans;
}
int invd[2005][2005];
int kd[2005],pkd;
void init_prod(int k){
pkd=1;
for(int i=1;i<=n;i++){
kd[i]=(k-x[i]);
if(kd[i]<MOD)
kd[i]+=MOD;
pkd=(ll)pkd*kd[i]%MOD;
for(int j=1;j<=n;j++){
if(i==j)
continue;
else{
int d=x[i]-x[j];
if(d<MOD)
d+=MOD;
invd[i][j]=qpow(d,MOD-2);;
}
}
}
}
int prod(int i,int k) {
ll ans=(ll)pkd*qpow(kd[i],MOD-2)%MOD;
for(int j=1;j<=n;j++){
if(j==i)
continue;
else{
ans=ans*(invd[i][j])%MOD;
}
}
return ans;
}
int lagrange(int k) {
ll ans=0;
init_prod(k);
for(int i=1;i<=n;i++){
ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
}
return ans;
}
int main() {
#ifdef Yinku
freopen("Yinku.in","r",stdin);
#endif // Yinku
int k;
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++){
scanf("%d%d",&x[i],&y[i]);
}
printf("%d\n",lagrange(k));
}
一般我们代入的点都是连续自然数,当 \(x_i\) 的取值是连续自然数时,插值法因为分母变得容易处理而变成n复杂度。
以从0开始取举例。
假设当前n为4,从1开始取。其实分子无论怎么取值都是可以用前缀积和后缀积预处理出来的,然后两段乘在一起。关键是分母,当分母是连续自然数的时候不必两两枚举。
\(y_1\frac{1*(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}+y_2\frac{(x-1)*(x-3)(x-4)}{(2-1)(2-3)(2-4)}+y_3\frac{(x-1)(x-2)*(x-4)}{(3-1)(3-2)(3-4)}+y_4\frac{(x-1)(x-2)(x-3)*1}{(4-1)(4-2)(4-3)}\)
分母是:\(1*(-1)(-2)(-3),(1)*(-1)*(-2),(2)(1)*(-1),(3)(2)(1)\)
也就是前面阶乘,后面是负数的阶乘,负数的阶乘也很好想,比如这里n=4,那么第一项有3个负数,阶乘为负数,第二项为正数,第三项为负数……也就是n-i为奇数的项是负数。
可不断插点的重心拉格朗日插值法。
原式
$f(x)=\sum\limits_{i=1}^{n} y_i \prod\limits_{i\ne j}\frac{x-x_j}{x_i-x_j} $
把分子补上缺的项,提出来,记 \(g=\prod\limits_{i=1}^{n} x-x_i\)
$f(x)=g\sum\limits_{i=1}^{n} \frac{1}{x-x_i} \prod\limits_{i\ne j}\frac{y_i}{x_i-x_j} $
再记 \(h_i= \prod\limits_{i\ne j}\frac{y_i}{x_i-x_j}\)
$f(x)=g\sum\limits_{i=1}^{n} \frac{h_i}{x-x_i} $
那么每次多加一个点的时候,更新 \(g\) 以及每一个 \(h_i\) 就可以了。