diff --git a/number_theory/two_square_sum/checker.cpp b/number_theory/two_square_sum/checker.cpp
new file mode 100644
index 000000000..5a2089f45
--- /dev/null
+++ b/number_theory/two_square_sum/checker.cpp
@@ -0,0 +1,40 @@
+#include "testlib.h"
+
+#include <vector>
+#include <tuple>
+#include <cassert>
+
+using namespace std;
+using ll = long long;
+
+ll read_ans_cnt(ll N, InStream& stream) {
+  ll n = stream.readInt();
+  vector<pair<ll, ll>> ANS(n);
+  for (int i = 0; i < n; ++i) {
+    ll x = stream.readLong(0LL, 1LL << 30);
+    ll y = stream.readLong(0LL, 1LL << 30);
+    stream.ensure(x * x + y * y == N);
+    ANS[i] = {x, y};
+  }
+  // check distinct
+  sort(ANS.begin(), ANS.end());
+  for (int i = 0; i < n - 1; ++i) stream.ensure(ANS[i] != ANS[i + 1]);
+  return n;
+}
+
+int main(int argc, char* argv[]) {
+  registerTestlibCmd(argc, argv);
+
+  int Q = inf.readInt();
+  for (int q = 0; q < Q; ++q) {
+    ll N = inf.readLong();
+    ll ans_cnt = read_ans_cnt(N, ans);
+    ll ouf_cnt = read_ans_cnt(N, ouf);
+    if (ans_cnt > ouf_cnt) {
+      quitf(_wa, "Not found all solutions. N=%lld, ans_cnt=%lld, ouf_cnt=%lld", N, ans_cnt, ouf_cnt);
+    } else if (ans_cnt < ouf_cnt) {
+      quitf(_fail, "Writer does not find all solutions! N=%lld, ans_cnt=%lld, ouf_cnt=%lld", N, ans_cnt, ouf_cnt);
+    }
+  }
+  quitf(_ok, "OK");
+}
diff --git a/number_theory/two_square_sum/gen/big2.cpp b/number_theory/two_square_sum/gen/big2.cpp
new file mode 100644
index 000000000..c47ef7883
--- /dev/null
+++ b/number_theory/two_square_sum/gen/big2.cpp
@@ -0,0 +1,41 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+vector<long long> enum_prime(long long st, long long ed) {
+    if (st == 1) st = 2;
+    vector<bool> is_prime(ed - st + 1, true);
+    for (long long i = 2; i * i <= ed; i++) {
+        for (long long j = (st + i - 1) / i * i; j <= ed; j += i) {
+            is_prime[j - st] = false;
+        }
+    }
+    vector<long long> primes;
+    for (long long i = st; i <= ed; i++) {
+        if (is_prime[i - st]) primes.push_back(i);
+    }
+    return primes;
+}
+
+int main(int, char* argv[]) {
+
+
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    auto primes = enum_prime(1'000'000'000 - 1'000'000, 1'000'000'000);
+    int k = int(primes.size());
+    int q = MAX_Q;
+    vector<long long> a(q);
+    for (int i = 0; i < q; i++) {
+        a[i] = primes[gen.uniform(0, k - 1)] * primes[gen.uniform(0, k - 1)];
+    }
+
+    printf("%d\n", q);
+    for (auto x: a) {
+        printf("%lld\n", x);
+    }
+    return 0;
+}
diff --git a/number_theory/two_square_sum/gen/big_semiprime_random.cpp b/number_theory/two_square_sum/gen/big_semiprime_random.cpp
new file mode 100644
index 000000000..32f88ead1
--- /dev/null
+++ b/number_theory/two_square_sum/gen/big_semiprime_random.cpp
@@ -0,0 +1,116 @@
+#include <iostream>
+#include <cmath>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using ll = long long;
+using ull = unsigned long long;
+template <class T> using V = vector<T>;
+
+// bit op
+int bsf(ull x) { return __builtin_ctzll(x); }
+
+// binary gcd
+ll gcd(ll _a, ll _b) {
+    ull a = abs(_a), b = abs(_b);
+    if (a == 0) return b;
+    if (b == 0) return a;
+    int shift = bsf(a | b);
+    a >>= bsf(a);
+    do {
+        b >>= bsf(b);
+        if (a > b) swap(a, b);
+        b -= a;
+    } while (b);
+    return (a << shift);
+}
+
+template <class T, class U> T pow_mod(T x, U n, T md) {
+    T r = 1 % md;
+    x %= md;
+    while (n) {
+        if (n & 1) r = (r * x) % md;
+        x = (x * x) % md;
+        n >>= 1;
+    }
+    return r;
+}
+
+bool is_prime(ll n) {
+    if (n <= 1) return false;
+    if (n == 2) return true;
+    if (n % 2 == 0) return false;
+    ll d = n - 1;
+    while (d % 2 == 0) d /= 2;
+    for (ll a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
+        if (n <= a) break;
+        ll t = d;
+        ll y = pow_mod<__int128_t>(a, t, n);  // over
+        while (t != n - 1 && y != 1 && y != n - 1) {
+            y = __int128_t(y) * y % n;  // flow
+            t <<= 1;
+        }
+        if (y != n - 1 && t % 2 == 0) {
+            return false;
+        }
+    }
+    return true;
+}
+
+ll pollard_single(ll n) {
+    if (is_prime(n)) return n;
+    if (n % 2 == 0) return 2;
+    ll st = 0;
+    auto f = [&](ll x) { return (__int128_t(x) * x + st) % n; };
+    while (true) {
+        st++;
+        ll x = st, y = f(x);
+        while (true) {
+            ll p = gcd((y - x + n), n);
+            if (p == 0 || p == n) break;
+            if (p != 1) return p;
+            x = f(x);
+            y = f(f(y));
+        }
+    }
+}
+
+V<ll> pollard(ll n) {
+    if (n == 1) return {};
+    ll x = pollard_single(n);
+    if (x == n) return {x};
+    V<ll> le = pollard(x);
+    V<ll> ri = pollard(n / x);
+    le.insert(le.end(), ri.begin(), ri.end());
+    return le;
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int q = MAX_Q;
+    V<ll> a(q);
+    for (int i = 0; i < q; ++i) {
+        while(true) {
+            ll x = gen.uniform(MAX_A - ll(sqrt(MAX_A)), MAX_A);
+            auto v = pollard(x);
+            if (v.size() == 2) {
+                ll p = max(v[0], v[1]);
+                ll q = min(v[0], v[1]);
+                if (p / q < 100) {
+                    a[i] = x;
+                    break;
+                }
+            }
+        }
+    }
+
+    printf("%d\n", q);
+    for (auto x: a) {
+        printf("%lld\n", x);
+    }
+    return 0;
+}
diff --git a/number_theory/two_square_sum/gen/example_00.in b/number_theory/two_square_sum/gen/example_00.in
new file mode 100644
index 000000000..96a552b52
--- /dev/null
+++ b/number_theory/two_square_sum/gen/example_00.in
@@ -0,0 +1,12 @@
+11
+0
+1
+2
+3
+4
+5
+6
+7
+8
+9
+10
diff --git a/number_theory/two_square_sum/gen/example_01.in b/number_theory/two_square_sum/gen/example_01.in
new file mode 100644
index 000000000..6525a53ff
--- /dev/null
+++ b/number_theory/two_square_sum/gen/example_01.in
@@ -0,0 +1,11 @@
+10
+5
+25
+125
+13
+65
+998244353
+1000000007
+1000000009
+536870912
+1073741824
diff --git a/number_theory/two_square_sum/gen/fixed_RNG_buster_00.in b/number_theory/two_square_sum/gen/fixed_RNG_buster_00.in
new file mode 100644
index 000000000..383392389
--- /dev/null
+++ b/number_theory/two_square_sum/gen/fixed_RNG_buster_00.in
@@ -0,0 +1,2 @@
+1
+124376107291
diff --git a/number_theory/two_square_sum/gen/max.cpp b/number_theory/two_square_sum/gen/max.cpp
new file mode 100644
index 000000000..54266058a
--- /dev/null
+++ b/number_theory/two_square_sum/gen/max.cpp
@@ -0,0 +1,13 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "../params.h"
+
+int main(int, char **argv) {
+  const long long seed = atoll(argv[1]);
+  printf("%lld\n", MAX_Q);
+  for (int i = 0; i < MAX_Q; ++i) {
+    const long long a = MAX_A - seed * MAX_Q - i;
+    printf("%lld\n", a);
+  }
+  return 0;
+}
diff --git a/number_theory/two_square_sum/gen/num_of_solution_max_00.in b/number_theory/two_square_sum/gen/num_of_solution_max_00.in
new file mode 100644
index 000000000..ca7bfcbb5
--- /dev/null
+++ b/number_theory/two_square_sum/gen/num_of_solution_max_00.in
@@ -0,0 +1,101 @@
+100
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
+331854559557386125
+663709119114772250
diff --git a/number_theory/two_square_sum/gen/num_of_solution_max_01.in b/number_theory/two_square_sum/gen/num_of_solution_max_01.in
new file mode 100644
index 000000000..8e03b57f7
--- /dev/null
+++ b/number_theory/two_square_sum/gen/num_of_solution_max_01.in
@@ -0,0 +1,101 @@
+100
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
+663709119114772250
diff --git a/number_theory/two_square_sum/gen/power_of_2_3_5_13.cpp b/number_theory/two_square_sum/gen/power_of_2_3_5_13.cpp
new file mode 100644
index 000000000..e0db92fe9
--- /dev/null
+++ b/number_theory/two_square_sum/gen/power_of_2_3_5_13.cpp
@@ -0,0 +1,31 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <vector>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+using ll = long long;
+
+int main(int, char **argv) {
+  const long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int primes[] = {2, 3, 5, 13};
+  ll p = primes[seed % 4];
+
+  vector<ll> A = {1};
+  while (1) {
+    ll a = A.back();
+    if (a > MAX_A / p) break;
+    A.emplace_back(a * p);
+  }
+
+  printf("%lld\n", MAX_Q);
+  for (int i = 0; i < MAX_Q; ++i) {
+    ll k = gen.uniform<int>(0, A.size() - 1);
+    ll a = A[k];
+    printf("%lld\n", a);
+  }
+  return 0;
+}
diff --git a/number_theory/two_square_sum/gen/prime_1_mod_4.cpp b/number_theory/two_square_sum/gen/prime_1_mod_4.cpp
new file mode 100644
index 000000000..29a54524e
--- /dev/null
+++ b/number_theory/two_square_sum/gen/prime_1_mod_4.cpp
@@ -0,0 +1,25 @@
+#include <iostream>
+#include <cinttypes>
+#include "random.h"
+#include "../params.h"
+#include "../lib/prime.hpp"
+
+using namespace std;
+using ll = long long;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+  PrimeGenerator pgen;
+
+  auto F = [&]() -> ll {
+    while (1) {
+      ll p = pgen.gen(MAX_A, gen);
+      if (p % 4 == 1) return p;
+    }
+  };
+
+  int Q = MAX_Q;
+  printf("%d\n", Q);
+  for (int i = 0; i < Q; i++) { printf("%lld\n", F()); }
+}
diff --git a/number_theory/two_square_sum/gen/prime_3_mod_4.cpp b/number_theory/two_square_sum/gen/prime_3_mod_4.cpp
new file mode 100644
index 000000000..3562826e6
--- /dev/null
+++ b/number_theory/two_square_sum/gen/prime_3_mod_4.cpp
@@ -0,0 +1,25 @@
+#include <iostream>
+#include <cinttypes>
+#include "random.h"
+#include "../params.h"
+#include "../lib/prime.hpp"
+
+using namespace std;
+using ll = long long;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+  PrimeGenerator pgen;
+
+  auto F = [&]() -> ll {
+    while (1) {
+      ll p = pgen.gen(MAX_A, gen);
+      if (p % 4 == 3) return p;
+    }
+  };
+
+  int Q = MAX_Q;
+  printf("%d\n", Q);
+  for (int i = 0; i < Q; i++) { printf("%lld\n", F()); }
+}
diff --git a/number_theory/two_square_sum/gen/random.cpp b/number_theory/two_square_sum/gen/random.cpp
new file mode 100644
index 000000000..4f1e53881
--- /dev/null
+++ b/number_theory/two_square_sum/gen/random.cpp
@@ -0,0 +1,18 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+  long long seed = atoll(argv[1]);
+  auto gen = Random(seed);
+
+  int q = gen.uniform(1LL, MAX_Q);
+  vector<long long> a(q);
+  for (int i = 0; i < q; i++) { a[i] = gen.uniform(1LL, MAX_A); }
+
+  printf("%d\n", q);
+  for (auto x: a) { printf("%lld\n", x); }
+  return 0;
+}
diff --git a/number_theory/two_square_sum/gen/small.cpp b/number_theory/two_square_sum/gen/small.cpp
new file mode 100644
index 000000000..88a817329
--- /dev/null
+++ b/number_theory/two_square_sum/gen/small.cpp
@@ -0,0 +1,23 @@
+#include <iostream>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+int main(int, char* argv[]) {
+
+
+    long long seed = atoll(argv[1]);
+
+    int q = MAX_Q;
+    vector<long long> a(q);
+    for (int i = 0; i < q; i++) {
+        a[i] = seed * 100 + 1 + i;
+    }
+
+    printf("%d\n", q);
+    for (auto x: a) {
+        printf("%lld\n", x);
+    }
+    return 0;
+}
diff --git a/number_theory/two_square_sum/gen/smooth.cpp b/number_theory/two_square_sum/gen/smooth.cpp
new file mode 100644
index 000000000..7c8fcfdcc
--- /dev/null
+++ b/number_theory/two_square_sum/gen/smooth.cpp
@@ -0,0 +1,33 @@
+#include <iostream>
+#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);
+
+  vector<ll> primes;
+  if (seed == 0) { primes = {2, 3, 5, 7, 11, 13, 17}; }
+  if (seed == 1) { primes = {2, 3 * 3, 5, 7 * 7, 11 * 11, 13, 17}; }
+  if (seed == 2) { primes = {2, 5, 13, 17, 29}; }
+  if (seed == 3) { primes = {5, 13, 17, 29, 37, 41}; }
+  if (seed == 4) { primes = {5, 13, 17}; }
+
+  auto F = [&]() -> ll {
+    ll x = 1;
+    while (x <= MAX_A / primes[0]) {
+      ll k = gen.uniform<int>(0, (primes.size()) - 1);
+      ll p = primes[k];
+      if (p > MAX_A / x) continue;
+      x *= p;
+    }
+    return x;
+  };
+
+  int Q = MAX_Q;
+  printf("%d\n", Q);
+  for (int i = 0; i < Q; i++) { printf("%lld\n", F()); }
+}
diff --git a/number_theory/two_square_sum/hash.json b/number_theory/two_square_sum/hash.json
new file mode 100644
index 000000000..e86767f01
--- /dev/null
+++ b/number_theory/two_square_sum/hash.json
@@ -0,0 +1,78 @@
+{
+  "big2_00.in": "b1ff3e437dcfa4fc71aeefc778abeb3d20817c49c6c062d22062d49f55900311",
+  "big2_00.out": "126a6079b6a399b42451ae4eb2e14ab92d730a409f756685faabd917c84aeb5f",
+  "big2_01.in": "15abd9777d2a819e2d77c38c3321ad8c958ebbb9fb280e4c0ef4cd28b9a5b7b5",
+  "big2_01.out": "612a652431529ee1b71350fcea821ddb10ab993e4539833f5f61ebce365fc18a",
+  "big2_02.in": "ef9b60b3978daf0901212d0609db7423f7d275cce5cc22347b00e6435cbf66f5",
+  "big2_02.out": "598a3ae7a390bf349e3057a6ee0138b0ac045df9f61bf813be394d9689be9213",
+  "big_semiprime_random_00.in": "f6c6a0cdecb4a51cd99ab5c58450c55d4ffb7a850c40a840aa52017b794207ad",
+  "big_semiprime_random_00.out": "ce6dfc38e2f0980238a13963430f479f6145e29d9ba4b4ad6fc2d6049dc1a4d0",
+  "big_semiprime_random_01.in": "020a9a3878409a88e481e407762db9776de31ca060f2bb66e4bffa8c7a8ab62e",
+  "big_semiprime_random_01.out": "0eb95512707bdb3ab52c940986a99de7a31bba1b7fa536676b71907ff52b4b40",
+  "example_00.in": "b9f1a59750f37363a2ce9860ac0fe18763030fb6b2ab726fe8e82208cbdd0b36",
+  "example_00.out": "5ce40c2573d0666dbc16832a343f93658a9f426b527c370ec41ac38cc31c2263",
+  "example_01.in": "029b9bc32892d2360c618ffd61d4d747a0babf90411d906713d184e1273935e1",
+  "example_01.out": "f325b189fa87322b4fdd0d80111834c0f34831210822e5fac17240d2ef8d7645",
+  "fixed_RNG_buster_00.in": "f4b9dc156a046a6ddbea461742d112bfd677212e5c9d0286d0c7914fa339467f",
+  "fixed_RNG_buster_00.out": "9a271f2a916b0b6ee6cecb2426f0b3206ef074578be55d9bc94f6f3fe3ab86aa",
+  "max_00.in": "677d6ee21d71b0d7109ce34deb51d8a4f76847589f36b4faddc24f1f9696435a",
+  "max_00.out": "32f4b7b8165ad0e8353524728527eed819558f247d8643c499c6e01d613c499c",
+  "max_01.in": "b5fcf44c6731b1658e29e7c720659ea6da9edc88f257c06f7612be25aa7ea23c",
+  "max_01.out": "21059dba01c7b4b583120d1f182c250798c9829706ac9240a20a5ebc13080506",
+  "max_02.in": "9909fa3b15e1a471dc8268d14e22d73647636540b105151341b035e9518bce36",
+  "max_02.out": "b83b860b6ac63a209ebe8f1bddf83df251bf83d9fba0631528d42e6b3a2db083",
+  "max_03.in": "5bccd055c9985ce42243d8ce5151198a136d15fc12dcf9b51df123aaf3fd99e1",
+  "max_03.out": "8be223a3702b05fb2ac4880696d734296b9d334ec3b391fff3873d67c5b0a3aa",
+  "max_04.in": "e6942922adfdd09f871dac7276871e0c0ed59e40f0be7c1f2e4e324e3f37278b",
+  "max_04.out": "34cdf840989b7b9da22eafe1c0431c2e4e9c90f7975f7118b0efd0b2cfc4cff4",
+  "num_of_solution_max_00.in": "2708bd767622db411cf3fd263a588f36a9d29e046fc0a70f502ca155112bd50c",
+  "num_of_solution_max_00.out": "ba6b433e539d4b28d3f6434400e5e70d6e7b1bc0039314fc40fe83b28f147caa",
+  "num_of_solution_max_01.in": "836331bdd91a0d0316311c2184b39b53a2fd54543d4d3b34c49badf88a7cb2da",
+  "num_of_solution_max_01.out": "daa915be573c4c5b77b5146b297233b964d9798c809e66d73b9cc27b173e02a1",
+  "power_of_2_3_5_13_00.in": "2948158d29a34b7ce9395cd2d593e36f2f93c0b6f3147aae2c05fb7bf3f244fc",
+  "power_of_2_3_5_13_00.out": "3d8117d54f3f4995bd988ca4f85128bd0228c8d2756a991003b3cde793e906d4",
+  "power_of_2_3_5_13_01.in": "ba3cc7d1aa2636fb47d8c5556ff7b6ac19c055d83eec1c45046db2a2882e3444",
+  "power_of_2_3_5_13_01.out": "733331b71bad0471d6e608f2ebb9ddbcb9a73271c9a19a70bed66b8734b8872c",
+  "power_of_2_3_5_13_02.in": "01f3c1f2787f3149d1bb8d1734fdc014e2d3a0cec670129c9bdc4743b5532645",
+  "power_of_2_3_5_13_02.out": "dda969dccd231e88c33d4f57502e0c3ba612339233262b0166c7f2303e011a4a",
+  "power_of_2_3_5_13_03.in": "60673faeb7efbd730581df5288e0cb1c210758be939bca76a479c309b0547c80",
+  "power_of_2_3_5_13_03.out": "f5450ea41aea45cc63f97cc213849497f96a7439a32c4561524d4a1971d47356",
+  "prime_1_mod_4_00.in": "7febb950922ee205d7845e55262e1f95e32960646176ae40fc1a6fe8f8049f2c",
+  "prime_1_mod_4_00.out": "34c71b86fa027b38d52cf25ca258ba3e7d2832bb1097220e1cf4a18eda573d1d",
+  "prime_1_mod_4_01.in": "5f7e846b3fb67ad68d83cb40cb215d3a94917609ada51016e45c1bd632d7e9d9",
+  "prime_1_mod_4_01.out": "fc3ca880012dc31d34a3ef7ab2449c42c0c7d84aa26be1db11180e7fc807ed42",
+  "prime_1_mod_4_02.in": "bccb22a3750bb791d306740b9f5033b5de0596ce70784d034050bc3aae760b18",
+  "prime_1_mod_4_02.out": "04d99e9a67b21b652d53ee426a8c027216a6aad9258100de9d5fa2fe96928311",
+  "prime_3_mod_4_00.in": "8513253b71035e3f6a2bebe163237ff0038b1ddf35ca7dc2574f085441ec788d",
+  "prime_3_mod_4_00.out": "56cf0eddf3379f6c97214bd16998261aecab2c19765ec2097cad997d4c54cd2b",
+  "random_00.in": "708fab146ee509592ee780ea08879ea067a95c40bba5cfa15c31a37e6d800b91",
+  "random_00.out": "0f91267e646fd8ed9cbceac84478fc725efd3ef5b0d561f6c04a241422bfc978",
+  "random_01.in": "89e49eda00d255047433bc1bd787a3ef4ba6625dec2b1fd59ada8664058b6a6e",
+  "random_01.out": "4aa6d11675d40f6bda28f0609fa0bfcfeff9fa1eaa0848141bbd3619f82131c8",
+  "random_02.in": "a812fff102f10bb3457014f82f3e71670febd767cb3e8dbdd15ffe69b790fdf8",
+  "random_02.out": "0371b728d1562e1131235bfead2e6d4070af56146c5673fca57e8f1243eb0e3e",
+  "random_03.in": "94fc703004ce27d94dc06eae1bd4c81dc0292ac33800639c261571d8552d5e76",
+  "random_03.out": "1dc1f3796664ceecd6a8acce51de22bf5dfef1247864b508754d7f19ee90ea2f",
+  "random_04.in": "c6aecb477eb159b64aa3991ca4e70996bd6b18f3ce0be6ec800aba8f3894a392",
+  "random_04.out": "b77fdfc650a0d5294de754eb0393f1ba50333514787e0e32488f00d8433728db",
+  "small_00.in": "3db8a2c06a16acef3f7f9a10d68c7f7b444f90fe303e0c495ff2c9246bdcf2e2",
+  "small_00.out": "c2d18418d995619440c541e8927152200dd707fa1c79e526755003ac04b144ff",
+  "small_01.in": "55f8894939cfc972bb977decdb03d9136d86e6bae836c801b89a6764e9e5b642",
+  "small_01.out": "67c8cc8a1af4d545c0848d88638df2c352d0bb907cb8e22ebe97873f9933faee",
+  "small_02.in": "84b4475e69d2debae9a454bcc58f5a5f36fddecce3f64a12bae015ae6e5163c7",
+  "small_02.out": "4cd18ab4bea303072f8b862451a1f57503383504b4c49cfe371772866d8f1b26",
+  "small_03.in": "690536a316f7dd7a530f47b7696edef88858aeaba0957580270d696b68eed57f",
+  "small_03.out": "4ad70ecb79bab1db696d875ec340616c56c675b35af64afc8ac81efc9278cfd1",
+  "small_04.in": "f0ffd843c4c00374a656d9e6078f6d830bfab3d9afa5c7b489d89dead0d388a7",
+  "small_04.out": "112908be42a641352e53706c61b8b20e58aaa226918228005b616215ee218f46",
+  "smooth_00.in": "30acabdbc75032daa4b896690d49b6ead90253c3427be98f3b199f1981daaa3d",
+  "smooth_00.out": "408d2c5c9ad6fa19c23ddc3c849f6284a3e59127d8f7bc7bb49923c3d2397cbf",
+  "smooth_01.in": "2231d085d65a8b113d848d6ce72c45ec7e1039e5e9c7b831d219cea3a4517794",
+  "smooth_01.out": "e5516ee2e81760f217e828d89e4342265a55eea819cb8101d2aa0161cea69060",
+  "smooth_02.in": "9efc4a5597b7e8746e52f850daa8f57cc9da48a9a1a7a9e85bc12b48338972f3",
+  "smooth_02.out": "11778914ef830d947e087b2e9ba82d1527a0cc885119d467920be4ede7d8b0cc",
+  "smooth_03.in": "d8ade3b676974636da4a8ada784738f6894b7955cdb389c32da6b77dfc38f34e",
+  "smooth_03.out": "fbd302944745e4e3058b7811ead0255eff6788bb5a3590c0183d5ef8fa446577",
+  "smooth_04.in": "ed8d1e98c1a861eff71832fcfc292af03f36a3ce7572ba5d340403687bb28d7e",
+  "smooth_04.out": "0016d80a8632d1b6597d4e9175286b24b1ab5f83e167a600b7ab3475e37856a8"
+}
\ No newline at end of file
diff --git a/number_theory/two_square_sum/info.toml b/number_theory/two_square_sum/info.toml
new file mode 100644
index 000000000..f70662efd
--- /dev/null
+++ b/number_theory/two_square_sum/info.toml
@@ -0,0 +1,48 @@
+title = 'Represent A Number As Two Square Sum'
+timelimit = 5.0
+forum = "https://github.com/yosupo06/library-checker-problems/issues/1241"
+
+[[tests]]
+    name = "example.in"
+    number = 2
+[[tests]]
+    name = "random.cpp"
+    number = 5
+[[tests]]
+    name = "small.cpp"
+    number = 5
+[[tests]]
+    name = "max.cpp"
+    number = 5
+[[tests]]
+    name = "big2.cpp"
+    number = 3
+[[tests]]
+    name = "big_semiprime_random.cpp"
+    number = 2
+[[tests]]
+    name = "fixed_RNG_buster.in"
+    number = 1
+
+[[tests]]
+    name = "num_of_solution_max.in"
+    number = 2
+
+[[tests]]
+    name = "power_of_2_3_5_13.cpp"
+    number = 4
+
+[[tests]]
+   name = "prime_1_mod_4.cpp"
+   number = 3
+[[tests]]
+    name = "prime_3_mod_4.cpp"
+    number = 1
+[[tests]]
+    name = "smooth.cpp"
+    number = 5
+
+
+[params]
+    MAX_Q = 100
+    MAX_A = 1_000_000_000_000_000_000
diff --git a/number_theory/two_square_sum/lib/montgomery.hpp b/number_theory/two_square_sum/lib/montgomery.hpp
new file mode 100644
index 000000000..9f18f6f64
--- /dev/null
+++ b/number_theory/two_square_sum/lib/montgomery.hpp
@@ -0,0 +1,120 @@
+#ifndef MONTGOMERY_HPP
+#define MONTGOMERY_HPP
+
+#include <utility>
+#include "util.hpp"
+
+class MontgomeryGenerator;
+
+class Montgomery {
+    private:
+
+    static const u64 MASK64 = ~0;
+    static const u64 MASK32 = (1LL << 32) - 1;
+
+    u64 MOD;
+    u64 INV; // = (-MOD)^(-1) mod 2^64
+    u64 v;
+
+    // a * b = res.first * 2^64 + res.second
+    std::pair<u64, u64> mul(u64 a, u64 b) const {
+        u64 a1  = a >> 32,
+            a2  = a & MASK32,
+            b1  = b >> 32,
+            b2  = b & MASK32,
+            c11 = a1 * b1,
+            c12 = a1 * b2,
+            c21 = a2 * b1,
+            c22 = a2 * b2;
+        u64 carry = ((c22 >> 32) + (c12 & MASK32) + (c21 & MASK32)) >> 32;
+        return std::make_pair(
+            c11 + (c12 >> 32) + (c21 >> 32) + carry,
+            (c12 << 32) + (c21 << 32) + c22
+        );
+    }
+
+    // x * (2^64)^-1 mod MOD
+    u64 reduce(u64 x2) const {
+        if(x2 == 0) return 0;
+        u64 b = x2 * INV;
+        // (x + MOD * b) / 2^64
+        auto y = mul(MOD, b);
+        return y.first + 1;
+    }
+
+    // (x1 * 2^64 + x2) * (2^64)^-1 mod MOD
+    u64 reduce(u64 x1, u64 x2) const {
+        if(x2 == 0) return x1;
+        u64 b = x2 * INV;
+        // (x + MOD * b) / 2^64
+        auto res = mul(MOD, b).first + x1 + 1;
+        if(MOD <= res) res -= MOD;
+        return res;
+    }
+
+    Montgomery make_montgomery(u64 x) {
+        return Montgomery(MOD, INV, x);
+    }
+
+    Montgomery(u64 mod, u64 inv, u64 x) : MOD(mod), INV(inv), v(modprod(x, (~MOD) + 1, MOD)) {
+        assert(mod < (1ULL << 63));
+    }
+
+    public:
+
+    u64 val() const {
+        return reduce(v);
+    }
+
+    Montgomery pow(u64 k) {
+        Montgomery now = *this, ret = make_montgomery(1);
+        while(k) {
+            if(k % 2 == 1) {
+                ret = ret * now;
+            }
+            k /= 2;
+            now = now * now;
+        }
+        return ret;
+    }
+
+    Montgomery operator*(const Montgomery &lhs) const {
+        assert(MOD == lhs.MOD);
+        auto tmp = mul(v, lhs.v);
+        auto res = *this;
+        res.v = reduce(tmp.first, tmp.second);
+        return res;
+    }
+
+    friend class MontgomeryGenerator;
+};
+
+class MontgomeryGenerator {
+    private:
+    const u64 MOD;
+    const u64 INV; // = (-MOD)^(-1) mod 2^64
+
+    u64 calc_inv() {
+        u64 res = 0,
+            t   = 0,
+            i   = 1;
+        for(int j = 0; j < 64; j++) {
+            if(t % 2 == 0) {
+                t += MOD;
+                res += i;
+            }
+            t /= 2;
+            i *= 2;
+        }
+        return res;
+    }
+
+    public:
+    MontgomeryGenerator(u64 m) : MOD(m), INV(calc_inv()) {}
+
+    Montgomery gen(u64 x) const {
+        return Montgomery(MOD, INV, x);
+    }
+};
+
+#endif // MONTGOMERY_HPP
diff --git a/number_theory/two_square_sum/lib/prime.hpp b/number_theory/two_square_sum/lib/prime.hpp
new file mode 100644
index 000000000..cca69e4b5
--- /dev/null
+++ b/number_theory/two_square_sum/lib/prime.hpp
@@ -0,0 +1,88 @@
+#ifndef PRIME_HPP
+#define PRIME_HPP
+
+#include <numeric>
+#include "util.hpp"
+#include "montgomery.hpp"
+
+// determine if x is prime.
+bool miller_rabin(u64 p) {
+    static const u64 witnesses[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
+
+    if(p == 1 || p % 2 == 0) {
+        return p == 2;
+    }
+
+    MontgomeryGenerator mg(p);
+    
+    // p - 1 = 2 ** r * d
+    u64 r = 0, d = p - 1;
+    while(d % 2 == 0) {
+        d /= 2;
+        r++;
+    }
+    
+    for(u64 wt: witnesses) {
+        if(wt % p == 0) {
+            continue;
+        }
+
+        Montgomery y = mg.gen(wt).pow(d);
+        u64 v = y.val();
+
+        if(v == 1 || v == p - 1) {
+            continue;
+        }
+        for(u64 c = 0; c < r; c++) {
+            y = y * y;
+            v = y.val();
+            if(v == p - 1) {
+                break;
+            }
+        }
+        if(v != p - 1) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+bool is_prime(u64 n) {
+    u64 primes[] = {2, 3, 5, 7, 11, 13, 17, 19};
+    for(u64 p : primes) {
+        if(p == n) return true;
+        if(n % p == 0) {
+            return false;
+        }
+    }
+    return miller_rabin(n);
+}
+
+class PrimeGenerator {
+    static const std::size_t M_LEN = 8;
+    static const std::size_t D     = 30;
+    static const int M[M_LEN];
+
+    size_t idx;
+
+    public:
+    PrimeGenerator() : idx(0) {}
+
+    // WARNIG: this method doesn't generate prime uniformly.
+    u64 gen(u64 upper, Random &gen) {
+        u64 gen_upper = upper / D;
+        assert(1 < gen_upper);
+        gen_upper--;
+        u64 res;
+        do {
+            res = gen.uniform<u64>(0, gen_upper) * D + M[idx];
+            idx = (idx == M_LEN - 1) ? 0 : idx + 1;
+        } while(!miller_rabin(res));
+        return res;
+    }
+};
+
+const int PrimeGenerator::M[8] = {1, 7, 11, 13, 17, 19, 23, 29};
+
+#endif // PRIME_HPP
diff --git a/number_theory/two_square_sum/lib/util.hpp b/number_theory/two_square_sum/lib/util.hpp
new file mode 100644
index 000000000..02dce2bc1
--- /dev/null
+++ b/number_theory/two_square_sum/lib/util.hpp
@@ -0,0 +1,40 @@
+#ifndef UTIL_HPP
+#define UTIL_HPP
+
+#include <cstdint>
+#include <cassert>
+#include <limits>
+#include "random.h"
+
+using u64 = std::uint64_t;
+
+// x * y mod m
+u64 modprod(u64 x, u64 y, u64 m) {
+    assert(m <= std::numeric_limits<u64>::max());
+
+    x %= m;
+    y %= m;
+
+    if(m <= std::numeric_limits<u64>::max() / m) {
+        return x * y % m;
+    }
+
+    u64 res = 0, crt = x;
+    while(0 < y) {
+        if(y & 1) {
+            res += crt;
+            if(m <= res) {
+                res -= m;
+            }
+        }
+        crt *= 2;
+        if(m <= crt) {
+            crt -= m;
+        }
+        y /= 2;
+    }
+    return res;
+}
+
+
+#endif // UTIL_HPP
diff --git a/number_theory/two_square_sum/sol/correct.cpp b/number_theory/two_square_sum/sol/correct.cpp
new file mode 100644
index 000000000..09ac266e3
--- /dev/null
+++ b/number_theory/two_square_sum/sol/correct.cpp
@@ -0,0 +1,549 @@
+#include <cstdlib>
+#include <cstdio>
+#include <vector>
+#include <numeric>
+#include <algorithm>
+#include <tuple>
+#include <cassert>
+#include <cstdint>
+#include <limits>
+
+using namespace std;
+
+using ll = long long;
+using u32 = unsigned int;
+using u64 = unsigned long long;
+using i128 = __int128;
+using u128 = unsigned __int128;
+
+using pi = pair<ll, ll>;
+using vi = vector<ll>;
+template <class T>
+using vc = vector<T>;
+
+// https://trap.jp/post/1224/
+#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 FOR_subset(t, s) for (ll t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s)))
+#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
+
+// (0, 1, 2, 3, 4) -> (-1, 0, 1, 0, 2)
+int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
+int lowbit(u32 x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
+int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
+int lowbit(u64 x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
+int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
+int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
+
+template <typename T>
+T floor(T a, T b) {
+  return a / b - (a % b && (a ^ b) < 0);
+}
+template <typename T, typename U>
+T SUM(const vector<U> &A) {
+  T sm = 0;
+  for (auto &&a: A) sm += a;
+  return sm;
+}
+
+#define MIN(v) *min_element(all(v))
+#define MAX(v) *max_element(all(v))
+#define LB(c, x) distance((c).begin(), lower_bound(all(c), (x)))
+#define UB(c, x) distance((c).begin(), upper_bound(all(c), (x)))
+#define UNIQUE(x) sort(all(x)), x.erase(unique(all(x)), x.end()), x.shrink_to_fit()
+
+template <typename T>
+T POP(vc<T> &que) {
+  T a = que.back();
+  que.pop_back();
+  return a;
+}
+
+template <class T, class S>
+inline bool chmax(T &a, const S &b) {
+  return (a < b ? a = b, 1 : 0);
+}
+template <class T, class S>
+inline bool chmin(T &a, const S &b) {
+  return (a > b ? a = b, 1 : 0);
+}
+
+u64 RNG_64() {
+  // library checker なので fixed seed
+  static uint64_t x_ = 10402452321532715448ULL;
+  x_ ^= x_ << 7;
+  return x_ ^= x_ >> 9;
+}
+
+u64 RNG(u64 lim) { return RNG_64() % lim; }
+
+ll RNG(ll l, ll r) { return l + RNG_64() % (r - l); }
+
+// odd mod.
+// x の代わりに rx を持つ
+template <int id, typename U1, typename U2>
+struct Mongomery_modint {
+  using mint = Mongomery_modint;
+  inline static U1 m, r, n2;
+  static constexpr int W = numeric_limits<U1>::digits;
+
+  static void set_mod(U1 mod) {
+    assert(mod & 1 && mod <= U1(1) << (W - 2));
+    m = mod, n2 = -U2(m) % m, r = m;
+    FOR(5) r *= 2 - m * r;
+    r = -r;
+    assert(r * m == U1(-1));
+  }
+  static U1 reduce(U2 b) { return (b + U2(U1(b) * r) * m) >> W; }
+
+  U1 x;
+  Mongomery_modint() : x(0) {}
+  Mongomery_modint(U1 x) : x(reduce(U2(x) * n2)){};
+  U1 val() const {
+    U1 y = reduce(x);
+    return y >= m ? y - m : y;
+  }
+  mint &operator+=(mint y) {
+    x = ((x += y.x) >= m ? x - m : x);
+    return *this;
+  }
+  mint &operator-=(mint y) {
+    x -= (x >= y.x ? y.x : y.x - m);
+    return *this;
+  }
+  mint &operator*=(mint y) {
+    x = reduce(U2(x) * y.x);
+    return *this;
+  }
+  mint operator+(mint y) const { return mint(*this) += y; }
+  mint operator-(mint y) const { return mint(*this) -= y; }
+  mint operator*(mint y) const { return mint(*this) *= y; }
+  bool operator==(mint y) const { return (x >= m ? x - m : x) == (y.x >= m ? y.x - m : y.x); }
+  bool operator!=(mint y) const { return not operator==(y); }
+  mint pow(ll n) const {
+    assert(n >= 0);
+    mint y = 1, z = *this;
+    for (; n; n >>= 1, z *= z)
+      if (n & 1) y *= z;
+    return y;
+  }
+};
+
+template <int id>
+using Mongomery_modint_32 = Mongomery_modint<id, u32, u64>;
+template <int id>
+using Mongomery_modint_64 = Mongomery_modint<id, u64, u128>;
+
+bool primetest(const u64 x) {
+  assert(x < u64(1) << 62);
+  if (x == 2 or x == 3 or x == 5 or x == 7) return true;
+  if (x % 2 == 0 or x % 3 == 0 or x % 5 == 0 or x % 7 == 0) return false;
+  if (x < 121) return x > 1;
+  const u64 d = (x - 1) >> lowbit(x - 1);
+
+  using mint = Mongomery_modint_64<202311020>;
+
+  mint::set_mod(x);
+  const mint one(u64(1)), minus_one(x - 1);
+  auto ok = [&](u64 a) -> bool {
+    auto y = mint(a).pow(d);
+    u64 t = d;
+    while (y != one && y != minus_one && t != x - 1) y *= y, t <<= 1;
+    if (y != minus_one && t % 2 == 0) return false;
+    return true;
+  };
+  if (x < (u64(1) << 32)) {
+    for (u64 a: {2, 7, 61})
+      if (!ok(a)) return false;
+  } else {
+    for (u64 a: {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
+      if (!ok(a)) return false;
+    }
+  }
+  return true;
+}
+
+template <typename mint>
+ll rho(ll n, ll c) {
+  assert(n > 1);
+  const mint cc(c);
+  auto f = [&](mint x) { return x * x + cc; };
+  mint x = 1, y = 2, z = 1, q = 1;
+  ll g = 1;
+  const ll m = 1LL << (topbit(n) / 5);
+  for (ll r = 1; g == 1; r <<= 1) {
+    x = y;
+    FOR(r) y = f(y);
+    for (ll k = 0; k < r && g == 1; k += m) {
+      z = y;
+      FOR(min(m, r - k)) y = f(y), q *= x - y;
+      g = gcd(q.val(), n);
+    }
+  }
+  if (g == n) do {
+      z = f(z);
+      g = gcd((x - z).val(), n);
+    } while (g == 1);
+  return g;
+}
+
+ll find_prime_factor(ll n) {
+  assert(n > 1);
+  if (primetest(n)) return n;
+  FOR(100) {
+    ll m = 0;
+    if (n < (1 << 30)) {
+      using mint = Mongomery_modint_32<20231025>;
+      mint::set_mod(n);
+      m = rho<mint>(n, RNG(0, n));
+    } else {
+      using mint = Mongomery_modint_64<20231025>;
+      mint::set_mod(n);
+      m = rho<mint>(n, RNG(0, n));
+    }
+    if (primetest(m)) return m;
+    n = m;
+  }
+  assert(0);
+  return -1;
+}
+
+// ソートしてくれる
+vc<pair<ll, int>> factor(ll n) {
+  assert(n >= 1);
+  vc<pair<ll, int>> pf;
+  FOR(p, 2, 100) {
+    if (p * p > n) break;
+    if (n % p == 0) {
+      ll e = 0;
+      do { n /= p, e += 1; } while (n % p == 0);
+      pf.eb(p, e);
+    }
+  }
+  while (n > 1) {
+    ll p = find_prime_factor(n);
+    ll e = 0;
+    do { n /= p, e += 1; } while (n % p == 0);
+    pf.eb(p, e);
+  }
+  sort(all(pf));
+  return pf;
+}
+
+vc<pair<ll, int>> factor_by_lpf(ll n, vc<int> &lpf) {
+  vc<pair<ll, int>> res;
+  while (n > 1) {
+    int p = lpf[n];
+    int e = 0;
+    while (n % p == 0) {
+      n /= p;
+      ++e;
+    }
+    res.eb(p, e);
+  }
+  return res;
+}
+
+// https://github.com/atcoder/ac-library/blob/master/atcoder/internal_math.hpp
+struct Barrett {
+  u32 m;
+  u64 im;
+  explicit Barrett(u32 m = 1) : m(m), im(u64(-1) / m + 1) {}
+  u32 umod() const { return m; }
+  u32 modulo(u64 z) {
+    if (m == 1) return 0;
+    u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
+    u64 y = x * m;
+    return (z - y + (z < y ? m : 0));
+  }
+  u64 floor(u64 z) {
+    if (m == 1) return z;
+    u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
+    u64 y = x * m;
+    return (z < y ? x - 1 : x);
+  }
+  pair<u64, u32> divmod(u64 z) {
+    if (m == 1) return {z, 0};
+    u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
+    u64 y = x * m;
+    if (z < y) return {x - 1, z - y + m};
+    return {x, z - y};
+  }
+  u32 mul(u32 a, u32 b) { return modulo(u64(a) * b); }
+};
+
+struct Barrett_64 {
+  u128 mod, mh, ml;
+
+  explicit Barrett_64(u64 mod = 1) : mod(mod) {
+    u128 m = u128(-1) / mod;
+    if (m * mod + mod == u128(0)) ++m;
+    mh = m >> 64;
+    ml = m & u64(-1);
+  }
+
+  u64 umod() const { return mod; }
+
+  u64 modulo(u128 x) {
+    u128 z = (x & u64(-1)) * ml;
+    z = (x & u64(-1)) * mh + (x >> 64) * ml + (z >> 64);
+    z = (x >> 64) * mh + (z >> 64);
+    x -= z * mod;
+    return x < mod ? x : x - mod;
+  }
+
+  u64 mul(u64 a, u64 b) { return modulo(u128(a) * b); }
+};
+
+u32 mod_pow(int a, ll n, int mod) {
+  assert(n >= 0);
+  a = ((a %= mod) < 0 ? a + mod : a);
+  if ((mod & 1) && (mod < (1 << 30))) {
+    using mint = Mongomery_modint_32<202311021>;
+    mint::set_mod(mod);
+    return mint(a).pow(n).val();
+  }
+  Barrett bt(mod);
+  int r = 1;
+  while (n) {
+    if (n & 1) r = bt.mul(r, a);
+    a = bt.mul(a, a), n >>= 1;
+  }
+  return r;
+}
+
+u64 mod_pow_64(ll a, ll n, u64 mod) {
+  assert(n >= 0);
+  a = ((a %= mod) < 0 ? a + mod : a);
+  if ((mod & 1) && (mod < (u64(1) << 62))) {
+    using mint = Mongomery_modint_64<202311021>;
+    mint::set_mod(mod);
+    return mint(a).pow(n).val();
+  }
+  Barrett_64 bt(mod);
+  ll r = 1;
+  while (n) {
+    if (n & 1) r = bt.mul(r, a);
+    a = bt.mul(a, a), n >>= 1;
+  }
+  return r;
+}
+
+template <typename T>
+struct Gaussian_Integer {
+  T x, y;
+  using G = Gaussian_Integer;
+
+  Gaussian_Integer(T x = 0, T y = 0) : x(x), y(y) {}
+  Gaussian_Integer(pair<T, T> p) : x(p.fi), y(p.se) {}
+
+  T norm() const { return x * x + y * y; }
+  G conjugate() const { return G(x, -y); }
+
+  G &operator+=(const G &g) {
+    x += g.x, y += g.y;
+    return *this;
+  }
+  G &operator-=(const G &g) {
+    x -= g.x, y -= g.y;
+    return *this;
+  }
+  G &operator*=(const G &g) {
+    tie(x, y) = mp(x * g.x - y * g.y, x * g.y + y * g.x);
+    return *this;
+  }
+  G &operator/=(const G &g) {
+    *this *= g.conjugate();
+    T n = g.norm();
+    x = floor(x + n / 2, n);
+    y = floor(y + n / 2, n);
+    return *this;
+  }
+  G &operator%=(const G &g) {
+    auto q = G(*this) / g;
+    q *= g;
+    (*this) -= q;
+    return *this;
+  }
+  G operator-() { return G(-x, -y); }
+  G operator+(const G &g) { return G(*this) += g; }
+  G operator-(const G &g) { return G(*this) -= g; }
+  G operator*(const G &g) { return G(*this) *= g; }
+  G operator/(const G &g) { return G(*this) /= g; }
+  G operator%(const G &g) { return G(*this) %= g; }
+  bool operator==(const G &g) { return (x == g.x && y == g.y); }
+
+  static G gcd(G a, G b) {
+    while (b.x != 0 || b.y != 0) {
+      a %= b;
+      swap(a, b);
+    }
+    return a;
+  }
+
+  G pow(ll n) const {
+    assert(n >= 0);
+    G ret(1), mul(*this);
+    while (n > 0) {
+      if (n & 1) ret *= mul;
+      mul *= mul;
+      n >>= 1;
+    }
+    return ret;
+  }
+
+  // (g,x,y) s.t ax+by=g
+  static tuple<G, G, G> extgcd(G a, G b) {
+    if (b.x != 0 || b.y != 0) {
+      G q = a / b;
+      auto [g, x, y] = extgcd(b, a - q * b);
+      return {g, y, x - q * y};
+    }
+    return {a, G{1, 0}, G{0, 0}};
+  }
+};
+
+pair<ll, ll> solve_norm_equation_prime(ll p) {
+  using G = Gaussian_Integer<i128>;
+  assert(p == 2 || p % 4 == 1);
+  if (p == 2) return {1, 1};
+  ll x = [&]() -> ll {
+    ll x = 1;
+    while (1) {
+      ++x;
+      ll pow_x = 1;
+      if (p < (1 << 30)) {
+        pow_x = mod_pow(x, (p - 1) / 4, p);
+        if (pow_x * pow_x % p == p - 1) return pow_x;
+      } else {
+        pow_x = mod_pow_64(x, (p - 1) / 4, p);
+        if (i128(pow_x) * pow_x % p == p - 1) return pow_x;
+      }
+    }
+    return -1;
+  }();
+  G a(p, 0), b(x, 1);
+  a = G::gcd(a, b);
+  assert(a.norm() == p);
+  return {a.x, a.y};
+}
+
+template <typename T>
+vc<Gaussian_Integer<T>> solve_norm_equation_factor(vc<pair<ll, int>> pfs) {
+  using G = Gaussian_Integer<T>;
+  vc<G> res;
+  for (auto &&[p, e]: pfs) {
+    if (p % 4 == 3 && e % 2 == 1) return {};
+  }
+
+  res.eb(G(1, 0));
+  for (auto &&[p, e]: pfs) {
+    if (p % 4 == 3) {
+      T pp = 1;
+      FOR(e / 2) pp *= p;
+      for (auto &&g: res) {
+        g.x *= pp;
+        g.y *= pp;
+      }
+      continue;
+    }
+    auto [pix, piy] = solve_norm_equation_prime(p);
+    G pi(pix, piy);
+    vc<G> pows(e + 1);
+    pows[0] = G(1, 0);
+    FOR(i, e) pows[i + 1] = pows[i] * pi;
+    if (p == 2) {
+      for (auto &&g: res) g *= pows[e];
+      continue;
+    }
+    vc<G> pis(e + 1);
+    FOR(j, e + 1) { pis[j] = pows[j] * (pows[e - j].conjugate()); }
+    vc<G> new_res;
+    new_res.reserve(len(res) * (e + 1));
+    for (auto &&g: res) {
+      for (auto &&a: pis) { new_res.eb(g * a); }
+    }
+    swap(res, new_res);
+  }
+
+  for (auto &&g: res) {
+    while (g.x <= 0 || g.y < 0) { g = G(-g.y, g.x); }
+  }
+  return res;
+}
+
+// i128 を使うと N <= 10^{18} もできる
+// ノルムがとれるように、2 乗してもオーバーフローしない型を使おう
+// 0 <= arg < 90 となるもののみ返す。
+// 単数倍は作らないので、使うときに気を付ける。
+template <typename T>
+vc<Gaussian_Integer<T>> solve_norm_equation(T N) {
+  using G = Gaussian_Integer<T>;
+  vc<G> res;
+  if (N < 0) return {};
+  if (N == 0) {
+    res.eb(G(0, 0));
+    return res;
+  }
+  auto pfs = factor(N);
+  return solve_norm_equation_factor<T>(pfs);
+}
+
+void solve() {
+  ll N;
+  scanf("%lld", &N);
+  if (N == 0) {
+    printf("1\n");
+    printf("0 0\n");
+    return;
+  }
+  auto pfs = factor(N);
+
+  vc<pair<ll, ll>> ANS;
+  for (auto &g: solve_norm_equation_factor<i128>(pfs)) {
+    ll x = g.x, y = g.y;
+    FOR(4) {
+      tie(x, y) = mp(-y, x);
+      if (x >= 0 && y >= 0) ANS.eb(x, y);
+    }
+  }
+
+  ll expected_ANS_cnt = 1;
+  bool square = 1;
+  for (auto &[p, e]: pfs) {
+    if (e & 1) square = 0;
+    if (p == 2) {};
+    if (p % 4 == 1) { expected_ANS_cnt *= e + 1; }
+    if (p % 4 == 3) { expected_ANS_cnt *= (e % 2 == 0 ? 1 : 0); }
+  }
+  if (square) expected_ANS_cnt++;
+  assert(expected_ANS_cnt == len(ANS));
+
+  printf("%lld\n", len(ANS));
+  for (auto &[a, b]: ANS) { printf("%lld %lld\n", a, b); }
+}
+
+signed main() {
+  int T;
+  scanf("%d", &T);
+  FOR(T) solve();
+  return 0;
+}
diff --git a/number_theory/two_square_sum/task.md b/number_theory/two_square_sum/task.md
new file mode 100644
index 000000000..c85f29a6e
--- /dev/null
+++ b/number_theory/two_square_sum/task.md
@@ -0,0 +1,51 @@
+## @{keyword.statement}
+
+@{lang.en}
+
+Given $Q$ cases. 
+
+For each case, you are given an integer $N_i$. Find all pairs of non-negative integers (a,b) statisfying $N_i=a^2+b^2$. 
+
+@{lang.ja}
+
+クエリが $Q$ 個与えられます．
+
+各クエリでは整数 $N_i$ が与えられるので，非負整数の組 $(a,b)$ であって $N_i=a^2+b^2$ を満たすものをすべて求めてください．
+
+@{lang.end}
+
+## @{keyword.constraints}
+
+- $1 \leq Q \leq @{param.MAX_Q}$
+- $0 \leq N_i \leq @{param.MAX_A}$
+
+## @{keyword.input}
+
+```
+$Q$
+$N_0$
+$N_1$
+$\vdots$
+$N_{Q - 1}$
+```
+
+## @{keyword.output}
+
+@{lang.en}
+For each test case, let $n$ be the number of pairs $(a,b)$ that satisfy the condition. Output a total of $n+1$ lines in the following format. The output order can be arbitrary.
+@{lang.ja}
+テストケースごとに，条件を満たす組 $(a,b)$ の個数を $n$ として次の形式で合計 $n+1$ 行出力してください．どのような順番で出力してもよいです．
+@{lang.end}
+
+```
+$n$
+$a$ $b$
+$\vdots$
+$a$ $b$
+```
+
+## @{keyword.sample}
+
+@{example.example_00}
+
+@{example.example_01}
diff --git a/number_theory/two_square_sum/verifier.cpp b/number_theory/two_square_sum/verifier.cpp
new file mode 100644
index 000000000..0e234a892
--- /dev/null
+++ b/number_theory/two_square_sum/verifier.cpp
@@ -0,0 +1,15 @@
+#include <iostream>
+#include "params.h"
+#include "testlib.h"
+
+int main() {
+  registerValidation();
+  int q = inf.readInt(1, MAX_Q);
+  inf.readChar('\n');
+  for (int i = 0; i < q; i++) {
+    inf.readLong(0, MAX_A);
+    inf.readChar('\n');
+  }
+  inf.readEof();
+  return 0;
+}