+## 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()``

+ Yield prime numbers 2, 3, 5, 7, 11, ...

+ >>> [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()``:

+ Yield pairs of (isprime(i), i) for i=0,1,2,3,4,5...

+ 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:

+ ------------------ ----------------------------------------------------

+ ------------------ ----------------------------------------------------

+ 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.

+ ------------------ ----------------------------------------------------

+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:

+ Return True if n is prime, otherwise return False.

+Exact primality tests are:

+ Naive and slow trial division test for n being prime.

+ A less naive trial division test for n being prime.

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

+ >>> miller_rabin(14008971) # composite 3*947*4931

+These functions return or yield the prime factors of an integer.

+ Return a list of the prime factors of n.

+ Yield tuples (factor, count) for n.

+The ``factors(n)`` function lists repeated factors:

+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))

+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.

+ -------------- --------------------------------------------------------

+ -------------- --------------------------------------------------------

+ 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()``.

+from __future__ import division

+from re import match as _re_match

+__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',

+# ============================

+# Python 2.x/3.x compatibility

+# ============================

+# This module should support 2.5+, including Python 3.

+ # No next() builtin, so we're probably running Python 2.5.

+ # Use a simplified version (without support for default).

+ # 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

+ from itertools import ifilter as filter, izip as zip

+ # 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

+ from itertools import compress

+ # 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)

+ from math import isfinite

+ from math import isnan, isinf

+ # Python 2.5. Quick and dirty substitutes.

+ return not (isnan(x) or isinf(x))

+ """Raise an exception if obj is not an integer."""

+ m = int(obj + 0) # May raise TypeError, or OverflowError.

+ raise ValueError('expected an integer but got %r' % obj)

+ """Raise an exception if obj is not a finite real number."""

+ m = obj + 0 # May raise TypeError.

+ raise ValueError('expected a finite real number but got %r' % obj)

+def _base_to_bases(base, n):

+ if isinstance(base, tuple):

+ # 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

+ raise ValueError('base %d out of range 1...%d' % (b, n-1))

+# =======================

+# 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.)

+ """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.

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ >>> erat(10000) == list(primes_below(10000))

+ # Generate a fixed array of integers.

+ arr = list(range(n+1)) # A list is faster than an array

+ # Cross out 0 and 1 since they aren't prime.

+ # Cross out all the multiples of i starting from i**2.

+ for p in range(i*i, n+1, i):

+ # Advance to the next number not crossed off.

+ while i <= n and arr[i] is None:

+ return list(filter(None, arr))

+ """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

+ # Note: this explicit test is slightly faster than using

+ # prime = table.pop(i, None) and testing for None.

+ """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

+ # 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.

+ p = table[q]; del table[q] # Faster than pop.

+ while x in table or not (x & 1):

+ """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

+ # 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.

+ 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)

+ # 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.

+ while x in roots or (x % 30) not in primeroots:

+ """Generate prime numbers using wheel factorisation modulo 210."""

+ for i in (2, 3, 5, 7, 11):

+ # 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

+ found = [(11, 121)] # Smallest prime we care about, and its square.

+ for incr in itertools.cycle(spokes):

+ if p2 > i: # i must be a prime.

+ elif i % p == 0: # i must be composite.

+ 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.

+# === Algorithms to avoid ===

+ """Awful and naive prime functions namespace.

+ A collection of prime-related algorithms which are supplied for

+ educational purposes, as toys, curios, or as terrible warnings on

+ None of these methods have acceptable performance; they are barely

+ tolerable even for the first 100 primes.

+ # === Prime number generators ===

+ """Generate prime numbers naively, and REALLY slowly.

+ >>> p = Awful.naive_primes1()

+ >>> [next(p) for _ in range(10)]

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ This is about as awful as a straight-forward algorithm to generate

+ primes can get without deliberate pessimation. This algorithm does

+ not make even the most trivial optimizations:

+ - it tests all numbers as potential primes, whether odd or even,

+ instead of skipping even numbers apart from 2;

+ - it checks for primality by dividing against every number less

+ than the candidate prime itself, instead of stopping at the

+ square root of the candidate;

+ - it fails to bail out early when it finds a factor, instead

+ pointlessly keeps testing.

+ The result is that this is horribly slow.

+ if not composite: # It must be a prime.

+ """Generate prime numbers naively, and very slowly.

+ >>> p = Awful.naive_primes2()

+ >>> [next(p) for _ in range(10)]

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ This is a little better than ``naive_primes1``, but still horribly

+ slow. It makes a single optimization by using a short-circuit test

+ for primality testing: as soon as a factor is found, the candidate

+ is rejected immediately.

+ if all(i%p != 0 for p in range(2, i)):

+ """Generate prime numbers naively, and very slowly.

+ >>> p = Awful.naive_primes3()

+ >>> [next(p) for _ in range(10)]

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ This is an incremental improvement over ``naive_primes2`` by only

+ testing odd numbers as potential primes and factors.

+ if all(i%p != 0 for p in range(3, i, 2)):

+ """Generate prime numbers using a simple trial division algorithm.

+ >>> p = Awful.trial_division()

+ >>> [next(p) for _ in range(10)]

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ This is the first non-naive algorithm. Due to its simplicity, it may

+ perform acceptably for the first hundred or so primes, if your needs

+ are not very demanding. However, it does not scale well for large

+ This uses three optimizations:

+ - only test odd numbers for primality;

+ - only check against the prime factors already seen;

+ - stop checking at the square root of the number being tested.

+ With these three optimizations, we get asymptotic behaviour of

+ O(N*sqrt(N)/(log N)**2) where N is the number of primes found.

+ Despite these , this is still unacceptably slow, especially

+ as the list of memorised primes grows.

+ it = itertools.takewhile(lambda p, i=i: p*p <= i, primes)

+ if all(i%p != 0 for p in it):

+ """Generate prime numbers very slowly using Euler's sieve.

+ >>> [next(p) for _ in range(10)]

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ 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

+ sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p > 0]

+ This algorithm is sometimes wrongly described as the Sieve of

+ Eratosthenes, but it is not, it is a version of Euler's Sieve.

+ Although simple, it is extremely 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_primes1``.

+ O'Neill calls this the "Sleight on Eratosthenes".

+ # 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)

+ nums = filter(lambda v, p=prime: (v % p) != 0, nums)

+ # === Prime number testing ===

+ """Naive primality test using naive and unoptimized trial division.

+ >>> Awful.isprime_naive(17)

+ >>> Awful.isprime_naive(18)

+ 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.

+ if n < 2 or n % 2 == 0: return False

+ for i in range(3, int(n**0.5)+1, 2):

+ """Slow primality test using a regular expression.

+ >>> Awful.isprime_regex(11)

+ >>> Awful.isprime_regex(15)

+ Unsurprisingly, 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.

+ return not _re_match(r'^1?$|^(11+?)\1+$', '1'*n)

+ # For a Perl or Ruby version of this, see here:

+ # http://montreal.pm.org/tech/neil_kandalgaonkar.shtml

+ # http://www.noulakaz.net/weblog/2007/03/18/a-regular-expression-to-check-for-prime-numbers/

+ """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 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)]

+ _ = next(odd_primes) # Skip 2.

+ # Yield the non-primes between the previous prime and

+ for i in itertools.islice(itertools.count(prev + 2), 0, None, 2):

+ # And yield the current prime.

+ """Convenience function that yields the first n primes.

+ [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

+ return itertools.islice(primes(), n)

+ """Convenience function that yields primes strictly greater than x.

+ >>> next(primes_above(200))

+ # Consume the primes below x as fast as possible, then yield the rest.

+ """Convenience function yielding primes less than or equal to x.

+ >>> list(primes_below(20))

+ [2, 3, 5, 7, 11, 13, 17, 19]

+ 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.

+ # http://www.research.att.com/~njas/sequences/A000040

+ raise ValueError('argument must be a positive integer')

+ return next(itertools.islice(primes(), n-1, None))

+ """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....)

+ 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

+ return sum(1 for p in primes_below(x))

+ primesum(n) returns the sum of the first n primes.

+ 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

+ """Yield the partial sums of the prime numbers.

+ >>> [next(p) for _ in range(5)] # primes 2, 3, 5, 7, 11, ...

+def isprime(n, trials=25, warn=False):

+ """Return True if n is a prime number, and False if it is not.

+ ========== =======================================================

+ ========== =======================================================

+ n Number being tested for primality.

+ trials Count of primality tests to perform (default 25).

+ warn If true, warn on inexact results. (Default is false.)

+ ========== =======================================================

+ For values of ``n`` under approximately 341 trillion, this function is

+ exact and the arguments ``trials`` and ``warn`` are ignored.

+ Above this cut-off value, this function may be probabilistic with a small

+ chance of wrongly reporting a composite (non-prime) number as prime. Such

+ composite numbers wrongly reported as prime are "false positive" errors.

+ The argument ``trials`` controls the risk of a false positive error. The

+ larger number of trials, the less the chance of an error (and the slower

+ the function). With the default value of 25, you can expect roughly one

+ such error every million trillion tests, which in practical terms is

+ ``isprime`` cannot give a false negative error: if it reports a number is

+ composite, it is certainly composite, but if it reports a number is prime,

+ it may be only probably prime. If you pass a true value for argument

+ ``warn``, then a warning will be raised if the result is probabilistic.

+ # Deal with trivial cases first.

+ is_probabilistic, bases = _choose_bases(n, trials)

+ is_prime = miller_rabin(n, bases)

+ if is_prime and is_probabilistic and warn:

+ warnings.warn("number is only probably prime not certainly prime")

+def _choose_bases(n, count):

+ """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 selection of

+ With k distinct Miller-Rabin tests, the probability of a false

+ positive result is no more than 1/(4**k).

+ # 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: # ~1.3 million

+ elif n < 9080191: # ~9.0 million

+ elif n < 4759123141: # ~4.7 billion

+ # Note to self: checked up to approximately 394 million in 9 hours.

+ 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)

+ # n is sufficiently large that we have to use a probabilistic test.

+ bases = tuple([random.randint(2, n-1) for _ in range(count)])

+ # FIXME Because bases are chosen at random, there may be duplicates

+ # although with extremely small probability given the size of n.

+ # FIXME Is it worthwhile to special case some of the lower, easier

+ # bases? bases = [2, 3, 5, 7, 11, 13, 17] + [random... ]?

+ # 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.

+def isprime_division(n):

+ """isprime_division(integer) -> True|False

+ Exact primality test returning True if the argument is a prime number,

+ >>> isprime_division(11)

+ >>> isprime_division(12)

+ This function uses trial division by the primes, skipping non-primes.

+ for divisor in primes():

+ if divisor > limit: break

+ if n % divisor == 0: return False

+# === Probabilistic primality tests ===

+ """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:

+ 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.

+ 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

+ >>> fermat(41041, (17, 23, 356, 359, 363))

+ 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

+ bases = _base_to_bases(base, n)

+ # Deal with the simple deterministic cases first.

+ # Now the Fermat test proper.

+ 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)

+ >>> miller_rabin(29, 7)

+ We can conclude from this that 99 is definitely composite, and that 29 is

+ 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.

+ 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``

+ >>> miller_rabin(41041, (16, 92, 100, 256)) # 41041 = 7*11*13*41

+ >>> miller_rabin(41041, (16, 92, 100, 256, 288))

+ 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

+ bases = _base_to_bases(base, n)

+ # Deal with the trivial cases.

+ # Now perform the Miller-Rabin test proper.

+ # Start by writing n-1 as 2**s * d.

+ 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.

+ """Factorise positive integer n as d*2**i, and return (d, i).

+ Private function used internally by ``miller_rabin``.

+ assert n > 0 and int(n) == n

+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)

+ >>> _is_composite(221, 3, 7, 385)

+ Private function used internally by ``miller_rabin``.

+ if pow(b, 2**i * d, n) == n-1:

+ # 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.

+ """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.

+ for p, count in factorise(n):

+ result.extend([p]*count)

+ # The following test only occurs if assertions are on.

+ assert prod == n, ('factors(%d) failed multiplication test' % 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

+ # The following test only occurs if assertions are on.

+ assert isprime(n), ('failed isprime test for %d' % n)

+if __name__ == '__main__':