[pypy-commit] benchmarks default: merge
fijal
noreply at buildbot.pypy.org
Fri Oct 24 15:55:40 CEST 2014
Author: Maciej Fijalkowski <fijall at gmail.com>
Branch:
Changeset: r278:0e250ed49b43
Date: 2014-10-24 15:55 +0200
http://bitbucket.org/pypy/benchmarks/changeset/0e250ed49b43/
Log: merge
diff --git a/benchmarks.py b/benchmarks.py
--- a/benchmarks.py
+++ b/benchmarks.py
@@ -91,7 +91,7 @@
else:
iteration_scaling = 1.0
_register_new_bm_twisted(name, 'twisted_' + name,
- globals(), bm_env={'PYTHONPATH': ':'.join(TWISTED)},
+ globals(), bm_env={'PYTHONPATH': os.pathsep.join(TWISTED)},
iteration_scaling=iteration_scaling)
_register_new_bm('spitfire', 'spitfire', globals(),
diff --git a/multithread/bench.py b/multithread/bench.py
--- a/multithread/bench.py
+++ b/multithread/bench.py
@@ -86,6 +86,7 @@
finally:
if not args.q:
print "== times ==\n", "\n".join(map(str, times))
+ print "== times short ==\n", str(times[-min(5, len(times)):])
print "== reported results ==\n", "\n".join(
map(str, filter(None, results)))
diff --git a/multithread/common/abstract_threading.py b/multithread/common/abstract_threading.py
--- a/multithread/common/abstract_threading.py
+++ b/multithread/common/abstract_threading.py
@@ -4,7 +4,7 @@
try:
from atomic import (atomic, getsegmentlimit, print_abort_info,
- hint_commit_soon)
+ hint_commit_soon, is_atomic)
except:
atomic = RLock()
def getsegmentlimit():
@@ -13,6 +13,8 @@
pass
def hint_commit_soon():
pass
+ def is_atomic():
+ return atomic._RLock__count > 0
class TLQueue_concurrent(object):
@@ -156,9 +158,10 @@
def __call__(self):
- with self._cond:
- while not self._done:
- self._cond.wait()
+ if not self._done:
+ with self._cond:
+ while not self._done:
+ self._cond.wait()
if self._exception:
raise self._exception
@@ -171,8 +174,10 @@
def _task(self, func, *args, **kwargs):
with self._cond:
try:
+ hint_commit_soon()
with atomic:
self._result = func(*args, **kwargs)
+ hint_commit_soon()
except Exception as e:
self._exception = e
finally:
diff --git a/multithread/mersenne/mersenne.py b/multithread/mersenne/mersenne.py
new file mode 100644
--- /dev/null
+++ b/multithread/mersenne/mersenne.py
@@ -0,0 +1,89 @@
+# -*- coding: utf-8 -*-
+
+
+
+import sys
+import time, random
+from common.abstract_threading import (
+ atomic, Future, set_thread_pool, ThreadPool,
+ hint_commit_soon, print_abort_info)
+
+import itertools
+from collections import deque
+
+
+from sys import stdout
+from math import sqrt, log
+
+
+def chunks(iterable, size):
+ it = iter(iterable)
+ item = list(itertools.islice(it, size))
+ while item:
+ yield item
+ item = list(itertools.islice(it, size))
+
+
+def is_prime ( p ):
+ if p == 2: return True # Lucas-Lehmer test only works on odd primes
+ elif p <= 1 or p % 2 == 0: return False
+ else:
+ for i in range(3, int(sqrt(p))+1, 2 ):
+ if p % i == 0: return False
+ return True
+
+def is_mersenne_prime ( p ):
+ if p == 2:
+ return True
+ else:
+ m_p = ( 1 << p ) - 1
+ s = 4
+ for i in range(3, p+1):
+ s = (s ** 2 - 2) % m_p
+ return s == 0
+
+
+def work(ps, counter, upb_count):
+ if counter[0] >= upb_count:
+ return
+
+ for p in ps:
+ with atomic:
+ if is_prime(p) and is_mersenne_prime(p):
+ #print p
+ counter[0] += 1
+ if counter[0] >= upb_count:
+ break
+
+
+def run(threads=2, n=2000):
+ threads = int(threads)
+ n = int(n)
+
+ set_thread_pool(ThreadPool(threads))
+
+
+ precision = n # maximum requested number of decimal places of 2 ** MP-1 #
+ long_bits_width = precision * log(10, 2)
+ upb_prime = int( long_bits_width - 1 ) / 2 # no unsigned #
+ upb_count = 45 # find 45 mprimes if int was given enough bits #
+
+ print " Finding Mersenne primes in M[2..%d]:" % upb_prime
+
+ counter = [0]
+ fs = []
+ for ps in chunks(xrange(2, upb_prime+1), 500):
+ fs.append(Future(work, ps, counter, upb_count))
+
+ [f() for f in fs]
+ print "found", counter[0]
+
+ # shutdown current pool
+ set_thread_pool(None)
+
+ return
+
+
+
+if __name__ == "__main__":
+ run()
diff --git a/multithread/nqueens/nqueens.py b/multithread/nqueens/nqueens.py
new file mode 100644
--- /dev/null
+++ b/multithread/nqueens/nqueens.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+# Based on code from http://code.activestate.com/recipes/576647/
+# by Raymond Hettinger
+#
+# It is a bit problematic because find_solutions spends half
+# of the execution time constructing the permutations. Thus
+# only half the execution is parallelised.
+
+
+import sys
+import time, random
+from common.abstract_threading import (
+ atomic, Future, set_thread_pool, ThreadPool,
+ hint_commit_soon)
+
+from itertools import permutations
+import itertools
+
+
+def chunks(iterable, size):
+ it = iter(iterable)
+ item = list(itertools.islice(it, size))
+ while item:
+ yield item
+ item = list(itertools.islice(it, size))
+
+
+
+def check_solutions(n, cols, perms):
+ sols = []
+ with atomic:
+ for vec in perms:
+ if n == len(set(vec[i]+i for i in cols)) \
+ == len(set(vec[i]-i for i in cols)):
+ sols.append(vec)
+ return sols
+
+
+def find_solutions(n):
+ solutions = []
+ fs = []
+ cols = range(n)
+ for perms in chunks(permutations(cols), 100000):
+ hint_commit_soon()
+ fs.append(Future(check_solutions, n, cols, perms))
+ print "Futures:", len(fs)
+ for f in fs:
+ solutions.extend(f())
+
+ print "found:", len(solutions)
+
+
+def run(threads=2, n=10):
+ threads = int(threads)
+ n = int(n)
+
+ set_thread_pool(ThreadPool(threads))
+
+ find_solutions(n)
+
+
+ # shutdown current pool
+ set_thread_pool(None)
+
+
+
+if __name__ == "__main__":
+ run()
diff --git a/multithread/perlin_noise/perlin_noise.py b/multithread/perlin_noise/perlin_noise.py
new file mode 100644
--- /dev/null
+++ b/multithread/perlin_noise/perlin_noise.py
@@ -0,0 +1,95 @@
+# -*- coding: utf-8 -*-
+# from http://rosettacode.org/wiki/Perlin_noise#Python
+
+import sys
+import time, random
+from common.abstract_threading import (AtomicFuture,
+ atomic, Future, set_thread_pool, ThreadPool,
+ hint_commit_soon, print_abort_info)
+
+import itertools
+from collections import deque
+
+p = [None] * 512
+permutation = [151,160,137,91,90,15,
+ 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
+ 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
+ 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
+ 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
+ 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
+ 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
+ 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
+ 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
+ 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
+ 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
+ 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
+ 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180]
+
+for i in range(256):
+ p[256+i] = p[i] = permutation[i]
+
+def perlin_noise(x, y, z):
+ X = int(x) & 255 # FIND UNIT CUBE THAT
+ Y = int(y) & 255 # CONTAINS POINT.
+ Z = int(z) & 255
+ x -= int(x) # FIND RELATIVE X,Y,Z
+ y -= int(y) # OF POINT IN CUBE.
+ z -= int(z)
+ u = fade(x) # COMPUTE FADE CURVES
+ v = fade(y) # FOR EACH OF X,Y,Z.
+ w = fade(z)
+ A = p[X ]+Y; AA = p[A]+Z; AB = p[A+1]+Z # HASH COORDINATES OF
+ B = p[X+1]+Y; BA = p[B]+Z; BB = p[B+1]+Z # THE 8 CUBE CORNERS,
+
+ return lerp(w, lerp(v, lerp(u, grad(p[AA ], x , y , z ), # AND ADD
+ grad(p[BA ], x-1, y , z )), # BLENDED
+ lerp(u, grad(p[AB ], x , y-1, z ), # RESULTS
+ grad(p[BB ], x-1, y-1, z ))),# FROM 8
+ lerp(v, lerp(u, grad(p[AA+1], x , y , z-1 ), # CORNERS
+ grad(p[BA+1], x-1, y , z-1 )), # OF CUBE
+ lerp(u, grad(p[AB+1], x , y-1, z-1 ),
+ grad(p[BB+1], x-1, y-1, z-1 ))))
+
+def fade(t):
+ return t ** 3 * (t * (t * 6 - 15) + 10)
+
+def lerp(t, a, b):
+ return a + t * (b - a)
+
+def grad(hash, x, y, z):
+ h = hash & 15 # CONVERT LO 4 BITS OF HASH CODE
+ u = x if h<8 else y # INTO 12 GRADIENT DIRECTIONS.
+ v = y if h<4 else (x if h in (12, 14) else z)
+ return (u if (h&1) == 0 else -u) + (v if (h&2) == 0 else -v)
+
+
+
+
+def work(n, x):
+ res = []
+ hint_commit_soon()
+ with atomic:
+ for y in xrange(n):
+ for z in xrange(n):
+ res.append(perlin_noise(x/3., y/3., z/3.))
+ hint_commit_soon()
+ return res
+
+
+def run(threads=2, n=50):
+ threads = int(threads)
+ n = int(n)
+
+ set_thread_pool(ThreadPool(threads))
+
+ res = []
+ for x in range(n):
+ res.append(Future(work, 400, x))
+ res = [f() for f in res]
+
+ # shutdown current pool
+ set_thread_pool(None)
+
+
+if __name__ == "__main__":
+ run()
diff --git a/multithread/primes/primes.py b/multithread/primes/primes.py
new file mode 100644
--- /dev/null
+++ b/multithread/primes/primes.py
@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+
+# from https://github.com/Tinche/stm-playground
+
+
+import sys
+import time, random
+from common.abstract_threading import (
+ atomic, Future, set_thread_pool, ThreadPool,
+ hint_commit_soon, print_abort_info)
+
+from itertools import izip, chain, repeat
+
+from Queue import Queue
+from pyprimes import isprime
+import threading
+
+def check_prime(num):
+ return isprime(num), num
+
+
+def grouper(n, iterable, padvalue=None):
+ "grouper(3, 'abcdefg', 'x') --> ('a','b','c'), ('d','e','f'), ('g','x','x')"
+ return izip(*[chain(iterable, repeat(padvalue, n-1))]*n)
+
+
+poison_pill = object()
+
+def worker(tasks, results):
+ while True:
+ batch = tasks.get()
+ if batch is poison_pill:
+ tasks.task_done()
+ return
+
+ result = []
+ for task in batch:
+ with atomic:
+ result.append(check_prime(task))
+ results.put(result)
+
+ tasks.task_done()
+
+
+
+def run(threads=2, n=2000000):
+ threads = int(threads)
+ n = int(n)
+
+ LIMIT = n
+ BATCH_SIZE = 1000
+
+ tasks = Queue()
+ results = Queue()
+ print("Starting...")
+
+ with atomic:
+ for batch in grouper(BATCH_SIZE, xrange(LIMIT), 1):
+ tasks.put(list(batch))
+ for _ in xrange(threads):
+ tasks.put(poison_pill)
+
+ for _ in xrange(threads):
+ t = threading.Thread(target=worker, args=(tasks, results))
+ t.start()
+ tasks.join()
+
+ count = 0
+ while not results.empty():
+ batch_results = results.get()
+ count += sum(1 for res in batch_results if res[0])
+
+ return count
+
+
+
+if __name__ == "__main__":
+ run()
diff --git a/multithread/primes/pyprimes.py b/multithread/primes/pyprimes.py
new file mode 100644
--- /dev/null
+++ b/multithread/primes/pyprimes.py
@@ -0,0 +1,1193 @@
+#!/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()``.
+
+
+"""
+
+
+from __future__ import division
+
+
+import functools
+import itertools
+import random
+
+from re import match as _re_match
+
+
+# Module metadata.
+__version__ = "0.1.2a"
+__date__ = "2012-08-25"
+__author__ = "Steven D'Aprano"
+__author_email__ = "steve+python at 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, or OverflowError.
+ 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)) # A list is faster than an array
+ # 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 ===
+
+class Awful:
+ """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
+ what **not** to do.
+
+ None of these methods have acceptable performance; they are barely
+ tolerable even for the first 100 primes.
+ """
+
+ # === Prime number generators ===
+
+ @staticmethod
+ def naive_primes1():
+ """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.
+ """
+ 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
+
+ @staticmethod
+ def naive_primes2():
+ """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.
+ """
+ i = 2
+ yield i
+ while True:
+ i += 1
+ if all(i%p != 0 for p in range(2, i)):
+ yield i
+
+ @staticmethod
+ def naive_primes3():
+ """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.
+ """
+ yield 2
+ i = 3
+ yield i
+ while True:
+ i += 2
+ if all(i%p != 0 for p in range(3, i, 2)):
+ yield i
+
+ @staticmethod
+ def trial_division():
+ """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
+ numbers of primes.
+
+ 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.
+ """
+ yield 2
+ primes = [2]
+ i = 3
+ while True:
+ 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
+
+ @staticmethod
+ def turner():
+ """Generate prime numbers very slowly using Euler's sieve.
+
+ >>> p = Awful.turner()
+ >>> [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
+ variation of::
+
+ primes = sieve [2..]
+ 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".
+ """
+ # 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)
+
+ # === Prime number testing ===
+
+ @staticmethod
+ def isprime_naive(n):
+ """Naive primality test using naive and unoptimized trial division.
+
+ >>> Awful.isprime_naive(17)
+ True
+ >>> Awful.isprime_naive(18)
+ 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.
+ """
+ _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
+
+ @staticmethod
+ def isprime_regex(n):
+ """Slow primality test using a regular expression.
+
+ >>> Awful.isprime_regex(11)
+ True
+ >>> Awful.isprime_regex(15)
+ False
+
+ 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.
+ """
+ _validate_int(n)
+ 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/
+
+
+
+# =====================
+# 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
+# =================
+
+def isprime(n, trials=25, warn=False):
+ """Return True if n is a prime number, and False if it is not.
+
+ >>> isprime(101)
+ True
+ >>> isprime(102)
+ False
+
+ ========== =======================================================
+ Argument Description
+ ========== =======================================================
+ 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
+ essentially "never".
+
+ ``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.
+ """
+ _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
+ is_probabilistic, bases = _choose_bases(n, trials)
+ is_prime = miller_rabin(n, bases)
+ if is_prime and is_probabilistic and warn:
+ import warnings
+ warnings.warn("number is only probably prime not certainly prime")
+ return is_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
+ possibly random bases.
+
+ 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
+ prob = False
+ if n < 1373653: # ~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 sufficiently large that we have to use a probabilistic test.
+ prob = True
+ 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.
+ return prob, bases
+
+
+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
+
+
+# === 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()
+
diff --git a/multithread/quick_sort/quick_sort.py b/multithread/quick_sort/quick_sort.py
new file mode 100644
--- /dev/null
+++ b/multithread/quick_sort/quick_sort.py
@@ -0,0 +1,129 @@
+# -*- coding: utf-8 -*-
+
+
+
+import sys
+import time, random
+from common.abstract_threading import (
+ atomic, Future, set_thread_pool, ThreadPool,
+ hint_commit_soon, print_abort_info)
+
+
+import itertools
+from collections import deque
+
+
+def chunks(iterable, size):
+ it = iter(iterable)
+ item = list(itertools.islice(it, size))
+ while item:
+ yield item
+ item = list(itertools.islice(it, size))
+
+
+
+def qsort(xs, l0, n):
+ if n < 2:
+ return
+ pivot = xs[l0 + n // 2]
+ l = l0
+ r = l + n - 1
+ while l <= r:
+ if xs[l] < pivot:
+ l += 1
+ continue
+ if xs[r] > pivot:
+ r -= 1
+ continue
+ xs[l], xs[r] = xs[r], xs[l]
+ l += 1
+ r -= 1
+ qsort(xs, l0, r - l0 + 1)
+ qsort(xs, l, l0 + n - l)
+
+
+
+def qsort_f(xs, l0, n, level):
+ if n < 2:
+ return []
+
+ pivot = xs[l0 + n // 2]
+ l = l0
+ r = l + n - 1
+ while l <= r:
+ xl = xs[l]
+ if xl < pivot:
+ l += 1
+ continue
+ xr = xs[r]
+ if xr > pivot:
+ r -= 1
+ continue
+ xs[l], xs[r] = xr, xl
+ l += 1
+ r -= 1
+
+ fs = []
+ do_futures = level >= 4
+ largs = (xs, l0, r - l0 + 1, level+1)
+ rargs = (xs, l, l0 + n - l, level+1)
+ leftf, rightf = False, False
+
+ if do_futures:
+ if largs[2] > 1500:
+ fs.append(Future(qsort_f, *largs))
+ leftf = True
+
+ if rargs[2] > 1500:
+ fs.append(Future(qsort_f, *rargs))
+ rightf = True
+
+ if not leftf:
+ fs.extend(qsort_f(*largs))
+
+ if not rightf:
+ fs.extend(qsort_f(*rargs))
+ #print_abort_info(0.0000001)
+
+ return fs
+
+
+def wait_for_futures(fs):
+ c = 0
+ while fs:
+ f = fs.pop()
+ fs.extend(f())
+ c += 1
+ print "Futures:", c
+
+def run(threads=2, n=60000):
+ threads = int(threads)
+ n = int(n)
+
+ set_thread_pool(ThreadPool(threads))
+
+ to_sort = range(n)
+ t = 0
+ for i in range(5):
+ random.seed(i+32)
+ random.shuffle(to_sort)
+ s = deque(to_sort)
+ # qsort(s, 0, len(s))
+
+ t -= time.time()
+ # start as future, otherwise we get more threads
+ # than we want (+1 for the main thread)
+ fs = Future(qsort_f, s, 0, len(s), 0)
+ wait_for_futures(fs())
+ #assert sorted(to_sort) == list(s)
+ t += time.time()
+
+ # shutdown current pool
+ set_thread_pool(None)
+
+ return t
+
+
+
+if __name__ == "__main__":
+ run()
diff --git a/multithread/raytrace/raytrace.py b/multithread/raytrace/raytrace.py
--- a/multithread/raytrace/raytrace.py
+++ b/multithread/raytrace/raytrace.py
@@ -7,15 +7,15 @@
print_abort_info, hint_commit_soon)
import time
-import platform
-if platform.python_implementation() == "Jython":
- # be fair to jython and don't use a lock where none is required:
- class fakeatomic:
- def __enter__(self):
- pass
- def __exit__(self,*args):
- pass
- atomic = fakeatomic()
+#import platform
+#if platform.python_implementation() == "Jython":
+# # be fair to jython and don't use a lock where none is required:
+# class fakeatomic:
+# def __enter__(self):
+# pass
+# def __exit__(self,*args):
+# pass
+# atomic = fakeatomic()
AMBIENT = 0.1
diff --git a/own/pyflate-fast.py b/own/pyflate-fast.py
--- a/own/pyflate-fast.py
+++ b/own/pyflate-fast.py
@@ -645,7 +645,7 @@
def _main():
filename = os.path.join(os.path.dirname(__file__), "interpreter.tar.bz2")
- input = open(filename)
+ input = open(filename, 'rb')
field = RBitfield(input)
magic = field.readbits(16)
diff --git a/unladen_swallow/perf.py b/unladen_swallow/perf.py
--- a/unladen_swallow/perf.py
+++ b/unladen_swallow/perf.py
@@ -969,7 +969,7 @@
def MeasureRietveld(python, options):
- PYTHONPATH = ":".join([DJANGO_DIR,
+ PYTHONPATH = os.pathsep.join([DJANGO_DIR,
# These paths are lifted from
# lib/google_appengine.appcfg.py. Note that we use
# our own version of Django instead of Appengine's.
@@ -1406,7 +1406,7 @@
time it took to run the benchmark loop once; mem_usage is a list
of memory usage samples in kilobytes.
"""
- pypath = ":".join([Relative("lib/spambayes"), Relative("lib/lockfile")])
+ pypath = os.pathsep.join([Relative("lib/spambayes"), Relative("lib/lockfile")])
bm_path = Relative("performance/bm_spambayes.py")
bm_env = {"PYTHONPATH": pypath}
return MeasureGeneric(python, options, bm_path, bm_env)
More information about the pypy-commit
mailing list