Paste
Of Code

#include <bits/stdc++.h>
#include <chrono>
#include <random>

using namespace std;

#define intL unsigned __int128

void print128(intL n) {
  if (n == 0) {
    cout << 0;
    return;
  }
  string s = "";
  while (n > 0) {
    s += (char)('0' + (n % 10));
    n /= 10;
  }
  reverse(s.begin(), s.end());
  cout << s;
}

intL gcd(intL a, intL b) {
  while (b) {
    intL temp = b;
    b = a % b;
    a = temp;
  }
  return a;
}

intL mul(intL a, intL b, intL MOD) {
  return (intL)(((unsigned __int128)a * b) % MOD);
}

intL power(intL base, intL exp, intL mod) {
  intL res = 1;
  base %= mod;
  while (exp > 0) {
    if (exp % 2 == 1)
      res = mul(res, base, mod);
    base = mul(base, base, mod);
    exp /= 2;
  }
  return res;
}

bool miller_rabin(intL n) {
  if (n < 2)
    return false;
  if (n == 2 || n == 3)
    return true;
  if (n % 2 == 0 || n % 3 == 0)
    return false;

  intL d = n - 1;
  int s = 0;
  while (d % 2 == 0) {
    d /= 2;
    s++;
  }

  static const unsigned long long bases[] = {2,  3,  5,  7,  11, 13,
                                             17, 19, 23, 29, 31, 37};
  for (unsigned long long a : bases) {
    if (n <= a)
      break;
    intL x = power(a, d, n);
    if (x == 1 || x == n - 1)
      continue;
    bool composite = true;
    for (int r = 1; r < s; r++) {
      x = mul(x, x, n);
      if (x == n - 1) {
        composite = false;
        break;
      }
    }
    if (composite)
      return false;
  }
  return true;
}

const int MAXS = 1000000;
int spf[MAXS >> 1];

int get_spf(int s) {
  if (!(s & 1))
    return 2;
  int idx = s >> 1;
  return spf[idx] == -1 ? s : spf[idx];
}

intL brent(intL N) {
  if (N % 2 == 0)
    return 2;
  if (N % 3 == 0)
    return 3;

  mt19937_64 rng(1337);
  uniform_int_distribution<long long> seed(1, N - 1);

  intL y = seed(rng), c = seed(rng), m = seed(rng) % (N - 1) + 1, g = 1, r = 1,
       q = 1, x = y, ys = y;

  int iterations = 0;
  while (g == 1 && iterations < 500000) {
    x = y;
    for (long long i = 0; i < r; ++i) {
      y = (mul(y, y, N) + c) % N;
    }
    intL k = 0;
    while (k < r && g == 1) {
      ys = y;
      intL lim = min(m, r - k);
      for (intL i = 0; i < lim; ++i) {
        y = (mul(y, y, N) + c) % N;
        intL diff = x > y ? x - y : y - x;
        q = mul(q, diff, N);
      }
      g = gcd(q, N);
      k += m;
      iterations += lim;
    }
    r <<= 1;
  }

  if (g == N || g == 1)
    return 0;
  return g;
}

intL squfof(intL N) {
  if (N % 2 == 0)
    return 2;

  unsigned long long n = (unsigned long long)N;
  unsigned long long root = sqrt(n);
  if (root * root == n)
    return root;

  auto scaled_squfof = [&](unsigned long long multiplier) {
    unsigned long long nn = n * multiplier;
    unsigned long long p0 = sqrt(nn);
    unsigned long long q0 = 1;
    unsigned long long p = p0;
    unsigned long long q = nn - p0 * p0;
    if (q == 0)
      return (unsigned long long)1;

    unsigned long long b0 = (p0 + p) / q;
    unsigned long long b = b0;
    unsigned long long l1 = 2 * b0, l2 = 1;
    unsigned long long i = 0;

    vector<unsigned long long> saved;

    while (i < 300000) {
      unsigned long long prev_q = q0;
      q0 = q;
      p = b * q0 - p;
      q = prev_q + b * (p0 - p);
      if (q == 0)
        break;
      b = (p0 + p) / q;

      if (i % 2 == 1) {
        unsigned long long root_q = sqrt(q);
        if (root_q * root_q == q) {
          saved.push_back(root_q);
        }
      }
      i++;
    }
    return (unsigned long long)0;
  };

  static const unsigned long long multipliers[] = {1,  3,  5,  7, 11,
                                                   13, 17, 19, 23};
  for (auto mult : multipliers) {
    if (__int128(n) * mult > ~((unsigned long long)0))
      continue;
    unsigned long long f = scaled_squfof(mult);
    if (f > 1 && f < n) {
      intL g = gcd(f, N);
      if (g > 1 && g < N)
        return g;
    }
  }
  return 0;
}

void factorize(intL N, map<intL, intL> &mp) {
  if (N == 1)
    return;

  if (N < MAXS) {
    while (N > 1) {
      int p = get_spf(N);
      while (N % p == 0) {
        ++mp[p];
        N /= p;
      }
    }
    return;
  }

  for (intL p = 2; p * p <= N && p < MAXS; ++p) {
    while (N % p == 0) {
      ++mp[p];
      N /= p;
    }
  }
  if (N == 1)
    return;

  if (miller_rabin(N)) {
    ++mp[N];
    return;
  }

  intL Q = brent(N);
  if (Q == 0 || Q == 1 || Q == N) {
    Q = squfof(N);
  }

  if (Q == 0 || Q == 1 || Q == N) {
    ++mp[N];
    return;
  }

  factorize(Q, mp);
  factorize(N / Q, mp);
}

signed main() {
  ios::sync_with_stdio(false);
  cin.tie(0);

  fill(spf, spf + (MAXS >> 1), -1);
  for (int i = 3; i * i <= MAXS; i += 2) {
    if (spf[i >> 1] == -1) {
      for (int j = i * i; j <= MAXS; j += (i << 1))
        if (spf[j >> 1] == -1)
          spf[j >> 1] = i;
    }
  }

  string s;
  while (cin >> s && s != "0") {
    intL N = 0;
    for (char ch : s)
      N = N * 10 + (ch - '0');

    map<intL, intL> mp;
    factorize(N, mp);

    for (auto const &[prime, exp] : mp) {
      print128(prime);
      cout << "^";
      print128(exp);
      cout << " ";
    }
    cout << "\n";
  }

  return 0;
}