```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 ``` ```# 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 fs+[n] 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, fs+[d]) return facts(n, c+1, fs) if -1 <= n <= 1: return [n] if n < -1: return factors(-n)+[-1] fs = [] while n%2 == 0: n = n//2; fs = fs+[2] 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) # divisors def numdiv(n): fs = factors(n) f = fs.pop(0); d = 1; x = 2 while fs: if f == fs[0]: x += 1 else: d *= x; x = 2 f = fs.pop(0) return d * x print numdiv(76576500) ```
 ```1 2 3 4 ``` ```[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] 576 ```