読者です 読者をやめる 読者になる 読者になる

mayoko’s diary

プロコンとかいろいろ。

ラグランジュ補間について

なんとなくメモ。

ラグランジュ補間のアイデア

 x に関する N 次多項式  P(x) を,  P(x_0) = p_0, P(x_1) = p_1, ..., P(x_N) = p_N というヒントから表したいです。この時,  Q_i(x) = \frac{(x-x_0)(x-x_1)...(x-x_N)}{(x-x_i)} とすれば,  P(x) = c_0Q_0(x) + c_1Q_1(x) + ... + c_NQ_N(x) と表せることに注目します。


このように書くとハッピーなのは,  x = x_i とすると,  i \neq j を満たす  j について,  Q_j(x_i) = 0 が成り立つことです。このことから,  c_i = \frac{P(x_i)}{Q_i(x_i)} が得られます。

実装

問題から  x_i, p_i がわかっています。多項式の係数  c_i を求めるフェイズと実際に値を突っ込むフェイズに分けるとわかりやすい気がします。ちなみに yukicoder 42 を見なおしているときに書きたくなった記事なので yukicoder 42 っぽい感じで書いてます。

for (int i = 0; i < 6; i++) {
    // p を求める
    ll p = dp[6][q+i*500];
    // q を求める
    ll q = 1;
    for (int j = 0; j < 6; j++) {
        if (i==j) continue;
        q *= i-j;
    }
    q = (q+MOD)%MOD;
    // c = p/q
    c[i] = p*mod_inverse(q, MOD);
    c[i] %= MOD;
}
  • 実際に値を突っ込むフェイズ
for (int i = 0; i < 6; i++) {
    // Q_i(m) の計算
    ll tmp = 1;
    for (int j = 0; j < 6; j++) {
        if (i==j) continue;
        (tmp *= m-j) %= MOD;
    }
    // ret += c_i * Q_i(m)
    ret += c[i]*tmp % MOD;
}