Typical DP Contest - M 家

M: 家 - Typical DP Contest | AtCoder

• bit dp + 繰り返し二乗法。 bit dp すらまともに書けなかった。orz

template<typename T>
using V = vector<T>;
template<typename T>
using VV = V<V<T>>;

// O(N^3)
VV<ll> mul(const VV<ll> &a, const VV<ll> &b, const int M) {
        int n = a.size();
        VV<ll> c(n,V<ll>(n,0));
        REP(i, n)
        REP(j, n)
        REP(k, n)
            (c[i][j] += a[i][k]*b[k][j]) %= M;
        return move(c);
}
// O(N^3*log(k))
VV<ll>  mpow(VV<ll> &&a, int k, int M) {
        int n = a.size();
        VV<ll> b(n,V<ll>(n,0));
        REP(i, n) b[i][i] = 1;

        while (k > 0)
        {
            if (k & 1)
                b = mul(b,a,M);
            a = mul(a,a,M);
            k >>= 1;
        }
        return move(b);
}
class m_house {
public:
    void solve(void) {
            const int MOD = (int)(1e+9)+7;
            ll  H;
            int R;
            cin>>H>>R;

            vector<vector<int>> g(R, vector<int>(R,0));
            REP(i, R)
            REP(j, R)
                cin>>g[i][j];

            //
            // 一つの階を i->j まで移動する方法は bit DP で計算できる
            //
            // A[i][j] := i->j まで移動する方法の総数
            //
            // とすると A*A*...*A = A^H が H 階から降りる方法の総数となる。
            // A^H を繰り返し二乗法で計算すればよい
            //
            VV<ll> A(R, V<ll>(R, 0));

            // dp[S][i][j] := i から出発して S の部屋を通過して j へ移動する移動方法の総数
            VV<V<ll>> dp((1<<R), VV<ll>(R, V<ll>(R, 0)));

            REP(i, R)
                dp[(1<<i)][i][i] = 1;

            // O(2^R*R^3) <= 268435456 < 10^9
            // 8秒が制限なんでなんとか間に合いそう
            REP(S,(1<<R)) // S が小さいものからやる ( S|(1<<j) を更新対象にするため)
            REP(i,R)
            REP(k,R)
            {
                //  i->k->j と k 経由で j まで到達するときを加算する
                // S をとおる i->k path がないなら飛ばす
                if (dp[S][i][k] == 0)
                    continue;
                REP(j,R)
                {
                    // 二回同じ部屋kを通過してしまったり、j->k の移動ができないときは飛ばす
                    if ( S&(1<<j) || !g[j][k] )
                        continue;
                    int S2 = (S | (1<<j));
                    (dp[S2][i][j] += dp[S][i][k]) %= MOD;
                }
            }

            REP(i, R)
            REP(j, R)
            REP(S,(1<<R))
                (A[i][j] += dp[S][i][j]) %= MOD;

            // O(log(H)*R^3) <= 86016
            A = mpow(move(A), H, MOD);
            cout<<A[0][0]<<endl;

            return;
    }
};