#!/usr/bin/python

'''
Pyecm factors large integers (up to 50 digits) using the Elliptic Curve Method (ECM), a fast factoring algorithm.

Copyright (C) 2006, 2007 Martin Kelly

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
'''

'''
We are using curves in the form y^2 = x^3 - 3x + (a^2 - 2)

Points are in Jacobian coordinates, unless they are immediate results to calls of the affine or small_multiples functions. In both of those cases, they are in affine coodinates.

The double and add functions have their first point in affine coorinates if the have _affine after their name.
'''

import math
import sys

PRODUCT = 111546435 # 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23
SMALLEST_COUNTEREXAMPLE_FASTPRIME = 2047
VERSION = '1.2.2'

try: # Use psyco if available
	import psyco
	psyco.full()
except ImportError:
	pass
try:# Use gmpy if available (otherwise define our own functions to replace its functionality)
	from gmpy import gcd, invert, mpz
except ImportError:
	def gcd(a, b):
		'''gcd(a, b) --> GCD of a and b

Computes the Greatest Common Divisor of a and b using the standard quadratic time improvement to the Euclidean Algorithm.

Returns the GCD of a and b.'''

		if b == 0:
			return a
		elif a == 0:
			return b

		count = 0

		if a < 0:
			a = -a
		if b < 0:
			b = -b

		while not ((a & 1) | (b & 1)):
			count += 1
			a >>= 1
			b >>= 1

		while not a & 1:
			a >>= 1

		while not b & 1:
			b >>= 1

		if b > a:
			b,a = a,b

		while b != 0 and a != b:
			a -= b
			while not (a & 1):
				a >>= 1

			if b > a:
				b, a = a, b

		return a << count

	def invert(a, b):
		'''invert(a, b) --> inverse of a (mod b)

Computes the inverse of a modulo b. b must be odd.

Returns the inverse of a (mod b).'''

		if a == 0 or b == 0:
			return 0

		truth = False
		if a < 0:
			t = True
			a = -a

		b_orig = b
		alpha = 1
		beta = 0

		while not a & 1:
			if alpha & 1:
				alpha += b_orig
			alpha >>= 1
			a >>= 1

		if b > a:
			a, b = b, a
			alpha, beta = beta, alpha

		while b != 0 and a != b:
			a -= b
			alpha -= beta

			while not a & 1:
				if alpha & 1:
					alpha += b_orig
				alpha >>= 1
				a >>= 1

			if b > a:
				a,b = b,a
				alpha, beta = beta, alpha

		if a == b:
			a -= b
			alpha -= beta
			a, b = b, a
			alpha, beta = beta, alpha

		if a != 1:
			return 0

		if truth:
			alpha = b_orig - alpha

		return alpha

	def mpz(n):
		return n

M_ONE = mpz(-1) # Make sure everything works with or without gmpy
ZERO = mpz(0)
ONE = mpz(1)

def add(p1, p2,  n):
	'''add(p1, p2,  n) --> list of addition result

Adds two points on an elliptic curve.

Returns a list of the addition result.'''

	u1, v1, w1 = p1[0], p1[1], p1[2]
	u2, v2, w2 = p2[0], p2[1], p2[2]

	w2_sq = (w2 * w2) % n
	w1_sq = (w1 * w1) % n
	w1_cube = (w1 * w1_sq) % n
	w2_cube = (w2 * w2_sq) % n
	x1 = (u1 * w2_sq) % n
	x2 = (u2 * w1_sq) % n
	y1 = (v1 * w2_cube) % n
	y2 = (v2 * w1_cube) % n

	u_diff_sq = ((x2 - x1) * (x2 - x1)) % n

	if x1 == x2:
		return double(p1,  n)

	u_diff_cube = ((x2 - x1) * u_diff_sq) % n
	prod = (x1 * u_diff_sq) % n
	u3 = ((y2 - y1)*(y2 - y1) - u_diff_cube - (prod << 1)) % n
	v3 = (y2 - y1) * (prod - u3) - y1 * u_diff_cube
	w3 = ((w1 * w2) % n) * (x2 - x1)
	v3, w3 = v3 % n, w3 % n

	return [u3, v3, w3]

def add_affine(p1, p2, n):
	'''add(p1, p2,  n) --> list of addition result

Adds two points on an elliptic curve. p1 is in Affine coordinates.

Returns a list of the addition result.'''

	u1, v1 = p1[0], p1[1]
	u2, v2, w2 = p2[0], p2[1], p2[2]

	w2_sq = (w2 * w2) % n
	w2_cube = (w2 * w2_sq) % n
	x1 = (u1 * w2_sq) % n
	y1 = (v1 * w2_cube) % n

	u_diff_sq = ((u2 - x1) * (u2 - x1)) % n

	if x1 == u2:
		return double_affine(p1,  n)

	u_diff_cube = ((u2 - x1) * u_diff_sq) % n
	prod = (x1 * u_diff_sq) % n
	u3 = ((v2 - y1)*(v2 - y1) - u_diff_cube - (prod << 1)) % n
	v3 = (v2 - y1) * (prod - u3) - y1 * u_diff_cube
	w3 = w2 * (u2 - x1)
	v3, w3 = v3 % n, w3 % n

	return [u3, v3, w3]

def atdn(a, d, n):
	'''atdn(a, d, n) --> result

Calculates a to the dth modulus n. Required sub-algorithm for ECM.

Returns the calculation's result.'''

	x = 1
	while d > 0:
		e = d & 1
		if e == 1:
			x = (a * x) % n

		a = (a * a) % n
		d >>= 1

	return x

def double(p1,  n):
	'''double(p1,  n) --> list of addition result

This is a specialized version of add, used when we know the two points are the same; therefore, we can apply some optimizations and remove some checks.

Returns a list of the addition result'''

	u1, v1, w1 = p1[0], p1[1], p1[2]
	w1_sq = (w1 * w1) % n
	sub_expression = (u1  - w1_sq) * (u1 + w1_sq) % n
	sub_expression += sub_expression << 1
	v1_sq = (v1 * v1) % n
	other_sub = (u1 * v1_sq << 2) % n

	u3 = (sub_expression * sub_expression - (other_sub << 1)) % n
	v3 = sub_expression * (other_sub - u3) - (v1_sq * v1_sq << 3)
	w3 = ((v1 * w1) << 1) % n

	v3 %= n

	return [u3, v3, w3]

def double_affine(p1, n):
	'''double_affine(p1, n) --> list of addition result

This is a specialized version of add, used when we know the two points are the same. Therefore, we can apply some optimizations and remove some checks. double_affine is specifically for Affine coordinates.

Returns a list of the addition result.'''

	u1, v1 = p1[0], p1[1]
	sub_expression = (u1 * u1 - 1) % n
	sub_expression += sub_expression << 1
	v1_sq = (v1 * v1) % n
	other_sub = (u1 * v1_sq << 2) % n

	u3 = (sub_expression * sub_expression - (other_sub << 1)) % n
	v3 = sub_expression * (other_sub - u3) - (v1_sq * v1_sq << 3)
	w3 = (v1 << 1)

	if w3 > n:
		w3 -= n

	v3 %= n

	return [u3, v3, w3]

def next_good_u(old_u, log_n):
	'''next_good_u(old_u, log_n) --> next u > old_u, next_u is prime

Computes the next u that is greater than old_u and prime.

Returns the next u that that is greater than old_u and prime.'''

	u = 1
	while u < old_u:
		for _ in xrange(log_n):
			u += 2
			while not isprime(u):
				u += 2
	return u

def prod(seq):
	'''prod(seq) --> seq's product

Calculates the product of the sequence by breaking it into smaller multiplications.

Returns the product of seq.'''

	jump = 1

	while jump < len(seq):
		for i in xrange(0, len(seq) - jump, jump << 1):
			seq[i] *= seq[i+jump]

		jump <<= 1

	return seq[0]

def small_multiples(p1, n):
	'''small_multiples --> (3 * p1, 5 * p1)

Computes three times p1 and five times p1 in Affine coordinates.

Returns (3 * pi, 5 * pi) in Affine coordinates.'''

	p2 = double_affine(p1, n)
	p3 = to_affine(add_affine(p1, p2, n), n)
	p5 = to_affine(add_affine(p3, p2, n), n)

	return (p3, p5)

def to_affine(p, n):
	'''to_affine(p, n) --> p in Affine coordinates

Computes in Affine coordinates. p is given in Jacobian coordinates. n is the number we are factoring.

Returns in Affine coordinates. p is given in Jacobian coordinates.'''

	u, v, w = p[0], p[1], p[2]
	w_inv = invert(w, n)
	w_inv_sq = (w_inv * w_inv) % n

	affine_coords = [(u * w_inv_sq) % n, (((v * w_inv_sq) % n) * w_inv) % n]

	return affine_coords

def ecm(n):
	'''ecm(n) --> factors (yielded via generator)

Factors an integer n using the Elliptic Curve Method (ECM).

Generates factors of n.'''

	def mainloop(n, b):
		if isprime(n):
			yield n
			return

		points = [] # A collection of points to avoid recalculation of primes
		points_new = []
		log_n = math.log(n)
		log_log_n = math.log(log_n)
		u = 4 + int(math.e ** math.sqrt (.25 * log_n * (log_log_n  - math.log(log_log_n))))
		log_n = int(log_n)
		u = next_good_u(u, log_n)
		u_sqrt = int(math.sqrt(u)) + 1
		log_u = math.log(u)
		u2 = int((94 * math.e / (99 * math.log(2) * (math.e + 1))) * u * log_u / math.log(log_u))
		bound = 1 + int(log_u**3) / log_n
		gcd_block = 1 + log_n / bound
		a = ZERO
		f = ONE
		cutoff = int(math.log(u) * math.log(math.log(u)))
		p = []

		while f in (ONE, n):
			a += ONE
			b += 1
			if b & 1:
				if not n % b:
					while not n % b:
						yield b
						n /= b

					del p
					del points
					del points_new
					if n == 1:
						return
					for factor in mainloop(n, b):
						yield factor
					return

			p1 = [M_ONE, a, ONE]
			prime = 3

			for _ in xrange(int(math.log(u, 2))):
				p1 = double(p1,  n)

			f = gcd(p1[-1], n)

			while prime <= u_sqrt and f == ONE:
				p = []
				for _ in xrange(log_n):
					while not isprime(prime):
						prime += 2

					p.append(mpz(prime ** max(1, int(math.log(u, prime)))))
					prime += 2

				w = p1[-1]
				p1 = to_affine(p1, n)

				if p1 == [0, 0]:
					f = gcd(w, n)
				else:
					p3, p5 = small_multiples(p1, n)
					p1 = multiply(p1, p3, p5, prod(p),  n)

			while prime < u and f == ONE:
				p = []
				for _ in xrange(log_n):
					while not isprime(prime):
						prime += 2

					p.append(prime)
					prime += 2

				w = p1[-1]
				p1 = to_affine(p1, n)

				if p1 == [0, 0]:
					f = gcd(w, n)
				else:
					p3, p5 = small_multiples(p1, n)
					p1 = multiply(p1, p3, p5, prod(p),  n)

			p = []

			if f == ONE:
				f = gcd(p1[-1], n)

			if f == ONE:
				p1_ = p1[:]

				while not isprime(prime):
					prime += 2
				w = p1[-1]
				p1 = to_affine(p1, n)
				if p1 == [0, 0]:
					f = gcd(w, n)
				else:
					p3, p5 = small_multiples(p1, n)
					p1 = multiply(p1, p3, p5, prime, n)

					f = gcd(p1[-1], n)
			if f == ONE:
				p1_ = to_affine(p1_, n)
				p3, p5 = small_multiples(p1_, n)
				points.append({1:p1_, 3:p3, 5:p5})
				points_new.append(p1[:])
			if len(points) >= bound and f in (ONE, n): # Stage 2, log_n curves at a time.
				old_prime = prime
				while prime < u2:
					product = ONE
					for _ in xrange(gcd_block):
						prime += 2
						while not isprime(prime):
							prime += 2
						for i in xrange(bound):
							already_done = points[i]
							to_do = prime - old_prime

							if to_do in already_done:
								to_add = already_done[to_do]
							else:
								p1 = already_done[1]
								p3 = already_done[3]
								p5 = already_done[5]

								to_add = multiply(p1, p3, p5, to_do, n)

								if to_do < cutoff:
									to_add = to_affine(to_add, n)
									already_done[to_do] = to_add

							if len(to_add) == 2:
								points_new[i] = add_affine(to_add, points_new[i], n)
							else:
								points_new[i] = add(to_add, points_new[i], n)

							product *= points_new[i][-1]
							product %= n

						old_prime = prime
					f = gcd(n, product)

					if f not in (ONE, n):
						break

				points = []
				points_new = []

		del p	# Get rid of all of the things that use significant memory before recursion
		del points
		del points_new

		for factor in mainloop(f, b):
			yield factor
		for factor in mainloop(n / f, b):
			yield factor

	if not isinstance(n, int) and not isinstance(n, long):
		raise ValueError, 'Number given must be integer or long.'

	if n == 0:
		yield '0 does not have a well-defined factorization.'
		return
	elif n < 0:
		yield -1
		n = -n

	if n > 840: # The formulas break for n < a constant
		log_n = math.log(n)
		log_log_n = math.log(log_n)

		u = int(math.e ** math.sqrt (.25 * log_n * (log_log_n - math.log(log_log_n)))) + 4
	else:
		u = 6 # 6 is arbitrary. The value must be in [3, 23] for this to work, however

	trial_division_bound = max(23, int(u / (math.log(u)**2)))

	while not n & 1:
		n >>= 1
		yield 2

	prime = 3
	while prime < trial_division_bound: # Trial division
		while not n % prime:
			n /= prime
			yield prime
		prime += 2

	if n == 1:
		return

	n = mpz(n)

	for factor in mainloop(n, trial_division_bound - 1):
		yield factor

def fastprime(n):
	'''fastprime(n) --> True or False

Tests for primality of n using an algorithm that is very fast, O(N**3 / log(N)) (assuming quadratic multiplication) where n has N digits, but ocasionally inaccurate for n >= 2047.

Returns the primality of n.'''

	if n <= 28:
		if n in (2, 3, 5, 7, 11, 13, 17, 19, 23):
			return True
		else:
			return False

	if not n & 1:
		return False

	m = int(n % PRODUCT ) # Make sure m is not a long (longs are slow)
	for i in (3, 5, 7, 11, 13, 17, 19, 23):
		if not m % i:
			return False
	if n < 841: # 29**2
		return True

	for i in xrange(4, int(math.log(n)/(math.log(n) ** 2))): # Trial divide up to 6 * int(math.log(n)/(math.log(n) ** 2)) - 1
		p1 = (i << 1) + (i << 2) + 5 # 6i + 5
		p2 = p1 + 2
		prod = p1 * p2
		m = n % prod

		if not m % p1:
			return False
		if not m % p2:
			return False

		if n < (prod + (p1 << 1) + (p1 << 3) + 36): # If n < (p1 + 6)**2
			return True

	j = 1
	d = n >> 1

	while not d & 1:
		d >>= 1
		j += 1

	count = 0

	g = 0	# The next few lines do a binary reversal of d
	while d:
		g <<= 1
		g += d & 1
		d >>= 1

	p = 1
	while g:
		p = p * p % n
		p <<= g & 1 # Multiply p by 2 if g is odd
		g >>= 1
	# Now p = atdn(a, d, n). Because a = 2, we can do some tricks to speed this up.
	if p > n:
		p -= n

	if p == 1 or p == n - 1:
		return True

	while count < j - 1:
		p = p * p % n
		count += 1

		if p == 1:
			return False
		elif p == n - 1:
			return True

	return False

def isprime(n):
	''' isprime(n) --> True or False

Tests for primality of n trying first fastprime and then a slower but accurate algorithm. Time complexity is O(N**3) (assuming quadratic multiplication), where n has N digits.

Returns the primality of n.
'''
	if not fastprime(n):
		return False
	elif n < SMALLEST_COUNTEREXAMPLE_FASTPRIME:
		return True

	j = 1
	d = n >> 1

	while not d & 1:
		d >>= 1
		j += 1

	for a in xrange(3, int(0.75 * math.log(math.log(n)) * math.log(n))):
		if not fastprime(a):
			continue

		count = 0
		do_loop = 0
		p = atdn(a,d,n)

		if p == 1 or p == n - 1:
			continue

		while count < j:
			p = (p * p) % n
			count += 1

			if p == 1:
				return False
			elif p == n - 1:
				do_loop = 1
				break

		if do_loop:
			continue

		return False

	return True

def multiply(p1, p3, p5, d,  n):
	'''multiply(p1, p3, p5, d,  n) --> p1 multiplied by d

Multiplies p1 by d. p3 and p5 are the triples and quintuples of the point p1, in Affine coordinates, as to avoid recomputation.

Returns p1 multiplied by d.
'''

	g = 0
	count = 0
	while not d & 1:
		d >>= 1
		count += 1

	while d:	# This loop constructs the non-adjacent form of d
		g <<= 2
		g ^= d & 3
		d += (d & 2) >> 1
		d >>= 1

	p = p1[:] + [1]
	g >>= 2

	mp1 = [p1[0], -p1[1]]
	mp3 = [p3[0], -p3[1]]
	mp5 = [p5[0], -p5[1]]

	while g > 0:
		p = double(p, n)
		if g & 1:
			if g & 16:
				p = double(p, n)
				p = double(p, n)

				if (g & 32) == ((g & 2) << 4):
					if g & 2:
						p = add_affine(mp5, p, n)
					else:
						p = add_affine(p5, p,  n)

				else:
					if g & 2:
						p = add_affine(mp3, p, n)
					else:
						p = add_affine(p3, p,  n)

				g >>= 4

			else:
				if g & 2:
					p = add_affine(mp1, p, n)
				else:
					p = add_affine(p1, p,  n)

		g >>= 2

	for _ in xrange(count):
		p = double(p, n)

	return p

def help():
	print	'''\
Usage: pyecm [factors]...
Factor numbers using the Elliptic Curve Method.

With no factors given via command-line, read standard input.

Report bugs to Martin Kelly <aomighty@gmail.com>.'''

	sys.exit()

def command_line():
	try:
		factors = sys.argv[1:]
		try:
			factors = [int(factor) for factor in factors]
		except ValueError: # If non-int given
			help()

		first_item = True

		for factor in factors:
			if not first_item:
				print # Formatting to look right
			print	'Factoring %d:' % factor
			for factor in ecm(factor):
				print	factor

			first_item = False

	except IndexError: # Incorrect parameters given
		help()

	sys.exit()

def interactive():
	print	'pyecm v. %s (interactive mode):' % VERSION
	print	'Type "exit" at any time to quit.'
	print
	try:
		response = raw_input().lower().strip()
		while response != 'exit' and response != 'quit':
			try:
				n = int(response)
				print	'Factoring number %d:' % n
				for factor in ecm(n):
					print	factor
				print
			except (IndexError, ValueError):
				print	'Please enter an integer.'
			response = raw_input().lower().strip()
	except (EOFError, KeyboardInterrupt):
		sys.exit()

def main():
	for item in sys.argv[1:]:
		if item == '-h' or item == '--help':
			help()

	try:
		sys.argv[1]
		command_line()
	except IndexError:
		interactive()

if __name__ == '__main__':
	main()
