modules:pyprimes-source-code

这里会显示出您选择的修订版和当前版本之间的差别。

— |
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> | ||

modules/pyprimes-source-code.txt · 最后更改: 2014/07/15 15:55 (外部编辑)

除额外注明的地方外，本维基上的内容按下列许可协议发布： CC Attribution-Share Alike 4.0 International