aoj2397 Three-way Branch

問題ページ
障害物がなければ行列累乗をすればよい．障害物に対応するため障害物が存在する行ごとに区切って考える．行列累乗で 列目に到達する方法を求めたあと，障害物が存在する位置に到達する方法を0通りに置き換える． で解けた．

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG_   
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
const ll INF = 1LL<<60;  
  
template<ll MOD>  
struct modint {  
    ll x;  
    modint(): x(0) {}  
    modint(ll y) : x(y>=0 ? y%MOD : y%MOD+MOD) {}  
    static constexpr ll mod() { return MOD; }  
    // e乗  
    modint pow(ll e) {  
        ll a = 1, p = x;  
        while(e > 0) {  
            if(e%2 == 0) {p = (p*p) % MOD; e /= 2;}  
            else {a = (a*p) % MOD; e--;}  
        }  
        return modint(a);  
    }  
    modint inv() const {  
        ll a=x, b=MOD, u=1, y=1, v=0, z=0;  
        while(a) {  
            ll q = b/a;  
            swap(z -= q*u, u);  
            swap(y -= q*v, v);  
            swap(b -= q*a, a);  
        }  
        return z;  
    }  
    // Comparators  
    bool operator <(modint b) { return x < b.x; }  
    bool operator >(modint b) { return x > b.x; }  
    bool operator<=(modint b) { return x <= b.x; }  
    bool operator>=(modint b) { return x >= b.x; }  
    bool operator!=(modint b) { return x != b.x; }  
    bool operator==(modint b) { return x == b.x; }  
    // Basic Operations  
    modint operator+(modint r) const { return modint(*this) += r; }  
    modint operator-(modint r) const { return modint(*this) -= r; }  
    modint operator*(modint r) const { return modint(*this) *= r; }  
    modint operator/(modint r) const { return modint(*this) /= r; }  
    modint &operator+=(modint r) {  
        if((x += r.x) >= MOD) x -= MOD;  
        return *this;  
    }  
    modint &operator-=(modint r) {  
        if((x -= r.x) < 0) x += MOD;  
        return *this;  
    }  
    modint &operator*=(modint r) {  
    #if !defined(_WIN32) || defined(_WIN64)  
        x = x * r.x % MOD; return *this;  
    #endif  
        unsigned long long y = x * r.x;  
        unsigned xh = (unsigned) (y >> 32), xl = (unsigned) y, d, m;  
        asm(  
            "divl %4; \n\t"  
            : "=a" (d), "=d" (m)  
            : "d" (xh), "a" (xl), "r" (MOD)  
        );  
        x = m;  
        return *this;  
    }  
    modint &operator/=(modint r) { return *this *= r.inv(); }  
    // increment, decrement  
    modint operator++() { x++; return *this; }  
    modint operator++(signed) { modint t = *this; x++; return t; }  
    modint operator--() { x--; return *this; }  
    modint operator--(signed) { modint t = *this; x--; return t; }  
    // 平方剰余のうち一つを返す なければ-1  
    friend modint sqrt(modint a) {  
        if(a == 0) return 0;  
        ll q = MOD-1, s = 0;  
        while((q&1)==0) q>>=1, s++;  
        modint z=2;  
        while(1) {  
            if(z.pow((MOD-1)/2) == MOD-1) break;  
            z++;  
        }  
        modint c = z.pow(q), r = a.pow((q+1)/2), t = a.pow(q);  
        ll m = s;  
        while(t.x>1) {  
            modint tp=t;  
            ll k=-1;  
            FOR(i, 1, m) {  
                tp *= tp;  
                if(tp == 1) { k=i; break; }  
            }  
            if(k==-1) return -1;  
            modint cp=c;  
            REP(i, m-k-1) cp *= cp;  
            c = cp*cp, t = c*t, r = cp*r, m = k;  
        }  
        return r.x;  
    }  
  
    template<class T>  
    friend modint operator*(T l, modint r) { return modint(l) *= r; }  
    template<class T>  
    friend modint operator+(T l, modint r) { return modint(l) += r; }  
    template<class T>  
    friend modint operator-(T l, modint r) { return modint(l) -= r; }  
    template<class T>  
    friend modint operator/(T l, modint r) { return modint(l) /= r; }  
    template<class T>  
    friend bool operator==(T l, modint r) { return modint(l) == r; }  
    template<class T>  
    friend bool operator!=(T l, modint r) { return modint(l) != r; }  
    // Input/Output  
    friend ostream &operator<<(ostream& os, modint a) { return os << a.x; }  
    friend istream &operator>>(istream& is, modint &a) {   
        is >> a.x;  
        a.x = ((a.x%MOD)+MOD)%MOD;  
        return is;  
    }  
    friend string to_frac(modint v) {  
        static map<ll, PII> mp;  
        if(mp.empty()) {  
            mp[0] = mp[MOD] = {0, 1};  
            FOR(i, 2, 1001) FOR(j, 1, i) if(__gcd(i, j) == 1) {  
                mp[(modint(i) / j).x] = {i, j};  
            }  
        }  
        auto itr = mp.lower_bound(v.x);  
        if(itr != mp.begin() && v.x - prev(itr)->first < itr->first - v.x) --itr;  
        string ret = to_string(itr->second.first + itr->second.second * ((int)v.x - itr->first));  
        if(itr->second.second > 1) {  
            ret += '/';  
            ret += to_string(itr->second.second);  
        }  
        return ret;  
    }  
};  
using mint = modint<1000000009>;  
  
// 有限体の行列  
struct matrix {  
    int h, w;  
    vector<mint> dat;  
    matrix() {}  
    matrix(int h, int w) : h(h), w(w), dat(h*w) {}  
    mint& get(int y, int x) { return dat[y*w+x]; }  
    mint get(int y, int x) const { return dat[y*w+x]; }  
  
    matrix& operator+=(const matrix& r) {  
        assert(h==r.h && w==r.w);  
        REP(i, h*w) dat[i] += r.dat[i];  
        return *this;  
    }  
    matrix& operator-=(const matrix& r) {  
        assert(h==r.h && w==r.w);  
        REP(i, h*w) dat[i] -= r.dat[i];  
        return *this;  
    }  
    matrix& operator*=(const matrix& r) {  
        assert(w==r.h);  
        matrix ret(h, w);  
        REP(i, h) REP(j, r.w) REP(k, w) {  
            ret.dat[i*r.w+j] += dat[i*w+k] * r.dat[k*r.w+j];  
        }  
        return (*this) = ret;  
    }  
    matrix operator+(const matrix& r) { return matrix(*this) += r; }  
    matrix operator-(const matrix& r) { return matrix(*this) -= r; }  
    matrix operator*(const matrix& r) { return matrix(*this) *= r; }  
    bool operator==(const matrix& a) { return dat==a.dat; }  
    bool operator!=(const matrix& a) { return dat!=a.dat; }  
  
    friend matrix pow(matrix p, ll n) {  
        matrix ret(p.h, p.w);  
        REP(i, p.h) ret.get(i, i) = 1;  
        while(n > 0) {  
            if(n&1) {ret *= p; n--;}  
            else {p *= p; n >>= 1;}  
        }  
        return ret;  
    }  
    // 階段行列を求める O(HW^2)  
    friend int gauss_jordan(matrix& a) {  
        int rank = 0;  
        REP(i, a.w) {  
            int pivot = -1;  
            FOR(j, rank, a.h) if(a.get(j,i) != 0) { pivot = j; break; }  
            if(pivot == -1) continue;  
            REP(j, a.w) swap(a.get(rank,j), a.get(pivot,j));  
            const mint inv = a.get(rank,i).inv();  
            REP(j, a.w) a.get(rank,j) *= inv;  
            REP(j, a.h) if(j != rank && a.get(j,i) != 0) {  
                const mint num = a.get(j,i);  
                REP(k, a.w) a.get(j,k) -= a.get(rank,k) * num;  
            }  
            rank++;  
        }  
        return rank;  
    }  
  
    friend ostream &operator<<(ostream& os, matrix a) {  
        os << endl;  
        REP(i, a.h) {  
            REP(j, a.w) os << a.get(i,j) << " ";  
            if(i+1!=a.h) os << endl;  
        }  
        return os;  
    }  
};  
  
ll testcase;  
void solve() {  
    ll w, h, n;  
    cin >> w >> h >> n;  
    if(w == 0) exit(0);  
    map<ll, vector<ll>> mp;  
    REP(i, n) {  
        ll x, y;  
        cin >> x >> y;  
        x--, y--;  
        mp[y].push_back(x);  
    }  
    mp[h-1].push_back(0);  
  
    matrix mat(w, w);  
    REP(i, w) for(ll j=max(0LL,i-1); j<min(i+2,w); ++j) mat.get(i,j) = 1;  
  
    ll pre = 0;  
    vector<mint> cnt(w);  
    cnt[0] = 1;  
    for(auto p: mp) {  
        ll y = p.first;  
        auto v = p.second;  
  
        matrix pw = pow(mat, y-pre);  
        vector<mint> ncnt(w);  
        REP(i, w) REP(j, w) ncnt[i] += pw.get(i,j) * cnt[j];  
        swap(cnt, ncnt);  
  
        for(auto i: v) cnt[i] = 0;  
        pre = y;  
    }  
  
    cout << "Case " << testcase+1 << ": " << cnt[w-1] << endl;  
    testcase++;  
}  
  
int main(void) {  
    while(1) solve();  
  
    return 0;  
}