[ create a new paste ] login | about

Link: http://codepad.org/SQNNuTvI    [ raw code | output | fork ]

programmingpraxis - Python, pasted on Jul 3:
# prime numbers

def primes(n): # sieve of eratosthenes
    i, p, ps, m = 0, 3, [2], n // 2
    sieve = [True] * m
    while p <= n:
        if sieve[i]:
            ps.append(p)
            for j in range((p*p-3)/2, m, p):
                sieve[j] = False
        i, p = i+1, p+2
    return ps

# from random import randint

seed = 17500728 # RIP j s bach

def random(): # float on range [0,1)
    global seed
    seed = (69069 * seed + 1234567) % 4294967296
    return seed / 4294967296.0

def randint(lo,hi): # int on range [lo,hi)
    return int((hi - lo) * random()) + lo

def isPrime(n, k=5): # miller-rabin
    if n < 2: return False
    for p in [2,3,5,7,11,13,17,19,23,29]:
        if n % p == 0: return n == p
    s, d = 0, n-1
    while d % 2 == 0:
        s, d = s+1, d/2
    for i in range(k):
        x = pow(randint(2, n-1), d, n)
        if x == 1 or x == n-1: continue
        for r in range(1, s):
            x = (x * x) % n
            if x == 1: return False
            if x == n-1: break
        else: return False
    return True

# from fractions import gcd

def gcd(a,b): # greatest common divisor
    if b == 0: return a
    return gcd(b, a % b)

def insertSorted(x, xs): # insert x in order
    i, ln = 0, len(xs)
    while i < ln and xs[i] < x: i += 1
    xs.insert(i,x)
    return xs

def factors(n, b2=-1, b1=10000): # 2,3,5-wheel, then rho
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + factors(-n)
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, f, fs = 0, 2, []
    while f*f <= n and f < b1:
        while n % f == 0:
            fs.append(f)
            n /= f
        f, w = f + wheel[w], w+1
        if w == 11: w = 3
    if n == 1: return fs
    h, t, g, c = 1, 1, 1, 1
    while not isPrime(n):
        while b2 <> 0 and g == 1:
            h = (h*h+c)%n # the hare runs
            h = (h*h+c)%n # twice as fast
            t = (t*t+c)%n # as the tortoise
            g = gcd(t-h, n); b2 -= 1
        if b2 == 0: return fs
        if isPrime(g):
            while n % g == 0:
                fs = insertSorted(g, fs)
                n /= g
        h, t, g, c = 1, 1, 1, c+1
    return insertSorted(n, fs)

def sigma(x, n, fs=[]): # sum of x'th powers of divisors of n
    def add(s, p, m):
    	if x == 0: return s * (m+1)
    	return s * (p**(x*(m+1))-1) / (p**x-1)
    if fs == []: fs = factors(n)
    prev, mult, sum = fs.pop(0), 1, 1
    while len(fs) > 0:
        fact = fs.pop(0)
        if fact <> prev:
            sum, prev, mult = add(sum, prev, mult), fact, 1
        else: mult += 1
    return add(sum, prev, mult)

def aliquot(n): # print aliquot sequence
    s, ss, k, fs = n, [n], 0, factors(n)
    print n, k, s, fs
    while s > 1:
        s, k = sigma(1,s,fs) - s, k + 1
        fs = factors(s)
        print n, k, s, fs
        if s in ss: return "cycle"
        ss.append(s)
    return ss.pop(-2)

# project euler 12
i, t = 1, 1
while sigma(0, t) < 500:
    i += 1; t += i
print t


Output:
1
76576500


Create a new paste based on this one


Comments: