#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
#include <cmath>


#include "gmp.h"
#include <gmpxx.h>

using namespace std;

struct factorBaseElement
{
	double logp;
	int prime;
	int start1, start2;
};

const int EXTRA = 5;
int factorBaseSize = 8;
int interValSize = 300;

vector<factorBaseElement> factorBase(factorBaseSize);
vector<int> interval(interValSize);
vector<int> intervalOrig(interValSize);


vector<double> intervalLog(interValSize, 0);

mpz_t matrix[8+5];
mpz_t tal[8+5];
mpz_t talOrig[8+5];


mpz_t N;
double THRESHOLD;
int STARVALUE;


void shanksTonelliAlg(mpz_t res, mpz_t prime, mpz_t n);
bool trialDivision2(mpz_t res, mpz_t tal);
void trialDivision(int tal);


void init()
{
	mpz_array_init(matrix[0], factorBaseSize + EXTRA, factorBaseSize*2 + EXTRA);
	mpz_array_init(tal[0], factorBaseSize + EXTRA, factorBaseSize*2 + EXTRA);
	mpz_array_init(talOrig[0], factorBaseSize + EXTRA, factorBaseSize*2 + EXTRA);


	mpz_t n, tmp, tmp2, tmp3;
	mpz_init(tmp); mpz_init(tmp2); mpz_init(tmp3);

	mpz_sqrtrem(tmp, tmp2, N);
	STARVALUE = mpz_get_ui(tmp);
	if(mpz_cmp_ui(tmp2, 0) != 0)
		STARVALUE++;
	cout << "STARVALUE: " << STARVALUE << endl;

	mpz_sqrt(tmp, N);
	THRESHOLD = log(mpz_get_ui(N))/log(2.0) - log(100);
	cout << "THRESHOLD: " << THRESHOLD << endl;


	///////////////////////////////// init interval values
	mpz_set_ui(tmp, STARVALUE);
	for(int i = 0; i < interValSize; i++)
	{

		mpz_add_ui(tmp2, tmp, i);
		intervalOrig[i] = mpz_get_ui(tmp2);
		mpz_mul(tmp2, tmp2, tmp2);
		mpz_sub(tmp2, tmp2, N);
		interval[i] = mpz_get_ui(tmp2);;
	}
	//////////////////////////////////////

}

void initFactorBase()
{
	mpz_t p;
	mpz_init_set_ui(p, 3);
	factorBase[0].prime = 2;
	factorBase[0].logp = log(2.0)/log(2.0);
	factorBase[0].start1 = factorBase[0].start2 = 0;

	for(int i = 1; i < factorBaseSize; mpz_nextprime(p, p))
	{
		if(mpz_legendre(N, p) != 1)
			continue;
		int prime = mpz_get_ui(p);
		factorBase[i].prime = prime;
		factorBase[i].logp = log(prime)/log(2.0);

		mpz_t res; mpz_init(res);
		shanksTonelliAlg(res, p, N);

		int start1 = mpz_get_ui(res);
		int start2 = prime - start1;
		start1 = (int)ceil((STARVALUE - start1)/(double)prime)*prime + start1;
		start2 = (int)ceil((STARVALUE - start2)/(double)prime)*prime + start2;

		factorBase[i].start1 = start1 - STARVALUE;
		factorBase[i].start2 = start2 - STARVALUE;
		i++;
	}
}

void shanksTonelliAlg(mpz_t res, mpz_t prime, mpz_t n)
{
	if (mpz_legendre(n, prime) == -1)
		return;
	int s;
	mpz_t tmp, Q, W, R, V;
	mpz_init(tmp); mpz_sub_ui(tmp, prime, 1);
	mpz_init(Q); mpz_init(W); mpz_init(R); mpz_init(V);

	s = mpz_scan1 (tmp, 0); // p-1 = Q2^s
	mpz_div_2exp(Q, tmp, s);


	while(true)
	{
		mpz_random(W, 2);
		mpz_mod(W, W, prime);
		if(mpz_legendre(W, prime) == -1)
		{
			mpz_set_ui(res, 0);
			break;
		}
	}


	mpz_add_ui(tmp, Q, 1);
	mpz_div_2exp(tmp, tmp, 1);
	mpz_powm(R, n, tmp, prime); // X = n^(Q-1)/2 mod p

	mpz_powm(V, W, Q, prime); // Y = w^Q mod p

	mpz_t inv, RR;
	mpz_init_set(inv, n); mpz_init(RR);
	mpz_invert(inv, inv, prime);

	while(true)
	{
		mpz_mul(RR, R, R);
		mpz_mul(RR, RR, inv);
		int i;
		for(i = 0; i < s; i++)
		{
			mpz_mod(RR, RR, prime);
			if(mpz_cmp_ui(RR, 1) == 0)
				break;
			mpz_mul(RR, RR, RR);
		}
		if(i == 0)
			break;

		mpz_set_ui(tmp, 1);
		mpz_mul_2exp(tmp, tmp, s-i-1);

		mpz_powm(tmp, V, tmp, prime);
		mpz_mul(R, R, tmp);
		mpz_mod(R, R, prime);
	}
	mpz_set(res, R);
}

void sieving()
{
	mpz_t tmp, tmp2;
	mpz_init(tmp); mpz_init(tmp2);
	for(int i = 1; i < factorBaseSize; i++)
	{
		for(int j = 0; j < 2; j++)
		{
			int index = (j == 0)? factorBase[i].start1: factorBase[i].start2;
			for(; index < interValSize; index += factorBase[i].prime)
			{
				int tmp = interval[index];
				while(tmp % factorBase[i].prime == 0)
				{
					tmp /= factorBase[i].prime;
					intervalLog[index] += factorBase[i].logp;
				}
			}

		}
	}
}

void gauss_eliminate()
{
	int size = factorBaseSize+EXTRA;
	bool found;
	/* Initialize the identity matrix */
	for(int i = 0; i < size; i++)
		mpz_setbit(matrix[i], factorBaseSize+i);
	int i;
	for(i = 0; i < size; i++)
	{
		/* Does this row contain only zeros */
		if(mpz_cmp_ui(matrix[i], 0) == 0)
		{
			/* If there is an all zero row below this one then swap the rows */
			found = false;
			for(int j = i+1; j < size; j++)
			{
				if (mpz_cmp_ui(matrix[j], 0) != 0)
				{
					found = true;
					break;
				}
				if(found == false)
					return;
				else
					mpz_swap(matrix[i], matrix[j]);
			}
		}
		/* Find the index of the first non-zero entry in this row and if there is a non-zero entry in the matrix below and to the left then swap */
		found = true;
		int j, l;
		while(found)
		{
			found = false;
			j = mpz_scan1(matrix[i], 0);
			if(j > 0)
			{
				for (int k = i+1; k < size; k++)
				{
					l = mpz_scan1(matrix[k], 0);
					if (l < j)
					{
						mpz_swap(matrix[i], matrix[k]);
						found = 1;
						break;
					}
				}
			}
		}
		/* Change every other rows j'th entry to 0 */
		for (int k = 0; k < size; k++)
		{
			if (i != k)
			{
				l = mpz_scan1(matrix[k], 0);
				if (l == j)
					mpz_xor(matrix[k], matrix[i], matrix[k]);
			}
		}
	}
}
void createMatrix()
{
	mpz_t res; mpz_init(res);
	mpz_t tmp; mpz_init(tmp);
	mpz_t tmp2; mpz_init(tmp2);
	int cols = 0;
	for(int i = 0; i < interValSize; i++)
	{
		if(intervalLog[i] >= THRESHOLD)
		{
			mpz_set_ui(tmp, interval[i]);

			if(trialDivision2(res, tmp))
			{
				mpz_set_ui(tmp2, interval[i]);
				mpz_set(tal[cols], tmp2);

				mpz_set_ui(tmp2, intervalOrig[i]);
				mpz_set(talOrig[cols], tmp2);

				mpz_set(matrix[cols], res);
				cols++;
				if(cols >= factorBaseSize+EXTRA)
				{
					gauss_eliminate();
					return;
				}
			}
		}
	}
}

void findFactors()
{


	int first, retries, success;
	mpz_t x, y;
	mpz_t tmp, tmp2;
	mpz_t q, p, gcd;

	mpz_init(x); mpz_init(y); mpz_init(tmp); mpz_init(tmp2);
	mpz_init(q); mpz_init(p); mpz_init(gcd);
	mpz_set_ui(y, 0);
	for(int i = 1; i < factorBaseSize; i++)
		mpz_setbit(y, i);

	first = retries = success = 0;
	for (int i = factorBaseSize + EXTRA - 1; i > 0; i--)
	{
		mpz_and(x, matrix[i], y);
		if (mpz_cmp_ui(x, 0) != 0)
		{
			first = i + 1;
			break;
		}
	}
	/* Loop through all possible solutions */
	for (int i = first; i < factorBaseSize + EXTRA; i++)
	{
		mpz_set_ui(y, 1);
		mpz_set_ui(x, 1);
		for (int j = 0; j < factorBaseSize + EXTRA; j++)
		{
			if (mpz_tstbit(matrix[i], factorBaseSize + j))
			{
				cout << "j: "<< j << endl;
				mpz_set(tmp, tal[j]);
				mpz_set(tmp2, talOrig[j]);
				cout << "tmp: " << tmp << endl;
				mpz_mul(y, y, tmp);
				mpz_mul(x, x, tmp2);
			}
		}
		cout << "y: " << y << endl;
		mpz_sqrt(y, y);
		mpz_sub(gcd, y, x);
		mpz_gcd(p, gcd, N);

		cout << "y: " << y << endl;
		if (mpz_cmp_ui(p, 1) == 0 || mpz_cmp(p, N) == 0)
			/* We got the trivial solution. Let's try again */
			retries++;
		else {
			/* We found the non-trivial factorization */
			mpz_add(gcd, y, x);
			mpz_gcd(q, gcd, N);
			success = 1;
			break;
		}
	}

	if(!success)
	{
		printf("Unable to find a factorization of ");
		cout << N << endl;
	}
	else
		cout << N << " = " << q << " * " << p << endl;
}

void trialDivision(int tal)
{
	for(int i = 2; i <= tal; i++)
	{
		while(tal%i == 0)
		{
			cout << i << " ";
			tal /= i;
		}
	}
	cout << endl;
}

bool trialDivision2(mpz_t res, mpz_t tal)
{
	mpz_t tmp; mpz_init(tmp);
	mpz_set_ui(res, 0);

	for(int i = 0; i < factorBase.size(); i++)
	{
		while(mpz_cdiv_q_ui(tmp, tal, factorBase[i].prime) == 0)
		{
			//cout << factorBase[i].prime << " ";
			mpz_set(tal, tmp);
			mpz_combit(res, i);
		}
		if(mpz_cmp_ui(tal, 1) == 0)
		{
			//cout << endl;
			return true;
		}
	}
	//cout << endl;
	return false;
}


int main()
{
	mpz_t tmp; mpz_init(tmp);
	mpz_t tmp2; mpz_init(tmp2);
	gmp_scanf("%Zd", N);

	init();
	initFactorBase();
	cout << "primes: " << endl;
	for(int i = 0; i < factorBase.size(); i++)
		cout << factorBase[i].prime << " ";
	cout << endl;
	sieving();
	createMatrix();
	findFactors();


	return 0;
}