diff --git a/math/factorize/gen/big_semiprime_gen.cpp b/math/factorize/gen/big_semiprime_gen.cpp
new file mode 100644
index 000000000..938ef512f
--- /dev/null
+++ b/math/factorize/gen/big_semiprime_gen.cpp
@@ -0,0 +1,70 @@
+#include <iostream>
+#include <cmath>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+using ll = long long;
+
+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 generate_semiprime(ll s) {
+    for (ll t = MAX_A / s; ; --t) {
+        if (is_prime(t)) return s * t;
+    }
+}
+
+int main(int, char* argv[]) {
+    ll seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    // 1 <= t/s <= 100
+    ll gen_max = sqrt(MAX_A);
+    ll gen_min = sqrt(MAX_A / 100);
+
+    int q = MAX_Q;
+    vector<ll> a(q);
+    for (int i = 0; i < q; i++) {
+        ll s;
+        do {
+            s = gen.uniform(gen_min, gen_max);
+        } while(!is_prime(s));
+        a[i] = generate_semiprime(s);
+    }
+
+    printf("%d\n", q);
+    for (auto x: a) {
+        printf("%lld\n", x);
+    }
+    return 0;
+}
diff --git a/math/factorize/gen/big_semiprime_random.cpp b/math/factorize/gen/big_semiprime_random.cpp
new file mode 100644
index 000000000..32f88ead1
--- /dev/null
+++ b/math/factorize/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/math/factorize/hash.json b/math/factorize/hash.json
index fa9afabda..634c6fa60 100644
--- a/math/factorize/hash.json
+++ b/math/factorize/hash.json
@@ -11,6 +11,14 @@
   "big2_02.out": "ebd3b9062066121868f812ab9910492865b8472fefd71d392df5f20718bbd331",
   "big2_worse_00.in": "46df5fc652c3cc15dc3d5277c90525d0437b37d8b9c1a3c90bb962f329a2b33b",
   "big2_worse_00.out": "c2b5fa215419e9175ae1454a8d79c72bed489ba6bbb21b41c5e27cc55b821342",
+  "big_semiprime_gen_00.in": "70c3f9ce57a5456f6a066b942cbbc44bf408e22a8cc8c32008dbff14efadb3e9",
+  "big_semiprime_gen_00.out": "86a075f0deb5092008a9275b48e77974ca8ef174d449d35210c5f394c46d4809",
+  "big_semiprime_gen_01.in": "1de2462f9291762ad848ac2e67a8b1c639a253ebc33f4ed9e72c8cf62474bb4e",
+  "big_semiprime_gen_01.out": "1e04525609206b2e023eac56c4725d2fc77069993ed22b48df62dcfbe022adfe",
+  "big_semiprime_random_00.in": "f6c6a0cdecb4a51cd99ab5c58450c55d4ffb7a850c40a840aa52017b794207ad",
+  "big_semiprime_random_00.out": "3cfec6469acf1a6111b420265acd520871b9880d13d460593885e25acdce14a7",
+  "big_semiprime_random_01.in": "020a9a3878409a88e481e407762db9776de31ca060f2bb66e4bffa8c7a8ab62e",
+  "big_semiprime_random_01.out": "2d042dd98197d97635464c04e6b0739693805006f04ecaf83009339d7d389bd2",
   "carmichael_00.in": "88928eedee7789093e2dddfa5a47bb5dd241851e6cfed7ee71b437afb9228daa",
   "carmichael_00.out": "65ebadfe1f3be538d67bb7f11c3591d3f65cb0b8406ef6812daf29c3e0e39f08",
   "example_00.in": "b4791ee88c8a679a0c7a4d038f84a9543bc4a7995ecbfc6db562e74245579dc6",
diff --git a/math/factorize/info.toml b/math/factorize/info.toml
index 837f380da..3fe5e4467 100644
--- a/math/factorize/info.toml
+++ b/math/factorize/info.toml
@@ -20,6 +20,12 @@ forum = "https://github.com/yosupo06/library-checker-problems/issues/41"
 [[tests]]
     name = "big2_worse.cpp"
     number = 1
+[[tests]]
+    name = "big_semiprime_gen.cpp"
+    number = 2
+[[tests]]
+    name = "big_semiprime_random.cpp"
+    number = 2
 [[tests]]
     name = "pow2.cpp"
     number = 3
@@ -48,6 +54,9 @@ forum = "https://github.com/yosupo06/library-checker-problems/issues/41"
 [[solutions]]
     name = "fixed_RNG.cpp"
     wrong = true
+[[solutions]]
+    name = "rho_with_difference_of_square.cpp"
+    wrong = false
 
 [params]
     MAX_Q = 100
diff --git a/math/factorize/sol/rho_with_difference_of_square.cpp b/math/factorize/sol/rho_with_difference_of_square.cpp
new file mode 100644
index 000000000..221779f7e
--- /dev/null
+++ b/math/factorize/sol/rho_with_difference_of_square.cpp
@@ -0,0 +1,129 @@
+#include <iostream>
+#include <vector>
+#include <algorithm>
+#include <cmath>
+#include <map>
+
+using namespace std;
+using uint = unsigned int;
+using ll = long long;
+using ull = unsigned long long;
+constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n - 1); }
+template <class T> using V = vector<T>;
+template <class T> using VV = V<V<T>>;
+
+// bit op
+int popcnt(uint x) { return __builtin_popcount(x); }
+int popcnt(ull x) { return __builtin_popcountll(x); }
+int bsr(uint x) { return 31 - __builtin_clz(x); }
+int bsr(ull x) { return 63 - __builtin_clzll(x); }
+int bsf(uint x) { return __builtin_ctz(x); }
+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 difference_of_square(ll n) {
+    ll sqrt_n = ll(sqrt(n));
+    for (ll i = 0; i < 300; ++i) {
+        ll t = (sqrt_n + i) * (sqrt_n + i) - n;
+        ll sqrt_t = ll(sqrt(t));
+        if (sqrt_t * sqrt_t == t) return sqrt_n + i - sqrt_t;
+    }
+    return 0;
+}
+
+ll pollard_single(ll n) {
+    if (is_prime(n)) return n;
+    if (n % 2 == 0) return 2;
+    ll ds = difference_of_square(n);
+    if (ds > 1) return ds;
+    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 q;
+    scanf("%d", &q);
+    map<ll, V<ll>> cache;
+    for (int i = 0; i < q; i++) {
+        ll a;
+        scanf("%lld", &a);
+        if (!cache.count(a)) {
+            auto v = pollard(a);
+            sort(v.begin(), v.end());
+            cache[a] = v;
+        }
+        auto v = cache[a];
+        printf("%d", int(v.size()));
+        for (auto d: v) printf(" %lld", d);
+        printf("\n");
+    }
+    return 0;
+}