###### tags: `HUPC2021`

# B - チェーンレター 解説

###### 原案: [Nerve](https://twitter.com/nrvkpr) | 解説: [Nerve](https://twitter.com/nrvkpr)
<br>

$d$ 日目に友人 $i$ の手紙が友人 $j$ に届けられるかどうかをbool値として考えて、加法をBitwise OR、乗法をBitwise ANDとした行列累乗をすることで解くことができます。

各行・列はある友人が送った手紙が別の友人に届くまでの日数の最大値で拡張します（制約から、行数・列数は1100以下になります）。
行列全体が $d$ 日経過後の状態、行rが(友人 $i$, $d_i$ 日経過後)、列cが(友人 $j$, $d_j$ 日経過後)をさすとき、r行c列目の要素は「$d+d_i$ 日経過した友人 $i$ が持っている文字列が $d+d_j$ 日経過後の友人 $j$ に渡されるかどうか」を表します。
$1$ 日経過後の状態で1を立てる要素は以下の通りです。

- 行rが(友人 $i$, $d_i$ 日経過後)、行cが(友人 $i$, $d_i+1$ 日経過後) を表す要素
- 行rが(友人 $i$, $d_{ij} - 1$ 日経過後)、行cが(友人 $v_{ij}$, $0$ 日経過後) を表す要素

この行列累乗を適切に実装すると、$K$ を行数として $\Theta(K^3\log D)$ で計算できますが、このままだと実行時間制限に間に合いません。

各要素がbool値で加法がBitwise OR、乗法がBitwise ANDであることに注目すると、行・列ごとにbitsetでまとめて計算することで高速化することができます。
bitsetの配列で行列を計算することで、行列の積が $\Theta(K^3/wordsize)$ で計算できます。
行列累乗全体の計算量は $\Theta(K^3\log D/wordsize)$ です。

以下実装例

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

const size_t MAT_ORDER = 1100;

struct BitsetSquareMatrix {
    using BSMatrix = BitsetSquareMatrix;
    using matrix = array<bitset<MAT_ORDER>, MAT_ORDER>;

    matrix rmat, cmat;

    BitsetSquareMatrix(){}

    inline void set(int r, int c) {
        rmat[r].set(c);
        cmat[c].set(r);
    }

    inline bool get(int r, int c) {
        return rmat[r][c];
    }

    inline BSMatrix &operator*=(const BSMatrix &A) {
        matrix rres, cres;
        for (int i = 0; i < MAT_ORDER; i++){
            for (int j = 0; j < MAT_ORDER; j++){
                rres[i][j] = (rmat[i] & A.cmat[j]).any();
                cres[j][i] = rres[i][j];
            }
        }

        rmat.swap(rres);
        cmat.swap(cres);

        return (*this);
    }

    inline BSMatrix &operator^=(long long p) {
        BSMatrix res = BSMatrix();
        for (int i = 0; i < MAT_ORDER; i++){
            res.rmat[i].set(i);
            res.cmat[i].set(i);
        }

        while(p > 0){
            if (p & 1) res *= *this;
            *this *= *this;
            p >>= 1;
        }

        rmat.swap(res.rmat);
        cmat.swap(res.cmat);

        return (*this);
    }

    BSMatrix operator*(const BSMatrix &B) const { return (BSMatrix(*this) *= B); }
    BSMatrix operator^(const long long k) const { return (BSMatrix(*this) ^= k); }

    friend ostream &operator<<(ostream &os, BSMatrix &p) {
        for (int i = 0; i < MAT_ORDER; i++) {
            os << "[";
            auto str = p.rmat[i].to_string();
            reverse(str.begin(), str.end());
            os << str;
            os << "]\n";
        }
        return (os);
    }
};


int main() {
    int N, D;
    cin >> N >> D;
    vector maxD(N, 1);

    vector<vector<int>> v(N), d(N);
    for (int i = 0; i < N; i++) {
        int K;
        cin >> K;
        v[i].resize(K), d[i].resize(K);
        for (int j = 0; j < K; j++) {
            cin >> v[i][j] >> d[i][j];
            v[i][j]--;
            maxD[i] = max(maxD[i], d[i][j]);
        }
    }

    vector dacc(N + 1, 0);
    for (int i = 0; i < N; i++)
        dacc[i + 1] = dacc[i] + maxD[i];

    // idxのday日後に対応する行列Aのインデックスを求める
    auto get_id = [&](int idx, int day) -> int {
        assert(dacc[idx] + day >= 0);
        return dacc[idx] + day;
    };


    BitsetSquareMatrix A;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < v[i].size(); j++) {
            A.set(get_id(i, d[i][j] - 1), get_id(v[i][j], 0));
        }
    }
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < maxD[i] - 1; j++){
            A.set(get_id(i, j), get_id(i, j + 1));
        }
    }

    A ^= D;

    vector<int> ans;
    for (int i = 0; i < N; i++) {
        bool isok = true;
        for (int j = 0; j < N; j++) {
            if (!A.get(get_id(j, 0), get_id(i, 0))) isok = false;
        }

        if (isok) ans.emplace_back(i);
    }

    cout << ans.size() << endl;
    for (auto a : ans) cout << a + 1 << " ";
    cout << endl;
}

<!-- </div> </details> -->