拉格朗日插值法

中国剩余定理和Lagrange插值法

引言

很巧，高代课讲到了Lagrange插值法，凭着“上课听不明白只好下课自己搞”的意志，我学学证明顺带来写个板子。

前置知识：中国剩余定理（CRT）

中国剩余定理是用来求解形如
$$\{\begin{array}{lllll}x& \equiv a_1(mod{m}_1)x& \equiv a_2(mod{m}_2)& \dots \dots x& \equiv a_k(mod{m}_k)\end{array}$$
的线性同余方程组的解的。其中，$\mathrm{\forall }i\ne j(a_i,a_j)=1$，我们需要找出最小的非负整数解$x$。

解法

设$M=\prod _{i=0}^k{m}_i$,$M_i=\frac{M}{m_i}$,$M_i{T}_i\equiv 1(mod{m}_i)$，其中$1\le i,j\le n$。

故可以构造一个解$x=\sum _{i=1}^n{a}_i{M}_i{t}_i$

由此，可以求出任意解$$

**最小正整数解$$**即为$$

证明

显然，在第$$个同余方程里，$$

$$

所以$$，满足题意。

乘法逆元

若$$，则称$$是$$在$$意义下的逆元。

用扩展欧几里得$$来求解

过程如下：

$$
因为ax+by=(a, b)\
且b*x’+(a%b)y’=(b,a%b)\
a%b=a-\lfloor\frac{a}{b}\rfloorb
$$

$$
有ay’+b(x’-\lfloor\frac{a}{b}\rfloory’)=(a,b)\
故x=y’,y=(x’-\lfloor\frac{a}{b}\rfloory’)
$$

显然，$$，故令$$，则$$，故$$，可得$$为$$在$$下的逆元，即为题设所求。

代码

ll n, a[maxn], m[maxn], Mi[maxn], M = 1, X;

void exgcd(ll a, ll b, ll &x, ll &y) {
    if(b == 0) {x = 1; y = 0; return;}
    exgcd(b, a % b, x, y);
    ll z = x;
    x = y, y = z - y * (a / b);
}

void CRT() {
    n = read();
    rep(i, 1, n) {
        m[i] = read(), a[i] = read();
        M *= m[i];
    }
    rep(i, 1, n) {
        Mi[i] = M / m[i];
        ll x = 0, y = 0;
        exgcd(Mi[i], m[i], x, y);
        X += a[i] * Mi[i] * (x < 0 ? (x + m[i]):x);
    }
    cout << X % M << endl;
}
Copy

正片：Lagrange插值

多项式储备知识（多项式余数）

已知多项式
$$

$$

下证：

已知
$$

可知，当$$时，该子项$$恒成立
故余数必为$$的项之和，即
$$

模意义下的推导过程

由上述可以得到一个关于$$的同余方程组，即
$$
根据中国剩余定理
$$
得$$在$$意义下的逆元$$为
$$
故有
$$

又$$唯一，故
$$
即为Lagrange插值的表达式

通常意义下的推导过程

已知$$过$$个点$$，以及其在$$轴上的投影$$

考虑构造$$个函数$$使得其过点，
$$
\left{
$$
\right.
$$
则可知所需函数$$.

则可设$$，将$$代入可得$$，所以
$$
由此导出Lagrange插值在通常意义下的公式，如下
$$
与上述在模意义下的Lagrange插值公式相同

代码实现（对$$时，即需要求$$的值的代码）

const ll p = 998244353;
ll n, y[maxn], x[maxn], k;
ll ans, s1, s2;

ll fastmod(ll a, ll k) {
    ll base = 1;
    while(k) {
        if(k & 1) base = base * a % p;
        a = a * a % p;
        k >>= 1;
    }
    return base % p;
}

ll inv(ll x) {
    return fastmod(x, p - 2);
}

void Lagrange() {
    n = read(), k = read();
    rep(i, 1, n) {
        x[i] = read(), y[i] = read();
    }
    rep(i, 1, n) {
        s1 = y[i] % p;
        s2 = 1LL;
        rep(j, 1, n) {
            if(i != j) {
                s1 = s1 * (k - x[j]) % p;
                s2 = s2 * (x[i]-x[j]) % p;
            }
        }
        ans += s1 * inv(s2) % p;
    }
    printf("%lld\n", (ans % p + p) % p);
}
Copy