SRM 527 div1 hard:P8XCoinChange

問題

解法

ここに書く解法は, topcoder の解説を参考にして書いてます。
http://apps.topcoder.com/wiki/display/tc/SRM+527
こっちのほうが証明とかいろいろ丁寧です。英語面倒ですがこっち読むほうがオススメです。

まず, 問題を言い換えてみます。与えられた問題は以下の問題と同じです。

「あるウサギが 0 からスタートして coins_sum にちょうど着くようにジャンプしながら進もうとしている。ただし, ジャンプして進
むことの出来る距離は values に書かれるいずれかの距離のみで, また, ジャンプして進む距離は非増加でなければならない(例えば, ジャンプして進む距離が 2 -> 2 -> 4 -> 1 とかなるのは (2 進んだ後 4 進んでいるので) NG)。ウサギが coins_sum にピッタリたどり着くようなジャンプの列は何通りあるか。」

ジャンプして進む距離が非増加である, というのが大事です。

で, この問題に一見関係ないようですが, 以下のようなアルゴリズムを考えます。
「ある地点にいる時, coins_sum を超えない最大の values の要素を選び, 進む。」

例えば values = {1, 3, 6} で, coins_sum = 11 の時, 6, 9, 10, 11 と進んでいきます。このようにして進んだ時に止まった座標(さっきの例では 6, 9, 10, 11) は, 実は上記の問題では必ず通る地点になります。ちゃんとした証明が topcoder の解説には書いてありますが, 直感的に明らかでしょう。ジャンプして進む距離が非増加なので, 必ず踏むことになりそうです。これがキーになります。

values が要素を n 個持っているとします。
solve(length) という関数を考えます。これは n*n の行列を返し, その行列の (i, j) 成分は, 「values[i] 以下のジャンプしかしないで, 最後のジャンプが values[j] であり, 合計で length だけ進んでいるような場合の数」を返します。
こんな関数作れるなら solve(coins_sum) を求めればええやんけｗと思いますが, この length には条件があって, length として取れるものは, values の要素だけです。

solve 関数のアルゴリズムはひとまず置いておいて, これがわかるとどうやって答えがわかるのかを考えます。
先ほど示したとおり, coins_sum までの道のりでは, 必ず踏む地点があります。しかも, その必ず踏む地点と次の地点の間は, values の要素で表すことが出来ます。
ということで, 例えば coins_sum までに上のアルゴリズムで p 回 values[i] というジャンプを使うのであれば, solve(values[i])^p を計算すると, values[i] のジャンプを使う間のジャンプの仕方の場合の数がわかります。これはよくある行列累乗テクニックなので OK でしょう(solve(A+B) = solve(A)*solve(B) という関係から, 累乗にしても同じになる)。

よって, solve 関数ができれば, 答えが得られることがわかりました。あとは solve 関数のアルゴリズムを考えましょう。これはそれほど難しくありません。
まず solve(values[0]) は, 0 列目がすべて 1 であるような行列になります。定義から明らかです。
また, k >= 1 を満たす k について, solve(values[k]) は, values[k] のジャンプを使って一気に values[k] まで行くか, values[k-1] 以下のジャンプを最初に使って, 小刻みに values[k] までたどり着くかのいずれかです。
前者は, solve(values[k]) の行列を m として, m[j][k] = 1 (j >= k) となります。
後者の場合は, values[k-1] 以下のジャンプのみを用いて values[k-1] 進むのを values[k]/values[k-1] 回やるのと同じになるので, solve(values[k-1])^(values[k]/values[k-1]) を求めれば OK です。

typedef long long number;
typedef vector<number> vec;
typedef vector<vec> matrix;

const ll MOD = 1e6+3;

matrix Zero(int n) {
    matrix A(n, vec(n));
    return A;
}

// O( n )
matrix identity(int n) {
    matrix A(n, vec(n));
    for (int i = 0; i < n; ++i) A[i][i] = 1;
    return A;
}
// O( n^3 )
matrix mul(const matrix A, const matrix B) {
    matrix C(A.size(), vec(B[0].size()));
    for (int i = 0; i < (int)C.size(); ++i)
        for (int j = 0; j < (int)C[i].size(); ++j)
            for (int k = 0; k < (int)A[i].size(); ++k) {
                C[i][j] += A[i][k] * B[k][j];
                C[i][j] %= MOD;
            }
    return C;
}
// O( n^3 log e )
matrix pow(const matrix &A, ll e) {
    if (e == 0) return identity(A.size());
    if (e == 1) return A;
    if (e % 2 == 0) {
        matrix tmp = pow(A, e/2);
        return mul(tmp, tmp);
    } else {
        matrix tmp = pow(A, e-1);
        return mul(A, tmp);
    }
}

class P8XCoinChange {
public:
    int solve(long long coins_sum, vector<long long> values) {
        int n = values.size();
        vector<matrix> M(n);
        M[0] = Zero(n);
        for (int i = 0; i < n; i++) {
            M[0][i][0] = 1;
        }
        for (int i = 1; i < n; i++) {
            ll e = values[i]/values[i-1];
            M[i] = pow(M[i-1], e);
            for (int j = i; j < n; j++) M[i][j][i]++;
        }
        matrix total = identity(n);
        for (int i = n-1; i >= 0; i--) {
            ll p = coins_sum/values[i];
            coins_sum %= values[i];
            total = mul(total, pow(M[i], p));
        }
        ll ans = 0;
        for (int i = 0; i < n; i++) (ans += total[n-1][i]) %= MOD;
        return ans;
    }
};

mayoko’s diary

プロコンとかいろいろ。

SRM 527 div1 hard:P8XCoinChange

問題

解法