# prime.py -- modest prime number library

def primes(n):
    """
        list of primes not exceeding n in ascending
        order; assumes n is an integer greater than
        1; uses Sieve of Eratosthenes
    """
    m = (n-1) // 2
    b = [True] * m
    i, p, ps = 0, 3, [2]
    while p*p < n:
        if b[i]:
            ps.append(p)
            j = 2*i*i + 6*i + 3
            while j < m:
                b[j] = False
                j = j + 2*i + 3
        i += 1; p += 2
    while i < m:
        if b[i]:
            ps.append(p)
        i += 1; p += 2
    return ps

def is_prime(n):
    """
        False if n is provably composite, else
        True if n is probably prime; assumes n
        is an integer greater than 1; uses
        Miller-Rabin test on prime bases < 100
    """
    ps = [2,3,5,7,11,13,17,19,23,29,31,37,41,
         43,47,53,59,61,67,71,73,79,83,89,97]
    def is_spsp(n, a):
        d, s = n-1, 0
        while d%2 == 0:
            d /= 2; s += 1
        if pow(a,d,n) == 1:
            return True
        for r in xrange(s):
            if pow(a, d*pow(2,r), n) == n-1:
                return True
        return False
    if n in ps: return True
    for p in ps:
        if not is_spsp(n,p):
            return False
    return True

def factors(n):
    """
        list of prime factors of n in ascending
        order; assumes n is an integer, may be
        positive, zero or negative; uses Pollard's
        rho algorithm with Floyd's cycle finder
    """
    def gcd(a,b):
        while b: a, b = b, a%b
        return abs(a)
    def facts(n,c,fs):
        f = lambda(x): (x*x+c) % n
        if is_prime(n): return [n]+fs
        t, h, d = 2, 2, 1
        while d == 1:
            t = f(t); h = f(f(h))
            d = gcd(t-h, n)
        if d == n:
            return facts(n, c+1, fs)
        if is_prime(d):
            return facts(n//d, c+1, [d]+fs)
        return facts(n, c+1, fs)
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + factors(-n)
    fs = []
    while n%2 == 0:
        n = n//2; fs = [2]+fs
    if n == 1: return fs
    return sorted(facts(n,1,fs))

# simple demonstrations
print primes(100)
print is_prime(pow(2,89)-1)
print factors(600851475143)

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
True
[71L, 839L, 1471L, 6857L]