用户工具

站点工具


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 (外部编辑)