# Python俱乐部

## 链接

modules:pyprimes-source-code

# pyprimes 源代码

```#!/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

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

>>> 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 = arr = 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 = 
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://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
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)

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).
"""
# 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()
```
modules/pyprimes-source-code.txt · 最后更改: 2014/07/15 15:55 (外部编辑)

### 页面工具 