mayoko’s diary

プロコンとかいろいろ。

yukicoder No.42 貯金箱の溜息

解法

持っている硬貨をそれぞれ C[0] = 1, C[1] = 5, ..., C[5] = 500 とします。

普通の dp で解こうとすると,
dp[x][M] = (C[x] 円以下の硬貨で M 円を支払う場合の数)となりますが, この dp は

dp[x][M] = dp[x-1][M] + dp[x][M-C[x]]

とすることで解くことが出来ます。もちろんこれでは計算時間的にムリですが解法の肝になってるアイデアは(多分)上式から導かれます。

解法の肝は, dp[x][m*C[x]+q] (m*C[x]+q = M で q < C[x])が実は m の x 次式で表せる, ということです。この証明(あってるのか)は後にするとして, これがわかると, dp[x][q], dp[x][C[x]+q], dp[x][2*C[x]+q], ..., dp[x][5*C[x]+q] の値をはっきりさせることで, ラグランジュ補間することが出来, それで一般の M についても答えがわかります。
最初に dp を前計算して dp[x][y](y は 3000 程度までで良い) がわかるようにしておけば 各テストにかかる時間は O(K^2) (K は硬貨の数)なので, 十分間に合います。

後は dp[x][m*C[x]+q] が m の x 次式で表せることの証明ですね。q はパラメーターで, 各 q について上で言った x 次式は異なることに注意です。さっき書くの忘れました。
数学的帰納法で示します。

まず, dp[0][*] = 1 で, 1 は 0 次式なので, 確かに x 次式です。

次に, dp[x-1][*] が x-1 次式で書けるとします。すると, 上の漸化式より,
dp[x][m*C[x]+q] = dp[x][(m-1)*C[x]+q] + dp[x-1][m*C[x]+q]
dp[x][(m-1)*C[x]+q] = dp[x][(m-2)*C[x]+q] + dp[x-1][(m-1)*C[x]+q]
dp[x][(m-2)*C[x]+q] = dp[x][(m-3)*C[x]+q] + dp[x-1][(m-2)*C[x]+q]
...
dp[x][C[x]+q] = dp[x][q] + dp[x-1][C[x]+q]

となるので,
dp[x][m*C[x]+q] = dp[x][q] + dp[x-1][m*C[x]+q] + dp[x-1][(m-1)*C[x]+q] + ... + dp[x-1][C[x]+q] …(☆)

ここで, 数学的帰納法の仮定より, dp[x-1][a*C[x-1]+q] という形の式は a の x-1 次式になります。また, C[x-1] は問題で与えられている数値上 C[x] の約数になっているので, dp[x-1][a*C[x]+q] という式は, dp[x-1][a*C[x-1]+q] と同じく a の x-1 次式になるはずです。

ここで (☆) 式に戻ると, 右辺は 定数 + (m に関する x-1 次式 m 個の和) という形になっています。よって, 右辺は x 次式です。
以上より dp[x][m*C[x]+q] は x 次式になります。

得た知見

上でやってる解法とか証明とか, 「どうしてこんなこと思いつくねん」という感じなので, もうちょっと掘り下げてみます。

問題の本質は, 「二項間漸化式」です。一番最初に書いた式
dp[x][M] = dp[x-1][M] + dp[x][M-C[x]]
ですが, これは (dp[x][M] というのが x 次の多項式になることを知っていれば) よく受験に出る有名な形です。例えば,
 p_n = p_{n-1} + n^2+n+1
のような感じです。( p_n が dp[x][M] に,  p_{n-1} が dp[x][M-C[x]]に,  n^2+n+1 が dp[x-1][M] に対応してる感じです)このような式は, 階差数列なので当然一般項は n についての 3 次式になります。

ということで, dp[x][M] は x が上がるに連れてどんどん次数が上がるようになります。

このような観察が本質なのかなぁ〜と思います。

const ll MOD = 1e9+9;
const int MAXN = 3100;
int C[6] = {1, 5, 10, 50, 100, 500};
ll dp[MAXN];

ll pow_mod(ll x, ll p) {
    if (p == 0) return 1;
    if (p == 1) return x;
    if (p%2) return (pow_mod(x, p-1)*x) % MOD;
    ll tmp = pow_mod(x, p/2);
    return (tmp*tmp)%MOD;
}

// mod_inverse
ll mod_inverse(ll a, ll m) {
    return (pow_mod(a, m-2));
}

void solve(ll M) {
    if (M < MAXN) {
        cout << dp[M] << endl;
        return;
    }
    ll ans = 0;
    ll q = M%500;
    ll m = (M/500) % MOD;
    for (int i = 0; i < 6; i++) {
        ll tmp = 1;
        for (int j = 0; j < 6; j++) {
            if (i == j) continue;
            (tmp *= m-j) %= MOD;
            (tmp *= mod_inverse(i-j, MOD)) %= MOD;
        }
        (tmp *= dp[i*500+q]) %= MOD;
        (ans += tmp) %= MOD;
    }
    ans = (ans+MOD)%MOD;
    cout << ans << endl;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    int T;
    cin >> T;
    for (int i = 0; i < MAXN; i++) dp[i] = 1;
    for (int i = 1; i < 6; i++) {
        for (int j = C[i]; j < MAXN; j++) (dp[j] += dp[j-C[i]]) %= MOD;
    }
    while (T--) {
        ll M;
        cin >> M;
        solve(M);
    }
    return 0;
}