2020-11-03

AOJ 2638 Hyperrectangle (JAG 夏合宿 2014 day2-J) (550 点)

りんごさんセット！！！

問題へのリンク

問題概要

$D$ 次元空間において

$0 \le x_{i} \le L_{i}$ 　( $i = 1, 2, \dots, D$ )
$x_{1} + \dots + x_{D} \le S$

を満たす部分の体積を $V$ とすると $VD!$ は整数値となる。この整数値を 1000000007 で割ったあまりを求めよ。

制約

$2 \le D \le 300$
$1 \le L_{i} \le 300$
$0 \le S \le \sum_{i=1}^{N}L_{i}$

考えたこと

最初は色々迷走していた。三次元の場合を考えていて思っていたのは

$S$ と $l_{1}$ の大小関係
$S$ と $l_{2}$ の大小関係
$S$ と $l_{3}$ の大小関係
$S$ と $l_{1} + l_{2}$ の大小関係
$S$ と $l_{2} + l_{3}$ の大小関係
$S$ と $l_{1} + l_{3}$ の大小関係

がそれぞれ異なるパターンが考えられるということ。これだけでも $2^{6} = 64$ 通りのパターンがある (対称性やあり得ない場合の除去で多少減らせるが...)。一般の場合を考えるとなると大変そう。。。

(以後、ぼやき)
最初に考えていたのは、1 × 1 × 1 の立方体に帰着させる方法だった。1 × 1 × 1 の立方体の切り取られ方は、 $D$ 次元なら $D-1$ 通りある。それぞれについて体積を求めておいて、各切り取られ方をする立方体が何個あるかを数えるという方針だ。でもそれをやろうとしたときに包除原理が必要になって、それで直接一般の場合も計算できるやんとなった。

包除原理へ

$x_{i} \le L$ という条件がイヤで、これがなければ簡単に体積を求められる。たとえば

$x_{1}, \dots, x_{D} \ge 0$
$x_{1} + \dots + x_{D} \le S$

を満たす部分の体積は $\frac{S^{D}}{D!}$ で求められる。というわけで包除原理が有効だと気づく (図形的には三次元でいえば、三角錐の切り貼りで表現できることに対応する)。

たとえば $x_{1}$ のみ条件を違反している部分の体積は次のようにして求められる。

$x_{1} - L_{i} \ge 0$
$(x_{1} - L_{i}) + x_{2} + \dots + x_{D} \le S - L_{i}$

を満たす部分の体積を求める問題に帰着されるので、 $\frac{S(S-L_{i})^{D}}{D!}$ となる。

同じようなノリで、 $(1, 2, \dots, D)$ のすべての部分集合 $A$ について、 $S$ に対応する部分の制約を違反している体積を $V(A)$ を求めて、

$\sum_{A} (-1)^{|A|}V(A)$

を計算すれば OK。ここで次の値を DP によって求める。

dp[ 0 or 1 ][ val ] := $L_{1}, \dots, L_{D}$ のうち (0: 偶数個, 1: 奇数個) を使って値 val にするような場合の数

これが求められれば包除原理で OK。計算量は、 $O(DS)$ 。

コード

#include <bits/stdc++.h>
using namespace std;

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { 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 istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return modpow(modinv(r), -n);
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& 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);
        }
        return Fp<MOD>(u);
    }
};

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

int main() {
    int D, S;
    cin >> D;
    vector<int> L(D);
    for (int i = 0; i < D; ++i) cin >> L[i];
    cin >> S;

    vector<vector<mint>> dp(2, vector<mint>(S+1, 0)), ndp = dp;
    dp[0][0] = 1;
    for (int i = 0; i < D; ++i) {
        ndp.assign(2, vector<mint>(S+1, 0));
        for (int p = 0; p < 2; ++p) {
            for (int s = 0; s <= S; ++s) {
                ndp[p][s] += dp[p][s];
                if (s + L[i] <= S) ndp[1-p][s+L[i]] += dp[p][s];
            }
        }
        swap(dp, ndp);
    }
    mint res = 0;
    for (int p = 0; p < 2; ++p) {
        for (int s = 0; s <= S; ++s) {
            mint tmp = modpow(mint(S-s), D) * dp[p][s];
            if (p % 2 == 0) res += tmp;
            else res -= tmp;
        }
    }
    cout << res << endl;
}