# 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