拉格朗日插值到中国剩余定理
目录
拉格朗日插值
我们现在给定 \((n+1)\) 个点 \((x_i,y_i)\) ,请求解出一个不超过 \(n\) 次的多项式函数 \(P(x)\)
使得 \(\forall i\in[1,n+1]\bigcap Z,P(x_i)=y_i\)
对于这种问题,暴力设出 \(P(x)\) 的 \((n+1)\) 个系数,然后代入这些点,联立方程求解
这些需要用到克拉默法则或者高斯消元法求解,复杂度挺高的,为 \(O(n^3)\),我们思考有没有更优的方法:
假设 \(\displaystyle P_i(x)=(\prod_{j=1}^{i-1}(x-x_j)\ )\cdot(\prod_{j=i+1}^{n+1}(x-x_j)\ )\)
那么显然: \(\forall j\neq i,P(x_j)=0\)
假设 \(P(x_i)=b_i\) 则记 \(a_i={y_i\over b_i}\) 即 \(a_i={y_i\over P(x_i)}\)
因此 \(a_iP(x_i)=y_i\)
这样有什么好处呢?
很显然,所求一定可以写为 \(\displaystyle P(x)=\sum_{i=1}^{n+1}a_iP(x_i)\)
例如我们代入 \(x_j\)
则 \(\displaystyle P(x_j)=\sum_{i=1}^{n+1}a_iP(x_j)=\sum_{i\neq j}a_iP(x_j)+a_jP(x_j)=0+y_j=y_j\)
虽然看起来并没有更优,但如果我们所求的不是 \(P(x)\) 本身,而是 \(P(x)\) 在 \(m\) 个点 \(x=k\) 处的取值时,速度会快上不少:
使用原方法,克拉默法则或是高斯消元法都是 \(O(n^3)\) 的,加上每一次秦九韶算法 \(O(n)\) ,总复杂度 \(O(n^3+nm)\)
而如果使用高斯消元法:
\(\displaystyle P(k)=\sum_{i=1}^{n+1}{y_i\over P_i(x_i)}P_i(k)\)
如果 \(\exist i=k\) 则 \(P(k)=y_i\) ,我们接下来讨论其余情况:
先花费 \(O(n^2)\) 的时间处理出所有的 \(P_i(x_i)\)
接下来,对于每个 \(k\) :
花费 \(O(n)\) 的时间预处理 \(\displaystyle L_i(k)=\prod_{i=1}^{i-1}(k-x_i),R_i(x)=\prod_{i+1}^n(k-x_i)\)
则,每次所需要求的 \(P_i(k)=L_i(k)\cdot R_i(k)\)
最后, \(P(k)=\displaystyle \sum_{i=1}^{n+1}{y_i\cdot L_i(k)\cdot R_i(k)\over P_i(x_i)}\) ,时间 \(O(n)\)
考虑需要代入的值的个数为 \(m\) 个,则总时间复杂度 \(O(n^2+nm)\)
不论 \(m\) 的取值为多少, \(O(n^2+nm)\) 一定比 \(O(n^3+nm)\) 更优
当然,当 \(m>n^2\) 时,两个复杂度都视为 \(O(nm)\) 了,就要考虑常数问题了
高斯消元法的实现
double px[MAXN],l[MAXN],r[MAXN],x[MAXN],y[MAXN],k[MAXN],ans[MAXN];
...
for(int i=1;i<=n+1;i++){
px[i]=1;
for(int j=1;j<=n+1;j++){
if(i!=j){
px[i]*=x[i]-x[j];
}
}
}
for(int j=1;j<=m;j++){
bool ext=0;
for(int i=1;i<=n+1;i++){
if(x[i]==k[j]){
ans[j]=y[i];
ext=1;
break;
}
}
if(ext) continue;
l[1]=r[n+1]=1;
ans[j]=0;
for(int i=2;i<=n+1;i++){
l[i]=l[i-1]*(k-x[i-1]);
}
for(int i=n;i>=1;i--){
r[i]=r[i+1]*(k-x[i+1])
}
for(int i=1;i<=n+1;i++){
ans[j]+=y[i]*r[i]*l[i]/x[i];
}
}
当然,如果保证输入是整数,且最后的答案对某数取模,就在前面求 \(P_i(x_i)\) 的预处理时,进行求逆元的操作,复杂度是 \(O(n(n+\log n)\ )=O(n^2)\) ,不变的
我们考虑一下,真的拉格朗日插值求出来的这个多项式 \(P(x)\) 是符合所有 \(P(x_i)=y_i\) 的唯一多项式吗?
显然不是,它只是次数不大于 \(n\) 的唯一多项式
我们考虑 \(\displaystyle p(x)=\prod_{i=1}^{n+1}(x-x_i)\) 以及任意常数 \(C\)
显然, \(P(x)+C\cdot p(x)\) 也符合上述式子,当然,它的次数只有当 \(C=0\) 时才能保证不大于 \(n\) 了
中国剩余定理(CRT)
又称孙子定理,由《孙子算经》中首次提到
有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
其实质为解一个同余式: \(\begin{cases} x\equiv a_1(\mod m_1) \\\ \\ x\equiv a_2(\mod m_2) \\\ \\ x\equiv a_3(\mod m_3) \\\ \\ ...... \\\ \\ x\equiv a_n(\mod m_n) \end{cases}\)
其中,保证 \(\forall i,j\in[1,n]\bigcap Z,i\neq j\Leftrightarrow gcd(m_i,m_j)=1\) 且对每个同余式都有解
当然,最后的解一定也为同余式 \(x\equiv a(\mod m)\) 的形式
一般会要求输出同余式、最小正数解或区间内解、区间内解的个数之类的
我们先注意到性质:任意两个模数都互质,那么我们直接考虑:
设 \(\displaystyle M_i=(\prod_{j=1}^{i-1}m_j)\cdot (\prod_{j=i+1}^n m_j)\)
那么也很显然, \(\forall j\neq i,M_i\mod m_j=0\)
而同时,我们假定 \(M_i\mod m_i=b_i\) ,则令 \(c_i=a_i\cdot (b_i)^{-1}\)
其中 \((b_i)^{-1}\) 为 \(b_i\) 在模 \(m_i\) 意义下的逆元
因此 \(c_iM_i\equiv a_i\cdot (b_i)^{-1}\cdot b_i\equiv a_i(\mod m_i)\)
那么,有没有可能不存在这个 \(c_i\) 呢?
显然是不可能的,因为 \(M_i\mod m_i\neq 0\) 且保证了一定存在 \(x\) 使得 \(x\equiv a_i(\mod m_i)\)
因此,一定存在 \(c_i\) 使得 \(c_iM_i\equiv a_i(\mod m_i)\) 的
同样的,我们参考拉格朗日插值,可以想到:令 \(\displaystyle a=\sum_{i=1}^n a_iM_i\)
即可保证当 \(x=a\) 时上述同余方程组成立
但是,我们知道,其不止一种解,若我们令 \(\displaystyle m=\prod_{i=1}^n m_i\)
则对于 \(\forall k\in Z,a+km\) 也可使上述方程组成立
因此最后的结果为 \(x\equiv a(\mod m)\)
求出这个后,接下来的思路就简单了
代码类似于拉格朗日插值,就不写了