拉格朗日插值法

中国剩余定理和Lagrange插值法

引言

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

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

中国剩余定理是用来求解形如 \[ \begin{cases} x & \equiv a_1 (mod \quad m_1) \\ x & \equiv a_2 (mod \quad m_2) \\ &...... \\ x & \equiv a_k (mod \quad m_k) \end{cases} \] 的线性同余方程组的解的。其中，\(\forall i\neq j \quad (a_i, a_j)=1\)，我们需要找出最小的非负整数解\(x\)。

解法

设\(M=\prod_{i=0}^k m_i\),\(M_i=\frac{M}{m_i}\),\(M_iT_i \equiv 1(mod \ m_i)\)，其中\(1\leq i,j \leq n\)。

故可以构造一个解\(x=\sum_{i=1}^na_iM_it_i\)

由此，可以求出任意解\(x_0=x+k*M\)

最小正整数解\(X\)即为\(x_0\%M\)

证明

显然，在第\(i\)个同余方程里，\(a_jM_jt_j \equiv 0(mod \ m_i)\)

\(a_iM_it_i \equiv a_i(mod \ m_i)\)

所以\(\forall i,X \equiv a_i (mod \ m_i)\)，满足题意。

乘法逆元

若\(a*a^{-1} \equiv 1 (mod \ p)\)，则称\(a^{-1}\)是\(a\)在\(mod \ p\)意义下的逆元。

用扩展欧几里得\((exgcd)\)来求解

过程如下：

\[ 因为a*x+b*y=(a, b)\\ 且b*x'+(a\%b)*y'=(b,a\%b)\\ a\%b=a-\lfloor\frac{a}{b}\rfloor*b \]

\[ 有a*y'+b*(x'-\lfloor\frac{a}{b}\rfloor*y')=(a,b)\\ 故x=y',y=(x'-\lfloor\frac{a}{b}\rfloor*y') \]

显然，\((M_i, m_i)=1\)，故令\(a=M_i, b=m_i\)，则\(M_ix+m_iy=1\)，故\(M_ix=1-m_iy\)，可得\(x\)为\(M_i\)在\(mod\ m_i\)下的逆元，即为题设所求。

代码

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;
}

正片：Lagrange插值

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

已知多项式 \[ f(x)\in P[x],f(x)=\sum_{i=0}^{n}a_ix^i \]

\[ 故f(x)\%(x-\alpha)=\sum_{i=0}^{n}a_i\alpha^i \]

下证：

已知 \[ f(x)=\sum_{i=0}^{n}a_ix^i \\ f(x)=\sum_{i=0}^{n}a_ix^i=\sum_{i=0}^{n}(x-\alpha+\alpha)^i=\sum_{i=0}^{n}a_i\sum_{j=0}^{i}{i \choose j}(x-\alpha)^j\alpha^{i-j} \]

可知，当\(j\neq0\)时，该子项\(p\%(x-\alpha)=0\)恒成立 故余数必为\(j=0\)的项之和，即 \[ \begin{aligned} r(x)=\sum_{i=0}^{n}a_i\alpha^i \end{aligned} \]

模意义下的推导过程

由上述可以得到一个关于\(f(x)\)的同余方程组，即 \[ \begin{cases} f(x) & \equiv y_1 (mod \ (x-x_1)) \\ f(x) & \equiv y_2 (mod \ (x-x_2)) \\ &...... \\ f(x) & \equiv y_n (mod \ (x-x_n)) \end{cases} \] 根据中国剩余定理 \[ \begin{aligned}M=\prod_{i=1}^{n}(x-x_i),M_i=\frac{M}{x-x_i}=\prod_{j\neq i}(x-x_j)\end{aligned} \] 得\(M_i\)在\(mod(x-x_i)\)意义下的逆元\(T_i\)为 \[ \begin{aligned}M_i^{-1}=T_i=\prod_{j\neq i}\frac{1}{x_i-x_j}\end{aligned} \] 故有 \[ \begin{aligned} f(x) \equiv \sum_{i=1}^{n}y_iM_iT_i\equiv \sum_{i=1}^{n}y_i\prod_{j\neq i}M_iT_i\equiv \sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{x-x_i}{x_i-x_j} \quad mod(M) \end{aligned} \]

又\(f(x)\)唯一，故 \[ \begin{aligned} f(x) = \sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{x-x_i}{x_i-x_j} \end{aligned} \] 即为Lagrange插值的表达式

通常意义下的推导过程

已知\(f(x)\)过\(n\)个点\(P_1(x_1,y_1),P_2(x_2,y_2),...,P_n(x_n,y_n)\)，以及其在\(x\)轴上的投影\(P_i^{'}(x_i,0)\)

考虑构造\(n\)个函数\(F_i(x)\)使得其过点， \[ \left\{ \begin{aligned} & P_j^{'} \quad j\neq i \\ & P_i \quad j=i \end{aligned} \right. \] 则可知所需函数\(f(x)=\sum_{i=1}^{n}F_i(x)\).

则可设\(F_i(x)=a\prod_{j\neq i}(x_i-x_j)\)，将\(P_i\)代入可得\(a=\frac{y_i}{\prod_{j\neq i}(x_i-x_j)}\)，所以 \[ \begin{aligned}F_i(x)=y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}\end{aligned} \] 由此导出Lagrange插值在通常意义下的公式，如下 \[ \begin{aligned}f(x)=\sum_{i=1}^{n}F_i(x)=\sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}\end{aligned} \] 与上述在模意义下的Lagrange插值公式相同

代码实现（对\(x=k\)时，即需要求\(f(k)\)的值的代码）

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);
}