けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

diverta 2019_2 E - Balanced Piles (橙色, 800 点)

最初は絶望感がヤバイタイプの問題。とても解ける問題には見えない。でも落ち着くと解ける。
こういう問題を作れるのはすごい。

問題へのリンク

問題概要

 N 個のマスに積み木を重ねていく。最初はどのマスも 0 段である。以下の操作を繰り返して、最終的にすべてのマスが  H 段となるようにしたい。

  • 全マスのうち、積み木の段数が最も小さいマスを選ぶ。それが複数あるときは 1 つ選ぶ
  • 全マスの積み木の段数の最大値を  M としたとき、選んだマスの積み木の段数を、 M 以上  M + D にする

操作列として考えられるものが何通りあるか、1000000007 で割ったあまりを求めよ。

制約

  •  1 \le N, H, D \le 10^{6}

考えたこと

一目見たときの絶望感がやばい。DP を立てたいと思っても、一見ありうる状態が  (H+1)^{N} 通りある。対称性とかでちょっと削ってもそう簡単には減らない。

なので考えるべきことは

  • 実は考えるべき状態空間はもっとずっと小さいのではないか
  • そもそも根本的にもっとわかりやすい特徴づけを与えて、それを数えればよいのではないか

といったことだと思う。そんなときに着目したくなるのが

  • 積み木の段数の最大値とか最小値とか
  • それらが何個ずつあるかとか

そういった値だけ持っておけば状態をまとめられるのではないかということ。で、頑張ればそれが上手く行って解けるという筋書きに。

解法 1: 僕の解法 (わかりやすいものに特徴付ける)

僕のやり方は結局出来上がる DP は editorial と同じなのだが、ちょっと違う考察を経た気がするのでメモ。

ほとんどの人は「各ステップの積み木の最大段数」に着目していたが、僕は「最小段数」に着目した。そして最小段数のなす数列が

  • 最初に  0 N
  • それ以降は、「隣り合う2要素の差が  D 以下」「広義単調増加」「同じ値は  N 個以下」という条件を満たす数列
  • 最終項は  H で、 H は 1 個だけ

という風になっていることはすぐにわかる。例えば  N = 3, H = 2, D = 1 のとき、

(0, 0, 0) -> (0, 0, 1) -> (1, 0, 1) -> (1, 2, 1) -> (1, 2, 2) -> (2, 2, 2)

に対応する数列は

0, 0, 0, 1, 1, 2

となる。ここでこの数列に対応する操作が何通りあるかというと実は

  • 数列中に登場する各整数値について、その個数の積をとったもの

になることが言える。まず数列の  i 項目に対応する操作では、最下段を一つ選んで  i+N 項目の値にする必要があって、そうすれば対応する操作が作れることが言える。例えば上の数列の 3 番目の 0 では、6 番目が 2 であることから、0 を 2 に引き上げる必要があって、実際上の例でも (1, 0, 1) -> (1, 2, 1) としている。

なぜなら、その時点で最下段だったやつというのは  N ターン後に再び最下段になるからだ。厳密にはタイがあるのでややこしいが、うまいこと順序を決めたりしてあげるとちゃんと示せる。

また、数列中に値 v が  m 個あるとき「初めて v が最下位になったときには  m 個の v がある」ということを意味していて、それらをどの順序で引き上げていくかには任意性があるので  m! 通りの係数がかかることになる。よって、数列に対応する操作の個数は「数列中に登場する各整数値について、その個数の積をとったもの」になることが示された。

DP へ

ここまでわかれば DP するのみ。

  • dp[ v ] := 数列で値 v までを完遂した段階での対応する操作の個数

とすると、v の手前が u = v-D, v-D+1, ..., v-1 の D 通りを場合分けして

  • dp[ v ] += dp[ u ] × (1! + 2! + ... + N!)

となることがわかる。ここで

  • 1! + 2! + ... + N! の値は予め計算しておける
  • この漸化式は累積和で高速化できる

ということから全体で  O(H) で解くことができる。

#include <iostream>
#include <vector>
#include <cstring>
using namespace std;

template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) v += MOD;
    }
    constexpr int getmod() { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b; swap(a, b);
            u -= t * v; swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr istream& operator >> (istream &is, Fp<MOD>& x) noexcept {
        return is >> x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept {
        if (n == 0) return 1;
        auto t = modpow(a, n / 2);
        t = t * t;
        if (n & 1) t = t * a;
        return t;
    }
};

// 二項係数ライブラリ
template<class T> struct BiCoef {
    vector<T> fact_, inv_, finv_;
    constexpr BiCoef() {}
    constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
        init(n);
    }
    constexpr void init(int n) noexcept {
        fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
        int MOD = fact_[0].getmod();
        for(int i = 2; i < n; i++){
            fact_[i] = fact_[i-1] * i;
            inv_[i] = -inv_[MOD%i] * (MOD/i);
            finv_[i] = finv_[i-1] * inv_[i];
        }
    }
    constexpr T com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fact_[n] * finv_[k] * finv_[n-k];
    }
    constexpr T fact(int n) const noexcept {
        if (n < 0) return 0;
        return fact_[n];
    }
    constexpr T inv(int n) const noexcept {
        if (n < 0) return 0;
        return inv_[n];
    }
    constexpr T finv(int n) const noexcept {
        if (n < 0) return 0;
        return finv_[n];
    }
};

const int MOD = 1000000007;
using mint = Fp<MOD>;
BiCoef<mint> bc;

const int MAX = 2100000;

int N, H, D;
mint dp[MAX], sdp[MAX];

mint solve() {
    bc.init(MAX);
    mint fac = 0;
    for (int i = 1; i <= N; ++i) fac += bc.fact(i);
    memset(dp, 0, sizeof(dp));
    memset(sdp, 0, sizeof(sdp));
    dp[0] = bc.fact(N);
    sdp[1] = bc.fact(N);
    for (int v = 1; v <= H; ++v) {
        int left = max(0, v - D);
        int right = v;
        dp[v] = (sdp[right] - sdp[left]) * fac;
        sdp[v+1] = sdp[v] + dp[v];
    }

    mint res = 0;
    for (int v = max(0, H - D); v <= H-1; ++v) {
        res += dp[v];
    }
    return res;
}

int main() {
    cin >> N >> H >> D;
    cout << solve() << endl;
}

解法 2: DP の状態量を削る考え方

最大要素に着目して解いてみる。まず操作を以下の読み替える。積み木は常にソートされた状態にしておく。そして

  • 積み木の最も左にあるやつの高さが最大になるようにして、後ろの方に配置する
    • 最大のやつと同点になるなら、最大のやつの個数が  a 個あるとしたら、 a + 1 通りの挿入箇所がある
    • 最大のやつより大きくなるなら、最後尾に配置する

という操作と思うことにする。このとき初期状態の並べ方が  N! 通りあることに予め注意しておく。ルールとして重要なのは「最も左にあるやつの高さを最大になるようにする」という部分。これを決めておくことで考えやすくなる。

以上の考察をもとに

  • dp[ v ][ a ] := 最大高さが v で v が a 個ある状態にするまでの操作の個数

とすると

  • dp[ v ][ 1 ] += dp[ v - k ][ b ] (k = 1, 2, ..., D, b = 1, 2, ..., N)
  • dp[ v ][ k ] += dp[ v ][ k - 1 ] × k

という風になることがわかる。このままだと  O(NH) だが式をよく見ると、

  • dp[ v ][ 1 ], dp[ v ][ 2 ], ... の比率は、v の値によらず一定っぽい

ということがわかるので、これをまとめあげることができる。最終的に

  • dp[ v ] := dp[ v ][ 1 ]

とすると、

  • dp[ v ] += (1! + 2! + ... + N!) × dp[ v - k ] (k = 1, 2, ..., D)

という風になる。あとは解法 1 と一緒。