#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <functional>
#include <ctime>
using namespace std;
namespace bitmatrix {
typedef unsigned long long ull;
struct mat {
int n, m;
vector<vector<ull>> x;
mat(int m, int n) : m(m), n(n), x(1+m/8, vector<ull>(1+n/8)) { }
bool get(int i, int j) const {
return x[i/8][j/8] & (1ull << (8*(i%8)+(j%8)));
}
void set(int i, int j, int b) {
if (b) x[i/8][j/8] |= (1ull << (8*(i%8)+(j%8)));
else x[i/8][j/8] &= ~(1ull << (8*(i%8)+(j%8)));
}
};
ostream &operator<<(ostream &os, const mat &A) {
for (int i = 0; i < A.m; ++i) {
for (int j = 0; j < A.n; ++j)
os << A.get(i, j);
os << endl;
}
return os;
}
mat eye(int n) {
mat I(n, n);
for (int i = 0; i < I.x.size(); ++i)
I.x[i][i] = 0x8040201008040201;
return I;
}
mat add(mat A, const mat &B) {
for (int i = 0; i < A.x.size(); ++i)
for (int j = 0; j < A.x[0].size(); ++j)
A.x[i][j] |= B.x[i][j];
return A;
}
ull mul(ull a, ull b) { // C[i][j] |= A[i][k] & B[k][j]
const ull u = 0xff, v = 0x101010101010101;
ull c = 0;
for (;a && b; a >>= 1, b >>= 8)
c |= (((a & v) * u) & ((b & u) * v));
return c;
}
mat mul(mat A, mat B) {
mat C(A.n, B.m);
for (int i = 0; i < A.x.size(); ++i)
for (int k = 0; k < B.x.size(); ++k)
for (int j = 0; j < B.x[0].size(); ++j)
C.x[i][j] |= mul(A.x[i][k], B.x[k][j]);
return C;
}
mat pow(mat A, int k) {
mat X = eye(A.n);
for (; k > 0; k >>= 1) {
if (k & 1) X = mul(X, A);
A = mul(A, A);
}
return X;
}
ull transpose(ull a) {
ull t = (a ^ (a >> 7)) & 0xaa00aa00aa00aa;
a = a ^ t ^ (t << 7);
t = (a ^ (a >> 14)) & 0xcccc0000cccc;
a = a ^ t ^ (t << 14);
t = (a ^ (a >> 28)) & 0xf0f0f0f0;
a = a ^ t ^ (t << 28);
return a;
}
mat transpose(mat A) {
mat B(A.m, A.n);
for (int i = 0; i < A.x.size(); ++i)
for (int j = 0; j < A.x[0].size(); ++j)
B.x[j][i] = transpose(A.x[i][j]);
return B;
}
}
namespace vector_bool {
typedef vector<bool> vec;
typedef vector<vec> mat;
mat eye(int n) {
mat I(n, vec(n));
for (int i = 0; i < n; ++i)
I[i][i] = 1;
return I;
}
mat add(mat A, const mat &B) {
for (int i = 0; i < A.size(); ++i)
for (int j = 0; j < A[0].size(); ++j)
A[i][j] = (A[i][j] | B[i][j]);
return A;
}
mat mul(mat A, const mat &B) {
for (int i = 0; i < A.size(); ++i) {
vec x(A[0].size());
for (int k = 0; k < B.size(); ++k)
for (int j = 0; j < B[0].size(); ++j)
x[j] = (x[j] | (A[i][k] & B[k][j]));
A[i].swap(x);
}
return A;
}
mat pow(mat A, int k) {
mat X = eye(A.size());
for (; k > 0; k >>= 1) {
if (k & 1) X = mul(X, A);
A = mul(A, A);
}
return X;
}
}
int main() {
for (int n = 1; n < 512; n *= 2) {
printf("%d\t", n);
{
using namespace bitmatrix;
mat A(n, n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
A.set(i, j, rand() % 2);
}
}
double _time = clock();
A = pow(add(A, eye(n)), n); // (A + I)^n is a transitive closure
_time = clock() - _time;
printf("%f\t", _time / CLOCKS_PER_SEC);
}
{
using namespace vector_bool;
mat A(n, vec(n));
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
A[i][j] = rand() % 2;
}
A[i][i] = 1;
}
double _time = clock();
A = pow(add(A, eye(n)), n); // (A + I)^n is a transitive closure
_time = clock() - _time;
printf("%f\n", _time / CLOCKS_PER_SEC);
}
}
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8dmVjdG9yPgojaW5jbHVkZSA8Y3N0ZGlvPgojaW5jbHVkZSA8YWxnb3JpdGhtPgojaW5jbHVkZSA8ZnVuY3Rpb25hbD4KI2luY2x1ZGUgPGN0aW1lPgoKdXNpbmcgbmFtZXNwYWNlIHN0ZDsKCm5hbWVzcGFjZSBiaXRtYXRyaXggewp0eXBlZGVmIHVuc2lnbmVkIGxvbmcgbG9uZyB1bGw7CnN0cnVjdCBtYXQgewogIGludCBuLCBtOwogIHZlY3Rvcjx2ZWN0b3I8dWxsPj4geDsKICBtYXQoaW50IG0sIGludCBuKSA6IG0obSksIG4obiksIHgoMSttLzgsIHZlY3Rvcjx1bGw+KDErbi84KSkgeyB9CiAgYm9vbCBnZXQoaW50IGksIGludCBqKSBjb25zdCB7IAogICAgcmV0dXJuIHhbaS84XVtqLzhdICYgKDF1bGwgPDwgKDgqKGklOCkrKGolOCkpKTsKICB9CiAgdm9pZCBzZXQoaW50IGksIGludCBqLCBpbnQgYikgewogICAgaWYgKGIpIHhbaS84XVtqLzhdIHw9ICAoMXVsbCA8PCAoOCooaSU4KSsoaiU4KSkpOyAKICAgIGVsc2UgICB4W2kvOF1bai84XSAmPSB+KDF1bGwgPDwgKDgqKGklOCkrKGolOCkpKTsKICB9IAp9Owpvc3RyZWFtICZvcGVyYXRvcjw8KG9zdHJlYW0gJm9zLCBjb25zdCBtYXQgJkEpIHsKICBmb3IgKGludCBpID0gMDsgaSA8IEEubTsgKytpKSB7CiAgICBmb3IgKGludCBqID0gMDsgaiA8IEEubjsgKytqKSAKICAgICAgb3MgPDwgQS5nZXQoaSwgaik7CiAgICBvcyA8PCBlbmRsOwogIH0KICByZXR1cm4gb3M7Cn0KbWF0IGV5ZShpbnQgbikgewogIG1hdCBJKG4sIG4pOwogIGZvciAoaW50IGkgPSAwOyBpIDwgSS54LnNpemUoKTsgKytpKSAKICAgIEkueFtpXVtpXSA9IDB4ODA0MDIwMTAwODA0MDIwMTsKICByZXR1cm4gSTsKfQptYXQgYWRkKG1hdCBBLCBjb25zdCBtYXQgJkIpIHsgCiAgZm9yIChpbnQgaSA9IDA7IGkgPCBBLnguc2l6ZSgpOyArK2kpCiAgICBmb3IgKGludCBqID0gMDsgaiA8IEEueFswXS5zaXplKCk7ICsraikgCiAgICAgIEEueFtpXVtqXSB8PSBCLnhbaV1bal07CiAgcmV0dXJuIEE7Cn0KdWxsIG11bCh1bGwgYSwgdWxsIGIpIHsgLy8gQ1tpXVtqXSB8PSBBW2ldW2tdICYgQltrXVtqXQogIGNvbnN0IHVsbCB1ID0gMHhmZiwgdiA9IDB4MTAxMDEwMTAxMDEwMTAxOwogIHVsbCBjID0gMDsKICBmb3IgKDthICYmIGI7IGEgPj49IDEsIGIgPj49IDgpCiAgICBjIHw9ICgoKGEgJiB2KSAqIHUpICYgKChiICYgdSkgKiB2KSk7CiAgcmV0dXJuIGM7Cn0KbWF0IG11bChtYXQgQSwgbWF0IEIpIHsKICBtYXQgQyhBLm4sIEIubSk7CiAgZm9yIChpbnQgaSA9IDA7IGkgPCBBLnguc2l6ZSgpOyArK2kpIAogICAgZm9yIChpbnQgayA9IDA7IGsgPCBCLnguc2l6ZSgpOyArK2spIAogICAgICBmb3IgKGludCBqID0gMDsgaiA8IEIueFswXS5zaXplKCk7ICsraikgCiAgICAgICAgQy54W2ldW2pdIHw9IG11bChBLnhbaV1ba10sIEIueFtrXVtqXSk7CiAgcmV0dXJuIEM7Cn0KbWF0IHBvdyhtYXQgQSwgaW50IGspIHsKICBtYXQgWCA9IGV5ZShBLm4pOwogIGZvciAoOyBrID4gMDsgayA+Pj0gMSkgewogICAgaWYgKGsgJiAxKSBYID0gbXVsKFgsIEEpOwogICAgQSA9IG11bChBLCBBKTsKICB9CiAgcmV0dXJuIFg7Cn0KCnVsbCB0cmFuc3Bvc2UodWxsIGEpIHsKICB1bGwgdCA9IChhIF4gKGEgPj4gNykpICYgMHhhYTAwYWEwMGFhMDBhYTsKICBhID0gYSBeIHQgXiAodCA8PCA3KTsKICB0ID0gKGEgXiAoYSA+PiAxNCkpICYgMHhjY2NjMDAwMGNjY2M7CiAgYSA9IGEgXiB0IF4gKHQgPDwgMTQpOwogIHQgPSAoYSBeIChhID4+IDI4KSkgJiAweGYwZjBmMGYwOwogIGEgPSBhIF4gdCBeICh0IDw8IDI4KTsKICByZXR1cm4gYTsKfQptYXQgdHJhbnNwb3NlKG1hdCBBKSB7CiAgbWF0IEIoQS5tLCBBLm4pOwogIGZvciAoaW50IGkgPSAwOyBpIDwgQS54LnNpemUoKTsgKytpKSAKICAgIGZvciAoaW50IGogPSAwOyBqIDwgQS54WzBdLnNpemUoKTsgKytqKSAKICAgICAgQi54W2pdW2ldID0gdHJhbnNwb3NlKEEueFtpXVtqXSk7CiAgcmV0dXJuIEI7Cn0KfQoKbmFtZXNwYWNlIHZlY3Rvcl9ib29sIHsKdHlwZWRlZiB2ZWN0b3I8Ym9vbD4gdmVjOwp0eXBlZGVmIHZlY3Rvcjx2ZWM+IG1hdDsKbWF0IGV5ZShpbnQgbikgewogIG1hdCBJKG4sIHZlYyhuKSk7CiAgZm9yIChpbnQgaSA9IDA7IGkgPCBuOyArK2kpCiAgICBJW2ldW2ldID0gMTsKICByZXR1cm4gSTsKfQptYXQgYWRkKG1hdCBBLCBjb25zdCBtYXQgJkIpIHsKICBmb3IgKGludCBpID0gMDsgaSA8IEEuc2l6ZSgpOyArK2kpCiAgICBmb3IgKGludCBqID0gMDsgaiA8IEFbMF0uc2l6ZSgpOyArK2opCiAgICAgIEFbaV1bal0gPSAoQVtpXVtqXSB8IEJbaV1bal0pOwogIHJldHVybiBBOwp9Cm1hdCBtdWwobWF0IEEsIGNvbnN0IG1hdCAmQikgewogIGZvciAoaW50IGkgPSAwOyBpIDwgQS5zaXplKCk7ICsraSkgewogICAgdmVjIHgoQVswXS5zaXplKCkpOwogICAgZm9yIChpbnQgayA9IDA7IGsgPCBCLnNpemUoKTsgKytrKSAKICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBCWzBdLnNpemUoKTsgKytqKSAKICAgICAgICB4W2pdID0gKHhbal0gfCAoQVtpXVtrXSAmIEJba11bal0pKTsKICAgIEFbaV0uc3dhcCh4KTsKICB9CiAgcmV0dXJuIEE7Cn0KbWF0IHBvdyhtYXQgQSwgaW50IGspIHsKICBtYXQgWCA9IGV5ZShBLnNpemUoKSk7CiAgZm9yICg7IGsgPiAwOyBrID4+PSAxKSB7CiAgICBpZiAoayAmIDEpIFggPSBtdWwoWCwgQSk7CiAgICBBID0gbXVsKEEsIEEpOwogIH0KICByZXR1cm4gWDsKfQp9CgppbnQgbWFpbigpIHsKICBmb3IgKGludCBuID0gMTsgbiA8IDUxMjsgbiAqPSAyKSB7CiAgICBwcmludGYoIiVkXHQiLCBuKTsKICAgIHsKICAgICAgdXNpbmcgbmFtZXNwYWNlIGJpdG1hdHJpeDsKICAgICAgbWF0IEEobiwgbik7CiAgICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbjsgKytpKSB7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBuOyArK2opIHsKICAgICAgICAgIEEuc2V0KGksIGosIHJhbmQoKSAlIDIpOwogICAgICAgIH0KICAgICAgfQogICAgICBkb3VibGUgX3RpbWUgPSBjbG9jaygpOwogICAgICBBID0gcG93KGFkZChBLCBleWUobikpLCBuKTsgLy8gKEEgKyBJKV5uIGlzIGEgdHJhbnNpdGl2ZSBjbG9zdXJlCiAgICAgIF90aW1lID0gY2xvY2soKSAtIF90aW1lOwogICAgICBwcmludGYoIiVmXHQiLCBfdGltZSAvIENMT0NLU19QRVJfU0VDKTsKICAgIH0KICAgIHsKICAgICAgdXNpbmcgbmFtZXNwYWNlIHZlY3Rvcl9ib29sOwogICAgICBtYXQgQShuLCB2ZWMobikpOwogICAgICBmb3IgKGludCBpID0gMDsgaSA8IG47ICsraSkgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgbjsgKytqKSB7CiAgICAgICAgICBBW2ldW2pdID0gcmFuZCgpICUgMjsKICAgICAgICB9CiAgICAgICAgQVtpXVtpXSA9IDE7CiAgICAgIH0KICAgICAgZG91YmxlIF90aW1lID0gY2xvY2soKTsKICAgICAgQSA9IHBvdyhhZGQoQSwgZXllKG4pKSwgbik7IC8vIChBICsgSSlebiBpcyBhIHRyYW5zaXRpdmUgY2xvc3VyZQogICAgICBfdGltZSA9IGNsb2NrKCkgLSBfdGltZTsKICAgICAgcHJpbnRmKCIlZlxuIiwgX3RpbWUgLyBDTE9DS1NfUEVSX1NFQyk7CiAgICB9CiAgfQp9