这里会显示出您选择的修订版和当前版本之间的差别。
— |
modules:pyprimes-source-code [2014/07/15 15:55] (当前版本) |
||
---|---|---|---|
行 1: | 行 1: | ||
+ | ====== pyprimes 源代码 ====== | ||
+ | |||
+ | <code python> | ||
+ | #!/usr/bin/env python | ||
+ | |||
+ | ## Module pyprimes.py | ||
+ | ## | ||
+ | ## Copyright (c) 2012 Steven D'Aprano. | ||
+ | ## | ||
+ | ## Permission is hereby granted, free of charge, to any person obtaining | ||
+ | ## a copy of this software and associated documentation files (the | ||
+ | ## "Software"), to deal in the Software without restriction, including | ||
+ | ## without limitation the rights to use, copy, modify, merge, publish, | ||
+ | ## distribute, sublicense, and/or sell copies of the Software, and to | ||
+ | ## permit persons to whom the Software is furnished to do so, subject to | ||
+ | ## the following conditions: | ||
+ | ## | ||
+ | ## The above copyright notice and this permission notice shall be | ||
+ | ## included in all copies or substantial portions of the Software. | ||
+ | ## | ||
+ | ## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, | ||
+ | ## EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF | ||
+ | ## MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. | ||
+ | ## IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY | ||
+ | ## CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, | ||
+ | ## TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE | ||
+ | ## SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. | ||
+ | |||
+ | |||
+ | """Generate and test for small primes using a variety of algorithms | ||
+ | implemented in pure Python. | ||
+ | |||
+ | This module includes functions for generating prime numbers, primality | ||
+ | testing, and factorising numbers into prime factors. Prime numbers are | ||
+ | positive integers with no factors other than themselves and 1. | ||
+ | |||
+ | |||
+ | Generating prime numbers | ||
+ | ======================== | ||
+ | |||
+ | To generate an unending stream of prime numbers, use the ``primes()`` | ||
+ | generator function: | ||
+ | |||
+ | primes(): | ||
+ | Yield prime numbers 2, 3, 5, 7, 11, ... | ||
+ | |||
+ | |||
+ | >>> p = primes() | ||
+ | >>> [next(p) for _ in range(10)] | ||
+ | [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] | ||
+ | |||
+ | |||
+ | To efficiently generate pairs of (isprime(i), i) for integers i, use the | ||
+ | generator functions ``checked_ints()`` and ``checked_oddints()``: | ||
+ | |||
+ | checked_ints() | ||
+ | Yield pairs of (isprime(i), i) for i=0,1,2,3,4,5... | ||
+ | |||
+ | checked_oddints() | ||
+ | Yield pairs of (isprime(i), i) for odd i=1,3,5,7... | ||
+ | |||
+ | |||
+ | >>> it = checked_ints() | ||
+ | >>> [next(it) for _ in range(5)] | ||
+ | [(False, 0), (False, 1), (True, 2), (True, 3), (False, 4)] | ||
+ | |||
+ | |||
+ | Other convenience functions wrapping ``primes()`` are: | ||
+ | |||
+ | ------------------ ---------------------------------------------------- | ||
+ | Function Description | ||
+ | ------------------ ---------------------------------------------------- | ||
+ | nprimes(n) Yield the first n primes, then stop. | ||
+ | nth_prime(n) Return the nth prime number. | ||
+ | prime_count(x) Return the number of primes less than or equal to x. | ||
+ | primes_below(x) Yield the primes less than or equal to x. | ||
+ | primes_above(x) Yield the primes strictly greater than x. | ||
+ | primesum(n) Return the sum of the first n primes. | ||
+ | primesums() Yield the partial sums of the prime numbers. | ||
+ | ------------------ ---------------------------------------------------- | ||
+ | |||
+ | |||
+ | Primality testing | ||
+ | ================= | ||
+ | |||
+ | These functions test whether numbers are prime or not. Primality tests fall | ||
+ | into two categories: exact tests, and probabilistic tests. | ||
+ | |||
+ | Exact tests are guaranteed to give the correct result, but may be slow, | ||
+ | particularly for large arguments. Probabilistic tests do not guarantee | ||
+ | correctness, but may be much faster for large arguments. | ||
+ | |||
+ | To test whether an integer is prime, use the ``isprime`` function: | ||
+ | |||
+ | isprime(n) | ||
+ | Return True if n is prime, otherwise return False. | ||
+ | |||
+ | |||
+ | >>> isprime(101) | ||
+ | True | ||
+ | >>> isprime(102) | ||
+ | False | ||
+ | |||
+ | |||
+ | Exact primality tests are: | ||
+ | |||
+ | isprime_naive(n) | ||
+ | Naive and slow trial division test for n being prime. | ||
+ | |||
+ | isprime_division(n) | ||
+ | A less naive trial division test for n being prime. | ||
+ | |||
+ | isprime_regex(n) | ||
+ | Uses a regex to test if n is a prime number. | ||
+ | |||
+ | .. NOTE:: ``isprime_regex`` should be considered a novelty | ||
+ | rather than a serious test, as it is very slow. | ||
+ | |||
+ | |||
+ | Probabilistic tests do not guarantee correctness, but can be faster for | ||
+ | large arguments. There are two probabilistic tests: | ||
+ | |||
+ | fermat(n [, base]) | ||
+ | Fermat primality test, returns True if n is a weak probable | ||
+ | prime to the given base, otherwise False. | ||
+ | |||
+ | miller_rabin(n [, base]) | ||
+ | Miller-Rabin primality test, returns True if n is a strong | ||
+ | probable prime to the given base, otherwise False. | ||
+ | |||
+ | |||
+ | Both guarantee no false negatives: if either function returns False, the | ||
+ | number being tested is certainly composite. However, both are subject to false | ||
+ | positives: if they return True, the number is only possibly prime. | ||
+ | |||
+ | |||
+ | >>> fermat(12400013) # composite 23*443*1217 | ||
+ | False | ||
+ | >>> miller_rabin(14008971) # composite 3*947*4931 | ||
+ | False | ||
+ | |||
+ | |||
+ | Prime factorisation | ||
+ | =================== | ||
+ | |||
+ | These functions return or yield the prime factors of an integer. | ||
+ | |||
+ | factors(n) | ||
+ | Return a list of the prime factors of n. | ||
+ | |||
+ | factorise(n) | ||
+ | Yield tuples (factor, count) for n. | ||
+ | |||
+ | |||
+ | The ``factors(n)`` function lists repeated factors: | ||
+ | |||
+ | |||
+ | >>> factors(37*37*109) | ||
+ | [37, 37, 109] | ||
+ | |||
+ | |||
+ | The ``factorise(n)`` generator yields a 2-tuple for each unique factor, giving | ||
+ | the factor itself and the number of times it is repeated: | ||
+ | |||
+ | >>> list(factorise(37*37*109)) | ||
+ | [(37, 2), (109, 1)] | ||
+ | |||
+ | |||
+ | Alternative and toy prime number generators | ||
+ | =========================================== | ||
+ | |||
+ | These functions are alternative methods of generating prime numbers. Unless | ||
+ | otherwise stated, they generate prime numbers lazily on demand. These are | ||
+ | supplied for educational purposes and are generally slower or less efficient | ||
+ | than the preferred ``primes()`` generator. | ||
+ | |||
+ | -------------- -------------------------------------------------------- | ||
+ | Function Description | ||
+ | -------------- -------------------------------------------------------- | ||
+ | croft() Yield prime numbers using the Croft Spiral sieve. | ||
+ | erat(n) Return primes up to n by the sieve of Eratosthenes. | ||
+ | sieve() Yield primes using the sieve of Eratosthenes. | ||
+ | cookbook() Yield primes using "Python Cookbook" algorithm. | ||
+ | wheel() Yield primes by wheel factorization. | ||
+ | -------------- -------------------------------------------------------- | ||
+ | |||
+ | .. TIP:: In the current implementation, the fastest of these | ||
+ | generators is aliased as ``primes()``. | ||
+ | |||
+ | |||
+ | These generators however are extremely slow, and should be considered as | ||
+ | examples of algorithms which should *not* be used, supplied as a horrible | ||
+ | warning of how *not* to calculate primes. | ||
+ | |||
+ | -------------- -------------------------------------------------------- | ||
+ | Function Description | ||
+ | -------------- -------------------------------------------------------- | ||
+ | awful_primes Yield primes very slowly by an awful algorithm. | ||
+ | naive_primes1 Yield primes slowly by a naive algorithm. | ||
+ | naive_primes2 Yield primes slowly by a less naive algorithm. | ||
+ | trial_division Yield primes slowly using trial division. | ||
+ | turner Yield primes very slowly using Turner's algorithm. | ||
+ | -------------- -------------------------------------------------------- | ||
+ | |||
+ | """ | ||
+ | |||
+ | |||
+ | from __future__ import division | ||
+ | |||
+ | |||
+ | import functools | ||
+ | import itertools | ||
+ | import random | ||
+ | |||
+ | |||
+ | # Module metadata. | ||
+ | __version__ = "0.1.1a" | ||
+ | __date__ = "2012-02-22" | ||
+ | __author__ = "Steven D'Aprano" | ||
+ | __author_email__ = "steve+python@pearwood.info" | ||
+ | |||
+ | __all__ = ['primes', 'checked_ints', 'checked_oddints', 'nprimes', | ||
+ | 'primes_above', 'primes_below', 'nth_prime', 'prime_count', | ||
+ | 'primesum', 'primesums', 'warn_probably', 'isprime', 'factors', | ||
+ | 'factorise', | ||
+ | ] | ||
+ | |||
+ | |||
+ | # ============================ | ||
+ | # Python 2.x/3.x compatibility | ||
+ | # ============================ | ||
+ | |||
+ | # This module should support 2.5+, including Python 3. | ||
+ | |||
+ | try: | ||
+ | next | ||
+ | except NameError: | ||
+ | # No next() builtin, so we're probably running Python 2.5. | ||
+ | # Use a simplified version (without support for default). | ||
+ | def next(iterator): | ||
+ | return iterator.next() | ||
+ | |||
+ | try: | ||
+ | range = xrange | ||
+ | except NameError: | ||
+ | # No xrange built-in, so we're almost certainly running Python3 | ||
+ | # and range is already a lazy iterator. | ||
+ | assert type(range(3)) is not list | ||
+ | |||
+ | try: | ||
+ | from itertools import ifilter as filter, izip as zip | ||
+ | except ImportError: | ||
+ | # Python 3, where filter and zip are already lazy. | ||
+ | assert type(filter(None, [1, 2])) is not list | ||
+ | assert type(zip("ab", [1, 2])) is not list | ||
+ | |||
+ | try: | ||
+ | from itertools import compress | ||
+ | except ImportError: | ||
+ | # Must be Python 2.x, so we need to roll our own. | ||
+ | def compress(data, selectors): | ||
+ | """compress('ABCDEF', [1,0,1,0,1,1]) --> A C E F""" | ||
+ | return (d for d, s in zip(data, selectors) if s) | ||
+ | |||
+ | try: | ||
+ | from math import isfinite | ||
+ | except ImportError: | ||
+ | # Python 2.6 or older. | ||
+ | try: | ||
+ | from math import isnan, isinf | ||
+ | except ImportError: | ||
+ | # Python 2.5. Quick and dirty substitutes. | ||
+ | def isnan(x): | ||
+ | return x != x | ||
+ | def isinf(x): | ||
+ | return x - x != 0 | ||
+ | def isfinite(x): | ||
+ | return not (isnan(x) or isinf(x)) | ||
+ | |||
+ | |||
+ | # ===================== | ||
+ | # Helpers and utilities | ||
+ | # ===================== | ||
+ | |||
+ | def _validate_int(obj): | ||
+ | """Raise an exception if obj is not an integer.""" | ||
+ | m = int(obj + 0) # May raise TypeError. | ||
+ | if obj != m: | ||
+ | raise ValueError('expected an integer but got %r' % obj) | ||
+ | |||
+ | |||
+ | def _validate_num(obj): | ||
+ | """Raise an exception if obj is not a finite real number.""" | ||
+ | m = obj + 0 # May raise TypeError. | ||
+ | if not isfinite(m): | ||
+ | raise ValueError('expected a finite real number but got %r' % obj) | ||
+ | |||
+ | |||
+ | def _base_to_bases(base, n): | ||
+ | if isinstance(base, tuple): | ||
+ | bases = base | ||
+ | else: | ||
+ | bases = (base,) | ||
+ | for b in bases: | ||
+ | _validate_int(b) | ||
+ | if not 1 <= b < n: | ||
+ | # Note that b=1 is a degenerate case which is always a prime | ||
+ | # witness for both the Fermat and Miller-Rabin tests. I mention | ||
+ | # this for completeness, not because we need to do anything | ||
+ | # about it. | ||
+ | raise ValueError('base %d out of range 1...%d' % (b, n-1)) | ||
+ | return bases | ||
+ | |||
+ | |||
+ | # ======================= | ||
+ | # Prime number generators | ||
+ | # ======================= | ||
+ | |||
+ | # The preferred generator to use is ``primes()``, which will be set to the | ||
+ | # "best" of these generators. (If you disagree with my judgement of best, | ||
+ | # feel free to use the generator of your choice.) | ||
+ | |||
+ | |||
+ | def erat(n): | ||
+ | """Return a list of primes up to and including n. | ||
+ | |||
+ | This is a fixed-size version of the Sieve of Eratosthenes, using an | ||
+ | adaptation of the traditional algorithm. | ||
+ | |||
+ | >>> erat(30) | ||
+ | [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] | ||
+ | >>> erat(10000) == list(primes_below(10000)) | ||
+ | True | ||
+ | |||
+ | """ | ||
+ | _validate_int(n) | ||
+ | # Generate a fixed array of integers. | ||
+ | arr = list(range(n+1)) | ||
+ | # Cross out 0 and 1 since they aren't prime. | ||
+ | arr[0] = arr[1] = None | ||
+ | i = 2 | ||
+ | while i*i <= n: | ||
+ | # Cross out all the multiples of i starting from i**2. | ||
+ | for p in range(i*i, n+1, i): | ||
+ | arr[p] = None | ||
+ | # Advance to the next number not crossed off. | ||
+ | i += 1 | ||
+ | while i <= n and arr[i] is None: | ||
+ | i += 1 | ||
+ | return list(filter(None, arr)) | ||
+ | |||
+ | |||
+ | def sieve(): | ||
+ | """Yield prime integers using the Sieve of Eratosthenes. | ||
+ | |||
+ | This algorithm is modified to generate the primes lazily rather than the | ||
+ | traditional version which operates on a fixed size array of integers. | ||
+ | """ | ||
+ | # This is based on a paper by Melissa E. O'Neill, with an implementation | ||
+ | # given by Gerald Britton: | ||
+ | # http://mail.python.org/pipermail/python-list/2009-January/1188529.html | ||
+ | innersieve = sieve() | ||
+ | prevsq = 1 | ||
+ | table = {} | ||
+ | i = 2 | ||
+ | while True: | ||
+ | # Note: this explicit test is slightly faster than using | ||
+ | # prime = table.pop(i, None) and testing for None. | ||
+ | if i in table: | ||
+ | prime = table[i] | ||
+ | del table[i] | ||
+ | nxt = i + prime | ||
+ | while nxt in table: | ||
+ | nxt += prime | ||
+ | table[nxt] = prime | ||
+ | else: | ||
+ | yield i | ||
+ | if i > prevsq: | ||
+ | j = next(innersieve) | ||
+ | prevsq = j**2 | ||
+ | table[prevsq] = j | ||
+ | i += 1 | ||
+ | |||
+ | |||
+ | def cookbook(): | ||
+ | """Yield prime integers lazily using the Sieve of Eratosthenes. | ||
+ | |||
+ | Another version of the algorithm, based on the Python Cookbook, | ||
+ | 2nd Edition, recipe 18.10, variant erat2. | ||
+ | """ | ||
+ | # http://onlamp.com/pub/a/python/excerpt/pythonckbk_chap1/index1.html?page=2 | ||
+ | table = {} | ||
+ | yield 2 | ||
+ | # Iterate over [3, 5, 7, 9, ...]. The following is equivalent to, but | ||
+ | # faster than, (2*i+1 for i in itertools.count(1)) | ||
+ | for q in itertools.islice(itertools.count(3), 0, None, 2): | ||
+ | # Note: this explicit test is marginally faster than using | ||
+ | # table.pop(i, None) and testing for None. | ||
+ | if q in table: | ||
+ | p = table[q]; del table[q] # Faster than pop. | ||
+ | x = p + q | ||
+ | while x in table or not (x & 1): | ||
+ | x += p | ||
+ | table[x] = p | ||
+ | else: | ||
+ | table[q*q] = q | ||
+ | yield q | ||
+ | |||
+ | |||
+ | def croft(): | ||
+ | """Yield prime integers using the Croft Spiral sieve. | ||
+ | |||
+ | This is a variant of wheel factorisation modulo 30. | ||
+ | """ | ||
+ | # Implementation is based on erat3 from here: | ||
+ | # http://stackoverflow.com/q/2211990 | ||
+ | # and this website: | ||
+ | # http://www.primesdemystified.com/ | ||
+ | # Memory usage increases roughly linearly with the number of primes seen. | ||
+ | # dict ``roots`` stores an entry x:p for every prime p. | ||
+ | for p in (2, 3, 5): | ||
+ | yield p | ||
+ | roots = {9: 3, 25: 5} # Map d**2 -> d. | ||
+ | primeroots = frozenset((1, 7, 11, 13, 17, 19, 23, 29)) | ||
+ | selectors = (1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0) | ||
+ | for q in compress( | ||
+ | # Iterate over prime candidates 7, 9, 11, 13, ... | ||
+ | itertools.islice(itertools.count(7), 0, None, 2), | ||
+ | # Mask out those that can't possibly be prime. | ||
+ | itertools.cycle(selectors) | ||
+ | ): | ||
+ | # Using dict membership testing instead of pop gives a | ||
+ | # 5-10% speedup over the first three million primes. | ||
+ | if q in roots: | ||
+ | p = roots[q] | ||
+ | del roots[q] | ||
+ | x = q + 2*p | ||
+ | while x in roots or (x % 30) not in primeroots: | ||
+ | x += 2*p | ||
+ | roots[x] = p | ||
+ | else: | ||
+ | roots[q*q] = q | ||
+ | yield q | ||
+ | |||
+ | |||
+ | def wheel(): | ||
+ | """Generate prime numbers using wheel factorisation modulo 210.""" | ||
+ | for i in (2, 3, 5, 7, 11): | ||
+ | yield i | ||
+ | # The following constants are taken from the paper by O'Neill. | ||
+ | spokes = (2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, | ||
+ | 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, | ||
+ | 6, 4, 2, 4, 2, 10, 2, 10) | ||
+ | assert len(spokes) == 48 | ||
+ | # This removes about 77% of the composites that we would otherwise | ||
+ | # need to divide by. | ||
+ | found = [(11, 121)] # Smallest prime we care about, and its square. | ||
+ | for incr in itertools.cycle(spokes): | ||
+ | i += incr | ||
+ | for p, p2 in found: | ||
+ | if p2 > i: # i must be a prime. | ||
+ | found.append((i, i*i)) | ||
+ | yield i | ||
+ | break | ||
+ | elif i % p == 0: # i must be composite. | ||
+ | break | ||
+ | else: # This should never happen. | ||
+ | raise RuntimeError("internal error: ran out of prime divisors") | ||
+ | |||
+ | |||
+ | # This is the preferred way of generating prime numbers. Set this to the | ||
+ | # fastest/best generator. | ||
+ | primes = croft | ||
+ | |||
+ | |||
+ | # === Algorithms to avoid === | ||
+ | |||
+ | # The following algorithms are supplied for educational purposes, as toys, | ||
+ | # curios, or as terrible warnings on what NOT to use. | ||
+ | # | ||
+ | # None of these have acceptable performance; they are barely tolerable even | ||
+ | # for the first 100 primes. | ||
+ | |||
+ | def awful_primes(): | ||
+ | """Generate prime numbers naively, and REALLY slowly. | ||
+ | |||
+ | This is about as awful as you can get while still being a straight-forward | ||
+ | and unobfuscated implementation. What makes this particularly awful is | ||
+ | that it doesn't stop testing for factors when it finds one, but | ||
+ | pointlessly keeps testing. | ||
+ | """ | ||
+ | i = 2 | ||
+ | yield i | ||
+ | while True: | ||
+ | i += 1 | ||
+ | composite = False | ||
+ | for p in range(2, i): | ||
+ | if i%p == 0: | ||
+ | composite = True | ||
+ | if not composite: # It must be a prime. | ||
+ | yield i | ||
+ | |||
+ | |||
+ | def naive_primes1(): | ||
+ | """Generate prime numbers naively and slowly. | ||
+ | |||
+ | This is a bit better than awful_primes, as it short-circuits testing for | ||
+ | composites. If it finds a single factor, it stops testing. Nevertheless, | ||
+ | this is still terribly slow. | ||
+ | """ | ||
+ | i = 2 | ||
+ | yield i | ||
+ | while True: | ||
+ | i += 1 | ||
+ | if all(i%p != 0 for p in range(2, i)): | ||
+ | yield i | ||
+ | |||
+ | |||
+ | def naive_primes2(): | ||
+ | """Generate prime numbers naively and slowly. | ||
+ | |||
+ | This is an incremental improvement over ``naive_primes1``, by only testing | ||
+ | for odd factors. | ||
+ | """ | ||
+ | yield 2 | ||
+ | i = 3 | ||
+ | yield i | ||
+ | while True: | ||
+ | i += 2 | ||
+ | if all(i%p != 0 for p in range(3, i, 2)): | ||
+ | yield i | ||
+ | |||
+ | |||
+ | def trial_division(): | ||
+ | """Generate prime numbers slowly using a simple trial division algorithm. | ||
+ | |||
+ | This uses three optimizations: we only test odd numbers for primality, | ||
+ | we only test with prime factors, and that only up to the square root of | ||
+ | the number being tested. This gives us asymptotic behaviour of | ||
+ | O(N*sqrt(N)/(log N)**2) where N is the number of primes found. | ||
+ | |||
+ | Despite these optimizations, this is still unacceptably slow, especially | ||
+ | as the list of memorised primes grows. | ||
+ | """ | ||
+ | yield 2 | ||
+ | primes = [2] | ||
+ | i = 3 | ||
+ | while 1: | ||
+ | it = itertools.takewhile(lambda p, i=i: p*p <= i, primes) | ||
+ | if all(i%p != 0 for p in it): | ||
+ | primes.append(i) | ||
+ | yield i | ||
+ | i += 2 | ||
+ | |||
+ | |||
+ | def turner(): | ||
+ | """Yield prime numbers very slowly using Euler's sieve. | ||
+ | |||
+ | The function is named for David Turner, who developed this implementation | ||
+ | in a paper in 1975. Due to its simplicity, it has become very popular, | ||
+ | particularly in Haskell circles where it is usually implemented as some | ||
+ | variation of: | ||
+ | |||
+ | primes = sieve [2..] | ||
+ | sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p > 0] | ||
+ | |||
+ | This algorithm is often wrongly described as the Sieve of Eratosthenes, | ||
+ | but it is not. Although simple, it is slow and inefficient, with | ||
+ | asymptotic behaviour of O(N**2/(log N)**2), which is even worse than | ||
+ | trial_division, and only marginally better than naive_primes. O'Neill | ||
+ | calls this the "Sleight on Eratosthenes". | ||
+ | """ | ||
+ | # References: | ||
+ | # http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes | ||
+ | # http://en.literateprograms.org/Sieve_of_Eratosthenes_(Haskell) | ||
+ | # http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf | ||
+ | # http://www.haskell.org/haskellwiki/Prime_numbers | ||
+ | nums = itertools.count(2) | ||
+ | while True: | ||
+ | prime = next(nums) | ||
+ | yield prime | ||
+ | nums = filter(lambda v, p=prime: (v % p) != 0, nums) | ||
+ | |||
+ | |||
+ | # ===================== | ||
+ | # Convenience functions | ||
+ | # ===================== | ||
+ | |||
+ | def checked_ints(): | ||
+ | """Yield tuples (isprime(i), i) for integers i=0, 1, 2, 3, 4, ... | ||
+ | |||
+ | >>> it = checked_ints() | ||
+ | >>> [next(it) for _ in range(6)] | ||
+ | [(False, 0), (False, 1), (True, 2), (True, 3), (False, 4), (True, 5)] | ||
+ | |||
+ | """ | ||
+ | oddnums = checked_oddints() | ||
+ | yield (False, 0) | ||
+ | yield next(oddnums) | ||
+ | yield (True, 2) | ||
+ | for t in oddnums: | ||
+ | yield t | ||
+ | yield (False, t[1]+1) | ||
+ | |||
+ | |||
+ | def checked_oddints(): | ||
+ | """Yield tuples (isprime(i), i) for odd integers i=1, 3, 5, 7, 9, ... | ||
+ | |||
+ | >>> it = checked_oddints() | ||
+ | >>> [next(it) for _ in range(6)] | ||
+ | [(False, 1), (True, 3), (True, 5), (True, 7), (False, 9), (True, 11)] | ||
+ | >>> [next(it) for _ in range(6)] | ||
+ | [(True, 13), (False, 15), (True, 17), (True, 19), (False, 21), (True, 23)] | ||
+ | |||
+ | """ | ||
+ | yield (False, 1) | ||
+ | odd_primes = primes() | ||
+ | _ = next(odd_primes) # Skip 2. | ||
+ | prev = 1 | ||
+ | for p in odd_primes: | ||
+ | # Yield the non-primes between the previous prime and | ||
+ | # the current one. | ||
+ | for i in itertools.islice(itertools.count(prev + 2), 0, None, 2): | ||
+ | if i >= p: break | ||
+ | yield (False, i) | ||
+ | # And yield the current prime. | ||
+ | yield (True, p) | ||
+ | prev = p | ||
+ | |||
+ | |||
+ | def nprimes(n): | ||
+ | """Convenience function that yields the first n primes. | ||
+ | |||
+ | >>> list(nprimes(10)) | ||
+ | [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] | ||
+ | |||
+ | """ | ||
+ | _validate_int(n) | ||
+ | return itertools.islice(primes(), n) | ||
+ | |||
+ | |||
+ | def primes_above(x): | ||
+ | """Convenience function that yields primes strictly greater than x. | ||
+ | |||
+ | >>> next(primes_above(200)) | ||
+ | 211 | ||
+ | |||
+ | """ | ||
+ | _validate_num(x) | ||
+ | it = primes() | ||
+ | # Consume the primes below x as fast as possible, then yield the rest. | ||
+ | p = next(it) | ||
+ | while p <= x: | ||
+ | p = next(it) | ||
+ | yield p | ||
+ | for p in it: | ||
+ | yield p | ||
+ | |||
+ | |||
+ | def primes_below(x): | ||
+ | """Convenience function yielding primes less than or equal to x. | ||
+ | |||
+ | >>> list(primes_below(20)) | ||
+ | [2, 3, 5, 7, 11, 13, 17, 19] | ||
+ | |||
+ | """ | ||
+ | _validate_num(x) | ||
+ | for p in primes(): | ||
+ | if p > x: | ||
+ | return | ||
+ | yield p | ||
+ | |||
+ | |||
+ | def nth_prime(n): | ||
+ | """nth_prime(n) -> int | ||
+ | |||
+ | Return the nth prime number, starting counting from 1. Equivalent to | ||
+ | p-subscript-n in standard maths notation. | ||
+ | |||
+ | >>> nth_prime(1) # First prime is 2. | ||
+ | 2 | ||
+ | >>> nth_prime(5) | ||
+ | 11 | ||
+ | >>> nth_prime(50) | ||
+ | 229 | ||
+ | |||
+ | """ | ||
+ | # http://www.research.att.com/~njas/sequences/A000040 | ||
+ | _validate_int(n) | ||
+ | if n < 1: | ||
+ | raise ValueError('argument must be a positive integer') | ||
+ | return next(itertools.islice(primes(), n-1, None)) | ||
+ | |||
+ | |||
+ | def prime_count(x): | ||
+ | """prime_count(x) -> int | ||
+ | |||
+ | Returns the number of prime numbers less than or equal to x. | ||
+ | It is also known as the Prime Counting Function, or pi(x). | ||
+ | (Not to be confused with the constant pi = 3.1415....) | ||
+ | |||
+ | >>> prime_count(20) | ||
+ | 8 | ||
+ | >>> prime_count(10000) | ||
+ | 1229 | ||
+ | |||
+ | The number of primes less than x is approximately x/(ln x - 1). | ||
+ | """ | ||
+ | # See also: http://primes.utm.edu/howmany.shtml | ||
+ | # http://mathworld.wolfram.com/PrimeCountingFunction.html | ||
+ | _validate_num(x) | ||
+ | return sum(1 for p in primes_below(x)) | ||
+ | |||
+ | |||
+ | def primesum(n): | ||
+ | """primesum(n) -> int | ||
+ | |||
+ | primesum(n) returns the sum of the first n primes. | ||
+ | |||
+ | >>> primesum(9) | ||
+ | 100 | ||
+ | >>> primesum(49) | ||
+ | 4888 | ||
+ | |||
+ | The sum of the first n primes is approximately n**2*ln(n)/2. | ||
+ | """ | ||
+ | # See: http://mathworld.wolfram.com/PrimeSums.html | ||
+ | # http://www.research.att.com/~njas/sequences/A007504 | ||
+ | _validate_int(n) | ||
+ | return sum(nprimes(n)) | ||
+ | |||
+ | |||
+ | def primesums(): | ||
+ | """Yield the partial sums of the prime numbers. | ||
+ | |||
+ | >>> p = primesums() | ||
+ | >>> [next(p) for _ in range(5)] # primes 2, 3, 5, 7, 11, ... | ||
+ | [2, 5, 10, 17, 28] | ||
+ | |||
+ | """ | ||
+ | n = 0 | ||
+ | for p in primes(): | ||
+ | n += p | ||
+ | yield n | ||
+ | |||
+ | |||
+ | # ================= | ||
+ | # Primality testing | ||
+ | # ================= | ||
+ | |||
+ | # Set this to a true value to have isprime(n) warn if the result is | ||
+ | # probabilistic; set it to a false value to skip the warning. | ||
+ | warn_probably = True | ||
+ | |||
+ | |||
+ | def isprime(n): | ||
+ | """isprime(n) -> True|False | ||
+ | |||
+ | Returns True if integer n is prime number, otherwise return False. | ||
+ | |||
+ | For n less than approximately 341 trillion, ``isprime(n)`` is exact. Above | ||
+ | that value, it is probabilistic with a vanishingly small chance of wrongly | ||
+ | reporting a composite number as being prime. (It will never report a prime | ||
+ | as composite.) The probability of a false positive is less than 1/10**24, | ||
+ | or fewer than 1 time in a million trillion trillion tests. | ||
+ | |||
+ | If the global variable ``warn_probably`` is true (the default), isprime | ||
+ | will raise a warning when n is probably prime rather than certainly prime. | ||
+ | """ | ||
+ | _validate_int(n) | ||
+ | # Deal with trivial cases first. | ||
+ | if n < 2: | ||
+ | return False | ||
+ | elif n == 2: | ||
+ | return True | ||
+ | elif n%2 == 0: | ||
+ | return False | ||
+ | elif n <= 7: # 3, 5, 7 | ||
+ | return True | ||
+ | bases = _choose_bases(n) | ||
+ | flag = miller_rabin(n, bases) | ||
+ | if flag and len(bases) > 7 and warn_probably: | ||
+ | import warnings | ||
+ | warnings.warn("number is only probably prime not certainly prime") | ||
+ | return flag | ||
+ | |||
+ | |||
+ | def _choose_bases(n): | ||
+ | """Choose appropriate bases for the Miller-Rabin primality test. | ||
+ | |||
+ | If n is small enough, returns a tuple of bases which are provably | ||
+ | deterministic for that n. If n is too large, return a mostly random | ||
+ | selection of bases such that the chances of an error is less than | ||
+ | 1/4**40 = 8.2e-25. | ||
+ | """ | ||
+ | # The Miller-Rabin test is deterministic and completely accurate for | ||
+ | # moderate sizes of n using a surprisingly tiny number of tests. | ||
+ | # See: Pomerance, Selfridge and Wagstaff (1980), and Jaeschke (1993) | ||
+ | # http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test | ||
+ | if n < 1373653: # Blah, it's too hard to read big ints at a glance. | ||
+ | # ~1.3 million | ||
+ | bases = (2, 3) | ||
+ | elif n < 9080191: # ~9.0 million | ||
+ | bases = (31, 73) | ||
+ | elif n < 4759123141: # ~4.7 billion | ||
+ | # Note to self: checked up to approximately 394 million in 9 hours. | ||
+ | bases = (2, 7, 61) | ||
+ | elif n < 2152302898747: # ~2.1 trillion | ||
+ | bases = (2, 3, 5, 7, 11) | ||
+ | elif n < 3474749660383: # ~3.4 trillion | ||
+ | bases = (2, 3, 5, 7, 11, 13) | ||
+ | elif n < 341550071728321: # ~341 trillion | ||
+ | bases = (2, 3, 5, 7, 11, 13, 17) | ||
+ | else: | ||
+ | # n is too large, so we have to use a probabilistic test. There's no | ||
+ | # harm in trying some of the lower values for base first. | ||
+ | bases = (2, 3, 5, 7, 11, 13, 17) + tuple( | ||
+ | [random.randint(18, n-1) for _ in range(40)] | ||
+ | ) | ||
+ | # Note: we can always be deterministic, no matter how large N is, by | ||
+ | # exhaustive testing against each i in the inclusive range | ||
+ | # 1 ... min(n-1, floor(2*(ln N)**2)). We don't do this, because it is | ||
+ | # expensive for large N, and of no real practical benefit. | ||
+ | return bases | ||
+ | |||
+ | |||
+ | def isprime_naive(n): | ||
+ | """Naive test for primes. Returns True if int n is prime, otherwise False. | ||
+ | |||
+ | >>> isprime_naive(7) | ||
+ | True | ||
+ | >>> isprime_naive(8) | ||
+ | False | ||
+ | |||
+ | Naive, slow but thorough test for primality using unoptimized trial | ||
+ | division. This function does far too much work, and consequently is very | ||
+ | slow, but it is simple enough to verify by eye and can be used to check | ||
+ | the results of faster algorithms. (At least for very small n.) | ||
+ | """ | ||
+ | _validate_int(n) | ||
+ | if n == 2: return True | ||
+ | if n < 2 or n % 2 == 0: return False | ||
+ | for i in range(3, int(n**0.5)+1, 2): | ||
+ | if n % i == 0: | ||
+ | return False | ||
+ | return True | ||
+ | |||
+ | |||
+ | def isprime_division(n): | ||
+ | """isprime_division(integer) -> True|False | ||
+ | |||
+ | Exact primality test returning True if the argument is a prime number, | ||
+ | otherwise False. | ||
+ | |||
+ | >>> isprime_division(11) | ||
+ | True | ||
+ | >>> isprime_division(12) | ||
+ | False | ||
+ | |||
+ | This function uses trial division by the primes, skipping non-primes. | ||
+ | """ | ||
+ | _validate_int(n) | ||
+ | if n < 2: | ||
+ | return False | ||
+ | limit = n**0.5 | ||
+ | for divisor in primes(): | ||
+ | if divisor > limit: break | ||
+ | if n % divisor == 0: return False | ||
+ | return True | ||
+ | |||
+ | |||
+ | from re import match as _re_match | ||
+ | |||
+ | def isprime_regex(n): | ||
+ | """isprime_regex(n) -> True|False | ||
+ | |||
+ | Astonishingly, you can test whether a number is prime using a regex. | ||
+ | It goes without saying that this is not efficient, and should be treated | ||
+ | as a novelty rather than a serious implementation. It is O(N^2) in time | ||
+ | and O(N) in memory: in other words, slow and expensive. | ||
+ | """ | ||
+ | _validate_int(n) | ||
+ | return not _re_match(r'^1?$|^(11+?)\1+$', '1'*n) | ||
+ | # For a Perl version, see here: | ||
+ | # http://montreal.pm.org/tech/neil_kandalgaonkar.shtml | ||
+ | # And for a Ruby version, here: | ||
+ | # http://www.noulakaz.net/weblog/2007/03/18/a-regular-expression-to-check-for-prime-numbers/ | ||
+ | |||
+ | |||
+ | # === Probabilistic primality tests === | ||
+ | |||
+ | def fermat(n, base=2): | ||
+ | """fermat(n [, base]) -> True|False | ||
+ | |||
+ | ``fermat(n, base)`` is a probabilistic test for primality which returns | ||
+ | True if integer n is a weak probable prime to the given integer base, | ||
+ | otherwise n is definitely composite and False is returned. | ||
+ | |||
+ | ``base`` must be a positive integer between 1 and n-1 inclusive, or a | ||
+ | tuple of such bases. By default, base=2. | ||
+ | |||
+ | If ``fermat`` returns False, that is definite proof that n is composite: | ||
+ | there are no false negatives. However, if it returns True, that is only | ||
+ | provisional evidence that n is prime. For example: | ||
+ | |||
+ | >>> fermat(99, 7) | ||
+ | False | ||
+ | >>> fermat(29, 7) | ||
+ | True | ||
+ | |||
+ | We can conclude that 99 is definitely composite, and state that 7 is a | ||
+ | witness that 29 may be prime. | ||
+ | |||
+ | As the Fermat test is probabilistic, composite numbers will sometimes | ||
+ | pass a test, or even repeated tests: | ||
+ | |||
+ | >>> fermat(3*11*17, 7) # A pseudoprime to base 7. | ||
+ | True | ||
+ | |||
+ | You can perform multiple tests with a single call by passing a tuple of | ||
+ | ints as ``base``. The number must pass the Fermat test for all the bases | ||
+ | in order to return True. If any test fails, ``fermat`` will return False. | ||
+ | |||
+ | >>> fermat(41041, (17, 23, 356, 359)) # 41041 = 7*11*13*41 | ||
+ | True | ||
+ | >>> fermat(41041, (17, 23, 356, 359, 363)) | ||
+ | False | ||
+ | |||
+ | If a number passes ``k`` Fermat tests, we can conclude that the | ||
+ | probability that it is either a prime number, or a particular type of | ||
+ | pseudoprime known as a Carmichael number, is at least ``1 - (1/2**k)``. | ||
+ | """ | ||
+ | # http://en.wikipedia.org/wiki/Fermat_primality_test | ||
+ | _validate_int(n) | ||
+ | bases = _base_to_bases(base, n) | ||
+ | # Deal with the simple deterministic cases first. | ||
+ | if n < 2: | ||
+ | return False | ||
+ | elif n == 2: | ||
+ | return True | ||
+ | elif n % 2 == 0: | ||
+ | return False | ||
+ | # Now the Fermat test proper. | ||
+ | for a in bases: | ||
+ | if pow(a, n-1, n) != 1: | ||
+ | return False # n is certainly composite. | ||
+ | return True # All of the bases are witnesses for n being prime. | ||
+ | |||
+ | |||
+ | def miller_rabin(n, base=2): | ||
+ | """miller_rabin(integer [, base]) -> True|False | ||
+ | |||
+ | ``miller_rabin(n, base)`` is a probabilistic test for primality which | ||
+ | returns True if integer n is a strong probable prime to the given integer | ||
+ | base, otherwise n is definitely composite and False is returned. | ||
+ | |||
+ | ``base`` must be a positive integer between 1 and n-1 inclusive, or a | ||
+ | tuple of such bases. By default, base=2. | ||
+ | |||
+ | If ``miller_rabin`` returns False, that is definite proof that n is | ||
+ | composite: there are no false negatives. However, if it returns True, | ||
+ | that is only provisional evidence that n is prime: | ||
+ | |||
+ | >>> miller_rabin(99, 7) | ||
+ | False | ||
+ | >>> miller_rabin(29, 7) | ||
+ | True | ||
+ | |||
+ | We can conclude from this that 99 is definitely composite, and that 29 is | ||
+ | possibly prime. | ||
+ | |||
+ | As the Miller-Rabin test is probabilistic, composite numbers will | ||
+ | sometimes pass one or more tests: | ||
+ | |||
+ | >>> miller_rabin(3*11*17, 103) # 3*11*17=561, the 1st Carmichael number. | ||
+ | True | ||
+ | |||
+ | You can perform multiple tests with a single call by passing a tuple of | ||
+ | ints as ``base``. The number must pass the Miller-Rabin test for each of | ||
+ | the bases before it will return True. If any test fails, ``miller_rabin`` | ||
+ | will return False. | ||
+ | |||
+ | >>> miller_rabin(41041, (16, 92, 100, 256)) # 41041 = 7*11*13*41 | ||
+ | True | ||
+ | >>> miller_rabin(41041, (16, 92, 100, 256, 288)) | ||
+ | False | ||
+ | |||
+ | If a number passes ``k`` Miller-Rabin tests, we can conclude that the | ||
+ | probability that it is a prime number is at least ``1 - (1/4**k)``. | ||
+ | """ | ||
+ | # http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test | ||
+ | _validate_int(n) | ||
+ | bases = _base_to_bases(base, n) | ||
+ | # Deal with the trivial cases. | ||
+ | if n < 2: | ||
+ | return False | ||
+ | if n == 2: | ||
+ | return True | ||
+ | elif n % 2 == 0: | ||
+ | return False | ||
+ | # Now perform the Miller-Rabin test proper. | ||
+ | # Start by writing n-1 as 2**s * d. | ||
+ | d, s = _factor2(n-1) | ||
+ | for a in bases: | ||
+ | if _is_composite(a, d, s, n): | ||
+ | return False # n is definitely composite. | ||
+ | # If we get here, all of the bases are witnesses for n being prime. | ||
+ | return True | ||
+ | |||
+ | |||
+ | def _factor2(n): | ||
+ | """Factorise positive integer n as d*2**i, and return (d, i). | ||
+ | |||
+ | >>> _factor2(768) | ||
+ | (3, 8) | ||
+ | >>> _factor2(18432) | ||
+ | (9, 11) | ||
+ | |||
+ | Private function used internally by ``miller_rabin``. | ||
+ | """ | ||
+ | assert n > 0 and int(n) == n | ||
+ | i = 0 | ||
+ | d = n | ||
+ | while 1: | ||
+ | q, r = divmod(d, 2) | ||
+ | if r == 1: | ||
+ | break | ||
+ | i += 1 | ||
+ | d = q | ||
+ | assert d%2 == 1 | ||
+ | assert d*2**i == n | ||
+ | return (d, i) | ||
+ | |||
+ | |||
+ | def _is_composite(b, d, s, n): | ||
+ | """_is_composite(b, d, s, n) -> True|False | ||
+ | |||
+ | Tests base b to see if it is a witness for n being composite. Returns | ||
+ | True if n is definitely composite, otherwise False if it *may* be prime. | ||
+ | |||
+ | >>> _is_composite(4, 3, 7, 385) | ||
+ | True | ||
+ | >>> _is_composite(221, 3, 7, 385) | ||
+ | False | ||
+ | |||
+ | Private function used internally by ``miller_rabin``. | ||
+ | """ | ||
+ | assert d*2**s == n-1 | ||
+ | if pow(b, d, n) == 1: | ||
+ | return False | ||
+ | for i in range(s): | ||
+ | if pow(b, 2**i * d, n) == n-1: | ||
+ | return False | ||
+ | return True | ||
+ | |||
+ | |||
+ | # =================== | ||
+ | # Prime factorisation | ||
+ | # =================== | ||
+ | |||
+ | if __debug__: | ||
+ | # Set _EXTRA_CHECKS to True to enable potentially expensive assertions | ||
+ | # in the factors() and factorise() functions. This is only defined or | ||
+ | # checked when assertions are enabled. | ||
+ | _EXTRA_CHECKS = False | ||
+ | |||
+ | |||
+ | def factors(n): | ||
+ | """factors(integer) -> [list of factors] | ||
+ | |||
+ | Returns a list of the (mostly) prime factors of integer n. For negative | ||
+ | integers, -1 is included as a factor. If n is 0 or 1, [n] is returned as | ||
+ | the only factor. Otherwise all the factors will be prime. | ||
+ | |||
+ | >>> factors(-693) | ||
+ | [-1, 3, 3, 7, 11] | ||
+ | >>> factors(55614) | ||
+ | [2, 3, 13, 23, 31] | ||
+ | |||
+ | """ | ||
+ | _validate_int(n) | ||
+ | result = [] | ||
+ | for p, count in factorise(n): | ||
+ | result.extend([p]*count) | ||
+ | if __debug__: | ||
+ | # The following test only occurs if assertions are on. | ||
+ | if _EXTRA_CHECKS: | ||
+ | prod = 1 | ||
+ | for x in result: | ||
+ | prod *= x | ||
+ | assert prod == n, ('factors(%d) failed multiplication test' % n) | ||
+ | return result | ||
+ | |||
+ | |||
+ | def factorise(n): | ||
+ | """factorise(integer) -> yield factors of integer lazily | ||
+ | |||
+ | >>> list(factorise(3*7*7*7*11)) | ||
+ | [(3, 1), (7, 3), (11, 1)] | ||
+ | |||
+ | Yields tuples of (factor, count) where each factor is unique and usually | ||
+ | prime, and count is an integer 1 or larger. | ||
+ | |||
+ | The factors are prime, except under the following circumstances: if the | ||
+ | argument n is negative, -1 is included as a factor; if n is 0 or 1, it | ||
+ | is given as the only factor. For all other integer n, all of the factors | ||
+ | returned are prime. | ||
+ | """ | ||
+ | _validate_int(n) | ||
+ | if n in (0, 1, -1): | ||
+ | yield (n, 1) | ||
+ | return | ||
+ | elif n < 0: | ||
+ | yield (-1, 1) | ||
+ | n = -n | ||
+ | assert n >= 2 | ||
+ | for p in primes(): | ||
+ | if p*p > n: break | ||
+ | count = 0 | ||
+ | while n % p == 0: | ||
+ | count += 1 | ||
+ | n //= p | ||
+ | if count: | ||
+ | yield (p, count) | ||
+ | if n != 1: | ||
+ | if __debug__: | ||
+ | # The following test only occurs if assertions are on. | ||
+ | if _EXTRA_CHECKS: | ||
+ | assert isprime(n), ('failed isprime test for %d' % n) | ||
+ | yield (n, 1) | ||
+ | |||
+ | |||
+ | |||
+ | if __name__ == '__main__': | ||
+ | import doctest | ||
+ | doctest.testmod() | ||
+ | | ||
+ | </code> | ||