pyprimes.py

#!/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()