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

目录

目录地址

上一篇

下一篇


拉格朗日插值

我们现在给定 \((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)\)

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

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

posted @ 2020-03-05 10:18  JustinRochester  阅读(663)  评论(0编辑  收藏  举报