diff --git a/palinodrome_product.py b/palinodrome_product.py
index ac77e59..af40787 100644
--- a/palinodrome_product.py
+++ b/palinodrome_product.py
@@ -1,3 +1,67 @@
+#!/usr/bin/env python
+
 """A palindromic number reads the same both ways. The largest palindrome made
 from the product of two 2-digit numbers is 9009 = 91 99.Find the largest
 palindrome made from the product of two 3-digit numbers."""
+
+# an unadorned 'digit' term implies base-ten
+
+# http://stackoverflow.com/questions/4643647/fast-prime-factorization-module
+from primefactorization import primefactors
+
+def check(prod):
+  '''
+  check if prod is a palindrome
+  '''
+  if str(prod) == str(prod)[::-1]:
+    return True # palindrome
+  else:
+    return False # not palindrome
+
+def find(n,m,found):
+  '''
+  Find largest base-ten palindrome in the set of integers {0,...,n**2} by
+  checking products of the form (n - i)*(n - j), where 0 <= i,j <= n are
+  integers.  Recurse over the product set from largest to smallest.
+  n:     starting number
+  m:     size of domain of input sets over which to search for a palindrome in
+         the product set
+  found: collect found palindromes
+  '''
+  if not isinstance(found,set):
+    raise TypeError('found argument must be of type set')
+  for i in xrange(m):              # iterate over i
+    for j in xrange(m):            # and j
+      if check((n - i)*(n - j)):
+        found.add((n - i)*(n - j)) # success
+  if len(found) > 0:               # found at least one between n,n - m
+    return
+  elif n - m > 0:                  # only check positive input integers
+    find(n - m,m,found)            # look at next part of input sets
+  else:
+    pass
+    '''
+    a single digit positive integer is a palindrome, and not all single digit
+    integers are prime, so the algorithm will always find a palindrome of the
+    form (n - i)*(n - j) for any nonnegative integers 0 <= i,j <= n.
+    '''
+
+def main(n=999,m=100):
+  '''
+  n: the largest positive integer used to find palindromes
+  m: estimated range of numbers less than n needed to find the largest
+     palindrome, [n,n - m]
+  '''
+  f = set()           # store found palindromes here
+  find(n,m,f)         # find palindromes
+  pal = sorted(f)[-1] # largest palindrome found given params n,m
+  facts = primefactors(pal,True)
+  length = len(facts)
+  output = 'Largest palindrome in range [0,%d**2]: %d\n' % (n,pal)
+  output += 'Its prime factorization is: %s\n' % facts
+  output += 'These can be combined in 2**%d = %d different ' % (length,2**length)
+  output += 'ways to find the palindrome as a product of two integers in the '
+  output += 'range [0,%d].' % n
+  print output
+
+if __name__ == '__main__' : main()
diff --git a/prime.py b/prime.py
index b54dead..38d9511 100644
--- a/prime.py
+++ b/prime.py
@@ -1,2 +1,61 @@
+#!/usr/bin/env python
+# -*- encoding: utf-8 -*-
+
 """Write a program to find the 10 001st prime number using the sieve of
 aristothanes."""
+
+'''
+sieve of Eratosthenes (κόσκινον Ἐρατοσθένους)?
+'''
+from math import sqrt
+
+start = 2        # beginning integer of search range
+stop = 8000      # ending integer of search range
+target = 10001   # number of primes to find
+integers = []    # list of integers
+primes = []      # list of primes
+
+def remove_factors(p,integers):
+  '''
+  removes all factors of prime p from integers, where 1 < p*i < integers[-1]
+  '''
+  imin = integers[0]       # smallest integer in range
+  imax = integers[-1]      # largest integer in range
+  smallest = max(imin/p,2) # minimum integer that will have factors that land in
+                           # the provided integer list
+  largest = imax/p + 1     # largest of these integers
+  for i in xrange(smallest,largest):
+    factor = p*i
+    if factor in integers:
+      integers.remove(factor)
+    elif factor > imax:
+      break
+
+def sieve(start=start,stop=stop):
+  '''
+  Given a set of positive integers, remove all integer products of known
+  primes, since they are thus not prime.  When p is less than sqrt(stop) + 1,
+  then all the remaining integers the current range of integers are prime.
+  Recurse until the range containing p_target is processed and return that
+  prime.
+  start: starting integer
+  stop:  ending integer
+  '''
+  integers = [i for i in xrange(start,stop)]
+  if len(primes) == 0:                  # base case
+    p = start                           # initial working prime
+  else:
+    for p in primes:                    # first check primes we know
+      if p < int(sqrt(stop) + 1):
+        remove_factors(p,integers)
+    p = integers[0]                     # first integer left is now prime
+  while p < int(sqrt(stop) + 1):        # then sift through the new primes
+    remove_factors(p,integers)          # remove all factors of p from integers
+    p = integers[integers.index(p) + 1] # next integer in integers is now prime
+  primes.extend(integers)               # all remaining integers are primes
+  if len(primes) < target:              # need to look for more primes
+    sieve(stop,2*stop)                  # look at the next range of integers
+  else:                                 # found p_target
+    print('The %dst prime is %d' % (target,primes[target + 1]))
+
+if __name__ == '__main__' : sieve()
diff --git a/primefactorization.py b/primefactorization.py
new file mode 100644
index 0000000..cf78485
--- /dev/null
+++ b/primefactorization.py
@@ -0,0 +1,140 @@
+import random
+from math import sqrt
+
+# http://stackoverflow.com/questions/4643647/fast-prime-factorization-module
+def primesbelow(N):
+    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
+    #""" Input N>=6, Returns a list of primes, 2 <= p < N """
+    correction = N % 6 > 1
+    N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6]
+    sieve = [True] * (N // 3)
+    sieve[0] = False
+    for i in xrange(int(sqrt(N)) // 3 + 1):
+        if sieve[i]:
+            k = (3 * i + 1) | 1
+            sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1)
+            sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1)
+    return [2, 3] + [(3 * i + 1) | 1 for i in xrange(1, N//3 - correction) if sieve[i]]
+
+smallprimeset = set(primesbelow(100000))
+_smallprimeset = 100000
+def isprime(n, precision=7):
+    # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time
+    if n == 1 or n % 2 == 0:
+        return False
+    elif n < 1:
+        raise ValueError("Out of bounds, first argument must be > 0")
+    elif n < _smallprimeset:
+        return n in smallprimeset
+
+
+    d = n - 1
+    s = 0
+    while d % 2 == 0:
+        d //= 2
+        s += 1
+
+    for repeat in xrange(precision):
+        a = random.randrange(2, n - 2)
+        x = pow(a, d, n)
+
+        if x == 1 or x == n - 1: continue
+
+        for r in xrange(s - 1):
+            x = pow(x, 2, n)
+            if x == 1: return False
+            if x == n - 1: break
+        else: return False
+
+    return True
+
+# https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
+def pollard_brent(n):
+    if n % 2 == 0: return 2
+    if n % 3 == 0: return 3
+
+    y, c, m = random.randint(1, n-1), random.randint(1, n-1), random.randint(1, n-1)
+    g, r, q = 1, 1, 1
+    while g == 1:
+        x = y
+        for i in xrange(r):
+            y = (pow(y, 2, n) + c) % n
+
+        k = 0
+        while k < r and g==1:
+            ys = y
+            for i in xrange(min(m, r-k)):
+                y = (pow(y, 2, n) + c) % n
+                q = q * abs(x-y) % n
+            g = gcd(q, n)
+            k += m
+        r *= 2
+    if g == n:
+        while True:
+            ys = (pow(ys, 2, n) + c) % n
+            g = gcd(abs(x - ys), n)
+            if g > 1:
+                break
+
+    return g
+
+smallprimes = primesbelow(1000) # might seem low, but 1000*1000 = 1000000, so this will fully factor every composite < 1000000
+def primefactors(n, sort=False):
+    factors = []
+
+    limit = int(sqrt(n)) + 1
+    for checker in smallprimes:
+        if checker > limit: break
+        while n % checker == 0:
+            factors.append(checker)
+            n //= checker
+
+
+    if n < 2: return factors
+    else :
+        factors.extend(bigfactors(n,sort))
+        return factors
+
+def bigfactors(n, sort = False):
+    factors = []
+    while n > 1:
+        if isprime(n):
+            factors.append(n)
+            break
+        factor = pollard_brent(n)
+        factors.extend(bigfactors(factor,sort)) # recurse to factor the not necessarily prime factor returned by pollard-brent
+        n //= factor
+
+    if sort: factors.sort()
+    return factors
+
+def factorization(n):
+    factors = {}
+    for p1 in primefactors(n):
+        try:
+            factors[p1] += 1
+        except KeyError:
+            factors[p1] = 1
+    return factors
+
+totients = {}
+def totient(n):
+    if n == 0: return 1
+
+    try: return totients[n]
+    except KeyError: pass
+
+    tot = 1
+    for p, exp in factorization(n).items():
+        tot *= (p - 1)  *  p ** (exp - 1)
+
+    totients[n] = tot
+    return tot
+
+def gcd(a, b):
+    if a == b: return a
+    while b > 0: a, b = b, a % b
+    return a
+
+def lcm(a, b):
+    return abs(a * b) // gcd(a, b)
diff --git a/series_product.py b/series_product.py
index 2d6a53a..b841835 100644
--- a/series_product.py
+++ b/series_product.py
@@ -1,3 +1,50 @@
+#!/usr/bin/env python
 """Complete this program to find the greatest produdct of five consecutive
 digits in the list"""
-in_file = open('array.txt')
+
+def product(iterable):
+  '''
+  returns the product of the numbers in iterable
+  '''
+  prod = 1
+  for i in iterable:
+    prod *= i
+  return prod
+
+def result(num,array,prods,skeys):
+  '''
+  prints results
+  num:   number of consecutive digits in products
+  array: array of digits
+  prods: dictionary of consecutive products keyed by array position of first
+         digit in the product
+  skeys: list of keys of prods rsorted by product
+  '''
+  out = 'The greatest product of %d consecutive digits in the list is ' % num
+  out += '%d.  ' % prods[skeys[0]]
+  maxes = [skeys[0]] # how many times greatest product occurs
+  for i in xrange(len(skeys) - 1): # don't overrun if all products are equal
+    if prods[skeys[i]] == prods[skeys[i + 1]]:
+      maxes.append(skeys[i])
+    else:
+      break # found all of them
+  out += 'It occurs %d time%s.  ' % (len(maxes),'' if len(maxes) == 1 else 's')
+  out += 'The consecutive digits are:\n\n'
+  for m in maxes:
+    digits = tuple([array[m + i] for i in xrange(num)])
+    out += ('  %s\n' % ('%d '*num).rstrip()) % digits
+  print(out)
+
+def consecutive_product(num=5):
+  array = [] # array of input digits
+  with open('array.txt') as in_file:
+    for line in in_file:
+      digits = line.rstrip('\n') # remove trailing newline
+      array.extend([int(d) for d in digits])
+  prods = {} # computed products keyed by array position of starting digit
+  for i in xrange(len(array) - num): # iterate over possible consecutive products
+    prods[i] = product(array[i:i + num]) # each slice of num digits
+  skeys = sorted(prods,key=prods.get,reverse=True) # sort by product
+  result(num,array,prods,skeys) # print results
+
+if __name__ == '__main__' : consecutive_product()