• API
• FAQ
• Tools
• Archive
daily pastebin goal
41%
SHARE
TWEET

# Untitled

a guest Jul 15th, 2018 71 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. """Prime number generation and number factorization."""
2.
3. import bisect, itertools, random, sys
4.
5. _primes = [2]
6.
7. _miller_rabin_limit = 48611    # 5000th prime
8. _miller_rabin_security = 7
9.
10. def modpow (a, b, c):
11.     """Efficiently compute (a^b)%c where a, b and c are positive integers."""
12.     if b == 1: return a % c
13.     d = b//2
14.     x = modpow (a, d, c)
15.     x = x*x%c
16.     if b % 2 == 1: x = x*a%c
17.     return x
18.
19. def miller_rabin (n, t = _miller_rabin_security):
20.     """Apply Miller-Rabin primality test to detect whether n is prime or not.
21.     The test is run t times."""
22.     if n % 2 == 0: return False
23.     r = (n-1)//2
24.     for s in itertools.count (1):
25.         if r % 2: break
26.         r //= 2
27.     for tt in xrange (t):
28.         a = int (random.uniform (2, n-1))
29.         y = modpow (a, r, n)
30.         if y != 1 and y != n-1:
31.             for j in xrange (1, s):
32.                 y = y**2 % n
33.                 if y == 1: return False
34.                 if y == n-1: break
35.             if y != n-1: return False
36.     return True
37.
38. def _is_prime (i):
39.     if i > _miller_rabin_limit: return miller_rabin (i)
40.     s = int (i**0.5)
41.     for j in _primes:
42.         if i % j == 0: return False
43.         if j > s: return True
44.     return True
45.
46. def _gen_primes ():
47.     for i in xrange (3, sys.maxint, 2):
48.         if _is_prime (i):
49.             _primes.append (i)
50.             yield i
51.
52. _primary = _gen_primes ()
53.
54. def _gen_primes (minval):
55.     i = 0
56.     if minval:
57.         l = _primes[-1]
58.         if minval > l:
59.             while minval > l: l =_primary.next ()
60.             i = len (_primes) - 1
61.         else:
62.             i = bisect.bisect_left (_primes, minval)
63.     for i in itertools.count (i):
64.         if i == len (_primes): yield _primary.next ()
65.         else: yield _primes[i]
66.
67. def gen_primes (maxval = None, minval = None):
68.     """Prime numbers generator. If minval is given, only primes greater or
69.     equal to minval are returned. If maxval is given, only primes smaller
70.     or equal to maxval are returned.
71.
72.     >>> for p in gen_primes(5): print p
73.     ...
74.     2
75.     3
76.     5
77.     7
78.     11"""
79.     if maxval:
80.         return itertools.takewhile (lambda x: x<=maxval, _gen_primes (minval))
81.     return _gen_primes (minval)
82.
83. def gen_factors (n, duplicates = True):
84.     """Generator for factors of n (n > 1). If duplicates is False, do not
85.     send the same factor more than once."""
86.     assert n > 1
87.     if n > _miller_rabin_limit and miller_rabin (n):
88.         yield n
89.         return
90.     s = int (n**0.5)
91.     for i in gen_primes ():
92.         if i > s:
93.             yield n
94.             return
95.         if n % i == 0:
96.             yield i
97.             n //= i
98.             while n % i == 0:
99.                 if duplicates: yield i
100.                 n //= i
101.             if n == 1: return
102.             if n > _miller_rabin_limit and miller_rabin (n):
103.                 yield n
104.                 return
105.             s = int (n**0.5)
106.
107. def factors (n):
108.     """Factorize n (n > 1) into its prime factors. Return a dictionary where
109.     keys are prime factors and values are powers.
110.
111.     >>> factors(18)
112.     {2: 1, 3: 2}"""
113.     l = {}
114.     for i in gen_factors (n):
115.         try: l[i] += 1
116.         except KeyError: l[i] = 1
117.     return l
118.
119. def factorslist (n):
120.     """Return a list of prime factors of n (n > 1).
121.
122.     >>> factorslist(18)
123.     [2, 3, 3]"""
124.     return list (gen_factors (n))
125.
126. def is_prime (n):
127.     """Check whether n is prime or not."""
128.     if n < 2: return False
129.     return gen_factors(n).next() == n
130.
131. def ufactors (n):
132.     """Return the list of unique prime factors of n (n > 1).
133.
134.     >>> ufactors(100)
135.     [2, 5]"""
136.     return list (gen_factors(n, duplicates = False))
137.
138. def nfactors (n):
139.     """Return the number of unique prime factors of n."""
140.     return len (ufactors (n))
141.
142. def totient (n):
143.     """Euler's totient function. Returns the number of integers between
144.     1 and n-1 relatively prime to n (n > 0).
145.
146.     >>> totient(6)
147.     2
148.     (6 is relatively prime to 1 and 5)"""
149.     if n == 1: return 0
150.     num = den = 1
151.     for p in gen_factors (n, duplicates = False):
152.         num *= (p-1)
153.         den *= p
154.     return n * num // den
155.
156. _s = {1: 1}
157. def s (n):
158.     """Sum of divisors of n (n > 0).
159.
160.     >>> s(6)
161.     12
162.     (divisors of 6 are 1, 2, 3 and 6, summing to 12)"""
163.     if not _s.has_key (n):
164.         t = 1
165.         for (p, k) in factors(n).items():
166.             t *= (p**(k+1)-1)//(p-1)
167.         _s[n] = t
168.     return _s[n]
169.
170. if __name__ == '__main__':
171.     # Quick test -- add primes up to 100,000 and compare to J result to:
172.     #   +/p:i.(p:^:_1)100000x
173.     assert sum (gen_primes (100000)) == 454396537
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top