# Programming With Prime Numbers
# http://programmingpraxis.files.wordpress.com/
# 2012/09/programmingwithprimenumbers.pdf

def sieve(n):
    b, p, ps = [True] * (n+1), 2, []
    for p in xrange(2, n+1):
        if b[p]:
            ps.append(p)
            for i in xrange(p, n+1, p):
                b[i] = False
    return ps

def primes(n):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    if n < 2:
        raise ValueError('must be greater than one')
    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 td_prime(n, limit=1000000):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    if n % 2 == 0:
        return n == 2
    d = 3
    while d * d <= n:
        if limit < d:
            raise OverflowError('limit exceeded')
        if n % d == 0:
            return False
        d += 2
    return True

def td_factors(n, limit=1000000):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    fs = []
    while n % 2 == 0:
        fs += [2]
        n /= 2
    if n == 1:
        return fs
    f = 3
    while f * f <= n:
        if limit < f:
            raise OverflowError('limit exceeded')
        if n % f == 0:
            fs += [f]
            n /= f
        else:
            f += 2
    return fs + [n]

def is_prime(n):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    if n < 2:
        return False
    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
        t = pow(a,d,n)
        if t == 1:
            return True
        while s > 0:
            if t == n-1:
                return True
            t = (t*t) % n
            s -= 1
        return False
    if n in ps: return True
    for p in ps:
        if not is_spsp(n,p):
            return False
    return True

def rho_factors(n, limit=1000000):
    if type(n) != int and type(n) != long:
        raise TypeError('must be integer')
    def gcd(a,b):
        while b: a, b = b, a%b
        return abs(a)
    def rho_factor(n, c, limit):
        f = lambda(x): (x*x+c) % n
        t, h, d = 2, 2, 1
        while d == 1:
            if limit == 0:
                raise OverflowError('limit exceeded')
            t = f(t); h = f(f(h)); d = gcd(t-h, n)
        if d == n:
            return rho_factor(n, c+1, limit)
        if is_prime(d):
            return d
        return rho_factor(d, c+1, limit)
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + rho_factors(-n, limit)
    fs = []
    while n % 2 == 0:
        n = n // 2; fs = fs + [2]
    if n == 1: return fs
    while not is_prime(n):
        f = rho_factor(n, 1, limit)
        n = n / f
        fs = fs + [f]
    return sorted(fs + [n])

print primes(100)
print len(primes(1000000))
print td_prime(600851475143)
print td_factors(600851475143)
print is_prime(600851475143)
print is_prime(2305843009213693951)
print rho_factors(600851475143)

Follow

Get every new post delivered to your Inbox.

Join 643 other followers

%d bloggers like this: