```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 ``` ```# 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 ```
 ```1 ``` ```76576500 ```