From c88209fdaa55213fe8a036ef6bd3f6fc78e680bf Mon Sep 17 00:00:00 2001
From: Reiner Herrmann <reiner@reiner-h.de>
Date: Mon, 21 Jun 2010 15:31:27 +0200
Subject: pe solution 55

---
 src/projecteuler/055.py    | 17 ++++++++++
 src/projecteuler/common.py | 77 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 94 insertions(+)
 create mode 100644 src/projecteuler/055.py

diff --git a/src/projecteuler/055.py b/src/projecteuler/055.py
new file mode 100644
index 0000000..d623d9b
--- /dev/null
+++ b/src/projecteuler/055.py
@@ -0,0 +1,17 @@
+
+def lychrel(x):
+	for i in xrange(50):
+		t = x + int(str(x)[::-1])
+		if t == int(str(t)[::-1]):
+			return False
+		x = t
+	return True
+
+count = 0
+
+for i in xrange(1, 10000):
+	if lychrel(i):
+		count += 1
+
+print count
+
diff --git a/src/projecteuler/common.py b/src/projecteuler/common.py
index c85af4b..8bc4953 100644
--- a/src/projecteuler/common.py
+++ b/src/projecteuler/common.py
@@ -1,3 +1,4 @@
+import math
 
 class sieve:
 	def __init__(self, limit):
@@ -54,3 +55,79 @@ class pandigital:
 			result |= new
 		return result
 
+def ggt(a, b): # Stein's algorithm
+	k = 0
+	t = 0
+	while a&1==0 and b&1==0:
+		a = a>>1
+		b = b>>1
+		k += 1
+	if a&1==1:
+		t = -b
+	else:
+		t = a
+	while t != 0:
+		while t&1==0:
+			t = t>>1
+		if t>0:
+			a = t
+		else:
+			b = -t
+		t = a - b
+	return a*(1<<k)
+
+def composite(n): # check if n = c^b
+	primes = sieve(int(math.sqrt(n))+1).primes()
+	for c in primes:
+		number = c
+		while number < n:
+			number *= c
+		if number == n:
+			return True
+	return False
+
+def order(n, r):
+	res = n % r
+	k = 1
+	while res != 1:
+		res = (res * n) % r
+		k += 1
+	return k
+
+def phi(x):
+	primes = sieve(x+1).primes()
+	product = x
+	for p in primes:
+		if x % p != 0:
+			continue
+		product *= (1 - 1.0/p)
+	return int(product)
+		
+
+def prime_test(n): # aks test
+	if composite(n):
+		return False
+
+	logn = math.log(n, 2)
+	o = 0
+	r = 0
+	while o <= logn:
+		r += 1
+		o = order(n, r)
+
+	a = 1
+	while a <= r:
+		g = ggt(a, n)
+		if g < n and g > 1:
+			return False
+		a += 1
+
+	if n <= r*r:
+		return True
+
+	for a in xrange(1, math.floor(math.sqrt(phi(r))*logn)):
+		# TODO if math.pow(x+a, n) != 
+	
+	return True
+
+
-- 
cgit v1.2.3