Issue #7652: Integrate the decimal floating point libmpdec library to speed

up the decimal module. Performance gains of the new C implementation are
between 12x and 80x, depending on the application.

8 files changed, 2654 insertions, 0 deletions

diff --git a/Modules/_decimal/tests/README.txt b/Modules/_decimal/tests/README.txt
new file mode 100644
index 0000000..97b6ff6
--- /dev/null
+++ b/Modules/_decimal/tests/README.txt
@@ -0,0 +1,15 @@
+
+
+This directory contains extended tests and a benchmark against decimal.py:
+
+  bench.py  ->  Benchmark for small and large precisions.
+  Usage: ../../../python bench.py
+
+  formathelper.py   ->
+  randdec.py        ->  Generate test cases for deccheck.py.
+  randfloat.py      ->
+
+  deccheck.py  ->  Run extended tests.
+  Usage: ../../../python deccheck.py [--short|--medium|--long|--all]
+
+
diff --git a/Modules/_decimal/tests/bench.py b/Modules/_decimal/tests/bench.py
new file mode 100644
index 0000000..013f0f6
--- /dev/null
+++ b/Modules/_decimal/tests/bench.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+
+#
+# Copyright (C) 2001-2012 Python Software Foundation. All Rights Reserved.
+# Modified and extended by Stefan Krah.
+#
+
+# Usage: ../../../python bench.py
+
+
+import time
+from math import log, ceil
+from test.support import import_fresh_module
+
+C = import_fresh_module('decimal', fresh=['_decimal'])
+P = import_fresh_module('decimal', blocked=['_decimal'])
+
+
+# Pi function from the decimal.py documentation
+def pi_float(prec):
+    """native float"""
+    lasts, t, s, n, na, d, da = 0, 3.0, 3, 1, 0, 0, 24
+    while s != lasts:
+        lasts = s
+        n, na = n+na, na+8
+        d, da = d+da, da+32
+        t = (t * n) / d
+        s += t
+    return s
+
+def pi_cdecimal(prec):
+    """cdecimal"""
+    C.getcontext().prec = prec
+    D = C.Decimal
+    lasts, t, s, n, na, d, da = D(0), D(3), D(3), D(1), D(0), D(0), D(24)
+    while s != lasts:
+        lasts = s
+        n, na = n+na, na+8
+        d, da = d+da, da+32
+        t = (t * n) / d
+        s += t
+    return s
+
+def pi_decimal(prec):
+    """decimal"""
+    P.getcontext().prec = prec
+    D = P.Decimal
+    lasts, t, s, n, na, d, da = D(0), D(3), D(3), D(1), D(0), D(0), D(24)
+    while s != lasts:
+        lasts = s
+        n, na = n+na, na+8
+        d, da = d+da, da+32
+        t = (t * n) / d
+        s += t
+    return s
+
+def factorial(n, m):
+    if (n > m):
+        return factorial(m, n)
+    elif m == 0:
+        return 1
+    elif n == m:
+        return n
+    else:
+        return factorial(n, (n+m)//2) * factorial((n+m)//2 + 1, m)
+
+
+print("\n# ======================================================================")
+print("#                   Calculating pi, 10000 iterations")
+print("# ======================================================================\n")
+
+for prec in [9, 19]:
+    print("\nPrecision: %d decimal digits\n" % prec)
+    for func in [pi_float, pi_cdecimal, pi_decimal]:
+        start = time.time()
+        for i in range(10000):
+            x = func(prec)
+        print("%s:" % func.__name__.replace("pi_", ""))
+        print("result: %s" % str(x))
+        print("time: %fs\n" % (time.time()-start))
+
+
+print("\n# ======================================================================")
+print("#                               Factorial")
+print("# ======================================================================\n")
+
+C.getcontext().prec = C.MAX_PREC
+
+for n in [100000, 1000000]:
+
+    print("n = %d\n" % n)
+
+    # C version of decimal
+    start_calc = time.time()
+    x = factorial(C.Decimal(n), 0)
+    end_calc = time.time()
+    start_conv = time.time()
+    sx = str(x)
+    end_conv = time.time()
+    print("cdecimal:")
+    print("calculation time: %fs" % (end_calc-start_calc))
+    print("conversion time: %fs\n" % (end_conv-start_conv))
+
+    # Python integers
+    start_calc = time.time()
+    y = factorial(n, 0)
+    end_calc = time.time()
+    start_conv = time.time()
+    sy = str(y)
+    end_conv =  time.time()
+
+    print("int:")
+    print("calculation time: %fs" % (end_calc-start_calc))
+    print("conversion time: %fs\n\n" % (end_conv-start_conv))
+
+    assert(sx == sy)
diff --git a/Modules/_decimal/tests/deccheck.py b/Modules/_decimal/tests/deccheck.py
new file mode 100644
index 0000000..aefbed6
--- /dev/null
+++ b/Modules/_decimal/tests/deccheck.py
@@ -0,0 +1,1074 @@
+#!/usr/bin/env python
+
+#
+# Copyright (c) 2008-2012 Stefan Krah. All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+
+#
+# Usage: python deccheck.py [--short|--medium|--long|--all]
+#
+
+import sys, random
+from copy import copy
+from collections import defaultdict
+from test.support import import_fresh_module
+from randdec import randfloat, all_unary, all_binary, all_ternary
+from formathelper import rand_format, rand_locale
+
+C = import_fresh_module('decimal', fresh=['_decimal'])
+P = import_fresh_module('decimal', blocked=['_decimal'])
+EXIT_STATUS = 0
+
+
+# Contains all categories of Decimal methods.
+Functions = {
+    # Plain unary:
+    'unary': (
+        '__abs__', '__bool__', '__ceil__', '__complex__', '__copy__',
+        '__floor__', '__float__', '__hash__', '__int__', '__neg__',
+        '__pos__', '__reduce__', '__repr__', '__str__', '__trunc__',
+        'adjusted', 'as_tuple', 'canonical', 'conjugate', 'copy_abs',
+        'copy_negate', 'is_canonical', 'is_finite', 'is_infinite',
+        'is_nan', 'is_qnan', 'is_signed', 'is_snan', 'is_zero', 'radix'
+    ),
+    # Unary with optional context:
+    'unary_ctx': (
+        'exp', 'is_normal', 'is_subnormal', 'ln', 'log10', 'logb',
+        'logical_invert', 'next_minus', 'next_plus', 'normalize',
+        'number_class', 'sqrt', 'to_eng_string'
+    ),
+    # Unary with optional rounding mode and context:
+    'unary_rnd_ctx': ('to_integral', 'to_integral_exact', 'to_integral_value'),
+    # Plain binary:
+    'binary': (
+        '__add__', '__divmod__', '__eq__', '__floordiv__', '__ge__', '__gt__',
+        '__le__', '__lt__', '__mod__', '__mul__', '__ne__', '__pow__',
+        '__radd__', '__rdivmod__', '__rfloordiv__', '__rmod__', '__rmul__',
+        '__rpow__', '__rsub__', '__rtruediv__', '__sub__', '__truediv__',
+        'compare_total', 'compare_total_mag', 'copy_sign', 'quantize',
+        'same_quantum'
+    ),
+    # Binary with optional context:
+    'binary_ctx': (
+        'compare', 'compare_signal', 'logical_and', 'logical_or', 'logical_xor',
+        'max', 'max_mag', 'min', 'min_mag', 'next_toward', 'remainder_near',
+        'rotate', 'scaleb', 'shift'
+    ),
+    # Plain ternary:
+    'ternary': ('__pow__',),
+    # Ternary with optional context:
+    'ternary_ctx': ('fma',),
+    # Special:
+    'special': ('__format__', '__reduce_ex__', '__round__', 'from_float',
+                'quantize'),
+    # Properties:
+    'property': ('real', 'imag')
+}
+
+# Contains all categories of Context methods. The n-ary classification
+# applies to the number of Decimal arguments.
+ContextFunctions = {
+    # Plain nullary:
+    'nullary': ('context.__hash__', 'context.__reduce__', 'context.radix'),
+    # Plain unary:
+    'unary': ('context.abs', 'context.canonical', 'context.copy_abs',
+              'context.copy_decimal', 'context.copy_negate',
+              'context.create_decimal', 'context.exp', 'context.is_canonical',
+              'context.is_finite', 'context.is_infinite', 'context.is_nan',
+              'context.is_normal', 'context.is_qnan', 'context.is_signed',
+              'context.is_snan', 'context.is_subnormal', 'context.is_zero',
+              'context.ln', 'context.log10', 'context.logb',
+              'context.logical_invert', 'context.minus', 'context.next_minus',
+              'context.next_plus', 'context.normalize', 'context.number_class',
+              'context.plus', 'context.sqrt', 'context.to_eng_string',
+              'context.to_integral', 'context.to_integral_exact',
+              'context.to_integral_value', 'context.to_sci_string'
+    ),
+    # Plain binary:
+    'binary': ('context.add', 'context.compare', 'context.compare_signal',
+               'context.compare_total', 'context.compare_total_mag',
+               'context.copy_sign', 'context.divide', 'context.divide_int',
+               'context.divmod', 'context.logical_and', 'context.logical_or',
+               'context.logical_xor', 'context.max', 'context.max_mag',
+               'context.min', 'context.min_mag', 'context.multiply',
+               'context.next_toward', 'context.power', 'context.quantize',
+               'context.remainder', 'context.remainder_near', 'context.rotate',
+               'context.same_quantum', 'context.scaleb', 'context.shift',
+               'context.subtract'
+    ),
+    # Plain ternary:
+    'ternary': ('context.fma', 'context.power'),
+    # Special:
+    'special': ('context.__reduce_ex__', 'context.create_decimal_from_float')
+}
+
+# Functions that require a restricted exponent range for reasonable runtimes.
+UnaryRestricted = [
+  '__ceil__', '__floor__', '__int__', '__long__', '__trunc__',
+  'to_integral', 'to_integral_value'
+]
+
+BinaryRestricted = ['__round__']
+
+TernaryRestricted = ['__pow__', 'context.power']
+
+
+# ======================================================================
+#                            Unified Context
+# ======================================================================
+
+# Translate symbols.
+CondMap = {
+        C.Clamped:             P.Clamped,
+        C.ConversionSyntax:    P.ConversionSyntax,
+        C.DivisionByZero:      P.DivisionByZero,
+        C.DivisionImpossible:  P.InvalidOperation,
+        C.DivisionUndefined:   P.DivisionUndefined,
+        C.Inexact:             P.Inexact,
+        C.InvalidContext:      P.InvalidContext,
+        C.InvalidOperation:    P.InvalidOperation,
+        C.Overflow:            P.Overflow,
+        C.Rounded:             P.Rounded,
+        C.Subnormal:           P.Subnormal,
+        C.Underflow:           P.Underflow,
+        C.FloatOperation:      P.FloatOperation,
+}
+
+RoundMap = {
+        C.ROUND_UP:           P.ROUND_UP,
+        C.ROUND_DOWN:         P.ROUND_DOWN,
+        C.ROUND_CEILING:      P.ROUND_CEILING,
+        C.ROUND_FLOOR:        P.ROUND_FLOOR,
+        C.ROUND_HALF_UP:      P.ROUND_HALF_UP,
+        C.ROUND_HALF_DOWN:    P.ROUND_HALF_DOWN,
+        C.ROUND_HALF_EVEN:    P.ROUND_HALF_EVEN,
+        C.ROUND_05UP:         P.ROUND_05UP
+}
+RoundModes = RoundMap.items()
+
+
+class Context(object):
+    """Provides a convenient way of syncing the C and P contexts"""
+
+    __slots__ = ['c', 'p']
+
+    def __init__(self, c_ctx=None, p_ctx=None):
+        """Initialization is from the C context"""
+        self.c = C.getcontext() if c_ctx is None else c_ctx
+        self.p = P.getcontext() if p_ctx is None else p_ctx
+        self.p.prec = self.c.prec
+        self.p.Emin = self.c.Emin
+        self.p.Emax = self.c.Emax
+        self.p.rounding = RoundMap[self.c.rounding]
+        self.p.capitals = self.c.capitals
+        self.settraps([sig for sig in self.c.traps if self.c.traps[sig]])
+        self.setstatus([sig for sig in self.c.flags if self.c.flags[sig]])
+        self.p.clamp = self.c.clamp
+
+    def __str__(self):
+        return str(self.c) + '\n' + str(self.p)
+
+    def getprec(self):
+        assert(self.c.prec == self.p.prec)
+        return self.c.prec
+
+    def setprec(self, val):
+        self.c.prec = val
+        self.p.prec = val
+
+    def getemin(self):
+        assert(self.c.Emin == self.p.Emin)
+        return self.c.Emin
+
+    def setemin(self, val):
+        self.c.Emin = val
+        self.p.Emin = val
+
+    def getemax(self):
+        assert(self.c.Emax == self.p.Emax)
+        return self.c.Emax
+
+    def setemax(self, val):
+        self.c.Emax = val
+        self.p.Emax = val
+
+    def getround(self):
+        assert(self.c.rounding == RoundMap[self.p.rounding])
+        return self.c.rounding
+
+    def setround(self, val):
+        self.c.rounding = val
+        self.p.rounding = RoundMap[val]
+
+    def getcapitals(self):
+        assert(self.c.capitals == self.p.capitals)
+        return self.c.capitals
+
+    def setcapitals(self, val):
+        self.c.capitals = val
+        self.p.capitals = val
+
+    def getclamp(self):
+        assert(self.c.clamp == self.p.clamp)
+        return self.c.clamp
+
+    def setclamp(self, val):
+        self.c.clamp = val
+        self.p.clamp = val
+
+    prec = property(getprec, setprec)
+    Emin = property(getemin, setemin)
+    Emax = property(getemax, setemax)
+    rounding = property(getround, setround)
+    clamp = property(getclamp, setclamp)
+    capitals = property(getcapitals, setcapitals)
+
+    def clear_traps(self):
+        self.c.clear_traps()
+        for trap in self.p.traps:
+            self.p.traps[trap] = False
+
+    def clear_status(self):
+        self.c.clear_flags()
+        self.p.clear_flags()
+
+    def settraps(self, lst):
+        """lst: C signal list"""
+        self.clear_traps()
+        for signal in lst:
+            self.c.traps[signal] = True
+            self.p.traps[CondMap[signal]] = True
+
+    def setstatus(self, lst):
+        """lst: C signal list"""
+        self.clear_status()
+        for signal in lst:
+            self.c.flags[signal] = True
+            self.p.flags[CondMap[signal]] = True
+
+    def assert_eq_status(self):
+        """assert equality of C and P status"""
+        for signal in self.c.flags:
+            if self.c.flags[signal] == (not self.p.flags[CondMap[signal]]):
+                return False
+        return True
+
+
+# We don't want exceptions so that we can compare the status flags.
+context = Context()
+context.Emin = C.MIN_EMIN
+context.Emax = C.MAX_EMAX
+context.clear_traps()
+
+# When creating decimals, _decimal is ultimately limited by the maximum
+# context values. We emulate this restriction for decimal.py.
+maxcontext = P.Context(
+    prec=C.MAX_PREC,
+    Emin=C.MIN_EMIN,
+    Emax=C.MAX_EMAX,
+    rounding=P.ROUND_HALF_UP,
+    capitals=1
+)
+maxcontext.clamp = 0
+
+def RestrictedDecimal(value):
+    maxcontext.traps = copy(context.p.traps)
+    maxcontext.clear_flags()
+    if isinstance(value, str):
+        value = value.strip()
+    dec = maxcontext.create_decimal(value)
+    if maxcontext.flags[P.Inexact] or \
+       maxcontext.flags[P.Rounded] or \
+       maxcontext.flags[P.InvalidOperation]:
+        return context.p._raise_error(P.InvalidOperation)
+    if maxcontext.flags[P.FloatOperation]:
+        context.p.flags[P.FloatOperation] = True
+    return dec
+
+
+# ======================================================================
+#      TestSet: Organize data and events during a single test case
+# ======================================================================
+
+class RestrictedList(list):
+    """List that can only be modified by appending items."""
+    def __getattribute__(self, name):
+        if name != 'append':
+            raise AttributeError("unsupported operation")
+        return list.__getattribute__(self, name)
+    def unsupported(self, *_):
+        raise AttributeError("unsupported operation")
+    __add__ = __delattr__ = __delitem__ = __iadd__ = __imul__ = unsupported
+    __mul__ = __reversed__ = __rmul__ = __setattr__ = __setitem__ = unsupported
+
+class TestSet(object):
+    """A TestSet contains the original input operands, converted operands,
+       Python exceptions that occurred either during conversion or during
+       execution of the actual function, and the final results.
+
+       For safety, most attributes are lists that only support the append
+       operation.
+
+       If a function name is prefixed with 'context.', the corresponding
+       context method is called.
+    """
+    def __init__(self, funcname, operands):
+        if funcname.startswith("context."):
+            self.funcname = funcname.replace("context.", "")
+            self.contextfunc = True
+        else:
+            self.funcname = funcname
+            self.contextfunc = False
+        self.op = operands               # raw operand tuple
+        self.context = context           # context used for the operation
+        self.cop = RestrictedList()      # converted C.Decimal operands
+        self.cex = RestrictedList()      # Python exceptions for C.Decimal
+        self.cresults = RestrictedList() # C.Decimal results
+        self.pop = RestrictedList()      # converted P.Decimal operands
+        self.pex = RestrictedList()      # Python exceptions for P.Decimal
+        self.presults = RestrictedList() # P.Decimal results
+
+
+# ======================================================================
+#                SkipHandler: skip known discrepancies
+# ======================================================================
+
+class SkipHandler:
+    """Handle known discrepancies between decimal.py and _decimal.so.
+       These are either ULP differences in the power function or
+       extremely minor issues."""
+
+    def __init__(self):
+        self.ulpdiff = 0
+        self.powmod_zeros = 0
+        self.maxctx = P.Context(Emax=10**18, Emin=-10**18)
+
+    def default(self, t):
+        return False
+    __ge__ =  __gt__ = __le__ = __lt__ = __ne__ = __eq__ = default
+    __reduce__ = __format__ = __repr__ = __str__ = default
+
+    def harrison_ulp(self, dec):
+        """ftp://ftp.inria.fr/INRIA/publication/publi-pdf/RR/RR-5504.pdf"""
+        a = dec.next_plus()
+        b = dec.next_minus()
+        return abs(a - b)
+
+    def standard_ulp(self, dec, prec):
+        return P._dec_from_triple(0, '1', dec._exp+len(dec._int)-prec)
+
+    def rounding_direction(self, x, mode):
+        """Determine the effective direction of the rounding when
+           the exact result x is rounded according to mode.
+           Return -1 for downwards, 0 for undirected, 1 for upwards,
+           2 for ROUND_05UP."""
+        cmp = 1 if x.compare_total(P.Decimal("+0")) >= 0 else -1
+
+        if mode in (P.ROUND_HALF_EVEN, P.ROUND_HALF_UP, P.ROUND_HALF_DOWN):
+            return 0
+        elif mode == P.ROUND_CEILING:
+            return 1
+        elif mode == P.ROUND_FLOOR:
+            return -1
+        elif mode == P.ROUND_UP:
+            return cmp
+        elif mode == P.ROUND_DOWN:
+            return -cmp
+        elif mode == P.ROUND_05UP:
+            return 2
+        else:
+            raise ValueError("Unexpected rounding mode: %s" % mode)
+
+    def check_ulpdiff(self, exact, rounded):
+        # current precision
+        p = context.p.prec
+
+        # Convert infinities to the largest representable number + 1.
+        x = exact
+        if exact.is_infinite():
+            x = P._dec_from_triple(exact._sign, '10', context.p.Emax)
+        y = rounded
+        if rounded.is_infinite():
+            y = P._dec_from_triple(rounded._sign, '10', context.p.Emax)
+
+        # err = (rounded - exact) / ulp(rounded)
+        self.maxctx.prec = p * 2
+        t = self.maxctx.subtract(y, x)
+        if context.c.flags[C.Clamped] or \
+           context.c.flags[C.Underflow]:
+            # The standard ulp does not work in Underflow territory.
+            ulp = self.harrison_ulp(y)
+        else:
+            ulp = self.standard_ulp(y, p)
+        # Error in ulps.
+        err = self.maxctx.divide(t, ulp)
+
+        dir = self.rounding_direction(x, context.p.rounding)
+        if dir == 0:
+            if P.Decimal("-0.6") < err < P.Decimal("0.6"):
+                return True
+        elif dir == 1: # directed, upwards
+            if P.Decimal("-0.1") < err < P.Decimal("1.1"):
+                return True
+        elif dir == -1: # directed, downwards
+            if P.Decimal("-1.1") < err < P.Decimal("0.1"):
+                return True
+        else: # ROUND_05UP
+            if P.Decimal("-1.1") < err < P.Decimal("1.1"):
+                return True
+
+        print("ulp: %s  error: %s  exact: %s  c_rounded: %s"
+              % (ulp, err, exact, rounded))
+        return False
+
+    def bin_resolve_ulp(self, t):
+        """Check if results of _decimal's power function are within the
+           allowed ulp ranges."""
+        # NaNs are beyond repair.
+        if t.rc.is_nan() or t.rp.is_nan():
+            return False
+
+        # "exact" result, double precision, half_even
+        self.maxctx.prec = context.p.prec * 2
+
+        op1, op2 = t.pop[0], t.pop[1]
+        if t.contextfunc:
+            exact = getattr(self.maxctx, t.funcname)(op1, op2)
+        else:
+            exact = getattr(op1, t.funcname)(op2, context=self.maxctx)
+
+        # _decimal's rounded result
+        rounded = P.Decimal(t.cresults[0])
+
+        self.ulpdiff += 1
+        return self.check_ulpdiff(exact, rounded)
+
+    ############################ Correct rounding #############################
+    def resolve_underflow(self, t):
+        """In extremely rare cases where the infinite precision result is just
+           below etiny, cdecimal does not set Subnormal/Underflow. Example:
+
+           setcontext(Context(prec=21, rounding=ROUND_UP, Emin=-55, Emax=85))
+           Decimal("1.00000000000000000000000000000000000000000000000"
+                   "0000000100000000000000000000000000000000000000000"
+                   "0000000000000025").ln()
+        """
+        if t.cresults != t.presults:
+            return False # Results must be identical.
+        if context.c.flags[C.Rounded] and \
+           context.c.flags[C.Inexact] and \
+           context.p.flags[P.Rounded] and \
+           context.p.flags[P.Inexact]:
+            return True # Subnormal/Underflow may be missing.
+        return False
+
+    def exp(self, t):
+        """Resolve Underflow or ULP difference."""
+        return self.resolve_underflow(t)
+
+    def log10(self, t):
+        """Resolve Underflow or ULP difference."""
+        return self.resolve_underflow(t)
+
+    def ln(self, t):
+        """Resolve Underflow or ULP difference."""
+        return self.resolve_underflow(t)
+
+    def __pow__(self, t):
+        """Always calls the resolve function. C.Decimal does not have correct
+           rounding for the power function."""
+        if context.c.flags[C.Rounded] and \
+           context.c.flags[C.Inexact] and \
+           context.p.flags[P.Rounded] and \
+           context.p.flags[P.Inexact]:
+            return self.bin_resolve_ulp(t)
+        else:
+            return False
+    power = __rpow__ = __pow__
+
+    ############################## Technicalities #############################
+    def __float__(self, t):
+        """NaN comparison in the verify() function obviously gives an
+           incorrect answer:  nan == nan -> False"""
+        if t.cop[0].is_nan() and t.pop[0].is_nan():
+            return True
+        return False
+    __complex__ = __float__
+
+    def __radd__(self, t):
+        """decimal.py gives precedence to the first NaN; this is
+           not important, as __radd__ will not be called for
+           two decimal arguments."""
+        if t.rc.is_nan() and t.rp.is_nan():
+            return True
+        return False
+    __rmul__ = __radd__
+
+    ################################ Various ##################################
+    def __round__(self, t):
+        """Exception: Decimal('1').__round__(-100000000000000000000000000)
+           Should it really be InvalidOperation?"""
+        if t.rc is None and t.rp.is_nan():
+            return True
+        return False
+
+shandler = SkipHandler()
+def skip_error(t):
+    return getattr(shandler, t.funcname, shandler.default)(t)
+
+
+# ======================================================================
+#                      Handling verification errors
+# ======================================================================
+
+class VerifyError(Exception):
+    """Verification failed."""
+    pass
+
+def function_as_string(t):
+    if t.contextfunc:
+        cargs = t.cop
+        pargs = t.pop
+        cfunc = "c_func: %s(" % t.funcname
+        pfunc = "p_func: %s(" % t.funcname
+    else:
+        cself, cargs = t.cop[0], t.cop[1:]
+        pself, pargs = t.pop[0], t.pop[1:]
+        cfunc = "c_func: %s.%s(" % (repr(cself), t.funcname)
+        pfunc = "p_func: %s.%s(" % (repr(pself), t.funcname)
+
+    err = cfunc
+    for arg in cargs:
+        err += "%s, " % repr(arg)
+    err = err.rstrip(", ")
+    err += ")\n"
+
+    err += pfunc
+    for arg in pargs:
+        err += "%s, " % repr(arg)
+    err = err.rstrip(", ")
+    err += ")"
+
+    return err
+
+def raise_error(t):
+    global EXIT_STATUS
+
+    if skip_error(t):
+        return
+    EXIT_STATUS = 1
+
+    err = "Error in %s:\n\n" % t.funcname
+    err += "input operands: %s\n\n" % (t.op,)
+    err += function_as_string(t)
+    err += "\n\nc_result: %s\np_result: %s\n\n" % (t.cresults, t.presults)
+    err += "c_exceptions: %s\np_exceptions: %s\n\n" % (t.cex, t.pex)
+    err += "%s\n\n" % str(t.context)
+
+    raise VerifyError(err)
+
+
+# ======================================================================
+#                        Main testing functions
+#
+#  The procedure is always (t is the TestSet):
+#
+#   convert(t) -> Initialize the TestSet as necessary.
+#
+#                 Return 0 for early abortion (e.g. if a TypeError
+#                 occurs during conversion, there is nothing to test).
+#
+#                 Return 1 for continuing with the test case.
+#
+#   callfuncs(t) -> Call the relevant function for each implementation
+#                   and record the results in the TestSet.
+#
+#   verify(t) -> Verify the results. If verification fails, details
+#                are printed to stdout.
+# ======================================================================
+
+def convert(t, convstr=True):
+    """ t is the testset. At this stage the testset contains a tuple of
+        operands t.op of various types. For decimal methods the first
+        operand (self) is always converted to Decimal. If 'convstr' is
+        true, string operands are converted as well.
+
+        Context operands are of type deccheck.Context, rounding mode
+        operands are given as a tuple (C.rounding, P.rounding).
+
+        Other types (float, int, etc.) are left unchanged.
+    """
+    for i, op in enumerate(t.op):
+
+        context.clear_status()
+
+        if not t.contextfunc and i == 0 or \
+           convstr and isinstance(op, str):
+            try:
+                c = C.Decimal(op)
+                cex = None
+            except (TypeError, ValueError, OverflowError) as e:
+                c = None
+                cex = e.__class__
+
+            try:
+                p = RestrictedDecimal(op)
+                pex = None
+            except (TypeError, ValueError, OverflowError) as e:
+                p = None
+                pex = e.__class__
+
+            t.cop.append(c)
+            t.cex.append(cex)
+            t.pop.append(p)
+            t.pex.append(pex)
+
+            if cex is pex:
+                if str(c) != str(p) or not context.assert_eq_status():
+                    raise_error(t)
+                if cex and pex:
+                    # nothing to test
+                    return 0
+            else:
+                raise_error(t)
+
+        elif isinstance(op, Context):
+            t.context = op
+            t.cop.append(op.c)
+            t.pop.append(op.p)
+
+        elif op in RoundModes:
+            t.cop.append(op[0])
+            t.pop.append(op[1])
+
+        else:
+            t.cop.append(op)
+            t.pop.append(op)
+
+    return 1
+
+def callfuncs(t):
+    """ t is the testset. At this stage the testset contains operand lists
+        t.cop and t.pop for the C and Python versions of decimal.
+        For Decimal methods, the first operands are of type C.Decimal and
+        P.Decimal respectively. The remaining operands can have various types.
+        For Context methods, all operands can have any type.
+
+        t.rc and t.rp are the results of the operation.
+    """
+    context.clear_status()
+
+    try:
+        if t.contextfunc:
+            cargs = t.cop
+            t.rc = getattr(context.c, t.funcname)(*cargs)
+        else:
+            cself = t.cop[0]
+            cargs = t.cop[1:]
+            t.rc = getattr(cself, t.funcname)(*cargs)
+        t.cex.append(None)
+    except (TypeError, ValueError, OverflowError, MemoryError) as e:
+        t.rc = None
+        t.cex.append(e.__class__)
+
+    try:
+        if t.contextfunc:
+            pargs = t.pop
+            t.rp = getattr(context.p, t.funcname)(*pargs)
+        else:
+            pself = t.pop[0]
+            pargs = t.pop[1:]
+            t.rp = getattr(pself, t.funcname)(*pargs)
+        t.pex.append(None)
+    except (TypeError, ValueError, OverflowError, MemoryError) as e:
+        t.rp = None
+        t.pex.append(e.__class__)
+
+def verify(t, stat):
+    """ t is the testset. At this stage the testset contains the following
+        tuples:
+
+            t.op: original operands
+            t.cop: C.Decimal operands (see convert for details)
+            t.pop: P.Decimal operands (see convert for details)
+            t.rc: C result
+            t.rp: Python result
+
+        t.rc and t.rp can have various types.
+    """
+    t.cresults.append(str(t.rc))
+    t.presults.append(str(t.rp))
+    if isinstance(t.rc, C.Decimal) and isinstance(t.rp, P.Decimal):
+        # General case: both results are Decimals.
+        t.cresults.append(t.rc.to_eng_string())
+        t.cresults.append(t.rc.as_tuple())
+        t.cresults.append(str(t.rc.imag))
+        t.cresults.append(str(t.rc.real))
+        t.presults.append(t.rp.to_eng_string())
+        t.presults.append(t.rp.as_tuple())
+        t.presults.append(str(t.rp.imag))
+        t.presults.append(str(t.rp.real))
+
+        nc = t.rc.number_class().lstrip('+-s')
+        stat[nc] += 1
+    else:
+        # Results from e.g. __divmod__ can only be compared as strings.
+        if not isinstance(t.rc, tuple) and not isinstance(t.rp, tuple):
+            if t.rc != t.rp:
+                raise_error(t)
+        stat[type(t.rc).__name__] += 1
+
+    # The return value lists must be equal.
+    if t.cresults != t.presults:
+        raise_error(t)
+    # The Python exception lists (TypeError, etc.) must be equal.
+    if t.cex != t.pex:
+        raise_error(t)
+    # The context flags must be equal.
+    if not t.context.assert_eq_status():
+        raise_error(t)
+
+
+# ======================================================================
+#                           Main test loops
+#
+#  test_method(method, testspecs, testfunc) ->
+#
+#     Loop through various context settings. The degree of
+#     thoroughness is determined by 'testspec'. For each
+#     setting, call 'testfunc'. Generally, 'testfunc' itself
+#     a loop, iterating through many test cases generated
+#     by the functions in randdec.py.
+#
+#  test_n-ary(method, prec, exp_range, restricted_range, itr, stat) ->
+#
+#     'test_unary', 'test_binary' and 'test_ternary' are the
+#     main test functions passed to 'test_method'. They deal
+#     with the regular cases. The thoroughness of testing is
+#     determined by 'itr'.
+#
+#     'prec', 'exp_range' and 'restricted_range' are passed
+#     to the test-generating functions and limit the generated
+#     values. In some cases, for reasonable run times a
+#     maximum exponent of 9999 is required.
+#
+#     The 'stat' parameter is passed down to the 'verify'
+#     function, which records statistics for the result values.
+# ======================================================================
+
+def log(fmt, args=None):
+    if args:
+        sys.stdout.write(''.join((fmt, '\n')) % args)
+    else:
+        sys.stdout.write(''.join((str(fmt), '\n')))
+    sys.stdout.flush()
+
+def test_method(method, testspecs, testfunc):
+    """Iterate a test function through many context settings."""
+    log("testing %s ...", method)
+    stat = defaultdict(int)
+    for spec in testspecs:
+        if 'samples' in spec:
+            spec['prec'] = sorted(random.sample(range(1, 101),
+                                  spec['samples']))
+        for prec in spec['prec']:
+            context.prec = prec
+            for expts in spec['expts']:
+                emin, emax = expts
+                if emin == 'rand':
+                    context.Emin = random.randrange(-1000, 0)
+                    context.Emax = random.randrange(prec, 1000)
+                else:
+                    context.Emin, context.Emax = emin, emax
+                if prec > context.Emax: continue
+                log("    prec: %d  emin: %d  emax: %d",
+                    (context.prec, context.Emin, context.Emax))
+                restr_range = 9999 if context.Emax > 9999 else context.Emax+99
+                for rounding in sorted(RoundMap):
+                    context.rounding = rounding
+                    context.capitals = random.randrange(2)
+                    if spec['clamp'] == 'rand':
+                        context.clamp = random.randrange(2)
+                    else:
+                        context.clamp = spec['clamp']
+                    exprange = context.c.Emax
+                    testfunc(method, prec, exprange, restr_range,
+                             spec['iter'], stat)
+    log("    result types: %s" % sorted([t for t in stat.items()]))
+
+def test_unary(method, prec, exp_range, restricted_range, itr, stat):
+    """Iterate a unary function through many test cases."""
+    if method in UnaryRestricted:
+        exp_range = restricted_range
+    for op in all_unary(prec, exp_range, itr):
+        t = TestSet(method, op)
+        try:
+            if not convert(t):
+                continue
+            callfuncs(t)
+            verify(t, stat)
+        except VerifyError as err:
+            log(err)
+
+def test_binary(method, prec, exp_range, restricted_range, itr, stat):
+    """Iterate a binary function through many test cases."""
+    if method in BinaryRestricted:
+        exp_range = restricted_range
+    for op in all_binary(prec, exp_range, itr):
+        t = TestSet(method, op)
+        try:
+            if not convert(t):
+                continue
+            callfuncs(t)
+            verify(t, stat)
+        except VerifyError as err:
+            log(err)
+
+def test_ternary(method, prec, exp_range, restricted_range, itr, stat):
+    """Iterate a ternary function through many test cases."""
+    if method in TernaryRestricted:
+        exp_range = restricted_range
+    for op in all_ternary(prec, exp_range, itr):
+        t = TestSet(method, op)
+        try:
+            if not convert(t):
+                continue
+            callfuncs(t)
+            verify(t, stat)
+        except VerifyError as err:
+            log(err)
+
+def test_format(method, prec, exp_range, restricted_range, itr, stat):
+    """Iterate the __format__ method through many test cases."""
+    for op in all_unary(prec, exp_range, itr):
+        fmt1 = rand_format(chr(random.randrange(32, 128)), 'EeGgn')
+        fmt2 = rand_locale()
+        for fmt in (fmt1, fmt2):
+            fmtop = (op[0], fmt)
+            t = TestSet(method, fmtop)
+            try:
+                if not convert(t, convstr=False):
+                    continue
+                callfuncs(t)
+                verify(t, stat)
+            except VerifyError as err:
+                log(err)
+    for op in all_unary(prec, 9999, itr):
+        fmt1 = rand_format(chr(random.randrange(32, 128)), 'Ff%')
+        fmt2 = rand_locale()
+        for fmt in (fmt1, fmt2):
+            fmtop = (op[0], fmt)
+            t = TestSet(method, fmtop)
+            try:
+                if not convert(t, convstr=False):
+                    continue
+                callfuncs(t)
+                verify(t, stat)
+            except VerifyError as err:
+                log(err)
+
+def test_round(method, prec, exprange, restricted_range, itr, stat):
+    """Iterate the __round__ method through many test cases."""
+    for op in all_unary(prec, 9999, itr):
+        n = random.randrange(10)
+        roundop = (op[0], n)
+        t = TestSet(method, roundop)
+        try:
+            if not convert(t):
+                continue
+            callfuncs(t)
+            verify(t, stat)
+        except VerifyError as err:
+            log(err)
+
+def test_from_float(method, prec, exprange, restricted_range, itr, stat):
+    """Iterate the __float__ method through many test cases."""
+    for rounding in sorted(RoundMap):
+        context.rounding = rounding
+        for i in range(1000):
+            f = randfloat()
+            op = (f,) if method.startswith("context.") else ("sNaN", f)
+            t = TestSet(method, op)
+            try:
+                if not convert(t):
+                    continue
+                callfuncs(t)
+                verify(t, stat)
+            except VerifyError as err:
+                log(err)
+
+def randcontext(exprange):
+    c = Context(C.Context(), P.Context())
+    c.Emax = random.randrange(1, exprange+1)
+    c.Emin = random.randrange(-exprange, 0)
+    maxprec = 100 if c.Emax >= 100 else c.Emax
+    c.prec = random.randrange(1, maxprec+1)
+    c.clamp = random.randrange(2)
+    c.clear_traps()
+    return c
+
+def test_quantize_api(method, prec, exprange, restricted_range, itr, stat):
+    """Iterate the 'quantize' method through many test cases, using
+       the optional arguments."""
+    for op in all_binary(prec, restricted_range, itr):
+        for rounding in RoundModes:
+            c = randcontext(exprange)
+            quantizeop = (op[0], op[1], rounding, c)
+            t = TestSet(method, quantizeop)
+            try:
+                if not convert(t):
+                    continue
+                callfuncs(t)
+                verify(t, stat)
+            except VerifyError as err:
+                log(err)
+
+
+def check_untested(funcdict, c_cls, p_cls):
+    """Determine untested, C-only and Python-only attributes.
+       Uncomment print lines for debugging."""
+    c_attr = set(dir(c_cls))
+    p_attr = set(dir(p_cls))
+    intersect = c_attr & p_attr
+
+    funcdict['c_only'] = tuple(sorted(c_attr-intersect))
+    funcdict['p_only'] = tuple(sorted(p_attr-intersect))
+
+    tested = set()
+    for lst in funcdict.values():
+        for v in lst:
+            v = v.replace("context.", "") if c_cls == C.Context else v
+            tested.add(v)
+
+    funcdict['untested'] = tuple(sorted(intersect-tested))
+
+    #for key in ('untested', 'c_only', 'p_only'):
+    #    s = 'Context' if c_cls == C.Context else 'Decimal'
+    #    print("\n%s %s:\n%s" % (s, key, funcdict[key]))
+
+
+if __name__ == '__main__':
+
+    import time
+
+    randseed = int(time.time())
+    random.seed(randseed)
+
+    # Set up the testspecs list. A testspec is simply a dictionary
+    # that determines the amount of different contexts that 'test_method'
+    # will generate.
+    base_expts = [(C.MIN_EMIN, C.MAX_EMAX)]
+    if C.MAX_EMAX == 999999999999999999:
+        base_expts.append((-999999999, 999999999))
+
+    # Basic contexts.
+    base = {
+        'expts': base_expts,
+        'prec': [],
+        'clamp': 'rand',
+        'iter': None,
+        'samples': None,
+    }
+    # Contexts with small values for prec, emin, emax.
+    small = {
+        'prec': [1, 2, 3, 4, 5],
+        'expts': [(-1, 1), (-2, 2), (-3, 3), (-4, 4), (-5, 5)],
+        'clamp': 'rand',
+        'iter': None
+    }
+    # IEEE interchange format.
+    ieee = [
+        # DECIMAL32
+        {'prec': [7], 'expts': [(-95, 96)], 'clamp': 1, 'iter': None},
+        # DECIMAL64
+        {'prec': [16], 'expts': [(-383, 384)], 'clamp': 1, 'iter': None},
+        # DECIMAL128
+        {'prec': [34], 'expts': [(-6143, 6144)], 'clamp': 1, 'iter': None}
+    ]
+
+    if '--medium' in sys.argv:
+        base['expts'].append(('rand', 'rand'))
+        # 5 random precisions
+        base['samples'] = 5
+        testspecs = [small] + ieee + [base]
+    if '--long' in sys.argv:
+        base['expts'].append(('rand', 'rand'))
+        # 10 random precisions
+        base['samples'] = 10
+        testspecs = [small] + ieee + [base]
+    elif '--all' in sys.argv:
+        base['expts'].append(('rand', 'rand'))
+        # All precisions in [1, 100]
+        base['samples'] = 100
+        testspecs = [small] + ieee + [base]
+    else: # --short
+        rand_ieee = random.choice(ieee)
+        base['iter'] = small['iter'] = rand_ieee['iter'] = 1
+        # 1 random precision and exponent pair
+        base['samples'] = 1
+        base['expts'] = [random.choice(base_expts)]
+        # 1 random precision and exponent pair
+        prec = random.randrange(1, 6)
+        small['prec'] = [prec]
+        small['expts'] = [(-prec, prec)]
+        testspecs = [small, rand_ieee, base]
+
+    check_untested(Functions, C.Decimal, P.Decimal)
+    check_untested(ContextFunctions, C.Context, P.Context)
+
+
+    log("\n\nRandom seed: %d\n\n", randseed)
+
+    # Decimal methods:
+    for method in Functions['unary'] + Functions['unary_ctx'] + \
+                  Functions['unary_rnd_ctx']:
+        test_method(method, testspecs, test_unary)
+
+    for method in Functions['binary'] + Functions['binary_ctx']:
+        test_method(method, testspecs, test_binary)
+
+    for method in Functions['ternary'] + Functions['ternary_ctx']:
+        test_method(method, testspecs, test_ternary)
+
+    test_method('__format__', testspecs, test_format)
+    test_method('__round__', testspecs, test_round)
+    test_method('from_float', testspecs, test_from_float)
+    test_method('quantize', testspecs, test_quantize_api)
+
+    # Context methods:
+    for method in ContextFunctions['unary']:
+        test_method(method, testspecs, test_unary)
+
+    for method in ContextFunctions['binary']:
+        test_method(method, testspecs, test_binary)
+
+    for method in ContextFunctions['ternary']:
+        test_method(method, testspecs, test_ternary)
+
+    test_method('context.create_decimal_from_float', testspecs, test_from_float)
+
+
+    sys.exit(EXIT_STATUS)
diff --git a/Modules/_decimal/tests/formathelper.py b/Modules/_decimal/tests/formathelper.py
new file mode 100644
index 0000000..9da2b37
--- /dev/null
+++ b/Modules/_decimal/tests/formathelper.py
@@ -0,0 +1,344 @@
+#
+# Copyright (c) 2008-2012 Stefan Krah. All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+
+
+# Generate PEP-3101 format strings.
+
+
+import os, sys, locale, random
+import platform, subprocess
+from test.support import import_fresh_module
+from distutils.spawn import find_executable
+
+C = import_fresh_module('decimal', fresh=['_decimal'])
+P = import_fresh_module('decimal', blocked=['_decimal'])
+
+
+windows_lang_strings = [
+  "chinese", "chinese-simplified", "chinese-traditional", "czech", "danish",
+  "dutch", "belgian", "english", "australian", "canadian", "english-nz",
+  "english-uk", "english-us", "finnish", "french", "french-belgian",
+  "french-canadian", "french-swiss", "german", "german-austrian",
+  "german-swiss", "greek", "hungarian", "icelandic", "italian", "italian-swiss",
+  "japanese", "korean", "norwegian", "norwegian-bokmal", "norwegian-nynorsk",
+  "polish", "portuguese", "portuguese-brazil", "russian", "slovak", "spanish",
+  "spanish-mexican", "spanish-modern", "swedish", "turkish",
+]
+
+preferred_encoding = {
+  'cs_CZ': 'ISO8859-2',
+  'cs_CZ.iso88592': 'ISO8859-2',
+  'czech': 'ISO8859-2',
+  'eesti': 'ISO8859-1',
+  'estonian': 'ISO8859-1',
+  'et_EE': 'ISO8859-15',
+  'et_EE.ISO-8859-15': 'ISO8859-15',
+  'et_EE.iso885915': 'ISO8859-15',
+  'et_EE.iso88591': 'ISO8859-1',
+  'fi_FI.iso88591': 'ISO8859-1',
+  'fi_FI': 'ISO8859-15',
+  'fi_FI@euro': 'ISO8859-15',
+  'fi_FI.iso885915@euro': 'ISO8859-15',
+  'finnish': 'ISO8859-1',
+  'lv_LV': 'ISO8859-13',
+  'lv_LV.iso885913': 'ISO8859-13',
+  'nb_NO': 'ISO8859-1',
+  'nb_NO.iso88591': 'ISO8859-1',
+  'bokmal': 'ISO8859-1',
+  'nn_NO': 'ISO8859-1',
+  'nn_NO.iso88591': 'ISO8859-1',
+  'no_NO': 'ISO8859-1',
+  'norwegian': 'ISO8859-1',
+  'nynorsk': 'ISO8859-1',
+  'ru_RU': 'ISO8859-5',
+  'ru_RU.iso88595': 'ISO8859-5',
+  'russian': 'ISO8859-5',
+  'ru_RU.KOI8-R': 'KOI8-R',
+  'ru_RU.koi8r': 'KOI8-R',
+  'ru_RU.CP1251': 'CP1251',
+  'ru_RU.cp1251': 'CP1251',
+  'sk_SK': 'ISO8859-2',
+  'sk_SK.iso88592': 'ISO8859-2',
+  'slovak': 'ISO8859-2',
+  'sv_FI': 'ISO8859-1',
+  'sv_FI.iso88591': 'ISO8859-1',
+  'sv_FI@euro': 'ISO8859-15',
+  'sv_FI.iso885915@euro': 'ISO8859-15',
+  'uk_UA': 'KOI8-U',
+  'uk_UA.koi8u': 'KOI8-U'
+}
+
+integers = [
+  "",
+  "1",
+  "12",
+  "123",
+  "1234",
+  "12345",
+  "123456",
+  "1234567",
+  "12345678",
+  "123456789",
+  "1234567890",
+  "12345678901",
+  "123456789012",
+  "1234567890123",
+  "12345678901234",
+  "123456789012345",
+  "1234567890123456",
+  "12345678901234567",
+  "123456789012345678",
+  "1234567890123456789",
+  "12345678901234567890",
+  "123456789012345678901",
+  "1234567890123456789012",
+]
+
+numbers = [
+  "0", "-0", "+0",
+  "0.0", "-0.0", "+0.0",
+  "0e0", "-0e0", "+0e0",
+  ".0", "-.0",
+  ".1", "-.1",
+  "1.1", "-1.1",
+  "1e1", "-1e1"
+]
+
+# Get the list of available locales.
+if platform.system() == 'Windows':
+    locale_list = windows_lang_strings
+else:
+    locale_list = ['C']
+    if os.path.isfile("/var/lib/locales/supported.d/local"):
+        # On Ubuntu, `locale -a` gives the wrong case for some locales,
+        # so we get the correct names directly:
+        with open("/var/lib/locales/supported.d/local") as f:
+            locale_list = [loc.split()[0] for loc in f.readlines() \
+                           if not loc.startswith('#')]
+    elif find_executable('locale'):
+        locale_list = subprocess.Popen(["locale", "-a"],
+                          stdout=subprocess.PIPE).communicate()[0]
+        try:
+            locale_list = locale_list.decode()
+        except UnicodeDecodeError:
+            # Some distributions insist on using latin-1 characters
+            # in their locale names.
+            locale_list = locale_list.decode('latin-1')
+        locale_list = locale_list.split('\n')
+try:
+    locale_list.remove('')
+except ValueError:
+    pass
+
+# Debian
+if os.path.isfile("/etc/locale.alias"):
+    with open("/etc/locale.alias") as f:
+        while 1:
+            try:
+                line = f.readline()
+            except UnicodeDecodeError:
+                continue
+            if line == "":
+                break
+            if line.startswith('#'):
+                continue
+            x = line.split()
+            if len(x) == 2:
+                if x[0] in locale_list:
+                    locale_list.remove(x[0])
+
+# FreeBSD
+if platform.system() == 'FreeBSD':
+    # http://www.freebsd.org/cgi/query-pr.cgi?pr=142173
+    # en_GB.US-ASCII has 163 as the currency symbol.
+    for loc in ['it_CH.ISO8859-1', 'it_CH.ISO8859-15', 'it_CH.UTF-8',
+                'it_IT.ISO8859-1', 'it_IT.ISO8859-15', 'it_IT.UTF-8',
+                'sl_SI.ISO8859-2', 'sl_SI.UTF-8',
+                'en_GB.US-ASCII']:
+        try:
+            locale_list.remove(loc)
+        except ValueError:
+            pass
+
+# Print a testcase in the format of the IBM tests (for runtest.c):
+def get_preferred_encoding():
+    loc = locale.setlocale(locale.LC_CTYPE)
+    if loc in preferred_encoding:
+        return preferred_encoding[loc]
+    else:
+        return locale.getpreferredencoding()
+
+def printit(testno, s, fmt, encoding=None):
+    if not encoding:
+        encoding = get_preferred_encoding()
+    try:
+        result = format(P.Decimal(s), fmt)
+        fmt = str(fmt.encode(encoding))[2:-1]
+        result = str(result.encode(encoding))[2:-1]
+        if "'" in result:
+            sys.stdout.write("xfmt%d  format  %s  '%s'  ->  \"%s\"\n"
+                             % (testno, s, fmt, result))
+        else:
+            sys.stdout.write("xfmt%d  format  %s  '%s'  ->  '%s'\n"
+                             % (testno, s, fmt, result))
+    except Exception as err:
+        sys.stderr.write("%s  %s  %s\n" % (err, s, fmt))
+
+
+# Check if an integer can be converted to a valid fill character.
+def check_fillchar(i):
+    try:
+        c = chr(i)
+        c.encode('utf-8').decode()
+        format(P.Decimal(0), c + '<19g')
+        if c in ("'", '"', '\\'):
+            return None
+        return c
+    except:
+        return None
+
+# Generate all unicode characters that are accepted as
+# fill characters by decimal.py.
+def all_fillchars():
+    for i in range(32, 0x110002):
+        c = check_fillchar(i)
+        if c: yield c
+
+# Return random fill character.
+def rand_fillchar():
+    while 1:
+        i = random.randrange(32, 0x110002)
+        c = check_fillchar(i)
+        if c: return c
+
+# Generate random format strings
+# [[fill]align][sign][#][0][width][.precision][type]
+def rand_format(fill, typespec='EeGgFfn%'):
+    active = sorted(random.sample(range(7), random.randrange(8)))
+    have_align = 0
+    s = ''
+    for elem in active:
+        if elem == 0: # fill+align
+            s += fill
+            s += random.choice('<>=^')
+            have_align = 1
+        elif elem == 1: # sign
+            s += random.choice('+- ')
+        elif elem == 2 and not have_align: # zeropad
+            s += '0'
+        elif elem == 3: # width
+            s += str(random.randrange(1, 100))
+        elif elem == 4: # thousands separator
+            s += ','
+        elif elem == 5: # prec
+            s += '.'
+            s += str(random.randrange(100))
+        elif elem == 6:
+            if 4 in active: c = typespec.replace('n', '')
+            else: c = typespec
+            s += random.choice(c)
+    return s
+
+# Partially brute force all possible format strings containing a thousands
+# separator. Fall back to random where the runtime would become excessive.
+# [[fill]align][sign][#][0][width][,][.precision][type]
+def all_format_sep():
+    for align in ('', '<', '>', '=', '^'):
+        for fill in ('', 'x'):
+            if align == '': fill = ''
+            for sign in ('', '+', '-', ' '):
+                for zeropad in ('', '0'):
+                    if align != '': zeropad = ''
+                    for width in ['']+[str(y) for y in range(1, 15)]+['101']:
+                        for prec in ['']+['.'+str(y) for y in range(15)]:
+                            # for type in ('', 'E', 'e', 'G', 'g', 'F', 'f', '%'):
+                            type = random.choice(('', 'E', 'e', 'G', 'g', 'F', 'f', '%'))
+                            yield ''.join((fill, align, sign, zeropad, width, ',', prec, type))
+
+# Partially brute force all possible format strings with an 'n' specifier.
+# [[fill]align][sign][#][0][width][,][.precision][type]
+def all_format_loc():
+    for align in ('', '<', '>', '=', '^'):
+        for fill in ('', 'x'):
+            if align == '': fill = ''
+            for sign in ('', '+', '-', ' '):
+                for zeropad in ('', '0'):
+                    if align != '': zeropad = ''
+                    for width in ['']+[str(y) for y in range(1, 20)]+['101']:
+                        for prec in ['']+['.'+str(y) for y in range(1, 20)]:
+                            yield ''.join((fill, align, sign, zeropad, width, prec, 'n'))
+
+# Generate random format strings with a unicode fill character
+# [[fill]align][sign][#][0][width][,][.precision][type]
+def randfill(fill):
+    active = sorted(random.sample(range(5), random.randrange(6)))
+    s = ''
+    s += str(fill)
+    s += random.choice('<>=^')
+    for elem in active:
+        if elem == 0: # sign
+            s += random.choice('+- ')
+        elif elem == 1: # width
+            s += str(random.randrange(1, 100))
+        elif elem == 2: # thousands separator
+            s += ','
+        elif elem == 3: # prec
+            s += '.'
+            s += str(random.randrange(100))
+        elif elem == 4:
+            if 2 in active: c = 'EeGgFf%'
+            else: c = 'EeGgFfn%'
+            s += random.choice(c)
+    return s
+
+# Generate random format strings with random locale setting
+# [[fill]align][sign][#][0][width][,][.precision][type]
+def rand_locale():
+    try:
+        loc = random.choice(locale_list)
+        locale.setlocale(locale.LC_ALL, loc)
+    except locale.Error as err:
+        pass
+    active = sorted(random.sample(range(5), random.randrange(6)))
+    s = ''
+    have_align = 0
+    for elem in active:
+        if elem == 0: # fill+align
+            s += chr(random.randrange(32, 128))
+            s += random.choice('<>=^')
+            have_align = 1
+        elif elem == 1: # sign
+            s += random.choice('+- ')
+        elif elem == 2 and not have_align: # zeropad
+            s += '0'
+        elif elem == 3: # width
+            s += str(random.randrange(1, 100))
+        elif elem == 4: # prec
+            s += '.'
+            s += str(random.randrange(100))
+    s += 'n'
+    return s
diff --git a/Modules/_decimal/tests/randdec.py b/Modules/_decimal/tests/randdec.py
new file mode 100644
index 0000000..80858bf
--- /dev/null
+++ b/Modules/_decimal/tests/randdec.py
@@ -0,0 +1,559 @@
+#!/usr/bin/env python
+
+#
+# Copyright (c) 2008-2012 Stefan Krah. All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+
+
+# Generate test cases for deccheck.py.
+
+
+#
+# Grammar from http://speleotrove.com/decimal/daconvs.html
+#
+# sign           ::=  '+' | '-'
+# digit          ::=  '0' | '1' | '2' | '3' | '4' | '5' | '6' | '7' |
+#                     '8' | '9'
+# indicator      ::=  'e' | 'E'
+# digits         ::=  digit [digit]...
+# decimal-part   ::=  digits '.' [digits] | ['.'] digits
+# exponent-part  ::=  indicator [sign] digits
+# infinity       ::=  'Infinity' | 'Inf'
+# nan            ::=  'NaN' [digits] | 'sNaN' [digits]
+# numeric-value  ::=  decimal-part [exponent-part] | infinity
+# numeric-string ::=  [sign] numeric-value | [sign] nan
+#
+
+
+from random import randrange, sample
+from fractions import Fraction
+from randfloat import un_randfloat, bin_randfloat, tern_randfloat
+
+
+def sign():
+    if randrange(2):
+        if randrange(2): return '+'
+        return ''
+    return '-'
+
+def indicator():
+    return "eE"[randrange(2)]
+
+def digits(maxprec):
+    if maxprec == 0: return ''
+    return str(randrange(10**maxprec))
+
+def dot():
+    if randrange(2): return '.'
+    return ''
+
+def decimal_part(maxprec):
+    if randrange(100) > 60: # integers
+        return digits(maxprec)
+    if randrange(2):
+        intlen = randrange(1, maxprec+1)
+        fraclen = maxprec-intlen
+        intpart = digits(intlen)
+        fracpart = digits(fraclen)
+        return ''.join((intpart, '.', fracpart))
+    else:
+        return ''.join((dot(), digits(maxprec)))
+
+def expdigits(maxexp):
+    return str(randrange(maxexp))
+
+def exponent_part(maxexp):
+    return ''.join((indicator(), sign(), expdigits(maxexp)))
+
+def infinity():
+    if randrange(2): return 'Infinity'
+    return 'Inf'
+
+def nan():
+    d = ''
+    if randrange(2):
+        d = digits(randrange(99))
+    if randrange(2):
+        return ''.join(('NaN', d))
+    else:
+        return ''.join(('sNaN', d))
+
+def numeric_value(maxprec, maxexp):
+    if randrange(100) > 90:
+        return infinity()
+    exp_part = ''
+    if randrange(100) > 60:
+        exp_part = exponent_part(maxexp)
+    return ''.join((decimal_part(maxprec), exp_part))
+
+def numeric_string(maxprec, maxexp):
+    if randrange(100) > 95:
+        return ''.join((sign(), nan()))
+    else:
+        return ''.join((sign(), numeric_value(maxprec, maxexp)))
+
+def randdec(maxprec, maxexp):
+    return numeric_string(maxprec, maxexp)
+
+def rand_adjexp(maxprec, maxadjexp):
+    d = digits(maxprec)
+    maxexp = maxadjexp-len(d)+1
+    if maxexp == 0: maxexp = 1
+    exp = str(randrange(maxexp-2*(abs(maxexp)), maxexp))
+    return ''.join((sign(), d, 'E', exp))
+
+
+def ndigits(n):
+    if n < 1: return 0
+    return randrange(10**(n-1), 10**n)
+
+def randtuple(maxprec, maxexp):
+    n = randrange(100)
+    sign = randrange(2)
+    coeff = ndigits(maxprec)
+    if n >= 95:
+        coeff = ()
+        exp = 'F'
+    elif n >= 85:
+        coeff = tuple(map(int, str(ndigits(maxprec))))
+        exp = "nN"[randrange(2)]
+    else:
+        coeff = tuple(map(int, str(ndigits(maxprec))))
+        exp = randrange(-maxexp, maxexp)
+    return (sign, coeff, exp)
+
+def from_triple(sign, coeff, exp):
+    return ''.join((str(sign*coeff), indicator(), str(exp)))
+
+
+# Close to 10**n
+def un_close_to_pow10(prec, maxexp, itr=None):
+    if itr is None:
+        lst = range(prec+30)
+    else:
+        lst = sample(range(prec+30), itr)
+    nines = [10**n - 1 for n in lst]
+    pow10 = [10**n for n in lst]
+    for coeff in nines:
+        yield coeff
+        yield -coeff
+        yield from_triple(1, coeff, randrange(2*maxexp))
+        yield from_triple(-1, coeff, randrange(2*maxexp))
+    for coeff in pow10:
+        yield coeff
+        yield -coeff
+
+# Close to 10**n
+def bin_close_to_pow10(prec, maxexp, itr=None):
+    if itr is None:
+        lst = range(prec+30)
+    else:
+        lst = sample(range(prec+30), itr)
+    nines = [10**n - 1 for n in lst]
+    pow10 = [10**n for n in lst]
+    for coeff in nines:
+        yield coeff, 1
+        yield -coeff, -1
+        yield 1, coeff
+        yield -1, -coeff
+        yield from_triple(1, coeff, randrange(2*maxexp)), 1
+        yield from_triple(-1, coeff, randrange(2*maxexp)), -1
+        yield 1, from_triple(1, coeff, -randrange(2*maxexp))
+        yield -1, from_triple(-1, coeff, -randrange(2*maxexp))
+    for coeff in pow10:
+        yield coeff, -1
+        yield -coeff, 1
+        yield 1, -coeff
+        yield -coeff, 1
+
+# Close to 1:
+def close_to_one_greater(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("1.", '0'*randrange(prec),
+                   str(randrange(rprec))))
+
+def close_to_one_less(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("0.9", '9'*randrange(prec),
+                   str(randrange(rprec))))
+
+# Close to 0:
+def close_to_zero_greater(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("0.", '0'*randrange(prec),
+                   str(randrange(rprec))))
+
+def close_to_zero_less(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("-0.", '0'*randrange(prec),
+                   str(randrange(rprec))))
+
+# Close to emax:
+def close_to_emax_less(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("9.", '9'*randrange(prec),
+                   str(randrange(rprec)), "E", str(emax)))
+
+def close_to_emax_greater(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("1.", '0'*randrange(prec),
+                   str(randrange(rprec)), "E", str(emax+1)))
+
+# Close to emin:
+def close_to_emin_greater(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("1.", '0'*randrange(prec),
+                   str(randrange(rprec)), "E", str(emin)))
+
+def close_to_emin_less(prec, emax, emin):
+    rprec = 10**prec
+    return ''.join(("9.", '9'*randrange(prec),
+                   str(randrange(rprec)), "E", str(emin-1)))
+
+# Close to etiny:
+def close_to_etiny_greater(prec, emax, emin):
+    rprec = 10**prec
+    etiny = emin - (prec - 1)
+    return ''.join(("1.", '0'*randrange(prec),
+                   str(randrange(rprec)), "E", str(etiny)))
+
+def close_to_etiny_less(prec, emax, emin):
+    rprec = 10**prec
+    etiny = emin - (prec - 1)
+    return ''.join(("9.", '9'*randrange(prec),
+                   str(randrange(rprec)), "E", str(etiny-1)))
+
+
+def close_to_min_etiny_greater(prec, max_prec, min_emin):
+    rprec = 10**prec
+    etiny = min_emin - (max_prec - 1)
+    return ''.join(("1.", '0'*randrange(prec),
+                   str(randrange(rprec)), "E", str(etiny)))
+
+def close_to_min_etiny_less(prec, max_prec, min_emin):
+    rprec = 10**prec
+    etiny = min_emin - (max_prec - 1)
+    return ''.join(("9.", '9'*randrange(prec),
+                   str(randrange(rprec)), "E", str(etiny-1)))
+
+
+close_funcs = [
+  close_to_one_greater, close_to_one_less, close_to_zero_greater,
+  close_to_zero_less, close_to_emax_less, close_to_emax_greater,
+  close_to_emin_greater, close_to_emin_less, close_to_etiny_greater,
+  close_to_etiny_less, close_to_min_etiny_greater, close_to_min_etiny_less
+]
+
+
+def un_close_numbers(prec, emax, emin, itr=None):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func in close_funcs:
+            yield func(prec, emax, emin)
+
+def bin_close_numbers(prec, emax, emin, itr=None):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func1 in close_funcs:
+            for func2 in close_funcs:
+                yield func1(prec, emax, emin), func2(prec, emax, emin)
+        for func in close_funcs:
+            yield randdec(prec, emax), func(prec, emax, emin)
+            yield func(prec, emax, emin), randdec(prec, emax)
+
+def tern_close_numbers(prec, emax, emin, itr):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func1 in close_funcs:
+            for func2 in close_funcs:
+                for func3 in close_funcs:
+                    yield (func1(prec, emax, emin), func2(prec, emax, emin),
+                           func3(prec, emax, emin))
+        for func in close_funcs:
+            yield (randdec(prec, emax), func(prec, emax, emin),
+                   func(prec, emax, emin))
+            yield (func(prec, emax, emin), randdec(prec, emax),
+                   func(prec, emax, emin))
+            yield (func(prec, emax, emin), func(prec, emax, emin),
+                   randdec(prec, emax))
+        for func in close_funcs:
+            yield (randdec(prec, emax), randdec(prec, emax),
+                   func(prec, emax, emin))
+            yield (randdec(prec, emax), func(prec, emax, emin),
+                   randdec(prec, emax))
+            yield (func(prec, emax, emin), randdec(prec, emax),
+                   randdec(prec, emax))
+
+
+# If itr == None, test all digit lengths up to prec + 30
+def un_incr_digits(prec, maxexp, itr):
+    if itr is None:
+        lst = range(prec+30)
+    else:
+        lst = sample(range(prec+30), itr)
+    for m in lst:
+        yield from_triple(1, ndigits(m), 0)
+        yield from_triple(-1, ndigits(m), 0)
+        yield from_triple(1, ndigits(m), randrange(maxexp))
+        yield from_triple(-1, ndigits(m), randrange(maxexp))
+
+# If itr == None, test all digit lengths up to prec + 30
+# Also output decimals im tuple form.
+def un_incr_digits_tuple(prec, maxexp, itr):
+    if itr is None:
+        lst = range(prec+30)
+    else:
+        lst = sample(range(prec+30), itr)
+    for m in lst:
+        yield from_triple(1, ndigits(m), 0)
+        yield from_triple(-1, ndigits(m), 0)
+        yield from_triple(1, ndigits(m), randrange(maxexp))
+        yield from_triple(-1, ndigits(m), randrange(maxexp))
+        # test from tuple
+        yield (0, tuple(map(int, str(ndigits(m)))), 0)
+        yield (1, tuple(map(int, str(ndigits(m)))), 0)
+        yield (0, tuple(map(int, str(ndigits(m)))), randrange(maxexp))
+        yield (1, tuple(map(int, str(ndigits(m)))), randrange(maxexp))
+
+# If itr == None, test all combinations of digit lengths up to prec + 30
+def bin_incr_digits(prec, maxexp, itr):
+    if itr is None:
+        lst1 = range(prec+30)
+        lst2 = range(prec+30)
+    else:
+        lst1 = sample(range(prec+30), itr)
+        lst2 = sample(range(prec+30), itr)
+    for m in lst1:
+        x = from_triple(1, ndigits(m), 0)
+        yield x, x
+        x = from_triple(-1, ndigits(m), 0)
+        yield x, x
+        x = from_triple(1, ndigits(m), randrange(maxexp))
+        yield x, x
+        x = from_triple(-1, ndigits(m), randrange(maxexp))
+        yield x, x
+    for m in lst1:
+        for n in lst2:
+            x = from_triple(1, ndigits(m), 0)
+            y = from_triple(1, ndigits(n), 0)
+            yield x, y
+            x = from_triple(-1, ndigits(m), 0)
+            y = from_triple(1, ndigits(n), 0)
+            yield x, y
+            x = from_triple(1, ndigits(m), 0)
+            y = from_triple(-1, ndigits(n), 0)
+            yield x, y
+            x = from_triple(-1, ndigits(m), 0)
+            y = from_triple(-1, ndigits(n), 0)
+            yield x, y
+            x = from_triple(1, ndigits(m), randrange(maxexp))
+            y = from_triple(1, ndigits(n), randrange(maxexp))
+            yield x, y
+            x = from_triple(-1, ndigits(m), randrange(maxexp))
+            y = from_triple(1, ndigits(n), randrange(maxexp))
+            yield x, y
+            x = from_triple(1, ndigits(m), randrange(maxexp))
+            y = from_triple(-1, ndigits(n), randrange(maxexp))
+            yield x, y
+            x = from_triple(-1, ndigits(m), randrange(maxexp))
+            y = from_triple(-1, ndigits(n), randrange(maxexp))
+            yield x, y
+
+
+def randsign():
+    return (1, -1)[randrange(2)]
+
+# If itr == None, test all combinations of digit lengths up to prec + 30
+def tern_incr_digits(prec, maxexp, itr):
+    if itr is None:
+        lst1 = range(prec+30)
+        lst2 = range(prec+30)
+        lst3 = range(prec+30)
+    else:
+        lst1 = sample(range(prec+30), itr)
+        lst2 = sample(range(prec+30), itr)
+        lst3 = sample(range(prec+30), itr)
+    for m in lst1:
+        for n in lst2:
+            for p in lst3:
+                x = from_triple(randsign(), ndigits(m), 0)
+                y = from_triple(randsign(), ndigits(n), 0)
+                z = from_triple(randsign(), ndigits(p), 0)
+                yield x, y, z
+
+
+# Tests for the 'logical' functions
+def bindigits(prec):
+    z = 0
+    for i in range(prec):
+        z += randrange(2) * 10**i
+    return z
+
+def logical_un_incr_digits(prec, itr):
+    if itr is None:
+        lst = range(prec+30)
+    else:
+        lst = sample(range(prec+30), itr)
+    for m in lst:
+        yield from_triple(1, bindigits(m), 0)
+
+def logical_bin_incr_digits(prec, itr):
+    if itr is None:
+        lst1 = range(prec+30)
+        lst2 = range(prec+30)
+    else:
+        lst1 = sample(range(prec+30), itr)
+        lst2 = sample(range(prec+30), itr)
+    for m in lst1:
+        x = from_triple(1, bindigits(m), 0)
+        yield x, x
+    for m in lst1:
+        for n in lst2:
+            x = from_triple(1, bindigits(m), 0)
+            y = from_triple(1, bindigits(n), 0)
+            yield x, y
+
+
+def randint():
+    p = randrange(1, 100)
+    return ndigits(p) * (1,-1)[randrange(2)]
+
+def randfloat():
+    p = randrange(1, 100)
+    s = numeric_value(p, 383)
+    try:
+        f = float(numeric_value(p, 383))
+    except ValueError:
+        f = 0.0
+    return f
+
+def randcomplex():
+    real = randfloat()
+    if randrange(100) > 30:
+        imag = 0.0
+    else:
+        imag = randfloat()
+    return complex(real, imag)
+
+def randfraction():
+    num = randint()
+    denom = randint()
+    if denom == 0:
+        denom = 1
+    return Fraction(num, denom)
+
+number_funcs = [randint, randfloat, randcomplex, randfraction]
+
+def un_random_mixed_op(itr=None):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func in number_funcs:
+            yield func()
+    # Test garbage input
+    for x in (['x'], ('y',), {'z'}, {1:'z'}):
+        yield x
+
+def bin_random_mixed_op(prec, emax, emin, itr=None):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func in number_funcs:
+            yield randdec(prec, emax), func()
+            yield func(), randdec(prec, emax)
+        for number in number_funcs:
+            for dec in close_funcs:
+                yield dec(prec, emax, emin), number()
+    # Test garbage input
+    for x in (['x'], ('y',), {'z'}, {1:'z'}):
+        for y in (['x'], ('y',), {'z'}, {1:'z'}):
+            yield x, y
+
+def tern_random_mixed_op(prec, emax, emin, itr):
+    if itr is None:
+        itr = 1000
+    for _ in range(itr):
+        for func in number_funcs:
+            yield randdec(prec, emax), randdec(prec, emax), func()
+            yield randdec(prec, emax), func(), func()
+            yield func(), func(), func()
+    # Test garbage input
+    for x in (['x'], ('y',), {'z'}, {1:'z'}):
+        for y in (['x'], ('y',), {'z'}, {1:'z'}):
+            for z in (['x'], ('y',), {'z'}, {1:'z'}):
+                yield x, y, z
+
+def all_unary(prec, exp_range, itr):
+    for a in un_close_to_pow10(prec, exp_range, itr):
+        yield (a,)
+    for a in un_close_numbers(prec, exp_range, -exp_range, itr):
+        yield (a,)
+    for a in un_incr_digits_tuple(prec, exp_range, itr):
+        yield (a,)
+    for a in un_randfloat():
+        yield (a,)
+    for a in un_random_mixed_op(itr):
+        yield (a,)
+    for a in logical_un_incr_digits(prec, itr):
+        yield (a,)
+    for _ in range(100):
+        yield (randdec(prec, exp_range),)
+    for _ in range(100):
+        yield (randtuple(prec, exp_range),)
+
+def all_binary(prec, exp_range, itr):
+    for a, b in bin_close_to_pow10(prec, exp_range, itr):
+        yield a, b
+    for a, b in bin_close_numbers(prec, exp_range, -exp_range, itr):
+        yield a, b
+    for a, b in bin_incr_digits(prec, exp_range, itr):
+        yield a, b
+    for a, b in bin_randfloat():
+        yield a, b
+    for a, b in bin_random_mixed_op(prec, exp_range, -exp_range, itr):
+        yield a, b
+    for a, b in logical_bin_incr_digits(prec, itr):
+        yield a, b
+    for _ in range(100):
+        yield randdec(prec, exp_range), randdec(prec, exp_range)
+
+def all_ternary(prec, exp_range, itr):
+    for a, b, c in tern_close_numbers(prec, exp_range, -exp_range, itr):
+        yield a, b, c
+    for a, b, c in tern_incr_digits(prec, exp_range, itr):
+        yield a, b, c
+    for a, b, c in tern_randfloat():
+        yield a, b, c
+    for a, b, c in tern_random_mixed_op(prec, exp_range, -exp_range, itr):
+        yield a, b, c
+    for _ in range(100):
+        a = randdec(prec, 2*exp_range)
+        b = randdec(prec, 2*exp_range)
+        c = randdec(prec, 2*exp_range)
+        yield a, b, c
diff --git a/Modules/_decimal/tests/randfloat.py b/Modules/_decimal/tests/randfloat.py
new file mode 100644
index 0000000..1687303
--- /dev/null
+++ b/Modules/_decimal/tests/randfloat.py
@@ -0,0 +1,250 @@
+# Copyright (c) 2010 Python Software Foundation. All Rights Reserved.
+# Adapted from Python's Lib/test/test_strtod.py (by Mark Dickinson)
+
+# More test cases for deccheck.py.
+
+import random
+
+TEST_SIZE = 2
+
+
+def test_short_halfway_cases():
+    # exact halfway cases with a small number of significant digits
+    for k in 0, 5, 10, 15, 20:
+        # upper = smallest integer >= 2**54/5**k
+        upper = -(-2**54//5**k)
+        # lower = smallest odd number >= 2**53/5**k
+        lower = -(-2**53//5**k)
+        if lower % 2 == 0:
+            lower += 1
+        for i in range(10 * TEST_SIZE):
+            # Select a random odd n in [2**53/5**k,
+            # 2**54/5**k). Then n * 10**k gives a halfway case
+            # with small number of significant digits.
+            n, e = random.randrange(lower, upper, 2), k
+
+            # Remove any additional powers of 5.
+            while n % 5 == 0:
+                n, e = n // 5, e + 1
+            assert n % 10 in (1, 3, 7, 9)
+
+            # Try numbers of the form n * 2**p2 * 10**e, p2 >= 0,
+            # until n * 2**p2 has more than 20 significant digits.
+            digits, exponent = n, e
+            while digits < 10**20:
+                s = '{}e{}'.format(digits, exponent)
+                yield s
+                # Same again, but with extra trailing zeros.
+                s = '{}e{}'.format(digits * 10**40, exponent - 40)
+                yield s
+                digits *= 2
+
+            # Try numbers of the form n * 5**p2 * 10**(e - p5), p5
+            # >= 0, with n * 5**p5 < 10**20.
+            digits, exponent = n, e
+            while digits < 10**20:
+                s = '{}e{}'.format(digits, exponent)
+                yield s
+                # Same again, but with extra trailing zeros.
+                s = '{}e{}'.format(digits * 10**40, exponent - 40)
+                yield s
+                digits *= 5
+                exponent -= 1
+
+def test_halfway_cases():
+    # test halfway cases for the round-half-to-even rule
+    for i in range(1000):
+        for j in range(TEST_SIZE):
+            # bit pattern for a random finite positive (or +0.0) float
+            bits = random.randrange(2047*2**52)
+
+            # convert bit pattern to a number of the form m * 2**e
+            e, m = divmod(bits, 2**52)
+            if e:
+                m, e = m + 2**52, e - 1
+            e -= 1074
+
+            # add 0.5 ulps
+            m, e = 2*m + 1, e - 1
+
+            # convert to a decimal string
+            if e >= 0:
+                digits = m << e
+                exponent = 0
+            else:
+                # m * 2**e = (m * 5**-e) * 10**e
+                digits = m * 5**-e
+                exponent = e
+            s = '{}e{}'.format(digits, exponent)
+            yield s
+
+def test_boundaries():
+    # boundaries expressed as triples (n, e, u), where
+    # n*10**e is an approximation to the boundary value and
+    # u*10**e is 1ulp
+    boundaries = [
+        (10000000000000000000, -19, 1110),   # a power of 2 boundary (1.0)
+        (17976931348623159077, 289, 1995),   # overflow boundary (2.**1024)
+        (22250738585072013831, -327, 4941),  # normal/subnormal (2.**-1022)
+        (0, -327, 4941),                     # zero
+        ]
+    for n, e, u in boundaries:
+        for j in range(1000):
+            for i in range(TEST_SIZE):
+                digits = n + random.randrange(-3*u, 3*u)
+                exponent = e
+                s = '{}e{}'.format(digits, exponent)
+                yield s
+            n *= 10
+            u *= 10
+            e -= 1
+
+def test_underflow_boundary():
+    # test values close to 2**-1075, the underflow boundary; similar
+    # to boundary_tests, except that the random error doesn't scale
+    # with n
+    for exponent in range(-400, -320):
+        base = 10**-exponent // 2**1075
+        for j in range(TEST_SIZE):
+            digits = base + random.randrange(-1000, 1000)
+            s = '{}e{}'.format(digits, exponent)
+            yield s
+
+def test_bigcomp():
+    for ndigs in 5, 10, 14, 15, 16, 17, 18, 19, 20, 40, 41, 50:
+        dig10 = 10**ndigs
+        for i in range(100 * TEST_SIZE):
+            digits = random.randrange(dig10)
+            exponent = random.randrange(-400, 400)
+            s = '{}e{}'.format(digits, exponent)
+            yield s
+
+def test_parsing():
+    # make '0' more likely to be chosen than other digits
+    digits = '000000123456789'
+    signs = ('+', '-', '')
+
+    # put together random short valid strings
+    # \d*[.\d*]?e
+    for i in range(1000):
+        for j in range(TEST_SIZE):
+            s = random.choice(signs)
+            intpart_len = random.randrange(5)
+            s += ''.join(random.choice(digits) for _ in range(intpart_len))
+            if random.choice([True, False]):
+                s += '.'
+                fracpart_len = random.randrange(5)
+                s += ''.join(random.choice(digits)
+                             for _ in range(fracpart_len))
+            else:
+                fracpart_len = 0
+            if random.choice([True, False]):
+                s += random.choice(['e', 'E'])
+                s += random.choice(signs)
+                exponent_len = random.randrange(1, 4)
+                s += ''.join(random.choice(digits)
+                             for _ in range(exponent_len))
+
+            if intpart_len + fracpart_len:
+                yield s
+
+test_particular = [
+     # squares
+    '1.00000000100000000025',
+    '1.0000000000000000000000000100000000000000000000000' #...
+    '00025',
+    '1.0000000000000000000000000000000000000000000010000' #...
+    '0000000000000000000000000000000000000000025',
+    '1.0000000000000000000000000000000000000000000000000' #...
+    '000001000000000000000000000000000000000000000000000' #...
+    '000000000025',
+    '0.99999999900000000025',
+    '0.9999999999999999999999999999999999999999999999999' #...
+    '999000000000000000000000000000000000000000000000000' #...
+    '000025',
+    '0.9999999999999999999999999999999999999999999999999' #...
+    '999999999999999999999999999999999999999999999999999' #...
+    '999999999999999999999999999999999999999990000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '0000000000000000000000000000025',
+
+    '1.0000000000000000000000000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '100000000000000000000000000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000001',
+    '1.0000000000000000000000000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '500000000000000000000000000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000005',
+    '1.0000000000000000000000000000000000000000000000000' #...
+    '000000000100000000000000000000000000000000000000000' #...
+    '000000000000000000250000000000000002000000000000000' #...
+    '000000000000000000000000000000000000000000010000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '0000000000000000001',
+    '1.0000000000000000000000000000000000000000000000000' #...
+    '000000000100000000000000000000000000000000000000000' #...
+    '000000000000000000249999999999999999999999999999999' #...
+    '999999999999979999999999999999999999999999999999999' #...
+    '999999999999999999999900000000000000000000000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '00000000000000000000000001',
+
+    '0.9999999999999999999999999999999999999999999999999' #...
+    '999999999900000000000000000000000000000000000000000' #...
+    '000000000000000000249999999999999998000000000000000' #...
+    '000000000000000000000000000000000000000000010000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '0000000000000000001',
+    '0.9999999999999999999999999999999999999999999999999' #...
+    '999999999900000000000000000000000000000000000000000' #...
+    '000000000000000000250000001999999999999999999999999' #...
+    '999999999999999999999999999999999990000000000000000' #...
+    '000000000000000000000000000000000000000000000000000' #...
+    '1',
+
+    # tough cases for ln etc.
+    '1.000000000000000000000000000000000000000000000000' #...
+    '00000000000000000000000000000000000000000000000000' #...
+    '00100000000000000000000000000000000000000000000000' #...
+    '00000000000000000000000000000000000000000000000000' #...
+    '0001',
+    '0.999999999999999999999999999999999999999999999999' #...
+    '99999999999999999999999999999999999999999999999999' #...
+    '99899999999999999999999999999999999999999999999999' #...
+    '99999999999999999999999999999999999999999999999999' #...
+    '99999999999999999999999999999999999999999999999999' #...
+    '9999'
+    ]
+
+
+TESTCASES = [
+      [x for x in test_short_halfway_cases()],
+      [x for x in test_halfway_cases()],
+      [x for x in test_boundaries()],
+      [x for x in test_underflow_boundary()],
+      [x for x in test_bigcomp()],
+      [x for x in test_parsing()],
+      test_particular
+]
+
+def un_randfloat():
+    for i in range(1000):
+        l = random.choice(TESTCASES[:6])
+        yield random.choice(l)
+    for v in test_particular:
+        yield v
+
+def bin_randfloat():
+    for i in range(1000):
+        l1 = random.choice(TESTCASES)
+        l2 = random.choice(TESTCASES)
+        yield random.choice(l1), random.choice(l2)
+
+def tern_randfloat():
+    for i in range(1000):
+        l1 = random.choice(TESTCASES)
+        l2 = random.choice(TESTCASES)
+        l3 = random.choice(TESTCASES)
+        yield random.choice(l1), random.choice(l2), random.choice(l3)
diff --git a/Modules/_decimal/tests/runall-memorydebugger.sh b/Modules/_decimal/tests/runall-memorydebugger.sh
new file mode 100644
index 0000000..44ff3b3
--- /dev/null
+++ b/Modules/_decimal/tests/runall-memorydebugger.sh
@@ -0,0 +1,175 @@
+#!/bin/sh
+
+#
+# Purpose: test with and without threads, all machine configurations, pydebug,
+#          refleaks, release build and release build with valgrind.
+#
+# Synopsis: ./runall-memorydebugger.sh [--all-configs64 | --all-configs32]
+#
+# Requirements: valgrind
+#
+
+# Set additional CFLAGS for ./configure
+ADD_CFLAGS=
+
+
+CONFIGS_64="x64 uint128 ansi64 universal"
+CONFIGS_32="ppro ansi32 ansi-legacy universal"
+
+VALGRIND="valgrind --tool=memcheck --leak-resolution=high \
+          --db-attach=yes --suppressions=Misc/valgrind-python.supp"
+
+# Get args
+case $@ in
+     *--all-configs64*)
+         CONFIGS=$CONFIGS_64
+         ;;
+     *--all-configs32*)
+         CONFIGS=$CONFIGS_32
+         ;;
+     *)
+         CONFIGS="auto"
+         ;;
+esac
+
+# gmake required
+GMAKE=`which gmake`
+if [ X"$GMAKE" = X"" ]; then
+    GMAKE=make
+fi
+
+# Pretty print configurations
+print_config ()
+{
+    len=`echo $@ | wc -c`
+    margin="#%"`expr \( 74 - $len \) / 2`"s"
+
+    echo ""
+    echo "# ========================================================================"
+    printf $margin ""
+    echo $@
+    echo "# ========================================================================"
+    echo ""
+}
+
+
+cd ..
+
+# test_decimal: refleak, regular and Valgrind tests
+for args in "--without-threads" ""; do
+    for config in $CONFIGS; do
+
+        unset PYTHON_DECIMAL_WITH_MACHINE
+        libmpdec_config=$config
+        if [ X"$config" != X"auto" ]; then
+            PYTHON_DECIMAL_WITH_MACHINE=$config
+            export PYTHON_DECIMAL_WITH_MACHINE
+        else
+            libmpdec_config=""
+        fi
+
+        ############ refleak tests ###########
+        print_config "refleak tests: config=$config" $args
+        printf "\nbuilding python ...\n\n"
+
+        cd ../../
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" --with-pydebug $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ======================== refleak tests ===========================\n\n"
+        ./python -m test -uall -R 2:2 test_decimal
+
+
+        ############ regular tests ###########
+        print_config "regular tests: config=$config" $args
+        printf "\nbuilding python ...\n\n"
+
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ======================== regular tests ===========================\n\n"
+        ./python -m test -uall test_decimal
+
+
+        ########### valgrind tests ###########
+        valgrind=$VALGRIND
+        case "$config" in
+            # Valgrind has no support for 80 bit long double arithmetic.
+            ppro) valgrind= ;;
+            auto) case `uname -m` in
+                      i386|i486|i586|i686) valgrind= ;;
+                  esac
+        esac
+
+        print_config "valgrind tests: config=$config" $args
+        printf "\nbuilding python ...\n\n"
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" --without-pymalloc $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ======================== valgrind tests ===========================\n\n"
+        $valgrind ./python -m test -uall test_decimal
+
+        cd Modules/_decimal
+    done
+done
+
+# deccheck
+cd ../../
+for config in $CONFIGS; do
+    for args in "--without-threads" ""; do
+
+        unset PYTHON_DECIMAL_WITH_MACHINE
+        if [ X"$config" != X"auto" ]; then
+            PYTHON_DECIMAL_WITH_MACHINE=$config
+            export PYTHON_DECIMAL_WITH_MACHINE
+        fi
+
+        ############ debug ############
+        print_config "deccheck: config=$config --with-pydebug" $args
+        printf "\nbuilding python ...\n\n"
+
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" --with-pydebug $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ========================== debug ===========================\n\n"
+        ./python Modules/_decimal/tests/deccheck.py
+
+        ########### regular ###########
+        print_config "deccheck: config=$config " $args
+        printf "\nbuilding python ...\n\n"
+
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ======================== regular ===========================\n\n"
+        ./python Modules/_decimal/tests/deccheck.py
+
+        ########### valgrind ###########
+        valgrind=$VALGRIND
+        case "$config" in
+            # Valgrind has no support for 80 bit long double arithmetic.
+            ppro) valgrind= ;;
+            auto) case `uname -m` in
+                      i386|i486|i586|i686) valgrind= ;;
+                  esac
+        esac
+
+        print_config "valgrind deccheck: config=$config " $args
+        printf "\nbuilding python ...\n\n"
+
+        $GMAKE distclean > /dev/null 2>&1
+        ./configure CFLAGS="$ADD_CFLAGS" --without-pymalloc $args > /dev/null 2>&1
+        $GMAKE | grep _decimal
+
+        printf "\n\n# ======================== valgrind ==========================\n\n"
+        $valgrind ./python Modules/_decimal/tests/deccheck.py
+    done
+done
+
+
+
diff --git a/Modules/_decimal/tests/runall.bat b/Modules/_decimal/tests/runall.bat
new file mode 100644
index 0000000..0d02e25
--- /dev/null
+++ b/Modules/_decimal/tests/runall.bat
@@ -0,0 +1,121 @@
+@ECHO OFF
+
+rem Test all machine configurations, pydebug, refleaks, release build.
+
+cd ..
+
+call vcvarsall x64
+echo.
+echo # ======================================================================
+echo #                       test_decimal: platform=x64
+echo # ======================================================================
+echo.
+
+cd ..\..\PCbuild
+echo # ==================== refleak tests =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Debug|x64" > NUL 2>&1
+amd64\python_d.exe -m test -uall -R 2:2 test_decimal
+echo.
+
+echo # ==================== regular tests =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Release|x64" > NUL 2>&1
+amd64\python.exe -m test -uall test_decimal
+echo.
+echo.
+
+
+call vcvarsall x86
+echo.
+echo # ======================================================================
+echo #                       test_decimal: platform=x86
+echo # ======================================================================
+echo.
+
+echo # ==================== refleak tests =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Debug|win32" > NUL 2>&1
+python_d.exe -m test -uall -R 2:2 test_decimal
+echo.
+
+echo # ==================== regular tests =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Release|win32" > NUL 2>&1
+python.exe -m test -uall test_decimal
+echo.
+echo.
+
+
+call vcvarsall x64
+echo.
+echo # ======================================================================
+echo #                         deccheck: platform=x64
+echo # ======================================================================
+echo.
+echo.
+echo # ==================== debug build =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Debug|x64" > NUL 2>&1
+amd64\python_d.exe ..\Modules\_decimal\tests\deccheck.py
+echo.
+echo.
+
+echo # =================== release build ======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Release|x64" > NUL 2>&1
+amd64\python.exe ..\Modules\_decimal\tests\deccheck.py
+echo.
+echo.
+
+
+call vcvarsall x86
+echo.
+echo # ======================================================================
+echo #                         deccheck: platform=x86
+echo # ======================================================================
+echo.
+echo.
+echo # ==================== debug build =======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Debug|win32" > NUL 2>&1
+python_d.exe ..\Modules\_decimal\tests\deccheck.py
+echo.
+echo.
+
+echo # =================== release build ======================
+echo.
+echo building python ...
+echo.
+vcbuild /clean pcbuild.sln > NUL 2>&1
+vcbuild pcbuild.sln "Release|win32" > NUL 2>&1
+python.exe ..\Modules\_decimal\tests\deccheck.py
+echo.
+echo.
+
+
+cd ..\Modules\_decimal\tests
+
+
+