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

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

AtCoder ABC 154 F - Many Many Paths (青色, 600 点)

すっごく色んな方法がありそう!!!

問題へのリンク

問題概要

正の整数  x1, y1, x2, y2 が与えられる。 x1 \le x \le x2,  y1 \le y \le y2 を満たすすべての整数  (x, y) についての

  •  C(x + y, y)

の総和を 1000000007 で割ったあまりを求めよ。

制約

  •  1 \le x1 \le x2 \le 10^{6}
  •  1 \le y1 \le y2 \le 10^{6}

考えたこと

とりあえず二項係数の計算自体は、適切に前処理をしておけば、 O(1) でできるようになる。詳しくはここに。

drken1215.hatenablog.com

さて、今回はまともにすべての二項係数を合計しようと思うと、最悪  10^{12} 個もの計算が必要になって間に合わない。よってこれらの二項係数に対して「いくつかグループに分けてグループごとに合計する」という方法を考えていくことになる。

y を固定する

つまり、 a = x1 + y,  b = x2 + y として、

  •  C(a, y) + C(a+1, y) + \dots + C(b, y)

の値を求める問題ということになる。そしてこの手の問題の定石として、

  1.  S(p) = C(1, y) + C(2, y) + \dots + C(p, y) として
  2.  S(b) - S(a-1) を計算する

という風にすることで、総和をとる添字のスタートを 1 とかにしてしまうことができる。でも  x \lt y のときは  C(x, y) = 0 であって意味がないので

  •  S(p) = C(y, y) + C(y+1, y) + \dots + C(p, y)

を求めることにする。二通りのやり方で求めてみる

  1. 経路数への帰着
  2. 流行りの多項式考察

経路数への帰着

 S(p) は、下図のような  p-y \times y+1 の格子を左上から右下まで最短で行く方法のうち、黒丸のうちのどこを通るかを場合分けして表したものとみなすことができる。よって

  •  S(p) = C(p+1, y+1)

が成立する。

多項式で求める

二項係数がきたら二項定理を考えるのは大アリだと思う。

  •  C(y, y) は、 (1 + x)^{y} x^{y} の係数
  •  C(y+1, y) は、 (1 + x)^{y+1} x^{y} の係数
  •  C(y+2, y) は、 (1 + x)^{y+2} x^{y} の係数
  • ...
  •  C(p, y) は、 (1 + x)^{p} x^{y} の係数

ということなので、 S(p) を求めるためには、これらの式を合計して  x^{y} の係数を考えれば OK。

 (1 + x)^{y} + (1 + x)^{y + 1} + \dots + (1 + x)^{p}
=  \frac{(1 +x)^{p+1} - (1 + x)^{y}}{x}

となる。この  x^{y} の係数を求めればよい。これは分子の右側は関係なくて  (1 + x)^{p+1} x^{y+1} の係数に一致する。よって

  •  S(p) = C(p+1, y+1)

が成立する。

まとめ

いずれにしても

  •  C(y, y) + C(y+1, y) + \dots + C(p, y) = C(p+1, y+1)

であることがわかったので、これを利用すれば各  y1 \le y \le y2 に対して、

  •  C(y, y) + C(y+1, y) + \dots + C(x2, y) = C(x2+y+1, y+1) から
  •  C(y, y) + C(y+1, y) + \dots + C(x1-1, y) = C(x1+1, y+1) を引いたもの

を合計すればよい。なおここからさらに式を簡単にもできるけど、ここまででも十分かなと。

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

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() { 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 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;

int main() {
    bc.init(2100000);
    long long x1, y1, x2, y2;
    cin >> x1 >> y1 >> x2 >> y2;
    mint res = 0;
    for (int r = y1; r <= y2; ++r) {
        mint tmp = bc.com(x2 + r + 1, r + 1) - bc.com(x1 + r, r + 1);
        res += tmp;
    }
    cout << res << endl;
}