大数素幂分解

using i128 = __int128_t;
const int mod = 1000000007;

i64 fpow(i64 a, i64 b, i64 mod) {
  i64 ans = 1;
  i64 base = a;
  while (b) {
    if (b & 1) ans = (i128)ans * base % mod;
    b >>= 1;
    base = (i128)base * base % mod;
  }
  return ans;
}

bool Miller_Rabbin(i64 a, i64 n) {
  i64 r = 0, d = n - 1;
  while (d % 2 == 0) {
    d >>= 1;
    r++;
  }

  i64 x = fpow(a, d, n);

  if (x == 1) return true;
  for (int i = 0; i < r; i++) {
    if (x == n - 1) return true;
    x = (i128)x * x % n;
  }

  return false;
}

bool isprime(i64 n) {
  if (n == 2) return true;
  if (n < 2 || n % 2 == 0) return false;
  vector<i64> prime = {2, 3, 5, 7, 9, 233, 331};
  for (auto t : prime) {
    if (n == t) return true;
    if (!Miller_Rabbin(t, n)) return false;
  }
  return true;
}

i64 gcd(i64 a, i64 b) {
  return b == 0 ? a : gcd(b, a % b);
}

i64 aabs(i64 w) {
  return w > 0 ? w : -w;
}

i64 Pollard_Rho(i64 n) {

  i64 c = 1ll * rand() % (n - 1) + 1;
  auto f = [&](i64 x) -> i64 {
    return ((i128)x * x + c) % n;
  };
  i64 l = 0, r = 0, cur = 1;

  for (i64 i = 1; ; i <<= 1, l = r, cur = 1) {
    for (i64 j = 1; j <= i; j++) {
      r = f(r);
      cur = (i128)cur * aabs(r - l) % n;
      if (j % 128 == 0) {
        i64 d = gcd(cur, n);
        if (d > 1) return d;
      }
    }
    i64 d = gcd(cur, n);
    if (d > 1) return d;
  }
}

set<i64> w;

void divide(i64 n) {
  if (n == 1) return;
  if (isprime(n)) {
    w.insert(n);
    return;
  }
  i64 p = n;
  while (p >= n) {
    p = Pollard_Rho(n);
  }
  while (n % p == 0) {
    n /= p;
  }
  divide(n);
  divide(p);
}

// divide(n) 后
// w 中会存储所有n的质因子，然后可以再 log 计算出每个质因子的指数