Add Pollard's Rho and Miller-Rabin

Signed-off-by: Bruno Freitas Tissei <bft15@inf.ufpr.br>

algorithms/graph/centroid_decomposition.cpp

@@ -18,42 +18,42 @@
 vector<int> graph[MAX];
 
 struct CentroidDecomposition {
-  vector<int> par, size, marked;
+  vector<int> par, size, vis;
 
   CentroidDecomposition(int N) :
-    par(N), size(N), marked(N) 
+    par(N), size(N), vis(N) 
   { init(); }
 
   void init() {
-    fill(all(marked), 0);
+    fill(all(vis), 0);
     build(0); // 0-indexed vertices
   }
 
   void build(int x, int p = -1) {
     int n = dfs(x);
-    int centroid = get_centroid(x, n);
+    int c = get_centroid(x, n);
 
-    marked[centroid] = 1;
-    par[centroid] = p;
+    vis[c] = 1;
+    par[c] = p;
 
-    for (auto i : graph[centroid])
-      if (!marked[i])
-        build(i, centroid);
+    for (auto i : graph[c])
+      if (!vis[i])
+        build(i, c);
   }
 
   // Calculates size of every subtree.
   int dfs(int x, int p = -1) {
     size[x] = 1;
     for (auto i : graph[x])
-      if (i != p && !marked[i])
+      if (i != p && !vis[i])
         size[x] += dfs(i, x);
     return size[x];
   }
 
   int get_centroid(int x, int n, int p = -1) {
     for (auto i : graph[x])
-      if (i != p && size[i] > n / 2 && !marked[i])
-        return get_centroid(i, x, n);
+      if (i != p && size[i] > n / 2 && !vis[i])
+        return get_centroid(i, n, x);
     return x;
   }
 

algorithms/math/binary_exponentiation.cpp

@@ -3,14 +3,23 @@
 /// Time: O(log n)
 /// Space: O(1)
 
-template <typename T>
-T bin_exp(T x, ll n) {
-  T ans = 1;
-  while (n) {
-    if (n & 1)
-      ans = ans * x;
-    n >>= 1;
-    x = x * x;
+ll bin_mul(ll a, ll b, ll M = MOD) {
+  ll x = 0;
+  a %= M;
+  while (b) {
+    if (b & 1) x = (x + a) % M;
+    b >>= 1;
+    a = (a * 2) % M;
   }
-  return ans;
+  return x % M;
+}
+
+ll bin_exp(ll a, ll e, ll M = MOD) {
+  ll x = 1;
+  while (e) {
+    if (e & 1) x = bin_mul(x, a, M);
+    e >>= 1;
+    a = bin_mul(a, a, M);
+  }
+  return x % M;
 }

algorithms/math/linear_recurrence.cpp

-/// Linear Recurrence
-///
-/// Time: O(log n)
-/// Space: O(1)
-/// 
-/// Include:
-///   - math/binary_exponentiation
-///   - math/matrix
-
-template <typename T>
-matrix<T> solve(int x, int y, int n) {
-  matrix<T> in(2, 2);
-
-  // Example
-  in[0][0] = x % MOD;
-  in[0][1] = y % MOD;
-  in[1][0] = 1;
-  in[1][1] = 0;
-
-  return fast_pow<T>(in, n);
-}

algorithms/math/miller_rabin.cpp

+/// Miller-Rabin primality test
+///
+/// Time: O(k * log^3 n)
+/// Space: O(1)
+
+const vector<ll> A = {
+  2, 325, 9375, 28178, 450775, 9780504, 1795265022
+};
+
+bool is_prime(ll n) {
+  if (n < 2 || n % 6 % 4 != 1)
+    return n - 2 < 2;
+
+  ll s = __builtin_ctzll(n - 1);
+  ll d = n >> s;
+
+  for (auto a : A) {
+    ll p = bin_exp(a, d, n), i = s;
+    while (p != 1 && p != n - 1 && a % n && i--)
+      p = bin_mul(p, p, n);
+    if (p != n - 1 && i != s) 
+      return 0;
+  }
+
+  return 1;
+}

algorithms/math/pollard_rho.cpp

+/// Pollard's Rho
+///
+/// Description: 
+///   Returns prime factorization of $n$.
+///
+/// Time: O(n^{1/4})
+/// Space: O(1)
+
+ll pollard(ll n) {
+  auto f = [n](ll x) {
+    return (bin_mul(x, x, n) + 1) % n;
+  };
+
+  if (n % 2 == 0) return 2;
+
+  for (ll i = 2;; ++i) {
+    ll x = i, y = f(x), p;
+    while ((p = __gcd(n + y - x, n)) == 1)
+      x = f(x), y = f(f(y));
+    if (p != n) return p;
+  }
+}
+
+vector<ll> factor(ll n) {
+  if (n == 1) return {};
+  if (is_prime(n)) // Use Miller-Rabin
+    return {n};
+
+  ll x = pollard(n);
+  auto l = factor(x);
+  auto r = factor(n/x);
+
+  l.insert(l.end(), all(r));
+  return l;
+}

contests/Cadernaveis/in

-3
-0 0
-3 0
-0 4
-4
-0 10
-10 0
-20 10
-10 20
-0

contests/Cadernaveis/out

-Case 1: 2.40
-Case 2: 14.15
-Case 3: 0.71
-Case 4: 0.71
-Case 5: 0.71
-Case 6: 0.71
-Case 7: 0.71
-Case 8: 0.71
-Case 9: 0.01
-Case 10: 0.01
-Case 11: 0.01
-Case 12: 0.01
-Case 13: 0.01
-Case 14: 0.01
-Case 15: 10000.00
-Case 16: 26.00
-Case 17: 26.00
-Case 18: 26.00
-Case 19: 26.00
-Case 20: 26.00
-Case 21: 26.00
-Case 22: 2.00
-Case 23: 2.00
-Case 24: 2.00
-Case 25: 2.00
-Case 26: 2.00
-Case 27: 2.00
-Case 28: 26.00
-Case 29: 26.00
-Case 30: 26.00
-Case 31: 26.00
-Case 32: 26.00
-Case 33: 26.00
-Case 34: 421.90
-Case 35: 421.90
-Case 36: 421.90
-Case 37: 421.90
-Case 38: 421.90
-Case 39: 421.90
-Case 40: 313.62
-Case 41: 313.62
-Case 42: 313.62
-Case 43: 313.62
-Case 44: 313.62
-Case 45: 313.62
-Case 46: 739.83
-Case 47: 739.83
-Case 48: 739.83
-Case 49: 739.83
-Case 50: 739.83
-Case 51: 739.83
-Case 52: 939.89
-Case 53: 939.89
-Case 54: 939.89
-Case 55: 939.89
-Case 56: 939.89
-Case 57: 939.89
-Case 58: 943.66
-Case 59: 943.66
-Case 60: 943.66
-Case 61: 943.66
-Case 62: 943.66
-Case 63: 943.66
-Case 64: 965.95
-Case 65: 965.95
-Case 66: 965.95
-Case 67: 965.95
-Case 68: 965.95
-Case 69: 965.95
-Case 70: 941.45
-Case 71: 941.45
-Case 72: 941.45
-Case 73: 941.45
-Case 74: 941.45
-Case 75: 941.45
-Case 76: 827.61
-Case 77: 827.61
-Case 78: 827.61
-Case 79: 827.61
-Case 80: 827.61
-Case 81: 827.61
-Case 82: 962.07
-Case 83: 962.07
-Case 84: 962.07
-Case 85: 962.07
-Case 86: 962.07
-Case 87: 962.07
-Case 88: 522.79
-Case 89: 522.79
-Case 90: 522.79
-Case 91: 522.79
-Case 92: 522.79
-Case 93: 522.79
-Case 94: 868.59
-Case 95: 868.59
-Case 96: 868.59
-Case 97: 868.59
-Case 98: 868.59
-Case 99: 868.59
-Case 100: 943.42
-Case 101: 943.42
-Case 102: 943.42
-Case 103: 943.42
-Case 104: 943.42
-Case 105: 943.42
-Case 106: 849.02
-Case 107: 849.02
-Case 108: 849.02
-Case 109: 849.02
-Case 110: 849.02
-Case 111: 849.02
-Case 112: 872.17
-Case 113: 872.17
-Case 114: 872.17
-Case 115: 872.17
-Case 116: 872.17
-Case 117: 872.17
-Case 118: 918.82
-Case 119: 918.82
-Case 120: 918.82
-Case 121: 918.82
-Case 122: 918.82
-Case 123: 918.82
-Case 124: 830.42
-Case 125: 830.42
-Case 126: 830.42
-Case 127: 830.42
-Case 128: 830.42
-Case 129: 830.42
-Case 130: 887.99
-Case 131: 887.99
-Case 132: 887.99
-Case 133: 887.99
-Case 134: 887.99
-Case 135: 887.99
-Case 136: 819.72
-Case 137: 819.72
-Case 138: 819.72
-Case 139: 819.72
-Case 140: 819.72
-Case 141: 819.72
-Case 142: 916.86
-Case 143: 916.86
-Case 144: 916.86
-Case 145: 916.86
-Case 146: 916.86
-Case 147: 916.86
-Case 148: 894.85
-Case 149: 894.85
-Case 150: 894.85
-Case 151: 894.85
-Case 152: 894.85
-Case 153: 894.85
-Case 154: 907.85
-Case 155: 907.85
-Case 156: 907.85
-Case 157: 907.85
-Case 158: 907.85
-Case 159: 907.85
-Case 160: 897.52
-Case 161: 897.52
-Case 162: 897.52
-Case 163: 897.52
-Case 164: 897.52
-Case 165: 897.52
-Case 166: 698.88
-Case 167: 698.88
-Case 168: 698.88
-Case 169: 698.88
-Case 170: 698.88
-Case 171: 698.88
-Case 172: 520.04
-Case 173: 520.04
-Case 174: 520.04
-Case 175: 520.04
-Case 176: 520.04
-Case 177: 520.04
-Case 178: 948.53
-Case 179: 948.53
-Case 180: 948.53
-Case 181: 948.53
-Case 182: 948.53
-Case 183: 948.53
-Case 184: 890.01
-Case 185: 890.01
-Case 186: 890.01
-Case 187: 890.01
-Case 188: 890.01
-Case 189: 890.01
-Case 190: 863.93
-Case 191: 863.93
-Case 192: 863.93
-Case 193: 863.93
-Case 194: 863.93
-Case 195: 863.93
-Case 196: 902.07
-Case 197: 902.07
-Case 198: 902.07
-Case 199: 902.07
-Case 200: 902.07
-Case 201: 902.07
-Case 202: 841.60
-Case 203: 841.60
-Case 204: 841.60
-Case 205: 841.60
-Case 206: 841.60
-Case 207: 841.60
-Case 208: 887.03
-Case 209: 887.03
-Case 210: 887.03
-Case 211: 887.03
-Case 212: 887.03
-Case 213: 887.03
-Case 214: 677.03
-Case 215: 677.03
-Case 216: 677.03
-Case 217: 677.03
-Case 218: 677.03
-Case 219: 677.03
-Case 220: 917.60
-Case 221: 917.60
-Case 222: 917.60
-Case 223: 917.60
-Case 224: 917.60
-Case 225: 917.60
-Case 226: 699.30
-Case 227: 699.30
-Case 228: 699.30
-Case 229: 699.30
-Case 230: 699.30
-Case 231: 699.30
-Case 232: 781.85
-Case 233: 781.85
-Case 234: 781.85
-Case 235: 781.85
-Case 236: 781.85
-Case 237: 781.85
-Case 238: 635.53
-Case 239: 635.53
-Case 240: 635.53
-Case 241: 635.53
-Case 242: 635.53
-Case 243: 635.53
-Case 244: 846.09
-Case 245: 846.09
-Case 246: 846.09
-Case 247: 846.09
-Case 248: 846.09
-Case 249: 846.09
-Case 250: 553.25
-Case 251: 553.25
-Case 252: 553.25
-Case 253: 553.25
-Case 254: 553.25
-Case 255: 553.25
-Case 256: 971.14
-Case 257: 971.14
-Case 258: 971.14
-Case 259: 971.14
-Case 260: 971.14
-Case 261: 971.14
-Case 262: 834.71
-Case 263: 834.71
-Case 264: 834.71
-Case 265: 834.71
-Case 266: 834.71
-Case 267: 834.71
-Case 268: 857.93
-Case 269: 857.93
-Case 270: 857.93
-Case 271: 857.93
-Case 272: 857.93
-Case 273: 857.93
-Case 274: 929.54
-Case 275: 929.54
-Case 276: 929.54
-Case 277: 929.54
-Case 278: 929.54
-Case 279: 929.54
-Case 280: 882.55
-Case 281: 882.55
-Case 282: 882.55
-Case 283: 882.55
-Case 284: 882.55
-Case 285: 882.55
-Case 286: 920.93
-Case 287: 920.93
-Case 288: 920.93
-Case 289: 920.93
-Case 290: 920.93
-Case 291: 920.93
-Case 292: 885.83
-Case 293: 885.83
-Case 294: 885.83
-Case 295: 885.83
-Case 296: 885.83
-Case 297: 885.83
-Case 298: 969.79
-Case 299: 969.79
-Case 300: 969.79
-Case 301: 969.79
-Case 302: 969.79
-Case 303: 969.79
-Case 304: 935.10
-Case 305: 935.10
-Case 306: 935.10
-Case 307: 935.10
-Case 308: 935.10
-Case 309: 935.10
-Case 310: 821.22
-Case 311: 821.22
-Case 312: 821.22
-Case 313: 821.22
-Case 314: 821.22
-Case 315: 821.22
-Case 316: 936.47
-Case 317: 936.47
-Case 318: 936.47
-Case 319: 936.47
-Case 320: 936.47
-Case 321: 936.47
-Case 322: 669.83
-Case 323: 669.83
-Case 324: 669.83
-Case 325: 669.83
-Case 326: 669.83
-Case 327: 669.83
-Case 328: 937.78
-Case 329: 937.78
-Case 330: 937.78
-Case 331: 937.78
-Case 332: 937.78
-Case 333: 937.78
-Case 334: 911.71
-Case 335: 911.71
-Case 336: 911.71
-Case 337: 911.71
-Case 338: 911.71
-Case 339: 911.71
-Case 340: 878.29
-Case 341: 878.29
-Case 342: 878.29
-Case 343: 878.29
-Case 344: 878.29
-Case 345: 878.29
-Case 346: 911.88
-Case 347: 911.88
-Case 348: 911.88
-Case 349: 911.88
-Case 350: 911.88
-Case 351: 911.88
-Case 352: 819.50
-Case 353: 819.50
-Case 354: 819.50
-Case 355: 819.50
-Case 356: 819.50
-Case 357: 819.50
-Case 358: 914.16
-Case 359: 914.16
-Case 360: 914.16
-Case 361: 914.16
-Case 362: 914.16
-Case 363: 914.16
-Case 364: 525.55
-Case 365: 525.55
-Case 366: 525.55
-Case 367: 525.55
-Case 368: 525.55
-Case 369: 525.55