diff --git a/math/pow_of_matrix/checker.cpp b/math/pow_of_matrix/checker.cpp
new file mode 100644
index 000000000..6a66d5330
--- /dev/null
+++ b/math/pow_of_matrix/checker.cpp
@@ -0,0 +1,62 @@
+// https://github.com/MikeMirzayanov/testlib/blob/master/checkers/wcmp.cpp
+
+// The MIT License (MIT)
+
+// Copyright (c) 2015 Mike Mirzayanov
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#include "testlib.h"
+
+using namespace std;
+
+int main(int argc, char * argv[])
+{
+    setName("compare sequences of tokens");
+    registerTestlibCmd(argc, argv);
+
+    int n = 0;
+    string j, p;
+
+    while (!ans.seekEof() && !ouf.seekEof()) 
+    {
+        n++;
+
+        ans.readWordTo(j);
+        ouf.readWordTo(p);
+        
+        if (j != p)
+            quitf(_wa, "%d%s words differ - expected: '%s', found: '%s'", n, englishEnding(n).c_str(), compress(j).c_str(), compress(p).c_str());
+    }
+
+    if (ans.seekEof() && ouf.seekEof())
+    {
+        if (n == 1)
+            quitf(_ok, "\"%s\"", compress(j).c_str());
+        else
+            quitf(_ok, "%d tokens", n);
+    }
+    else
+    {
+        if (ans.seekEof())
+            quitf(_wa, "Participant output contains extra tokens");
+        else
+            quitf(_wa, "Unexpected EOF in the participants output");
+    }
+}
diff --git a/math/pow_of_matrix/gen/example_00.in b/math/pow_of_matrix/gen/example_00.in
new file mode 100644
index 000000000..8223d21ce
--- /dev/null
+++ b/math/pow_of_matrix/gen/example_00.in
@@ -0,0 +1,3 @@
+2 7
+0 1
+1 1
diff --git a/math/pow_of_matrix/gen/example_01.in b/math/pow_of_matrix/gen/example_01.in
new file mode 100644
index 000000000..28cb328c1
--- /dev/null
+++ b/math/pow_of_matrix/gen/example_01.in
@@ -0,0 +1,4 @@
+3 0
+0 0 0
+0 0 0
+0 0 0
diff --git a/math/pow_of_matrix/gen/example_02.in b/math/pow_of_matrix/gen/example_02.in
new file mode 100644
index 000000000..586769b06
--- /dev/null
+++ b/math/pow_of_matrix/gen/example_02.in
@@ -0,0 +1,4 @@
+3 1000000000000000000
+1 2 3
+4 5 6
+7 8 9
diff --git a/math/pow_of_matrix/gen/lowrank_max_random.cpp b/math/pow_of_matrix/gen/lowrank_max_random.cpp
new file mode 100644
index 000000000..ad5cffc55
--- /dev/null
+++ b/math/pow_of_matrix/gen/lowrank_max_random.cpp
@@ -0,0 +1,59 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform(N_MAX - 10, N_MAX);
+  int rank = gen.uniform(1, n - 1);
+
+  vector<vector<ll>> a(n, vector<ll>(n));
+  for (int i = 0; i < rank; i++) {
+    a[i][i] = gen.uniform(1LL, MOD - 1);
+    for (int j = i + 1; j < n; j++) { a[i][j] = gen.uniform(0LL, MOD - 1); }
+  }
+  for (int ph = 0; ph < 10000; ph++) {
+    int x = gen.uniform(0, n - 1);
+    int y = gen.uniform(0, n - 1);
+    ll freq = gen.uniform(1LL, MOD - 1);
+    if (x == y) continue;
+    if (gen.uniform_bool()) {
+      for (int i = 0; i < n; i++) {
+        a[x][i] += freq * a[y][i];
+        a[x][i] %= MOD;
+      }
+    } else {
+      for (int i = 0; i < n; i++) {
+        a[i][x] += freq * a[i][y];
+        a[i][x] %= MOD;
+      }
+    }
+  }
+
+  gen.shuffle(a.begin(), a.end());
+  for (int i = 0; i < n; i++) {
+    for (int j = i + 1; j < n; j++) { swap(a[i][j], a[j][i]); }
+  }
+  gen.shuffle(a.begin(), a.end());
+
+  ll K = gen.uniform<ll>(0, K_MAX);
+  out(n, K, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/max_random.cpp b/math/pow_of_matrix/gen/max_random.cpp
new file mode 100644
index 000000000..72bb494c8
--- /dev/null
+++ b/math/pow_of_matrix/gen/max_random.cpp
@@ -0,0 +1,33 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = N_MAX;
+  vector<vector<ll>> a(n, vector<ll>(n));
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) { a[i][j] = gen.uniform<int>(0, MOD - 1); }
+  }
+
+  ll K = gen.uniform<ll>(0, K_MAX);
+  out(n, K, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/max_random_worst.cpp b/math/pow_of_matrix/gen/max_random_worst.cpp
new file mode 100644
index 000000000..57d07de93
--- /dev/null
+++ b/math/pow_of_matrix/gen/max_random_worst.cpp
@@ -0,0 +1,33 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = N_MAX;
+  vector<vector<ll>> a(n, vector<ll>(n));
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) { a[i][j] = gen.uniform<int>(0, MOD - 1); }
+  }
+
+  ll K = (1LL << 59) - 1;
+  out(n, K, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/perm_max_random.cpp b/math/pow_of_matrix/gen/perm_max_random.cpp
new file mode 100644
index 000000000..281702ad5
--- /dev/null
+++ b/math/pow_of_matrix/gen/perm_max_random.cpp
@@ -0,0 +1,31 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform(N_MAX - 10, N_MAX);
+  vector<vector<ll>> a(n, vector<ll>(n));
+  vector<int> p = gen.perm<int>(n);
+  for (int i = 0; i < n; i++) { a[i][p[i]] = gen.uniform<int>(0, MOD - 1); }
+
+  ll K = gen.uniform<ll>(0, K_MAX);
+  out(n, K, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/random.cpp b/math/pow_of_matrix/gen/random.cpp
new file mode 100644
index 000000000..04ddcb86a
--- /dev/null
+++ b/math/pow_of_matrix/gen/random.cpp
@@ -0,0 +1,33 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int n = gen.uniform<int>(1, N_MAX);
+  vector<vector<ll>> a(n, vector<ll>(n));
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) { a[i][j] = gen.uniform<int>(0, MOD - 1); }
+  }
+
+  ll K = gen.uniform<ll>(0, K_MAX);
+  out(n, K, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/signed_overflow.cpp b/math/pow_of_matrix/gen/signed_overflow.cpp
new file mode 100644
index 000000000..122b6fbf0
--- /dev/null
+++ b/math/pow_of_matrix/gen/signed_overflow.cpp
@@ -0,0 +1,35 @@
+#include <iostream>
+#include <set>
+#include <limits>
+
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+  int val = MOD / 2 - 1;
+  unsigned long long max_ll = numeric_limits<long long>::max();
+
+  // Find smallest size n where n * val * val > max_ll
+  int n = 0;
+  for (unsigned long long res = 0; res <= max_ll;
+       res += (long long)val * val, ++n)
+    ;
+  ll K = gen.uniform<ll>(0, K_MAX);
+
+  printf("%d %lld\n", n, K);
+
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (j) printf(" ");
+      printf("%d", val);
+    }
+    printf("\n");
+  }
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/small.cpp b/math/pow_of_matrix/gen/small.cpp
new file mode 100644
index 000000000..1e1cfa6bc
--- /dev/null
+++ b/math/pow_of_matrix/gen/small.cpp
@@ -0,0 +1,37 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+void out(int n, ll K, vector<vector<ll>> a) {
+  printf("%d %lld\n", n, K);
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      printf("%lld", a[i][j]);
+      if (j + 1 != n) printf(" ");
+    }
+    printf("\n");
+  }
+}
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  vector<int> ns = {1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3};
+  vector<int> ks = {0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3};
+
+  int n = ns[seed % 12];
+  int k = ks[seed % 12];
+
+  vector<vector<ll>> a(n, vector<ll>(n));
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) { a[i][j] = gen.uniform<int>(0, MOD - 1); }
+  }
+
+  out(n, k, a);
+  return 0;
+}
diff --git a/math/pow_of_matrix/gen/unsigned_overflow.cpp b/math/pow_of_matrix/gen/unsigned_overflow.cpp
new file mode 100644
index 000000000..67a6113cf
--- /dev/null
+++ b/math/pow_of_matrix/gen/unsigned_overflow.cpp
@@ -0,0 +1,39 @@
+#include <iostream>
+#include <set>
+#include <limits>
+
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+  int val = MOD - 2;
+
+  // Find smallest size n where n * val * val > numeric_limits<unsigned long
+  // long>::max()
+  int n = 0;
+  unsigned long long res = 0;
+  long long val2 = (long long)val * val;
+  while (res + val2 >= res) // Check if adding val2 to res will overflow
+  {
+    res += val2;
+    ++n;
+  }
+
+  ll K = gen.uniform<ll>(0, K_MAX);
+  printf("%d %lld\n", n, K);
+
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (j) printf(" ");
+      printf("%d", val);
+    }
+    printf("\n");
+  }
+  return 0;
+}
diff --git a/math/pow_of_matrix/hash.json b/math/pow_of_matrix/hash.json
new file mode 100644
index 000000000..c2d6bce4c
--- /dev/null
+++ b/math/pow_of_matrix/hash.json
@@ -0,0 +1,66 @@
+{
+  "example_00.in": "8c96d051dc1252449299b5cfb3b06c466f5cb0194f97172b33f5078642c1fb61",
+  "example_00.out": "b014e72211237f7f5001f372abe2082ea85fa3c7bb6ddda3778000eb5739a381",
+  "example_01.in": "81f6b4247d80a09621f0747ab246f4e707c0fc3926b0c1f522ac493e5b5eae97",
+  "example_01.out": "bf284252c00be268976cb5f85fe21220797083556c200aa90e9f159123e71873",
+  "example_02.in": "679e0c7dd8b57266890588b436e127e6bed407dd7cc92d7eac0b9f0914f3ef1e",
+  "example_02.out": "2f8dbf0dd57f04a042edb45647bbe4d0d50862d5086dcc4dacefacf8ea8e7a00",
+  "lowrank_max_random_00.in": "ea3ae3425946931d851a1db4d661468b6fb08fd92b8efe0ac01d60d60a1886ab",
+  "lowrank_max_random_00.out": "06f2597a17232e81af8683f6cd9989112e0e4a23fb82143f13ec4dfc753abd6e",
+  "lowrank_max_random_01.in": "5e3518e4a510bdb890e411179fa01f8bd41d418a929175a5609ce32c142bf878",
+  "lowrank_max_random_01.out": "82ba96ded3c793f7a454675dd7816298aa5c4657d5d67e85fa18b9e0613c3b36",
+  "max_random_00.in": "54f7d64aa5c5ee891ccaade66a0626be24a40017e7454be7e1b1e12524d3ce36",
+  "max_random_00.out": "d0e392f85f64d16be4d29d2b9cdb5e8ff65bf1e5bdc248c8b986989c2dbe8fe7",
+  "max_random_01.in": "077a541e659a94478f56fd388874a516956c31d22464bd24fd63085b15177c86",
+  "max_random_01.out": "ec7ed797d3a222d08345a18ddf4ba0301e8222c49aff33db5e67bd18616300d3",
+  "max_random_02.in": "bb7f0569218717fcf5dbc09ed54d599afbb66cc5ec38a8b4f7decd90b7f60820",
+  "max_random_02.out": "d8a902b636c27d77e3e43cb99d29a542012d5718de2b5d1b24a4245de4cc99b5",
+  "max_random_03.in": "9c48eb180419c4cadb21882448d7cebc789b1c2d86e27d0f026be1fe7abd2a8e",
+  "max_random_03.out": "0a85424e818527de34b61a68d92b106d0baeee2923f8a5e45e748ff4f41f3bf5",
+  "max_random_worst_00.in": "8bab459f1937906aa7d1ab4731b8ceb8fd24cbb13398785d6b64c57ca3f96c64",
+  "max_random_worst_00.out": "5b425c3c5d12e9c120f29c0f76c1ad582a9dad09a6d9099f2a6455195e8e5cd2",
+  "max_random_worst_01.in": "d99d119f3a52ca6066a99d42ba363e184d10b0d73cf6a1a752b44fa171b446fa",
+  "max_random_worst_01.out": "cb2c717e217aa1057dd92c047959ec9c6c42688a2f96f8c27f2a65eb35a8f6ca",
+  "max_random_worst_02.in": "49d58655ab64c821bfcec7450aa1ebecd6d9da5baac7b14dc29e948aa60ef858",
+  "max_random_worst_02.out": "76af4dc33e2d5c3feca7675bb8ded9affc5611d8775250fcac64aecb5912f82d",
+  "max_random_worst_03.in": "6ce66c9e829521fa1d416295348188d90e91e7d09b17d8a071565faf1a4618ab",
+  "max_random_worst_03.out": "0dfd9d4f2a6bdce3d9a8d45fe41cc60d560f7fd9f2150e0872a2810a1aa4861c",
+  "perm_max_random_00.in": "1a4c32f645aa0b2243ddff51630762e9ec352c58f40b2a3df1a99d9d8c176aa0",
+  "perm_max_random_00.out": "21a6b93cfa8aa616d60cd65ccbfecb7c98dc11dad8f23634a80055661608f3ab",
+  "perm_max_random_01.in": "866a0a66eb3d683e7150ccacc7c1d4885b555773459b4d45678654190c68414a",
+  "perm_max_random_01.out": "8e736a59917e0cf4b7744cef2d05ba6807c5b81ce9728da66ffe3e0e049c0652",
+  "random_00.in": "d9edb3b8efa26b8c40c708f37813a2ff44fffdc56a02af4ef0594e01dbc1b064",
+  "random_00.out": "af7a1ed00b25420cd595b792e945f79b4d149ae7fac87fbff3b4916f6cbbbd66",
+  "random_01.in": "18ff31aba7fd51754c0a37dad9260ccbfd9154603e948cd15afe9fee0426bfe9",
+  "random_01.out": "b33e04f85d2776d80130df76a1904921c3a714d38b46f47d09f78243ee2efa7c",
+  "random_02.in": "1bd68099f35fbde28c7a157c23be1743385903faaf927578b21320bbf6434e1b",
+  "random_02.out": "a6c239f0ec88f89641f4097afb7b83e365e6a76f5c60ac5142df8ddb90593fc5",
+  "signed_overflow_00.in": "cf53f3409bc01d23426bd52be940d4f5d02c1c30edad154a78e7c941cb7d347d",
+  "signed_overflow_00.out": "6811813e34ad291d293d8ec1807056845cc5f896ec3dd96877cfd74e02bd7d53",
+  "small_00.in": "3c45b1eda0cb7bcadd7f25fbcccb38d63b60f0a6aee84b54f8bce38d1a6e8c07",
+  "small_00.out": "4355a46b19d348dc2f57c046f8ef63d4538ebb936000f3c9ee954a27460dd865",
+  "small_01.in": "30844f59c667ed9b2501ec8e962f73f5ff242e931eb0ae9e94500ba395396dea",
+  "small_01.out": "70b2e20a859896a187c305eb3029848afd081fda4405f3ceb20285893410b217",
+  "small_02.in": "889ac6326529e21ac65f92e51846473fbe52ae6cb10a4644e55ced3015130db9",
+  "small_02.out": "bf284252c00be268976cb5f85fe21220797083556c200aa90e9f159123e71873",
+  "small_03.in": "ef7e985afdac0c67aa8a69681da734108b76f7c28e4dd59e40674593478677ea",
+  "small_03.out": "af3f76179322d5f9409d5fc655ae7f5e6ec6a572a0df43e580f552cbea13a9f2",
+  "small_04.in": "7c97dd267f6eb89f1b6313a78fc9998b3bfb475e45947479c178e0bb7c9e81d0",
+  "small_04.out": "cb76b5777ca834f966699ad212c68c3a056f7e749753736d80634f733ce1f7de",
+  "small_05.in": "de0a02afd7d1c85e76e5acc85b6be04c594e454a1b7a1e38ef49670490628b45",
+  "small_05.out": "4c9c3181f5a709a3af94947957f802409023802aa93c09698b4ac527c252c1dc",
+  "small_06.in": "5b60e41f367aea1ce1b1e74a9207470ac98dfd3524416fb2e6b6510c359a82bd",
+  "small_06.out": "b478b81b481e9c31746e6bb8761fa88553aa2fab0a8daa52dd51794699af5442",
+  "small_07.in": "e3031bda06a6d9d0c23fa208c25935626561a411b3616d0ea3a765285b58fcac",
+  "small_07.out": "2144ca19e9da783620e89cf7c0fa1c7b1aeb796b14ca9cf6c7025437a3920575",
+  "small_08.in": "78a81435da48a049189333778318a0fc4c6893364948133e3dd8545dec1c53f5",
+  "small_08.out": "0cf040b24cd023ef08ff924761ffcc76cb41b940d0390351d97c2aaee346d853",
+  "small_09.in": "1b6b05ebb30916894e7e44b8e0b729bb418870fabc254a120fa3f45b94ed2262",
+  "small_09.out": "f13b831ceaa02fc1eb3c79061ac9b34de17677df88448cc9d8994d6b28dd60c5",
+  "small_10.in": "44026854bfdb4ae1ca269bc6f3dcdfe728a64dbd1e2852f0e354d9f514988615",
+  "small_10.out": "50c0344d39ce80562bd1abef1a234745245de007af37752a72f526772f88e353",
+  "small_11.in": "fcd810efd07cd65684a9175aa2c160f02827a6f20e6a5a185e87d336bfdd6ab1",
+  "small_11.out": "6ed9021c9cc90afc948cfc9d8e1dd3b56f38bfb7d7566e6876a4aa7aa8f4f804",
+  "unsigned_overflow_00.in": "52d993215eb7fa8f0a370bd8cf58cfe9d49b99d9a12bdc59db247e6e74631e3b",
+  "unsigned_overflow_00.out": "7eadc006d6ec87ec0908f0f1a9093fdeddd57514ac0315c4c8655d4c215b5be5"
+}
\ No newline at end of file
diff --git a/math/pow_of_matrix/info.toml b/math/pow_of_matrix/info.toml
new file mode 100644
index 000000000..fc950b65c
--- /dev/null
+++ b/math/pow_of_matrix/info.toml
@@ -0,0 +1,36 @@
+title = 'Pow of Matrix'
+timelimit = 5.0
+forum = "https://github.com/yosupo06/library-checker-problems/issues/567"
+
+[[tests]]
+    name = "example.in"
+    number = 3
+[[tests]]
+    name = "random.cpp"
+   number = 3
+[[tests]]
+    name = "small.cpp"
+    number = 12
+[[tests]]
+    name = "max_random.cpp"
+    number = 4
+[[tests]]
+    name = "max_random_worst.cpp"
+    number = 4
+[[tests]]
+    name = "lowrank_max_random.cpp"
+    number = 2
+[[tests]]
+    name = "perm_max_random.cpp"
+    number = 2
+[[tests]]
+    name = "signed_overflow.cpp"
+    number = 1
+[[tests]]
+    name = "unsigned_overflow.cpp"
+    number = 1
+
+[params]
+    N_MAX = 200
+    K_MAX = 1_000_000_000_000_000_000
+    MOD = 998244353
diff --git a/math/pow_of_matrix/sol/correct.cpp b/math/pow_of_matrix/sol/correct.cpp
new file mode 100644
index 000000000..d03fe707d
--- /dev/null
+++ b/math/pow_of_matrix/sol/correct.cpp
@@ -0,0 +1,172 @@
+#include <cstdio>
+#include <cassert>
+#include <vector>
+
+using namespace std;
+
+using ll = long long;
+using u32 = unsigned int;
+using u64 = unsigned long long;
+using u128 = unsigned __int128;
+
+using vi = vector<ll>;
+template <class T>
+using vc = vector<T>;
+template <class T>
+using vvc = vector<vc<T>>;
+
+#define vv(type, name, h, ...) \
+  vector<vector<type>> name(h, vector<type>(__VA_ARGS__))
+
+#define FOR1(a) for (ll _ = 0; _ < ll(a); ++_)
+#define FOR2(i, a) for (ll i = 0; i < ll(a); ++i)
+#define FOR3(i, a, b) for (ll i = a; i < ll(b); ++i)
+#define FOR4(i, a, b, c) for (ll i = a; i < ll(b); i += (c))
+#define FOR1_R(a) for (ll i = (a)-1; i >= ll(0); --i)
+#define FOR2_R(i, a) for (ll i = (a)-1; i >= ll(0); --i)
+#define FOR3_R(i, a, b) for (ll i = (b)-1; i >= ll(a); --i)
+#define overload4(a, b, c, d, e, ...) e
+#define overload3(a, b, c, d, ...) d
+#define FOR(...) overload4(__VA_ARGS__, FOR4, FOR3, FOR2, FOR1)(__VA_ARGS__)
+#define FOR_R(...) overload3(__VA_ARGS__, FOR3_R, FOR2_R, FOR1_R)(__VA_ARGS__)
+
+#define all(x) x.begin(), x.end()
+#define len(x) ll(x.size())
+#define elif else if
+
+#define eb emplace_back
+#define mp make_pair
+#define mt make_tuple
+#define fi first
+#define se second
+
+template <int mod>
+struct modint {
+  static constexpr u32 umod = u32(mod);
+  static_assert(umod < u32(1) << 31);
+  u32 val;
+
+  static modint raw(u32 v) {
+    modint x;
+    x.val = v;
+    return x;
+  }
+  constexpr modint() : val(0) {}
+  constexpr modint(u32 x) : val(x % umod) {}
+  constexpr modint(u64 x) : val(x % umod) {}
+  constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){};
+  constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){};
+  bool operator<(const modint &other) const { return val < other.val; }
+  modint &operator+=(const modint &p) {
+    if ((val += p.val) >= umod) val -= umod;
+    return *this;
+  }
+  modint &operator-=(const modint &p) {
+    if ((val += umod - p.val) >= umod) val -= umod;
+    return *this;
+  }
+  modint &operator*=(const modint &p) {
+    val = u64(val) * p.val % umod;
+    return *this;
+  }
+  modint &operator/=(const modint &p) {
+    *this *= p.inverse();
+    return *this;
+  }
+  modint operator-() const { return modint::raw(val ? mod - val : u32(0)); }
+  modint operator+(const modint &p) const { return modint(*this) += p; }
+  modint operator-(const modint &p) const { return modint(*this) -= p; }
+  modint operator*(const modint &p) const { return modint(*this) *= p; }
+  modint operator/(const modint &p) const { return modint(*this) /= p; }
+  bool operator==(const modint &p) const { return val == p.val; }
+  bool operator!=(const modint &p) const { return val != p.val; }
+  modint inverse() const {
+    int a = val, b = mod, u = 1, v = 0, t;
+    while (b > 0) {
+      t = a / b;
+      swap(a -= t * b, b), swap(u -= t * v, v);
+    }
+    return modint(u);
+  }
+  modint pow(ll n) const {
+    assert(n >= 0);
+    modint ret(1), mul(val);
+    while (n > 0) {
+      if (n & 1) ret *= mul;
+      mul *= mul;
+      n >>= 1;
+    }
+    return ret;
+  }
+  static constexpr int get_mod() { return mod; }
+  static constexpr pair<int, int> ntt_info() {
+    if (mod == 167772161) return {25, 17};
+    if (mod == 469762049) return {26, 30};
+    if (mod == 754974721) return {24, 362};
+    if (mod == 880803841) return {23, 211};
+    if (mod == 943718401) return {22, 663003469};
+    if (mod == 998244353) return {23, 31};
+    if (mod == 1045430273) return {20, 363};
+    if (mod == 1051721729) return {20, 330};
+    if (mod == 1053818881) return {20, 2789};
+    return {-1, -1};
+  }
+  static constexpr bool can_ntt() { return ntt_info().fi != -1; }
+};
+
+using modint998 = modint<998244353>;
+
+template <typename T>
+vc<vc<T>> matrix_mul(const vc<vc<T>> &A, const vc<vc<T>> &B, int N1 = -1,
+                     int N2 = -1, int N3 = -1) {
+  if (N1 == -1) { N1 = len(A), N2 = len(B), N3 = len(B[0]); }
+  vv(u32, b, N3, N2);
+  FOR(i, N2) FOR(j, N3) b[j][i] = B[i][j].val;
+  vv(T, C, N1, N3);
+
+  FOR(i, N1) FOR(j, N3) {
+    u128 sm = 0;
+    FOR(m, N2) sm += u64(A[i][m].val) * b[j][m];
+    C[i][j] = T::raw(sm % (T::get_mod()));
+  }
+  return C;
+}
+
+template <typename T>
+vc<vc<T>> matrix_pow(vc<vc<T>> A, ll n) {
+  int N = len(A);
+  vv(T, ret, N, N);
+  FOR(i, N) ret[i][i] = T(1);
+  while (n) {
+    if (n & 1) ret = matrix_mul(ret, A, N, N, N);
+    n /= 2;
+    if (n) A = matrix_mul(A, A, N, N, N);
+  }
+  return ret;
+}
+
+using mint = modint998;
+
+void solve() {
+  ll N, K;
+  scanf("%lld %lld", &N, &K);
+  vv(mint, A, N, N);
+  FOR(i, N) FOR(j, N) {
+    int x;
+    scanf("%d", &x);
+    A[i][j] = x;
+  }
+  A = matrix_pow(A, K);
+  FOR(i, N) {
+    FOR(j, N) {
+      if (j) printf(" ");
+      printf("%d", A[i][j].val);
+    }
+    printf("\n");
+  }
+}
+
+signed main() {
+  solve();
+  return 0;
+}
diff --git a/math/pow_of_matrix/task.md b/math/pow_of_matrix/task.md
new file mode 100644
index 000000000..fe96f35cb
--- /dev/null
+++ b/math/pow_of_matrix/task.md
@@ -0,0 +1,39 @@
+## @{keyword.statement}
+
+@{lang.en}
+Given $N \times N$ matrix $A = \lbrace a_{ij} \rbrace$ with entries in $\mathbb{Z}/@{param.MOD}\mathbb{Z}$, and a non-negative integer $K$. 
+print $A^K = (b_{ij})_{i,j}$. 
+@{lang.ja}
+$\mathbb{Z}/@{param.MOD}\mathbb{Z}$ 成分の $N \times N$ 正方行列 $A=a_{ij}$ および非負整数 $K$ が与えられます．$A^K = (b_{ij})_{i,j}$ を求めてください．
+@{lang.end}
+
+## @{keyword.constraints}
+
+- $1 \leq N \leq @{param.N_MAX}$
+- $0 \leq K \leq @{param.K_MAX}$
+- $0 \leq a_{ij} < @{param.MOD}$
+
+## @{keyword.input}
+
+```
+$N$ $K$
+$a_{11}$ $a_{12}$ ... $a_{1N}$
+$a_{21}$ $a_{22}$ ... $a_{2N}$
+:
+$a_{N1}$ $a_{N2}$ ... $a_{NN}$
+```
+
+## @{keyword.output}
+
+```
+$b_{11}$ $b_{12}$ ... $b_{1N}$
+$b_{21}$ $b_{22}$ ... $b_{2N}$
+:
+$b_{N1}$ $b_{N2}$ ... $b_{NN}$
+```
+
+## @{keyword.sample}
+
+@{example.example_00}
+
+@{example.example_01}
diff --git a/math/pow_of_matrix/verifier.cpp b/math/pow_of_matrix/verifier.cpp
new file mode 100644
index 000000000..4b5aabe38
--- /dev/null
+++ b/math/pow_of_matrix/verifier.cpp
@@ -0,0 +1,22 @@
+#include <iostream>
+
+#include "testlib.h"
+#include "params.h"
+
+int main() {
+  registerValidation();
+
+  int n = inf.readInt(1, N_MAX);
+  inf.readSpace();
+  inf.readLong(0LL, K_MAX);
+  inf.readChar('\n');
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      inf.readInt(0, MOD - 1);
+      if (j + 1 != n) inf.readSpace();
+    }
+    inf.readChar('\n');
+  }
+  inf.readEof();
+  return 0;
+}