#!/bin/python3
import sys
import random
import gmpy2

PRODUCT_PRIMES_UNDER_100 = 2305567963945518424753102147331756070

def gcd(a, b):
	a, b = (max(a, b), min(a, b))

	while b != 0:
		a, b = (b, a % b)

	return a

# Modern implementations use x^2 + 1 rather than x^2 - 1 like in the original algorithm
def g(x, c, n):
	return (x * x + c) % n

def pollards_rho(n, starting_x, c):
	x_i = g(starting_x, c, n)
	x_j = g(g(starting_x, c, n), c, n)

	d = gcd(x_i - x_j, n)

	while (d <= 1):
		x_i = g(x_i, c, n)
		x_j = g(g(x_j, c, n), c, n)

		d = gcd(x_i - x_j, n)

	# At this point, we have now found an x_i and x_j that have created a cycle

	if d == 1 or d == n:
		return pollards_rho(n, starting_x, c + 1)

	return d

# The following is an implementation of the Miller-Rabin primality test with 40 iterations
def miller_rabin(n: int) -> int:
	# Repeat the test 40 times, thus the maximum error rate is 4^(-40) (basically negligible)
	for _ in range(40):
		a = random.randint(2, n - 2)

		exponent = n - 1

		while True:
			# If "exponent" is odd, we've reached the end of the factorization, so we may now check the last two factors for proper divisibility.  If they work, "n" passed the test for "a".  If they don't, we can be 100% positive that "n" is composite).
			if exponent & 1:
				if gmpy2.powmod(a, exponent, n) == 1 or gmpy2.powmod(a, exponent, n) == n - 1:
					break
				else:
					return False
			elif gmpy2.powmod(a, exponent, n) == n - 1:
				break

			# This simply divides the exponent by 2
			exponent >>= 1

	return True

def is_prime(n: int) -> int:
	# Basic test to eliminate trivial numbers and numbers that are easily not prime
	if 1 < n and n < 4:
		return True
	elif not n & 1 or n == 1 or (gcd(PRODUCT_PRIMES_UNDER_100, n) != 1 and n > 97):
		return False

	# At this point, use the Miller-Rabin test
	return miller_rabin(n)

# The final procedure, combining everything together
def factor(n: int) -> list[int]:
	factors = []

	def factor_iter(n):
		if is_prime(n):
			factors.append(n)
			return

		if gmpy2.is_square(n):
			f = gmpy2.isqrt(n)
		elif (d := gcd(PRODUCT_PRIMES_UNDER_100, n)) > 1 and d != n:
			f = n // d
		else:
			f = pollards_rho(n, 2, 1)

		# Using // is okay here (we know that everything is divisible), and it is necessary so as to prevent Python from casting the two integers to floats
		f_ = n // f

		factor_iter(f)
		factor_iter(f_)

	factor_iter(n)

	factors.sort()
	return factors

number_to_factor = int(sys.argv[1]) if len(sys.argv) > 1 else int(input("n: "))

print(f"{number_to_factor}:", end='')
for f in factor(number_to_factor):
	print(f" {f}", end='')
print()

"""for n in range(2, 100000):
	factor(n)"""