// Proof-of-concept implementation of factoring into coprimes
// 2004.04.07, Daniel J. Bernstein
//
// To quote the paper: ``Warning: The algorithms in this paper are _not_
// designed to be as fast as possible. They are designed to be as _simple_
// as possible, under the constraints of (1) running in essentially linear
// time and (2) using a limited set of arithmetic operations. There are
// many ways to save more time; these speedups are important but are beyond
// the scope of this paper.''
//
// Sample use: echo '[91 119 221 1547 6898073]' | ./dcba-concept

#include <iostream>
#include <LiDIA/bigint.h>
#include <LiDIA/base_vector.h>
using namespace LiDIA;
using namespace std;

/* a^2^n */
bigint twopower(bigint a,unsigned int n)
{
  /* Algorithm 10.1 */
  for (;;) {
    if (!n) return a;
    a = a * a;
    n = n - 1;
  }
}

/* gout = gcd{a,c}; x = ppi(a,c); y = ppo(a,c) */
void gcdppippo(bigint &gout,bigint &x,bigint &y,bigint a,bigint c)
{
  /* Algorithm 11.3 */
  bigint g;
  x = gcd(a,c);
  gout = x;
  y = a / x;
  for (;;) {
    g = gcd(x,y);
    if (g == 1) return;
    x = x * g;
    y = y / g;
  }
}

/* x = ppi(a,c); y = ppo(a,c) */
void ppippo(bigint &x,bigint &y,bigint a,bigint c)
{
  bigint g;
  gcdppippo(g,x,y,a,c);
}

/* ppi(a,c) */
bigint ppi(bigint a,bigint c)
{
  bigint g;
  bigint x;
  bigint y;
  gcdppippo(g,x,y,a,c);
  return x;
}

/* gout = gcd{a,b}; x = ppg(a,b); y = pple(a,b) */
void gcdppgpple(bigint &gout,bigint &x,bigint &y,bigint a,bigint b)
{
  /* Algorithm 11.4 */
  bigint g;
  y = gcd(a,b);
  gout = y;
  x = a / y;
  for (;;) {
    g = gcd(x,y);
    if (g == 1) return;
    x = x * g;
    y = y / g;
  }
}

/* out = out + cb{a,b} */
void append_cb(base_vector<bigint> &out,bigint a,bigint b)
{
  /* Algorithm 13.2 */
  if (b == 1) {
    if (a != 1) out[out.size()] = a;
    return;
  }
  bigint r;
  ppippo(a,r,a,b);
  if (r != 1) out[out.size()] = r;
  bigint g;
  bigint h;
  bigint c;
  gcdppgpple(g,h,c,a,b);
  bigint c0 = c;
  bigint x = c0;
  unsigned int n = 1;
  for (;;) {
    gcdppgpple(g,h,c,h,g * g);
    bigint d = gcd(c,b);
    x = x * d;
    bigint y = twopower(d,n - 1);
    append_cb(out,c / y,d);
    if (h == 1) break;
    n = n + 1;
  }
  append_cb(out,b / x,c0);
}

/* prod S */
bigint prod(base_vector<bigint> S)
{
  /* Algorithm 14.1 */
  unsigned int n = S.size();
  if (!n) return 1;
  if (n == 1) return S[0];
  base_vector<bigint> T(0,EXPAND);
  T.assign(0,S,0,n/2 - 1);
  bigint X = prod(T);
  T.assign(0,S,n/2,n - 1);
  bigint Y = prod(T);
  return X * Y;
}

/* append each ppi(a,P[i]) to out */
void splitappend(base_vector<bigint> &out,bigint a,base_vector<bigint> P)
{
  /* Algorithm 15.3 */
  unsigned int n = P.size();
  if (!n) return;
  bigint b = ppi(a,prod(P));
  if (n == 1) { out[out.size()] = b; return; }
  base_vector<bigint> Q(0,EXPAND);
  Q.assign(0,P,0,n/2 - 1);
  splitappend(out,b,Q);
  Q.assign(0,P,n/2,n - 1);
  splitappend(out,b,Q);
}

/* split(a,P) */
base_vector<bigint> split(bigint a,base_vector<bigint> P)
{
  base_vector<bigint> result(0,EXPAND);
  splitappend(result,a,P);
  return result;
}

/* cb(P union {b}) */
base_vector<bigint> cbextend(base_vector<bigint> P,bigint b)
{
  /* Algorithm 16.2 */
  base_vector<bigint> result(0,EXPAND);
  if (!P.size()) {
    if (b != 1) result[0] = b;
    return result;
  }
  bigint x = prod(P);
  bigint a;
  bigint r;
  ppippo(a,r,b,x);
  if (r != 1) result[0] = r;
  base_vector<bigint> S(0,EXPAND);
  S = split(a,P);
  for (unsigned int i = 0;i < P.size();++i)
    append_cb(result,P[i],S[i]);
  return result;
}

/* bit_i k */
int bit(unsigned int i,unsigned int k)
{
  if (k & (1 << i)) return 1;
  return 0;
}

/* cb(P union Q) */
base_vector<bigint> cbmerge(base_vector<bigint> P,base_vector<bigint> Q)
{
  /* Algorithm 17.3 */
  base_vector<bigint> S(0,EXPAND);
  base_vector<bigint> T(0,EXPAND);
  base_vector<bigint> R(0,EXPAND);
  unsigned int n = Q.size();
  unsigned int b = 1;
  while ((1 << b) < n) ++b;
  S = P;
  unsigned int i = 0;
  for (;;) {
    if (i == b) return S;
    R.set_size(0);
    for (unsigned int k = 0;k < n;++k)
      if (bit(i,k) == 0) R[R.size()] = Q[k];
    T = cbextend(S,prod(R));
    R.set_size(0);
    for (unsigned int k = 0;k < n;++k)
      if (bit(i,k) == 1) R[R.size()] = Q[k];
    S = cbextend(T,prod(R));
    ++i;
  }
}

/* cb S */
base_vector<bigint> cb(base_vector<bigint> S)
{
  /* Algorithm 18.1 */
  base_vector<bigint> result(0,EXPAND);
  unsigned int n = S.size();
  if (n == 0) return result;
  if (n == 1) {
    if (S[0] != 1) result[0] = S[0];
    return result;
  }
  base_vector<bigint> T(0,EXPAND);
  T.assign(0,S,0,n/2 - 1);
  base_vector<bigint> P = cb(T);
  T.assign(0,S,n/2,n - 1);
  base_vector<bigint> Q = cb(T);
  return cbmerge(P,Q);
}

/* (i,pai) = reduce(p,a) */
void reduce(bigint &i,bigint &pai,bigint p,bigint a)
{
  /* Algorithm 19.2 */
  if (!!(a % p)) { i = 0; pai = a; return; }
  bigint j;
  bigint b;
  reduce(j,b,p * p,a / p);
  if (!(b % p)) { i = 2 * j + 2; pai = b / p; return; }
  i = 2 * j + 1; pai = b;
}

int printfactors(bigint a,base_vector<bigint> P)
{
  /* Algorithm 20.1 */
  unsigned int n = P.size();
  if (!n) {
    if (a != 1) { cout << "fail"; return 0; }
    return 1;
  }
  if (n == 1) {
    bigint m;
    bigint c;
    reduce(m,c,P[0],a);
    if (c != 1) { cout << "fail"; return 0; }
    cout << P[0] << "^" << m << " ";
    return 1;
  }
  base_vector<bigint> Q(0,EXPAND);
  Q.assign(0,P,0,n/2 - 1);
  bigint y = prod(Q);
  bigint b;
  bigint c;
  ppippo(b,c,a,y);
  if (!printfactors(b,Q)) return 0;
  Q.assign(0,P,n/2,n - 1);
  if (!printfactors(c,Q)) return 0;
  return 1;
}

void printfactors(base_vector<bigint> S,base_vector<bigint> P)
{
  /* Algorithm 21.2 */
  unsigned int n = S.size();
  if (!n) return;
  bigint x = prod(P);
  bigint y = prod(S);
  bigint z = ppi(x,y);
  base_vector<bigint> D = split(z,P);
  base_vector<bigint> Q(0,EXPAND);
  for (unsigned int i = 0;i < P.size();++i)
    if (D[i] == P[i]) Q[Q.size()] = P[i];
  if (n == 1) {
    cout << y << " = ";
    printfactors(y,Q);
    cout << "\n";
    return;
  }
  base_vector<bigint> T(0,EXPAND);
  T.assign(0,S,0,n/2 - 1);
  printfactors(T,Q);
  T.assign(0,S,n/2,n - 1);
  printfactors(T,Q);
}

int main()
{
  base_vector<bigint> S;
  base_vector<bigint> P;

  cin >> S;
  P = cb(S);
  cout << "S = " << S << "\n";
  cout << "cb S = " << P << "\n";
  printfactors(S,P);
  return 0;
}