拉格朗日插值到中国剩余定理

目录

目录地址

上一篇

下一篇

---

拉格朗日插值

我们现在给定 ((n+1)) 个点 ((x_i,y_i)) ，请求解出一个不超过 (n) 次的多项式函数 (P(x))

使得 (forall iin[1,n+1]igcap 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 eq i,P(x_j)=0)

假设 (P(x_i)=b_i) 则记 (a_i={y_iover b_i}) 即 (a_i={y_iover 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 eq 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_iover 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_icdot 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)+Ccdot p(x)) 也符合上述式子，当然，它的次数只有当 (C=0) 时才能保证不大于 (n) 了

---

中国剩余定理(CRT)

又称孙子定理，由《孙子算经》中首次提到

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

其实质为解一个同余式： (egin{cases} xequiv a_1(mod m_1) \ \ xequiv a_2(mod m_2) \ \ xequiv a_3(mod m_3) \ \ ...... \ \ xequiv a_n(mod m_n) end{cases})

其中，保证 (forall i,jin[1,n]igcap Z,i eq jLeftrightarrow gcd(m_i,m_j)=1) 且对每个同余式都有解

当然，最后的解一定也为同余式 (xequiv a(mod m)) 的形式

一般会要求输出同余式、最小正数解或区间内解、区间内解的个数之类的

我们先注意到性质：任意两个模数都互质，那么我们直接考虑：

设 (displaystyle M_i=(prod_{j=1}^{i-1}m_j)cdot (prod_{j=i+1}^n m_j))

那么也很显然， (forall j eq i,M_imod m_j=0)

而同时，我们假定 (M_imod m_i=b_i) ，则令 (c_i=a_icdot (b_i)^{-1})

其中 ((b_i)^{-1}) 为 (b_i) 在模 (m_i) 意义下的逆元

因此 (c_iM_iequiv a_icdot (b_i)^{-1}cdot b_iequiv a_i(mod m_i))

那么，有没有可能不存在这个 (c_i) 呢？

显然是不可能的，因为 (M_imod m_i eq 0) 且保证了一定存在 (x) 使得 (xequiv a_i(mod m_i))

因此，一定存在 (c_i) 使得 (c_iM_iequiv 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 kin Z,a+km) 也可使上述方程组成立

因此最后的结果为 (xequiv a(mod m))

求出这个后，接下来的思路就简单了

代码类似于拉格朗日插值，就不写了