ラグランジュ補間 - ferinの競プロ帳

次の問題を解くアルゴリズム
$n$ 次以下の多項式 $f$ が存在する． $f(x_i) = y_i$ となる $(x_i, y_i)$ の組が $n+1$ 個存在するときに $f$ を復元しろ．

求める $f$ の形式で場合分けする．

$f(T)$ を一つ求める ( $T$ はgiven)

$f_i(x) = \prod_{j \neq i} (x-x_j) = \prod_{j} (x-x_j) / (x-x_i)$ とすると $f(x) = \sum_{i=0}^n c_i f_i(x)$ (ただし $c_i$ は未知の係数) と書ける．このように $f_i(x)$ を定めると $f_i(x_j) = 0\ \ (i \neq j)$ となるのが都合がいい． $f(x_i) = c_i f_i(x_i)$ となるため $c_i = f(x_i) / f_i(x_i)$ と係数 $c_i$ を求められる．
$f_i(x_i)$ は $O(n)$ で求められるので $O(n^2)$ で $c_i\ \ (0 \le i \le n)$ を求められる．

$f_i(T)$ は $O(n)$ で求められるので $f(T) = c_0 f_0(T) + c_1 f_1(T) + \ldots + c_n f_n(T)$ は $O(n^2)$ で求められる．

// x座標が相異なるn+1点(x_i,y_i)を通るn次以下の多項式f(T)の値を返す
// O(n^2)
mint lagrange_interpolation(vector<mint> x, vector<mint> y, mint T) {
    const ll n = x.size() - 1;
    mint ret = 0;
    REP(i, n+1) {
        mint t = 1;
        REP(j, n+1) {
            if(i == j) continue;
            t *= T-x[j];
            t /= x[i]-x[j];
        }
        ret += t * y[i];
    }
    return ret;
}

ARC033 D - 見たことのない多項式の部分点「提出」

$f(x)$ の係数を求める

$c_i$ を求めるところまでは同じようにやる．

$f(x) = \sum_{i=0}^n c_i f_i(x)$ を求めるために $f_i(x)$ を $O(n)$ で求める．これを戻すDPを使って行う． $(x-x_0)(x-x_1)\cdots(x-x_n)$ の $x^j$ の係数をDPを使って計算する． $dp(i,j) = (x-a_i$ まで見たときの $x^j$ の係数 $)$ とする．遷移は $dp(i,j) = dp(i-1,j) \times -x_i + dp(i-1,j-1)$ となる．このDPの計算は $O(n^2)$ でできる．

$dp(i,j)$ から $dp(i-1,j)$ を復元することを考える． $j$ が小さい方から順番に見ていく． $dp(i,0) = dp(i-1,0) \times -x_i$ であるから $dp(i-1,0) = -dp(i,0) / x_i$ と求められる． $dp(i,j) = dp(i-1,j)\times -x_i + dp(i-1,j-1)$ であるから $dp(i-1,j) = -(dp(i,j)-dp(i-1,j-1))/x_i$ と求められる．このようにdpを戻す操作は $(x-x_0)(x-x_1)\cdots(x-x_n)$ を $(x-x_i)$ で割る操作に対応する．

よって $f_i(x)$ を $O(n)$ で求めることができたので $f(x) = \sum_{i=0}^n c_i f_i(x)$ を $O(n^2)$ で求められた．

// x座標が相異なるn+1点(x_i,y_i)を通るn次以下の多項式f(x)を返す
// O(n^2)
vector<mint> lagrange_interpolation(vector<mint> x, vector<mint> y) {
    const ll n = x.size() - 1;
    REP(i, n+1) {
        mint t = 1;
        REP(j, n+1) if(i != j) t *= x[i]-x[j];
        y[i] /= t;
    }
    ll cur = 0, nxt = 1;
    vector<vector<mint>> dp(2, vector<mint>(n+2));
    dp[0][0] = -1 * x[0], dp[0][1] = 1;
    FOR(i, 1, n+1) {
        REP(j, n+2) {
            dp[nxt][j] = dp[cur][j] * -1 * x[i];
            if(j >= 1) dp[nxt][j] += dp[cur][j-1];
        }
        swap(nxt, cur);
    }
    vector<mint> xinv(n+1);
    REP(i, n+1) xinv[i] = x[i].inv();
    vector<mint> ret(n+1);  // f(x)
    REP(i, n+1) {
        if(y[i]==0) continue;
        // 0割り対策の場合分け
        if(x[i] == 0) {
            REP(j, n+1) ret[j] += dp[cur][j+1] * y[i];
        } else {
            ret[0] -= dp[cur][0] * xinv[i] * y[i];
            mint pre = -1 * dp[cur][0] * xinv[i];
            FOR(j, 1, n+1) {
                ret[j] -= (dp[cur][j] - pre) * xinv[i] * y[i];
                pre = -1 * (dp[cur][j] - pre) * xinv[i];
            }
        }
    }
    return ret;
}

ABC137 F - Polynomial Construction 「提出」

与えられる点の座標が等差数列

$x_i$ が等差数列になってるケースではもっと高速に解ける．
初項 $a$ で交差 $d$ の等差数列とする．

$f(T)$ を一つ求める (Tはgiven)

求めるのに $O(n^2)$ かけていた $f_i(x_i), f_i(T)$ を高速に求めることができればうれしい．

$f_i(x_i)$ と $f_{i+1}(x_{i+1})$ の形が似ていることを使う．
$f_i(x_i) = (x_i-x_0)\times(x_i-x_1)\times\cdots\times d\times -d \times\cdots\times(x_i-x_{n-1})\times(x_i-x_n)$
$f_{i+1}(x_{i+1}) = (x_i+d-x_0)\times(x_i-x_0)\times\cdots\times d\times -d \times\cdots\times(x_i-x_{n-2})\times(x_i-x_{n-1})$
比較すると $(x_i+d-x_0)$ が増えて $(x_i-x_n)$ が消えている．したがって $O(1)$ で $f_{i}(x_i)$ から $f_{i+1}(x_{i+1})$ を求められる．

$f_i(T) = (T-x_0)(T-x_1)\cdots(T-x_n)/(T-x_i)$ を計算するには，予め $\prod_i (T-x_i)$ を計算しておいて， $T-x_i$ で割ることで高速に計算できる．

modの割り算が必須なので計算量は全体で $O(n\log \text{mod})$ となる．

// x座標が相異なるn+1点(x_i,y_i)を通るn次以下の多項式f(T)の値を返す
// x_i = a + i*d 0<=i<=n (等差数列)
// O(nlog(MOD))
mint lagrange_interpolation_arithmetic(mint a, mint d, vector<mint> y, mint T) {
    const ll n = y.size() - 1;
    mint ret = 0, ft = 1;
    REP(i, n+1) ft *= T-(a+d*i);
    // f_0(x_0)
    mint f = 1;
    FOR(i, 1, n+1) f *= -1*i*d;
    ret += c[0] * ft / (T-a);
    // f_i(x_i) → f_{i+1}(x_{i+1})
    REP(i, n) {
        f *= d*(i+1) / (d*(i-n));
        ret += y[i+1] / f * ft / (T-a-d*(i+1));
    }
    return ret;
}

ARC033 D - 見たことのない多項式「提出」

例題

ネタバレ注意
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

codeforces Educational Codeforces Round 7 F. The Sum of the k-th Powers

問題ページ
$S(x) = \sum_{i=0}^x i^k \bmod 10^9+7$ とすると答えは $S(n)$ になる． $S(x)$ はk+1次式っぽいので $S(0), S(1), \ldots, S(k+1)$ のk+2点を使ってラグランジュ補間をすれば $S(n)$ が求まる．

答えが多項式の値なことはわかっているが具体的に多項式を構成するのがつらいときに使う

Submission #58554396 - Codeforces

yukicoder No.42 貯金箱の溜息

問題ページ
dp[i番目のコインまで][金額j] = (支払う方法の数) としたdpテーブルでf(m) = dp[5][500m+q] q<500 はmについて5次式になっている（らしい）．f(0)=dp[5][500×0+q], f(1)=dp[5][500×1+q], … , f(5)=dp[5][500×5+q] の6点を用いて多項式fを復元し，f(M/500)を求める．

答えが多項式の値なことはわかっているが具体的に多項式を構成するのがつらいときに使う
(これが多項式だって気づく難易度高すぎない？)

提出 #367653 No.42 貯金箱の溜息 - yukicoder

ライブラリ作ったはいいけど今後使うことはあるのだろうか

を一つ求める (はgiven)

の係数を求める