tests/perf_bench: Add some benchmarks from python-performance.

From https://github.com/python/pyperformance commit
6690642ddeda46fc5ee6e97c3ef4b2f292348ab8

6 files changed, 1182 insertions, 0 deletions

diff --git a/tests/perf_bench/bm_chaos.py b/tests/perf_bench/bm_chaos.py
new file mode 100644
index 000000000..04e04531c
--- /dev/null
+++ b/tests/perf_bench/bm_chaos.py
@@ -0,0 +1,274 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# create chaosgame-like fractals
+# Copyright (C) 2005 Carl Friedrich Bolz
+
+import math
+import random
+
+
+class GVector(object):
+
+    def __init__(self, x=0, y=0, z=0):
+        self.x = x
+        self.y = y
+        self.z = z
+
+    def Mag(self):
+        return math.sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2)
+
+    def dist(self, other):
+        return math.sqrt((self.x - other.x) ** 2
+                         + (self.y - other.y) ** 2
+                         + (self.z - other.z) ** 2)
+
+    def __add__(self, other):
+        if not isinstance(other, GVector):
+            raise ValueError("Can't add GVector to " + str(type(other)))
+        v = GVector(self.x + other.x, self.y + other.y, self.z + other.z)
+        return v
+
+    def __sub__(self, other):
+        return self + other * -1
+
+    def __mul__(self, other):
+        v = GVector(self.x * other, self.y * other, self.z * other)
+        return v
+    __rmul__ = __mul__
+
+    def linear_combination(self, other, l1, l2=None):
+        if l2 is None:
+            l2 = 1 - l1
+        v = GVector(self.x * l1 + other.x * l2,
+                    self.y * l1 + other.y * l2,
+                    self.z * l1 + other.z * l2)
+        return v
+
+    def __str__(self):
+        return "<%f, %f, %f>" % (self.x, self.y, self.z)
+
+    def __repr__(self):
+        return "GVector(%f, %f, %f)" % (self.x, self.y, self.z)
+
+
+class Spline(object):
+    """Class for representing B-Splines and NURBS of arbitrary degree"""
+
+    def __init__(self, points, degree, knots):
+        """Creates a Spline.
+
+        points is a list of GVector, degree is the degree of the Spline.
+        """
+        if len(points) > len(knots) - degree + 1:
+            raise ValueError("too many control points")
+        elif len(points) < len(knots) - degree + 1:
+            raise ValueError("not enough control points")
+        last = knots[0]
+        for cur in knots[1:]:
+            if cur < last:
+                raise ValueError("knots not strictly increasing")
+            last = cur
+        self.knots = knots
+        self.points = points
+        self.degree = degree
+
+    def GetDomain(self):
+        """Returns the domain of the B-Spline"""
+        return (self.knots[self.degree - 1],
+                self.knots[len(self.knots) - self.degree])
+
+    def __call__(self, u):
+        """Calculates a point of the B-Spline using de Boors Algorithm"""
+        dom = self.GetDomain()
+        if u < dom[0] or u > dom[1]:
+            raise ValueError("Function value not in domain")
+        if u == dom[0]:
+            return self.points[0]
+        if u == dom[1]:
+            return self.points[-1]
+        I = self.GetIndex(u)
+        d = [self.points[I - self.degree + 1 + ii]
+             for ii in range(self.degree + 1)]
+        U = self.knots
+        for ik in range(1, self.degree + 1):
+            for ii in range(I - self.degree + ik + 1, I + 2):
+                ua = U[ii + self.degree - ik]
+                ub = U[ii - 1]
+                co1 = (ua - u) / (ua - ub)
+                co2 = (u - ub) / (ua - ub)
+                index = ii - I + self.degree - ik - 1
+                d[index] = d[index].linear_combination(d[index + 1], co1, co2)
+        return d[0]
+
+    def GetIndex(self, u):
+        dom = self.GetDomain()
+        for ii in range(self.degree - 1, len(self.knots) - self.degree):
+            if u >= self.knots[ii] and u < self.knots[ii + 1]:
+                I = ii
+                break
+        else:
+            I = dom[1] - 1
+        return I
+
+    def __len__(self):
+        return len(self.points)
+
+    def __repr__(self):
+        return "Spline(%r, %r, %r)" % (self.points, self.degree, self.knots)
+
+
+def write_ppm(im, w, h, filename):
+    with open(filename, "wb") as f:
+        f.write(b'P6\n%i %i\n255\n' % (w, h))
+        for j in range(h):
+            for i in range(w):
+                val = im[j * w + i]
+                c = val * 255
+                f.write(b'%c%c%c' % (c, c, c))
+
+
+class Chaosgame(object):
+
+    def __init__(self, splines, thickness, subdivs):
+        self.splines = splines
+        self.thickness = thickness
+        self.minx = min([p.x for spl in splines for p in spl.points])
+        self.miny = min([p.y for spl in splines for p in spl.points])
+        self.maxx = max([p.x for spl in splines for p in spl.points])
+        self.maxy = max([p.y for spl in splines for p in spl.points])
+        self.height = self.maxy - self.miny
+        self.width = self.maxx - self.minx
+        self.num_trafos = []
+        maxlength = thickness * self.width / self.height
+        for spl in splines:
+            length = 0
+            curr = spl(0)
+            for i in range(1, subdivs + 1):
+                last = curr
+                t = 1 / subdivs * i
+                curr = spl(t)
+                length += curr.dist(last)
+            self.num_trafos.append(max(1, int(length / maxlength * 1.5)))
+        self.num_total = sum(self.num_trafos)
+
+    def get_random_trafo(self):
+        r = random.randrange(int(self.num_total) + 1)
+        l = 0
+        for i in range(len(self.num_trafos)):
+            if r >= l and r < l + self.num_trafos[i]:
+                return i, random.randrange(self.num_trafos[i])
+            l += self.num_trafos[i]
+        return len(self.num_trafos) - 1, random.randrange(self.num_trafos[-1])
+
+    def transform_point(self, point, trafo=None):
+        x = (point.x - self.minx) / self.width
+        y = (point.y - self.miny) / self.height
+        if trafo is None:
+            trafo = self.get_random_trafo()
+        start, end = self.splines[trafo[0]].GetDomain()
+        length = end - start
+        seg_length = length / self.num_trafos[trafo[0]]
+        t = start + seg_length * trafo[1] + seg_length * x
+        basepoint = self.splines[trafo[0]](t)
+        if t + 1 / 50000 > end:
+            neighbour = self.splines[trafo[0]](t - 1 / 50000)
+            derivative = neighbour - basepoint
+        else:
+            neighbour = self.splines[trafo[0]](t + 1 / 50000)
+            derivative = basepoint - neighbour
+        if derivative.Mag() != 0:
+            basepoint.x += derivative.y / derivative.Mag() * (y - 0.5) * \
+                self.thickness
+            basepoint.y += -derivative.x / derivative.Mag() * (y - 0.5) * \
+                self.thickness
+        else:
+            # can happen, especially with single precision float
+            pass
+        self.truncate(basepoint)
+        return basepoint
+
+    def truncate(self, point):
+        if point.x >= self.maxx:
+            point.x = self.maxx
+        if point.y >= self.maxy:
+            point.y = self.maxy
+        if point.x < self.minx:
+            point.x = self.minx
+        if point.y < self.miny:
+            point.y = self.miny
+
+    def create_image_chaos(self, w, h, iterations, rng_seed):
+        # Always use the same sequence of random numbers
+        # to get reproductible benchmark
+        random.seed(rng_seed)
+
+        im = bytearray(w * h)
+        point = GVector((self.maxx + self.minx) / 2,
+                        (self.maxy + self.miny) / 2, 0)
+        for _ in range(iterations):
+            point = self.transform_point(point)
+            x = (point.x - self.minx) / self.width * w
+            y = (point.y - self.miny) / self.height * h
+            x = int(x)
+            y = int(y)
+            if x == w:
+                x -= 1
+            if y == h:
+                y -= 1
+            im[(h - y - 1) * w + x] = 1
+
+        return im
+
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (100, 50): (0.25, 100, 50, 50, 50, 1234),
+    (1000, 1000): (0.25, 200, 400, 400, 1000, 1234),
+    (5000, 1000): (0.25, 400, 500, 500, 7000, 1234),
+}
+
+def bm_setup(params):
+    splines = [
+        Spline([
+            GVector(1.597, 3.304, 0.0),
+            GVector(1.576, 4.123, 0.0),
+            GVector(1.313, 5.288, 0.0),
+            GVector(1.619, 5.330, 0.0),
+            GVector(2.890, 5.503, 0.0),
+            GVector(2.373, 4.382, 0.0),
+            GVector(1.662, 4.360, 0.0)],
+            3, [0, 0, 0, 1, 1, 1, 2, 2, 2]),
+        Spline([
+            GVector(2.805, 4.017, 0.0),
+            GVector(2.551, 3.525, 0.0),
+            GVector(1.979, 2.620, 0.0),
+            GVector(1.979, 2.620, 0.0)],
+            3, [0, 0, 0, 1, 1, 1]),
+        Spline([
+            GVector(2.002, 4.011, 0.0),
+            GVector(2.335, 3.313, 0.0),
+            GVector(2.367, 3.233, 0.0),
+            GVector(2.367, 3.233, 0.0)],
+            3, [0, 0, 0, 1, 1, 1])
+    ]
+
+    chaos = Chaosgame(splines, params[0], params[1])
+    image = None
+
+    def run():
+        nonlocal image
+        _, _, width, height, iter, rng_seed = params
+        image = chaos.create_image_chaos(width, height, iter, rng_seed)
+
+    def result():
+        norm = params[4]
+        # Images are not the same when floating point behaviour is different,
+        # so return percentage of pixels that are set (rounded to int).
+        #write_ppm(image, params[2], params[3], 'out-.ppm')
+        pix = int(100 * sum(image) / len(image))
+        return norm, pix
+
+    return run, result
diff --git a/tests/perf_bench/bm_fannkuch.py b/tests/perf_bench/bm_fannkuch.py
new file mode 100644
index 000000000..c9782e3e9
--- /dev/null
+++ b/tests/perf_bench/bm_fannkuch.py
@@ -0,0 +1,67 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# The Computer Language Benchmarks Game
+# http://benchmarksgame.alioth.debian.org/
+# Contributed by Sokolov Yura, modified by Tupteq.
+
+def fannkuch(n):
+    count = list(range(1, n + 1))
+    max_flips = 0
+    m = n - 1
+    r = n
+    check = 0
+    perm1 = list(range(n))
+    perm = list(range(n))
+    perm1_ins = perm1.insert
+    perm1_pop = perm1.pop
+
+    while 1:
+        if check < 30:
+            check += 1
+
+        while r != 1:
+            count[r - 1] = r
+            r -= 1
+
+        if perm1[0] != 0 and perm1[m] != m:
+            perm = perm1[:]
+            flips_count = 0
+            k = perm[0]
+            while k:
+                perm[:k + 1] = perm[k::-1]
+                flips_count += 1
+                k = perm[0]
+
+            if flips_count > max_flips:
+                max_flips = flips_count
+
+        while r != n:
+            perm1_ins(r, perm1_pop(0))
+            count[r] -= 1
+            if count[r] > 0:
+                break
+            r += 1
+        else:
+            return max_flips
+
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (50, 10): (5,),
+    (100, 10): (6,),
+    (500, 10): (7,),
+    (1000, 10): (8,),
+    (5000, 10): (9,),
+}
+
+def bm_setup(params):
+    state = None
+    def run():
+        nonlocal state
+        state = fannkuch(params[0])
+    def result():
+        return params[0], state
+    return run, result
diff --git a/tests/perf_bench/bm_float.py b/tests/perf_bench/bm_float.py
new file mode 100644
index 000000000..5a66b9bb3
--- /dev/null
+++ b/tests/perf_bench/bm_float.py
@@ -0,0 +1,70 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# Artificial, floating point-heavy benchmark originally used by Factor.
+
+from math import sin, cos, sqrt
+
+
+class Point(object):
+    __slots__ = ('x', 'y', 'z')
+
+    def __init__(self, i):
+        self.x = x = sin(i)
+        self.y = cos(i) * 3
+        self.z = (x * x) / 2
+
+    def __repr__(self):
+        return "<Point: x=%s, y=%s, z=%s>" % (self.x, self.y, self.z)
+
+    def normalize(self):
+        x = self.x
+        y = self.y
+        z = self.z
+        norm = sqrt(x * x + y * y + z * z)
+        self.x /= norm
+        self.y /= norm
+        self.z /= norm
+
+    def maximize(self, other):
+        self.x = self.x if self.x > other.x else other.x
+        self.y = self.y if self.y > other.y else other.y
+        self.z = self.z if self.z > other.z else other.z
+        return self
+
+
+def maximize(points):
+    next = points[0]
+    for p in points[1:]:
+        next = next.maximize(p)
+    return next
+
+
+def benchmark(n):
+    points = [None] * n
+    for i in range(n):
+        points[i] = Point(i)
+    for p in points:
+        p.normalize()
+    return maximize(points)
+
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (50, 25): (1, 150),
+    (100, 100): (1, 250),
+    (1000, 1000): (10, 1500),
+    (5000, 1000): (20, 3000),
+}
+
+def bm_setup(params):
+    state = None
+    def run():
+        nonlocal state
+        for _ in range(params[0]):
+            state = benchmark(params[1])
+    def result():
+        return params[0] * params[1], 'Point(%.4f, %.4f, %.4f)' % (state.x, state.y, state.z)
+    return run, result
diff --git a/tests/perf_bench/bm_hexiom.py b/tests/perf_bench/bm_hexiom.py
new file mode 100644
index 000000000..3a6f1f6c4
--- /dev/null
+++ b/tests/perf_bench/bm_hexiom.py
@@ -0,0 +1,647 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# Solver of Hexiom board game.
+# Benchmark from Laurent Vaucher.
+# Source: https://github.com/slowfrog/hexiom : hexiom2.py, level36.txt
+# (Main function tweaked by Armin Rigo.)
+
+
+##################################
+class Dir(object):
+
+    def __init__(self, x, y):
+        self.x = x
+        self.y = y
+
+
+DIRS = [Dir(1, 0),
+        Dir(-1, 0),
+        Dir(0, 1),
+        Dir(0, -1),
+        Dir(1, 1),
+        Dir(-1, -1)]
+
+EMPTY = 7
+
+##################################
+
+
+class Done(object):
+    MIN_CHOICE_STRATEGY = 0
+    MAX_CHOICE_STRATEGY = 1
+    HIGHEST_VALUE_STRATEGY = 2
+    FIRST_STRATEGY = 3
+    MAX_NEIGHBORS_STRATEGY = 4
+    MIN_NEIGHBORS_STRATEGY = 5
+
+    def __init__(self, count, empty=False):
+        self.count = count
+        self.cells = None if empty else [
+            [0, 1, 2, 3, 4, 5, 6, EMPTY] for i in range(count)]
+
+    def clone(self):
+        ret = Done(self.count, True)
+        ret.cells = [self.cells[i][:] for i in range(self.count)]
+        return ret
+
+    def __getitem__(self, i):
+        return self.cells[i]
+
+    def set_done(self, i, v):
+        self.cells[i] = [v]
+
+    def already_done(self, i):
+        return len(self.cells[i]) == 1
+
+    def remove(self, i, v):
+        if v in self.cells[i]:
+            self.cells[i].remove(v)
+            return True
+        else:
+            return False
+
+    def remove_all(self, v):
+        for i in range(self.count):
+            self.remove(i, v)
+
+    def remove_unfixed(self, v):
+        changed = False
+        for i in range(self.count):
+            if not self.already_done(i):
+                if self.remove(i, v):
+                    changed = True
+        return changed
+
+    def filter_tiles(self, tiles):
+        for v in range(8):
+            if tiles[v] == 0:
+                self.remove_all(v)
+
+    def next_cell_min_choice(self):
+        minlen = 10
+        mini = -1
+        for i in range(self.count):
+            if 1 < len(self.cells[i]) < minlen:
+                minlen = len(self.cells[i])
+                mini = i
+        return mini
+
+    def next_cell_max_choice(self):
+        maxlen = 1
+        maxi = -1
+        for i in range(self.count):
+            if maxlen < len(self.cells[i]):
+                maxlen = len(self.cells[i])
+                maxi = i
+        return maxi
+
+    def next_cell_highest_value(self):
+        maxval = -1
+        maxi = -1
+        for i in range(self.count):
+            if (not self.already_done(i)):
+                maxvali = max(k for k in self.cells[i] if k != EMPTY)
+                if maxval < maxvali:
+                    maxval = maxvali
+                    maxi = i
+        return maxi
+
+    def next_cell_first(self):
+        for i in range(self.count):
+            if (not self.already_done(i)):
+                return i
+        return -1
+
+    def next_cell_max_neighbors(self, pos):
+        maxn = -1
+        maxi = -1
+        for i in range(self.count):
+            if not self.already_done(i):
+                cells_around = pos.hex.get_by_id(i).links
+                n = sum(1 if (self.already_done(nid) and (self[nid][0] != EMPTY)) else 0
+                        for nid in cells_around)
+                if n > maxn:
+                    maxn = n
+                    maxi = i
+        return maxi
+
+    def next_cell_min_neighbors(self, pos):
+        minn = 7
+        mini = -1
+        for i in range(self.count):
+            if not self.already_done(i):
+                cells_around = pos.hex.get_by_id(i).links
+                n = sum(1 if (self.already_done(nid) and (self[nid][0] != EMPTY)) else 0
+                        for nid in cells_around)
+                if n < minn:
+                    minn = n
+                    mini = i
+        return mini
+
+    def next_cell(self, pos, strategy=HIGHEST_VALUE_STRATEGY):
+        if strategy == Done.HIGHEST_VALUE_STRATEGY:
+            return self.next_cell_highest_value()
+        elif strategy == Done.MIN_CHOICE_STRATEGY:
+            return self.next_cell_min_choice()
+        elif strategy == Done.MAX_CHOICE_STRATEGY:
+            return self.next_cell_max_choice()
+        elif strategy == Done.FIRST_STRATEGY:
+            return self.next_cell_first()
+        elif strategy == Done.MAX_NEIGHBORS_STRATEGY:
+            return self.next_cell_max_neighbors(pos)
+        elif strategy == Done.MIN_NEIGHBORS_STRATEGY:
+            return self.next_cell_min_neighbors(pos)
+        else:
+            raise Exception("Wrong strategy: %d" % strategy)
+
+##################################
+
+
+class Node(object):
+
+    def __init__(self, pos, id, links):
+        self.pos = pos
+        self.id = id
+        self.links = links
+
+##################################
+
+
+class Hex(object):
+
+    def __init__(self, size):
+        self.size = size
+        self.count = 3 * size * (size - 1) + 1
+        self.nodes_by_id = self.count * [None]
+        self.nodes_by_pos = {}
+        id = 0
+        for y in range(size):
+            for x in range(size + y):
+                pos = (x, y)
+                node = Node(pos, id, [])
+                self.nodes_by_pos[pos] = node
+                self.nodes_by_id[node.id] = node
+                id += 1
+        for y in range(1, size):
+            for x in range(y, size * 2 - 1):
+                ry = size + y - 1
+                pos = (x, ry)
+                node = Node(pos, id, [])
+                self.nodes_by_pos[pos] = node
+                self.nodes_by_id[node.id] = node
+                id += 1
+
+    def link_nodes(self):
+        for node in self.nodes_by_id:
+            (x, y) = node.pos
+            for dir in DIRS:
+                nx = x + dir.x
+                ny = y + dir.y
+                if self.contains_pos((nx, ny)):
+                    node.links.append(self.nodes_by_pos[(nx, ny)].id)
+
+    def contains_pos(self, pos):
+        return pos in self.nodes_by_pos
+
+    def get_by_pos(self, pos):
+        return self.nodes_by_pos[pos]
+
+    def get_by_id(self, id):
+        return self.nodes_by_id[id]
+
+
+##################################
+class Pos(object):
+
+    def __init__(self, hex, tiles, done=None):
+        self.hex = hex
+        self.tiles = tiles
+        self.done = Done(hex.count) if done is None else done
+
+    def clone(self):
+        return Pos(self.hex, self.tiles, self.done.clone())
+
+##################################
+
+
+def constraint_pass(pos, last_move=None):
+    changed = False
+    left = pos.tiles[:]
+    done = pos.done
+
+    # Remove impossible values from free cells
+    free_cells = (range(done.count) if last_move is None
+                  else pos.hex.get_by_id(last_move).links)
+    for i in free_cells:
+        if not done.already_done(i):
+            vmax = 0
+            vmin = 0
+            cells_around = pos.hex.get_by_id(i).links
+            for nid in cells_around:
+                if done.already_done(nid):
+                    if done[nid][0] != EMPTY:
+                        vmin += 1
+                        vmax += 1
+                else:
+                    vmax += 1
+
+            for num in range(7):
+                if (num < vmin) or (num > vmax):
+                    if done.remove(i, num):
+                        changed = True
+
+    # Computes how many of each value is still free
+    for cell in done.cells:
+        if len(cell) == 1:
+            left[cell[0]] -= 1
+
+    for v in range(8):
+        # If there is none, remove the possibility from all tiles
+        if (pos.tiles[v] > 0) and (left[v] == 0):
+            if done.remove_unfixed(v):
+                changed = True
+        else:
+            possible = sum((1 if v in cell else 0) for cell in done.cells)
+            # If the number of possible cells for a value is exactly the number of available tiles
+            # put a tile in each cell
+            if pos.tiles[v] == possible:
+                for i in range(done.count):
+                    cell = done.cells[i]
+                    if (not done.already_done(i)) and (v in cell):
+                        done.set_done(i, v)
+                        changed = True
+
+    # Force empty or non-empty around filled cells
+    filled_cells = (range(done.count) if last_move is None
+                    else [last_move])
+    for i in filled_cells:
+        if done.already_done(i):
+            num = done[i][0]
+            empties = 0
+            filled = 0
+            unknown = []
+            cells_around = pos.hex.get_by_id(i).links
+            for nid in cells_around:
+                if done.already_done(nid):
+                    if done[nid][0] == EMPTY:
+                        empties += 1
+                    else:
+                        filled += 1
+                else:
+                    unknown.append(nid)
+            if len(unknown) > 0:
+                if num == filled:
+                    for u in unknown:
+                        if EMPTY in done[u]:
+                            done.set_done(u, EMPTY)
+                            changed = True
+                        # else:
+                        #    raise Exception("Houston, we've got a problem")
+                elif num == filled + len(unknown):
+                    for u in unknown:
+                        if done.remove(u, EMPTY):
+                            changed = True
+
+    return changed
+
+
+ASCENDING = 1
+DESCENDING = -1
+
+
+def find_moves(pos, strategy, order):
+    done = pos.done
+    cell_id = done.next_cell(pos, strategy)
+    if cell_id < 0:
+        return []
+
+    if order == ASCENDING:
+        return [(cell_id, v) for v in done[cell_id]]
+    else:
+        # Try higher values first and EMPTY last
+        moves = list(reversed([(cell_id, v)
+                               for v in done[cell_id] if v != EMPTY]))
+        if EMPTY in done[cell_id]:
+            moves.append((cell_id, EMPTY))
+        return moves
+
+
+def play_move(pos, move):
+    (cell_id, i) = move
+    pos.done.set_done(cell_id, i)
+
+
+def print_pos(pos, output):
+    hex = pos.hex
+    done = pos.done
+    size = hex.size
+    for y in range(size):
+        print(" " * (size - y - 1), end="", file=output)
+        for x in range(size + y):
+            pos2 = (x, y)
+            id = hex.get_by_pos(pos2).id
+            if done.already_done(id):
+                c = done[id][0] if done[id][0] != EMPTY else "."
+            else:
+                c = "?"
+            print("%s " % c, end="", file=output)
+        print(end="\n", file=output)
+    for y in range(1, size):
+        print(" " * y, end="", file=output)
+        for x in range(y, size * 2 - 1):
+            ry = size + y - 1
+            pos2 = (x, ry)
+            id = hex.get_by_pos(pos2).id
+            if done.already_done(id):
+                c = done[id][0] if done[id][0] != EMPTY else "."
+            else:
+                c = "?"
+            print("%s " % c, end="", file=output)
+        print(end="\n", file=output)
+
+
+OPEN = 0
+SOLVED = 1
+IMPOSSIBLE = -1
+
+
+def solved(pos, output, verbose=False):
+    hex = pos.hex
+    tiles = pos.tiles[:]
+    done = pos.done
+    exact = True
+    all_done = True
+    for i in range(hex.count):
+        if len(done[i]) == 0:
+            return IMPOSSIBLE
+        elif done.already_done(i):
+            num = done[i][0]
+            tiles[num] -= 1
+            if (tiles[num] < 0):
+                return IMPOSSIBLE
+            vmax = 0
+            vmin = 0
+            if num != EMPTY:
+                cells_around = hex.get_by_id(i).links
+                for nid in cells_around:
+                    if done.already_done(nid):
+                        if done[nid][0] != EMPTY:
+                            vmin += 1
+                            vmax += 1
+                    else:
+                        vmax += 1
+
+                if (num < vmin) or (num > vmax):
+                    return IMPOSSIBLE
+                if num != vmin:
+                    exact = False
+        else:
+            all_done = False
+
+    if (not all_done) or (not exact):
+        return OPEN
+
+    print_pos(pos, output)
+    return SOLVED
+
+
+def solve_step(prev, strategy, order, output, first=False):
+    if first:
+        pos = prev.clone()
+        while constraint_pass(pos):
+            pass
+    else:
+        pos = prev
+
+    moves = find_moves(pos, strategy, order)
+    if len(moves) == 0:
+        return solved(pos, output)
+    else:
+        for move in moves:
+            # print("Trying (%d, %d)" % (move[0], move[1]))
+            ret = OPEN
+            new_pos = pos.clone()
+            play_move(new_pos, move)
+            # print_pos(new_pos)
+            while constraint_pass(new_pos, move[0]):
+                pass
+            cur_status = solved(new_pos, output)
+            if cur_status != OPEN:
+                ret = cur_status
+            else:
+                ret = solve_step(new_pos, strategy, order, output)
+            if ret == SOLVED:
+                return SOLVED
+    return IMPOSSIBLE
+
+
+def check_valid(pos):
+    hex = pos.hex
+    tiles = pos.tiles
+    # fill missing entries in tiles
+    tot = 0
+    for i in range(8):
+        if tiles[i] > 0:
+            tot += tiles[i]
+        else:
+            tiles[i] = 0
+    # check total
+    if tot != hex.count:
+        raise Exception(
+            "Invalid input. Expected %d tiles, got %d." % (hex.count, tot))
+
+
+def solve(pos, strategy, order, output):
+    check_valid(pos)
+    return solve_step(pos, strategy, order, output, first=True)
+
+
+# TODO Write an 'iterator' to go over all x,y positions
+
+def read_file(file):
+    lines = [line.strip("\r\n") for line in file.splitlines()]
+    size = int(lines[0])
+    hex = Hex(size)
+    linei = 1
+    tiles = 8 * [0]
+    done = Done(hex.count)
+    for y in range(size):
+        line = lines[linei][size - y - 1:]
+        p = 0
+        for x in range(size + y):
+            tile = line[p:p + 2]
+            p += 2
+            if tile[1] == ".":
+                inctile = EMPTY
+            else:
+                inctile = int(tile)
+            tiles[inctile] += 1
+            # Look for locked tiles
+            if tile[0] == "+":
+                # print("Adding locked tile: %d at pos %d, %d, id=%d" %
+                #      (inctile, x, y, hex.get_by_pos((x, y)).id))
+                done.set_done(hex.get_by_pos((x, y)).id, inctile)
+
+        linei += 1
+    for y in range(1, size):
+        ry = size - 1 + y
+        line = lines[linei][y:]
+        p = 0
+        for x in range(y, size * 2 - 1):
+            tile = line[p:p + 2]
+            p += 2
+            if tile[1] == ".":
+                inctile = EMPTY
+            else:
+                inctile = int(tile)
+            tiles[inctile] += 1
+            # Look for locked tiles
+            if tile[0] == "+":
+                # print("Adding locked tile: %d at pos %d, %d, id=%d" %
+                #      (inctile, x, ry, hex.get_by_pos((x, ry)).id))
+                done.set_done(hex.get_by_pos((x, ry)).id, inctile)
+        linei += 1
+    hex.link_nodes()
+    done.filter_tiles(tiles)
+    return Pos(hex, tiles, done)
+
+
+def solve_file(file, strategy, order, output):
+    pos = read_file(file)
+    solve(pos, strategy, order, output)
+
+
+LEVELS = {}
+
+LEVELS[2] = ("""
+2
+  . 1
+ . 1 1
+  1 .
+""", """\
+ 1 1
+. . .
+ 1 1
+""")
+
+LEVELS[10] = ("""
+3
+  +.+. .
+ +. 0 . 2
+ . 1+2 1 .
+  2 . 0+.
+   .+.+.
+""", """\
+  . . 1
+ . 1 . 2
+0 . 2 2 .
+ . . . .
+  0 . .
+""")
+
+LEVELS[20] = ("""
+3
+   . 5 4
+  . 2+.+1
+ . 3+2 3 .
+ +2+. 5 .
+   . 3 .
+""", """\
+  3 3 2
+ 4 5 . 1
+3 5 2 . .
+ 2 . . .
+  . . .
+""")
+
+LEVELS[25] = ("""
+3
+   4 . .
+  . . 2 .
+ 4 3 2 . 4
+  2 2 3 .
+   4 2 4
+""", """\
+  3 4 2
+ 2 4 4 .
+. . . 4 2
+ . 2 4 3
+  . 2 .
+""")
+
+LEVELS[30] = ("""
+4
+    5 5 . .
+   3 . 2+2 6
+  3 . 2 . 5 .
+ . 3 3+4 4 . 3
+  4 5 4 . 5 4
+   5+2 . . 3
+    4 . . .
+""", """\
+   3 4 3 .
+  4 6 5 2 .
+ 2 5 5 . . 2
+. . 5 4 . 4 3
+ . 3 5 4 5 4
+  . 2 . 3 3
+   . . . .
+""")
+
+LEVELS[36] = ("""
+4
+    2 1 1 2
+   3 3 3 . .
+  2 3 3 . 4 .
+ . 2 . 2 4 3 2
+  2 2 . . . 2
+   4 3 4 . .
+    3 2 3 3
+""", """\
+   3 4 3 2
+  3 4 4 . 3
+ 2 . . 3 4 3
+2 . 1 . 3 . 2
+ 3 3 . 2 . 2
+  3 . 2 . 2
+   2 2 . 1
+""")
+
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (100, 100): (1, 10, DESCENDING, Done.FIRST_STRATEGY),
+    (1000, 1000): (1, 25, DESCENDING, Done.FIRST_STRATEGY),
+    (5000, 1000): (10, 25, DESCENDING, Done.FIRST_STRATEGY),
+}
+
+def bm_setup(params):
+    try:
+        import uio as io
+    except ImportError:
+        import io
+
+    loops, level, order, strategy = params
+
+    board, solution = LEVELS[level]
+    board = board.strip()
+    expected = solution.rstrip()
+    output = None
+
+    def run():
+        nonlocal output
+        for _ in range(loops):
+            stream = io.StringIO()
+            solve_file(board, strategy, order, stream)
+            output = stream.getvalue()
+            stream = None
+
+    def result():
+        norm = params[0] * params[1]
+        out = '\n'.join(line.rstrip() for line in output.splitlines())
+        return norm, ((out == expected), out)
+
+    return run, result
diff --git a/tests/perf_bench/bm_nqueens.py b/tests/perf_bench/bm_nqueens.py
new file mode 100644
index 000000000..87e7245c0
--- /dev/null
+++ b/tests/perf_bench/bm_nqueens.py
@@ -0,0 +1,62 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# Simple, brute-force N-Queens solver.
+# author: collinwinter@google.com (Collin Winter)
+# n_queens function: Copyright 2009 Raymond Hettinger
+
+# Pure-Python implementation of itertools.permutations().
+def permutations(iterable, r=None):
+    """permutations(range(3), 2) --> (0,1) (0,2) (1,0) (1,2) (2,0) (2,1)"""
+    pool = tuple(iterable)
+    n = len(pool)
+    if r is None:
+        r = n
+    indices = list(range(n))
+    cycles = list(range(n - r + 1, n + 1))[::-1]
+    yield tuple(pool[i] for i in indices[:r])
+    while n:
+        for i in reversed(range(r)):
+            cycles[i] -= 1
+            if cycles[i] == 0:
+                indices[i:] = indices[i + 1:] + indices[i:i + 1]
+                cycles[i] = n - i
+            else:
+                j = cycles[i]
+                indices[i], indices[-j] = indices[-j], indices[i]
+                yield tuple(pool[i] for i in indices[:r])
+                break
+        else:
+            return
+
+# From http://code.activestate.com/recipes/576647/
+def n_queens(queen_count):
+    """N-Queens solver.
+    Args: queen_count: the number of queens to solve for, same as board size.
+    Yields: Solutions to the problem, each yielded value is a N-tuple.
+    """
+    cols = range(queen_count)
+    for vec in permutations(cols):
+        if (queen_count == len(set(vec[i] + i for i in cols))
+                        == len(set(vec[i] - i for i in cols))):
+            yield vec
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (50, 25): (1, 5),
+    (100, 25): (1, 6),
+    (1000, 100): (1, 7),
+    (5000, 100): (1, 8),
+}
+
+def bm_setup(params):
+    res = None
+    def run():
+        nonlocal res
+        for _ in range(params[0]):
+            res = len(list(n_queens(params[1])))
+    def result():
+        return params[0] * 10 ** (params[1] - 3), res
+    return run, result
diff --git a/tests/perf_bench/bm_pidigits.py b/tests/perf_bench/bm_pidigits.py
new file mode 100644
index 000000000..ca2e5297d
--- /dev/null
+++ b/tests/perf_bench/bm_pidigits.py
@@ -0,0 +1,62 @@
+# Source: https://github.com/python/pyperformance
+# License: MIT
+
+# Calculating some of the digits of π.
+# This benchmark stresses big integer arithmetic.
+# Adapted from code on: http://benchmarksgame.alioth.debian.org/
+
+
+def compose(a, b):
+    aq, ar, as_, at = a
+    bq, br, bs, bt = b
+    return (aq * bq,
+            aq * br + ar * bt,
+            as_ * bq + at * bs,
+            as_ * br + at * bt)
+
+
+def extract(z, j):
+    q, r, s, t = z
+    return (q * j + r) // (s * j + t)
+
+
+def gen_pi_digits(n):
+    z = (1, 0, 0, 1)
+    k = 1
+    digs = []
+    for _ in range(n):
+        y = extract(z, 3)
+        while y != extract(z, 4):
+            z = compose(z, (k, 4 * k + 2, 0, 2 * k + 1))
+            k += 1
+            y = extract(z, 3)
+        z = compose((10, -10 * y, 0, 1), z)
+        digs.append(y)
+    return digs
+
+
+###########################################################################
+# Benchmark interface
+
+bm_params = {
+    (50, 25): (1, 35),
+    (100, 100): (1, 65),
+    (1000, 1000): (2, 250),
+    (5000, 1000): (3, 350),
+}
+
+def bm_setup(params):
+    state = None
+
+    def run():
+        nonlocal state
+        nloop, ndig = params
+        ndig = params[1]
+        for _ in range(nloop):
+            state = None # free previous result
+            state = gen_pi_digits(ndig)
+
+    def result():
+        return params[0] * params[1], ''.join(str(d) for d in state)
+
+    return run, result