Overflow benchmarks.

2 files changed, 240 insertions, 0 deletions

diff --git a/Demos/overflow_perf.pyx b/Demos/overflow_perf.pyx
new file mode 100644
index 000000000..d8dc28990
--- /dev/null
+++ b/Demos/overflow_perf.pyx
@@ -0,0 +1,187 @@
+# distutils: extra_compile_args = -O3
+
+import cython
+
+ctypedef fused INT:
+    int
+    long long
+    unsigned int
+    unsigned long long
+    object
+
+ctypedef fused C_INT:
+    int
+    long long
+    unsigned int
+    unsigned long long
+
+@cython.overflowcheck(False)
+def fib(INT n):
+    """
+    >>> [fib(k) for k in range(10)]
+    [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
+    """
+    cdef INT a, b, k
+    a, b = 0, 1
+    for k in range(n):
+        a, b = b, a + b
+    return int(b)
+
+
+@cython.overflowcheck(True)
+def fib_overflow(INT n):
+    """
+    >>> [fib_overflow(k) for k in range(10)]
+    [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
+    """
+    cdef INT a, b, k
+    a, b = 0, 1
+    for k in range(n):
+        a, b = b, a + b
+    return int(b)
+
+
+@cython.overflowcheck(False)
+def collatz(INT n):
+    """
+    >>> collatz(1)
+    0
+    >>> collatz(5)
+    5
+    >>> collatz(10)
+    6
+    """
+    cdef INT k = 0
+    while n != 1:
+        if n % 2 == 0:
+            n /= 2
+        else:
+            n = 3*n + 1
+        k += 1
+    return int(k)
+
+@cython.overflowcheck(True)
+@cython.overflowcheck.fold(False)
+def collatz_overflow(INT n):
+    """
+    >>> collatz_overflow(1)
+    0
+    >>> collatz_overflow(5)
+    5
+    >>> collatz_overflow(10)
+    6
+    """
+    cdef INT k = 0
+    while n != 1:
+        if n % 2 == 0:
+            n /= 2
+        else:
+            n = 3*n + 1
+        k += 1
+    return int(k)
+
+@cython.overflowcheck(True)
+@cython.overflowcheck.fold(True)
+def collatz_overflow_fold(INT n):
+    """
+    >>> collatz_overflow_fold(1)
+    0
+    >>> collatz_overflow_fold(5)
+    5
+    >>> collatz_overflow_fold(10)
+    6
+    """
+    cdef INT k = 0
+    while n != 1:
+        if n % 2 == 0:
+            n /= 2
+        else:
+            n = 3*n + 1
+        k += 1
+    return int(k)
+
+
+
+@cython.overflowcheck(False)
+def factorial(INT n):
+    """
+    >>> factorial(2)
+    2
+    >>> factorial(5)
+    120
+    """
+    cdef INT k, res = 1
+    for k in range(2, n+1):
+        res = res * k
+    return int(res)
+
+@cython.overflowcheck(True)
+def factorial_overflow(INT n):
+    """
+    >>> factorial_overflow(2)
+    2
+    >>> factorial_overflow(5)
+    120
+    """
+    cdef INT k, res = 1
+    for k in range(2, n+1):
+        res = res * k
+    return int(res)
+
+
+
+@cython.overflowcheck(False)
+def most_orthogonal(C_INT[:,::1] vectors):
+    cdef C_INT n = vectors.shape[0]
+    cdef C_INT* a
+    cdef C_INT* b
+    cdef double min_dot = 2 # actual max is 1
+    for i in range(n):
+        for j in range(i):
+            a = &vectors[i, 0]
+            b = &vectors[j, 0]
+            # A highly nested arithmetic expression...
+            normalized_dot = (1.0 * (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) /
+                ((a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) * (b[0]*b[0] + b[1]*b[1]+b[2]*b[2])))
+            if normalized_dot < min_dot:
+                min_dot = normalized_dot
+                min_pair = i, j
+    return vectors[i], vectors[j]
+
+@cython.overflowcheck(True)
+@cython.overflowcheck.fold(False)
+def most_orthogonal_overflow(C_INT[:,::1] vectors):
+    cdef C_INT n = vectors.shape[0]
+    cdef C_INT* a
+    cdef C_INT* b
+    cdef double min_dot = 2 # actual max is 1
+    for i in range(n):
+        for j in range(i):
+            a = &vectors[i, 0]
+            b = &vectors[j, 0]
+            # A highly nested arithmetic expression...
+            normalized_dot = ((a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) /
+                (1.0 * (a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) * (b[0]*b[0] + b[1]*b[1]+b[2]*b[2])))
+            if normalized_dot < min_dot:
+                min_dot = normalized_dot
+                min_pair = i, j
+    return vectors[i], vectors[j]
+
+@cython.overflowcheck(True)
+@cython.overflowcheck.fold(True)
+def most_orthogonal_overflow_fold(C_INT[:,::1] vectors):
+    cdef C_INT n = vectors.shape[0]
+    cdef C_INT* a
+    cdef C_INT* b
+    cdef double min_dot = 2 # actual max is 1
+    for i in range(n):
+        for j in range(i):
+            a = &vectors[i, 0]
+            b = &vectors[j, 0]
+            # A highly nested arithmetic expression...
+            normalized_dot = ((a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) /
+                (1.0 * (a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) * (b[0]*b[0] + b[1]*b[1]+b[2]*b[2])))
+            if normalized_dot < min_dot:
+                min_dot = normalized_dot
+                min_pair = i, j
+    return vectors[i], vectors[j]
diff --git a/Demos/overflow_perf_run.py b/Demos/overflow_perf_run.py
new file mode 100644
index 000000000..ec84947c2
--- /dev/null
+++ b/Demos/overflow_perf_run.py
@@ -0,0 +1,53 @@
+from overflow_perf import *
+
+import sys
+import timeit
+try:
+    import numpy as np
+except ImportError:
+    np = None
+
+
+def run_tests(N):
+    global f
+    for func in most_orthogonal, fib, collatz, factorial:
+        print func.__name__
+        for type in ['int', 'unsigned int', 'long long', 'unsigned long long', 'object']:
+            if func == most_orthogonal:
+                if type == 'object' or np == None:
+                    continue
+                type_map = {'int': 'int32', 'unsigned int': 'uint32', 'long long': 'int64', 'unsigned long long': 'uint64'}
+                shape = N, 3
+                arg = np.ndarray(shape, dtype=type_map[type])
+                arg[:] = 1000 * np.random.random(shape)
+            else:
+                arg = N
+            try:
+                print "%s[%s](%s)" % (func.__name__, type, N)
+                with_overflow = my_timeit(globals()[func.__name__ + "_overflow"][type], arg)
+                no_overflow = my_timeit(func[type], arg)
+                print "\t%0.04e\t%0.04e\t%0.04f" % (no_overflow, with_overflow, with_overflow / no_overflow)
+                if func.__name__ + "_overflow_fold" in globals():
+                    with_overflow = my_timeit(globals()[func.__name__ + "_overflow_fold"][type], arg)
+                    print  "\t%0.04e\t%0.04e\t%0.04f" % (no_overflow, with_overflow, with_overflow / no_overflow), "(folded)"
+            except OverflowError:
+                print "    ", "Overflow"
+
+def my_timeit(func, N):
+    global f, arg
+    f = func
+    arg = N
+    for exponent in range(10, 30):
+        times = 2 ** exponent
+        res = min(timeit.repeat("f(arg)", setup="from __main__ import f, arg", repeat=5, number=times))
+        if res > .25:
+            break
+    return res / times
+
+params = sys.argv[1:]
+if not params:
+    params = [129, 9, 97]
+for arg in params:
+    print
+    print "N", arg
+    run_tests(int(arg))