fork(1) download
copy 
  1. #include <iostream>
  2. #include <vector>
  3. #include <cstdio>
  4. #include <algorithm>
  5. #include <functional>
  6. #include <ctime>
  7.  
  8. using namespace std;
  9.  
  10. namespace bitmatrix {
  11. typedef unsigned long long ull;
  12. struct mat {
  13.   int n, m;
  14.   vector<vector<ull>> x;
  15.   mat(int m, int n) : m(m), n(n), x(1+m/8, vector<ull>(1+n/8)) { }
  16.   bool get(int i, int j) const { 
  17.     return x[i/8][j/8] & (1ull << (8*(i%8)+(j%8)));
  18.   }
  19.   void set(int i, int j, int b) {
  20.     if (b) x[i/8][j/8] |=  (1ull << (8*(i%8)+(j%8))); 
  21.     else   x[i/8][j/8] &= ~(1ull << (8*(i%8)+(j%8)));
  22.   } 
  23. };
  24. ostream &operator<<(ostream &os, const mat &A) {
  25.   for (int i = 0; i < A.m; ++i) {
  26.     for (int j = 0; j < A.n; ++j) 
  27.       os << A.get(i, j);
  28.     os << endl;
  29.   }
  30.   return os;
  31. }
  32. mat eye(int n) {
  33.   mat I(n, n);
  34.   for (int i = 0; i < I.x.size(); ++i) 
  35.     I.x[i][i] = 0x8040201008040201;
  36.   return I;
  37. }
  38. mat add(mat A, const mat &B) { 
  39.   for (int i = 0; i < A.x.size(); ++i)
  40.     for (int j = 0; j < A.x[0].size(); ++j) 
  41.       A.x[i][j] |= B.x[i][j];
  42.   return A;
  43. }
  44. ull mul(ull a, ull b) { // C[i][j] |= A[i][k] & B[k][j]
  45.   const ull u = 0xff, v = 0x101010101010101;
  46.   ull c = 0;
  47.   for (;a && b; a >>= 1, b >>= 8)
  48.     c |= (((a & v) * u) & ((b & u) * v));
  49.   return c;
  50. }
  51. mat mul(mat A, mat B) {
  52.   mat C(A.n, B.m);
  53.   for (int i = 0; i < A.x.size(); ++i) 
  54.     for (int k = 0; k < B.x.size(); ++k) 
  55.       for (int j = 0; j < B.x[0].size(); ++j) 
  56.         C.x[i][j] |= mul(A.x[i][k], B.x[k][j]);
  57.   return C;
  58. }
  59. mat pow(mat A, int k) {
  60.   mat X = eye(A.n);
  61.   for (; k > 0; k >>= 1) {
  62.     if (k & 1) X = mul(X, A);
  63.     A = mul(A, A);
  64.   }
  65.   return X;
  66. }
  67.  
  68. ull transpose(ull a) {
  69.   ull t = (a ^ (a >> 7)) & 0xaa00aa00aa00aa;
  70.   a = a ^ t ^ (t << 7);
  71.   t = (a ^ (a >> 14)) & 0xcccc0000cccc;
  72.   a = a ^ t ^ (t << 14);
  73.   t = (a ^ (a >> 28)) & 0xf0f0f0f0;
  74.   a = a ^ t ^ (t << 28);
  75.   return a;
  76. }
  77. mat transpose(mat A) {
  78.   mat B(A.m, A.n);
  79.   for (int i = 0; i < A.x.size(); ++i) 
  80.     for (int j = 0; j < A.x[0].size(); ++j) 
  81.       B.x[j][i] = transpose(A.x[i][j]);
  82.   return B;
  83. }
  84. }
  85.  
  86. namespace vector_bool {
  87. typedef vector<bool> vec;
  88. typedef vector<vec> mat;
  89. mat eye(int n) {
  90.   mat I(n, vec(n));
  91.   for (int i = 0; i < n; ++i)
  92.     I[i][i] = 1;
  93.   return I;
  94. }
  95. mat add(mat A, const mat &B) {
  96.   for (int i = 0; i < A.size(); ++i)
  97.     for (int j = 0; j < A[0].size(); ++j)
  98.       A[i][j] = (A[i][j] | B[i][j]);
  99.   return A;
  100. }
  101. mat mul(mat A, const mat &B) {
  102.   for (int i = 0; i < A.size(); ++i) {
  103.     vec x(A[0].size());
  104.     for (int k = 0; k < B.size(); ++k) 
  105.       for (int j = 0; j < B[0].size(); ++j) 
  106.         x[j] = (x[j] | (A[i][k] & B[k][j]));
  107.     A[i].swap(x);
  108.   }
  109.   return A;
  110. }
  111. mat pow(mat A, int k) {
  112.   mat X = eye(A.size());
  113.   for (; k > 0; k >>= 1) {
  114.     if (k & 1) X = mul(X, A);
  115.     A = mul(A, A);
  116.   }
  117.   return X;
  118. }
  119. }
  120.  
  121. int main() {
  122.   for (int n = 1; n < 512; n *= 2) {
  123.     printf("%d\t", n);
  124.     {
  125.       using namespace bitmatrix;
  126.       mat A(n, n);
  127.       for (int i = 0; i < n; ++i) {
  128.         for (int j = 0; j < n; ++j) {
  129.           A.set(i, j, rand() % 2);
  130.         }
  131.       }
  132.       double _time = clock();
  133.       A = pow(add(A, eye(n)), n); // (A + I)^n is a transitive closure
  134.       _time = clock() - _time;
  135.       printf("%f\t", _time / CLOCKS_PER_SEC);
  136.     }
  137.     {
  138.       using namespace vector_bool;
  139.       mat A(n, vec(n));
  140.       for (int i = 0; i < n; ++i) {
  141.         for (int j = 0; j < n; ++j) {
  142.           A[i][j] = rand() % 2;
  143.         }
  144.         A[i][i] = 1;
  145.       }
  146.       double _time = clock();
  147.       A = pow(add(A, eye(n)), n); // (A + I)^n is a transitive closure
  148.       _time = clock() - _time;
  149.       printf("%f\n", _time / CLOCKS_PER_SEC);
  150.     }
  151.   }
  152. }
Success #stdin #stdout 1.2s 3472KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1	0.000009	0.000004
2	0.000003	0.000004
4	0.000004	0.000012
8	0.000009	0.000040
16	0.000017	0.000210
32	0.000063	0.001596
64	0.000429	0.013526
128	0.003398	0.117845
256	0.029370	1.028783
https://ideone.com/8AsuI2
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!