Merged revisions 58095-58132,58136-58148,58151-58197 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk

........
  r58096 | brett.cannon | 2007-09-10 23:38:27 +0200 (Mon, 10 Sep 2007) | 4 lines

  Fix a possible segfault from recursing too deep to get the repr of a list.

  Closes issue #1096.
........
  r58097 | bill.janssen | 2007-09-10 23:51:02 +0200 (Mon, 10 Sep 2007) | 33 lines

  More work on SSL support.

  * Much expanded test suite:

    All protocols tested against all other protocols.
    All protocols tested with all certificate options.
    Tests for bad key and bad cert.
    Test of STARTTLS functionality.
    Test of RAND_* functions.

  * Fixes for threading/malloc bug.

  * Issue 1065 fixed:

    sslsocket class renamed to SSLSocket.
    sslerror class renamed to SSLError.
    Function "wrap_socket" now used to wrap an existing socket.

  * Issue 1583946 finally fixed:

    Support for subjectAltName added.
    Subject name now returned as proper DN list of RDNs.

  * SSLError exported from socket as "sslerror".

  * RAND_* functions properly exported from ssl.py.

  * Documentation improved:

    Example of how to create a self-signed certificate.
    Better indexing.
........
  r58098 | guido.van.rossum | 2007-09-11 00:02:25 +0200 (Tue, 11 Sep 2007) | 9 lines

  Patch # 1140 (my code, approved by Effbot).

  Make sure the type of the return value of re.sub(x, y, z) is the type
  of y+x (i.e. unicode if either is unicode, str if they are both str)
  even if there are no substitutions or if x==z (which triggered various
  special cases in join_list()).

  Could be backported to 2.5; no need to port to 3.0.
........
  r58099 | guido.van.rossum | 2007-09-11 00:36:02 +0200 (Tue, 11 Sep 2007) | 8 lines

  Patch # 1026 by Benjamin Aranguren (with Alex Martelli):
  Backport abc.py and isinstance/issubclass overloading to 2.6.

  I had to backport test_typechecks.py myself, and make one small change
  to abc.py to avoid duplicate work when x.__class__ and type(x) are the
  same.
........
  r58100 | bill.janssen | 2007-09-11 01:41:24 +0200 (Tue, 11 Sep 2007) | 3 lines

  A better way of finding an open port to test with.
........
  r58101 | bill.janssen | 2007-09-11 03:09:19 +0200 (Tue, 11 Sep 2007) | 4 lines

  Make sure test_ssl doesn't reference the ssl module in a
  context where it can't be imported.
........
  r58102 | bill.janssen | 2007-09-11 04:42:07 +0200 (Tue, 11 Sep 2007) | 3 lines

  Fix some documentation bugs.
........
  r58103 | nick.coghlan | 2007-09-11 16:01:18 +0200 (Tue, 11 Sep 2007) | 1 line

  Always use the -E flag when spawning subprocesses in test_cmd_line (Issue 1056)
........
  r58106 | thomas.heller | 2007-09-11 21:17:48 +0200 (Tue, 11 Sep 2007) | 3 lines

  Disable some tests that fail on the 'ppc Debian unstable' buildbot to
  find out if they cause the segfault on the 'alpha Debian' machine.
........
  r58108 | brett.cannon | 2007-09-11 23:02:28 +0200 (Tue, 11 Sep 2007) | 6 lines

  Generators had their throw() method allowing string exceptions.  That's a
  no-no.

  Fixes issue #1147.  Need to fix 2.5 to raise a proper warning if a string
  exception is passed in.
........
  r58112 | georg.brandl | 2007-09-12 20:03:51 +0200 (Wed, 12 Sep 2007) | 3 lines

  New documentation page for the bdb module.
  (This doesn't need to be merged to Py3k.)
........
  r58114 | georg.brandl | 2007-09-12 20:05:57 +0200 (Wed, 12 Sep 2007) | 2 lines

  Bug #1152: use non-deprecated name in example.
........
  r58115 | georg.brandl | 2007-09-12 20:08:33 +0200 (Wed, 12 Sep 2007) | 2 lines

  Fix #1122: wrong return type documented for various _Size() functions.
........
  r58117 | georg.brandl | 2007-09-12 20:10:56 +0200 (Wed, 12 Sep 2007) | 2 lines

  Fix #1139: PyFile_Encoding really is PyFile_SetEncoding.
........
  r58119 | georg.brandl | 2007-09-12 20:29:18 +0200 (Wed, 12 Sep 2007) | 2 lines

  bug #1154: release memory allocated by "es" PyArg_ParseTuple format specifier.
........
  r58121 | bill.janssen | 2007-09-12 20:52:05 +0200 (Wed, 12 Sep 2007) | 1 line

  root certificate for https://svn.python.org/, used in test_ssl
........
  r58122 | georg.brandl | 2007-09-12 21:00:07 +0200 (Wed, 12 Sep 2007) | 3 lines

  Bug #1153: repr.repr() now doesn't require set and dictionary items
  to be orderable to properly represent them.
........
  r58125 | georg.brandl | 2007-09-12 21:29:28 +0200 (Wed, 12 Sep 2007) | 4 lines

  #1120: put explicit version in the shebang lines of pydoc, idle
  and smtpd.py scripts that are installed by setup.py. That way, they
  work when only "make altinstall" is used.
........
  r58139 | mark.summerfield | 2007-09-13 16:54:30 +0200 (Thu, 13 Sep 2007) | 9 lines

  Replaced variable o with obj in operator.rst because o is easy to
  confuse.

  Added a note about Python 3's collections.Mapping etc., above section
  that describes isMappingType() etc.

  Added xrefs between os, os.path, fileinput, and open().
........
  r58143 | facundo.batista | 2007-09-13 20:13:15 +0200 (Thu, 13 Sep 2007) | 7 lines


  Merged the decimal-branch (revisions 54886 to 58140). Decimal is now
  fully updated to the latests Decimal Specification (v1.66) and the
  latests test cases (v2.56).

  Thanks to Mark Dickinson for all his help during this process.
........
  r58145 | facundo.batista | 2007-09-13 20:42:09 +0200 (Thu, 13 Sep 2007) | 7 lines


  Put the parameter watchexp back in (changed watchexp from an int
  to a bool).  Also second argument to watchexp is now converted
  to Decimal, just as with all the other two-argument operations.

  Thanks Mark Dickinson.
........
  r58147 | andrew.kuchling | 2007-09-14 00:49:34 +0200 (Fri, 14 Sep 2007) | 1 line

  Add various items
........
  r58148 | andrew.kuchling | 2007-09-14 00:50:10 +0200 (Fri, 14 Sep 2007) | 1 line

  Make target unique
........
  r58154 | facundo.batista | 2007-09-14 20:58:34 +0200 (Fri, 14 Sep 2007) | 3 lines


  Included the new functions, and new descriptions.
........
  r58155 | thomas.heller | 2007-09-14 21:40:35 +0200 (Fri, 14 Sep 2007) | 2 lines

  ctypes.util.find_library uses dump(1) instead of objdump(1) on Solaris.
  Fixes issue #1777530; will backport to release25-maint.
........
  r58159 | facundo.batista | 2007-09-14 23:29:52 +0200 (Fri, 14 Sep 2007) | 3 lines


  Some additions (examples and a bit on the tutorial).
........
  r58160 | georg.brandl | 2007-09-15 18:53:36 +0200 (Sat, 15 Sep 2007) | 2 lines

  Remove bdb from the "undocumented modules" list.
........
  r58164 | bill.janssen | 2007-09-17 00:06:00 +0200 (Mon, 17 Sep 2007) | 15 lines

  Add support for asyncore server-side SSL support.  This requires
  adding the 'makefile' method to ssl.SSLSocket, and importing the
  requisite fakefile class from socket.py, and making the appropriate
  changes to it to make it use the SSL connection.

  Added sample HTTPS server to test_ssl.py, and test that uses it.

  Change SSL tests to use https://svn.python.org/, instead of
  www.sf.net and pop.gmail.com.

  Added utility function to ssl module, get_server_certificate,
  to wrap up the several things to be done to pull a certificate
  from a remote server.
........
  r58173 | bill.janssen | 2007-09-17 01:16:46 +0200 (Mon, 17 Sep 2007) | 1 line

  use binary mode when reading files for testAsyncore to make Windows happy
........
  r58175 | raymond.hettinger | 2007-09-17 02:55:00 +0200 (Mon, 17 Sep 2007) | 7 lines

  Sync-up named tuples with the latest version of the ASPN recipe.
  Allows optional commas in the field-name spec (help when named tuples are used in conjuction with sql queries).
  Adds the __fields__ attribute for introspection and to support conversion to dictionary form.
  Adds a  __replace__() method similar to str.replace() but using a named field as a target.
  Clean-up spelling and presentation in doc-strings.
........
  r58176 | brett.cannon | 2007-09-17 05:28:34 +0200 (Mon, 17 Sep 2007) | 5 lines

  Add a bunch of GIL release/acquire points in tp_print implementations and for
  PyObject_Print().

  Closes issue #1164.
........
  r58177 | sean.reifschneider | 2007-09-17 07:45:04 +0200 (Mon, 17 Sep 2007) | 2 lines

  issue1597011: Fix for bz2 module corner-case error due to error checking bug.
........
  r58180 | facundo.batista | 2007-09-17 18:26:50 +0200 (Mon, 17 Sep 2007) | 3 lines


  Decimal is updated, :)
........
  r58181 | facundo.batista | 2007-09-17 19:30:13 +0200 (Mon, 17 Sep 2007) | 5 lines


  The methods always return Decimal classes, even if they're
  executed through a subclass (thanks Mark Dickinson).
  Added a bit of testing for this.
........
  r58183 | sean.reifschneider | 2007-09-17 22:53:21 +0200 (Mon, 17 Sep 2007) | 2 lines

  issue1082: Fixing platform and system for Vista.
........
  r58185 | andrew.kuchling | 2007-09-18 03:36:16 +0200 (Tue, 18 Sep 2007) | 1 line

  Add item; sort properly
........
  r58186 | raymond.hettinger | 2007-09-18 05:33:19 +0200 (Tue, 18 Sep 2007) | 1 line

  Handle corner cased on 0-tuples and 1-tuples.  Add verbose option so people can see how it works.
........
  r58192 | georg.brandl | 2007-09-18 09:24:40 +0200 (Tue, 18 Sep 2007) | 2 lines

  A bit of reordering, also show more subheadings in the lang ref index.
........
  r58193 | facundo.batista | 2007-09-18 18:53:18 +0200 (Tue, 18 Sep 2007) | 4 lines


  Speed up of the various division operations (remainder, divide,
  divideint and divmod). Thanks Mark Dickinson.
........
  r58197 | raymond.hettinger | 2007-09-19 00:18:02 +0200 (Wed, 19 Sep 2007) | 1 line

  Cleanup docs for NamedTuple.
........
diff --git a/Lib/decimal.py b/Lib/decimal.py
index 8afeb4b..d7bd127 100644
--- a/Lib/decimal.py
+++ b/Lib/decimal.py
@@ -128,7 +128,7 @@
 
     # Constants for use in setting up contexts
     'ROUND_DOWN', 'ROUND_HALF_UP', 'ROUND_HALF_EVEN', 'ROUND_CEILING',
-    'ROUND_FLOOR', 'ROUND_UP', 'ROUND_HALF_DOWN',
+    'ROUND_FLOOR', 'ROUND_UP', 'ROUND_HALF_DOWN', 'ROUND_05UP',
 
     # Functions for manipulating contexts
     'setcontext', 'getcontext', 'localcontext'
@@ -144,6 +144,7 @@
 ROUND_FLOOR = 'ROUND_FLOOR'
 ROUND_UP = 'ROUND_UP'
 ROUND_HALF_DOWN = 'ROUND_HALF_DOWN'
+ROUND_05UP = 'ROUND_05UP'
 
 # Rounding decision (not part of the public API)
 NEVER_ROUND = 'NEVER_ROUND'    # Round in division (non-divmod), sqrt ONLY
@@ -204,13 +205,22 @@
     x ** (non-integer)
     x ** (+-)INF
     An operand is invalid
+
+    The result of the operation after these is a quiet positive NaN,
+    except when the cause is a signaling NaN, in which case the result is
+    also a quiet NaN, but with the original sign, and an optional
+    diagnostic information.
     """
     def handle(self, context, *args):
         if args:
             if args[0] == 1:  # sNaN, must drop 's' but keep diagnostics
-                return Decimal( (args[1]._sign, args[1]._int, 'n') )
+                ans = Decimal((args[1]._sign, args[1]._int, 'n'))
+                return ans._fix_nan(context)
+            elif args[0] == 2:
+                return Decimal( (args[1], args[2], 'n') )
         return NaN
 
+
 class ConversionSyntax(InvalidOperation):
     """Trying to convert badly formed string.
 
@@ -218,9 +228,8 @@
     converted to a number and it does not conform to the numeric string
     syntax.  The result is [0,qNaN].
     """
-
     def handle(self, context, *args):
-        return (0, (0,), 'n')  # Passed to something which uses a tuple.
+        return NaN
 
 class DivisionByZero(DecimalException, ZeroDivisionError):
     """Division by 0.
@@ -235,9 +244,7 @@
     -0, for power.
     """
 
-    def handle(self, context, sign, double = None, *args):
-        if double is not None:
-            return (Infsign[sign],)*2
+    def handle(self, context, sign, *args):
         return Infsign[sign]
 
 class DivisionImpossible(InvalidOperation):
@@ -249,7 +256,7 @@
     """
 
     def handle(self, context, *args):
-        return (NaN, NaN)
+        return NaN
 
 class DivisionUndefined(InvalidOperation, ZeroDivisionError):
     """Undefined result of division.
@@ -259,9 +266,7 @@
     the dividend is also zero.  The result is [0,qNaN].
     """
 
-    def handle(self, context, tup=None, *args):
-        if tup is not None:
-            return (NaN, NaN)  # for 0 %0, 0 // 0
+    def handle(self, context, *args):
         return NaN
 
 class Inexact(DecimalException):
@@ -340,7 +345,7 @@
 
     def handle(self, context, sign, *args):
         if context.rounding in (ROUND_HALF_UP, ROUND_HALF_EVEN,
-                                     ROUND_HALF_DOWN, ROUND_UP):
+                                ROUND_HALF_DOWN, ROUND_UP):
             return Infsign[sign]
         if sign == 0:
             if context.rounding == ROUND_CEILING:
@@ -520,7 +525,7 @@
         Decimal("3.14")
         >>> Decimal((0, (3, 1, 4), -2))  # tuple (sign, digit_tuple, exponent)
         Decimal("3.14")
-        >>> Decimal(314)                 # int or long
+        >>> Decimal(314)                 # int
         Decimal("314")
         >>> Decimal(Decimal(314))        # another decimal instance
         Decimal("314")
@@ -545,7 +550,7 @@
             return self
 
         # From an integer
-        if isinstance(value, (int,int)):
+        if isinstance(value, int):
             if value >= 0:
                 self._sign = 0
             else:
@@ -561,8 +566,8 @@
             if value[0] not in (0,1):
                 raise ValueError('Invalid sign')
             for digit in value[1]:
-                if not isinstance(digit, (int,int)) or digit < 0:
-                    raise ValueError("The second value in the tuple must be"
+                if not isinstance(digit, int) or digit < 0:
+                    raise ValueError("The second value in the tuple must be "
                                 "composed of non negative integer elements.")
             self._sign = value[0]
             self._int  = tuple(value[1])
@@ -596,10 +601,6 @@
             if _isnan(value):
                 sig, sign, diag = _isnan(value)
                 self._is_special = True
-                if len(diag) > context.prec:  # Diagnostic info too long
-                    self._sign, self._int, self._exp = \
-                                context._raise_error(ConversionSyntax)
-                    return self
                 if sig == 1:
                     self._exp = 'n'  # qNaN
                 else:  # sig == 2
@@ -611,9 +612,8 @@
                 self._sign, self._int, self._exp = _string2exact(value)
             except ValueError:
                 self._is_special = True
-                self._sign, self._int, self._exp = \
-                              context._raise_error(ConversionSyntax,
-                                  "Invalid literal for Decimal: %r" % value)
+                return context._raise_error(ConversionSyntax,
+                                   "Invalid literal for Decimal: %r" % value)
             return self
 
         raise TypeError("Cannot convert %r to Decimal" % value)
@@ -622,7 +622,7 @@
         """Returns whether the number is not actually one.
 
         0 if a number
-        1 if NaN
+        1 if NaN  (it could be a normal quiet NaN or a phantom one)
         2 if sNaN
         """
         if self._is_special:
@@ -646,7 +646,7 @@
             return 1
         return 0
 
-    def _check_nans(self, other = None, context=None):
+    def _check_nans(self, other=None, context=None):
         """Returns whether the number is not actually one.
 
         if self, other are sNaN, signal
@@ -673,9 +673,9 @@
                 return context._raise_error(InvalidOperation, 'sNaN',
                                         1, other)
             if self_is_nan:
-                return self
+                return self._fix_nan(context)
 
-            return other
+            return other._fix_nan(context)
         return 0
 
     def __bool__(self):
@@ -688,21 +688,28 @@
             return True
         return sum(self._int) != 0
 
-    def __cmp__(self, other, context=None):
+    def __cmp__(self, other):
         other = _convert_other(other)
         if other is NotImplemented:
-            return other
+            # Never return NotImplemented
+            return 1
 
         if self._is_special or other._is_special:
-            ans = self._check_nans(other, context)
-            if ans:
+            # check for nans, without raising on a signaling nan
+            if self._isnan() or other._isnan():
                 return 1  # Comparison involving NaN's always reports self > other
 
             # INF = INF
             return cmp(self._isinfinity(), other._isinfinity())
 
-        if not self and not other:
-            return 0  # If both 0, sign comparison isn't certain.
+        # check for zeros;  note that cmp(0, -0) should return 0
+        if not self:
+            if not other:
+                return 0
+            else:
+                return -((-1)**other._sign)
+        if not other:
+            return (-1)**self._sign
 
         # If different signs, neg one is less
         if other._sign < self._sign:
@@ -712,35 +719,15 @@
 
         self_adjusted = self.adjusted()
         other_adjusted = other.adjusted()
-        if self_adjusted == other_adjusted and \
-           self._int + (0,)*(self._exp - other._exp) == \
-           other._int + (0,)*(other._exp - self._exp):
-            return 0  # equal, except in precision. ([0]*(-x) = [])
-        elif self_adjusted > other_adjusted and self._int[0] != 0:
+        if self_adjusted == other_adjusted:
+            self_padded = self._int + (0,)*(self._exp - other._exp)
+            other_padded = other._int + (0,)*(other._exp - self._exp)
+            return cmp(self_padded, other_padded) * (-1)**self._sign
+        elif self_adjusted > other_adjusted:
             return (-1)**self._sign
-        elif self_adjusted < other_adjusted and other._int[0] != 0:
+        else: # self_adjusted < other_adjusted
             return -((-1)**self._sign)
 
-        # Need to round, so make sure we have a valid context
-        if context is None:
-            context = getcontext()
-
-        context = context._shallow_copy()
-        rounding = context._set_rounding(ROUND_UP)  # round away from 0
-
-        flags = context._ignore_all_flags()
-        res = self.__sub__(other, context=context)
-
-        context._regard_flags(*flags)
-
-        context.rounding = rounding
-
-        if not res:
-            return 0
-        elif res._sign:
-            return -1
-        return 1
-
     def __eq__(self, other):
         if not isinstance(other, (Decimal, int)):
             return NotImplemented
@@ -780,9 +767,7 @@
         NaN => one is NaN
         Like __cmp__, but returns Decimal instances.
         """
-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
+        other = _convert_other(other, raiseit=True)
 
         # Compare(NaN, NaN) = NaN
         if (self._is_special or other and other._is_special):
@@ -790,7 +775,7 @@
             if ans:
                 return ans
 
-        return Decimal(self.__cmp__(other, context))
+        return Decimal(self.__cmp__(other))
 
     def __hash__(self):
         """x.__hash__() <==> hash(x)"""
@@ -819,7 +804,7 @@
         # Invariant:  eval(repr(d)) == d
         return 'Decimal("%s")' % str(self)
 
-    def __str__(self, eng = 0, context=None):
+    def __str__(self, eng=False, context=None):
         """Return string representation of the number in scientific notation.
 
         Captures all of the information in the underlying representation.
@@ -909,7 +894,7 @@
 
         Same rules for when in exponential and when as a value as in __str__.
         """
-        return self.__str__(eng=1, context=context)
+        return self.__str__(eng=True, context=context)
 
     def __neg__(self, context=None):
         """Returns a copy with the sign switched.
@@ -923,17 +908,15 @@
 
         if not self:
             # -Decimal('0') is Decimal('0'), not Decimal('-0')
-            sign = 0
-        elif self._sign:
-            sign = 0
+            ans = self.copy_sign(Dec_0)
         else:
-            sign = 1
+            ans = self.copy_negate()
 
         if context is None:
             context = getcontext()
         if context._rounding_decision == ALWAYS_ROUND:
-            return Decimal((sign, self._int, self._exp))._fix(context)
-        return Decimal( (sign, self._int, self._exp))
+            return ans._fix(context)
+        return ans
 
     def __pos__(self, context=None):
         """Returns a copy, unless it is a sNaN.
@@ -945,19 +928,16 @@
             if ans:
                 return ans
 
-        sign = self._sign
         if not self:
             # + (-0) = 0
-            sign = 0
+            ans = self.copy_sign(Dec_0)
+        else:
+            ans = Decimal(self)
 
         if context is None:
             context = getcontext()
-
         if context._rounding_decision == ALWAYS_ROUND:
-            ans = self._fix(context)
-        else:
-            ans = Decimal(self)
-        ans._sign = sign
+            return ans._fix(context)
         return ans
 
     def __abs__(self, round=1, context=None):
@@ -1020,16 +1000,19 @@
             sign = min(self._sign, other._sign)
             if negativezero:
                 sign = 1
-            return Decimal( (sign, (0,), exp))
+            ans = Decimal( (sign, (0,), exp))
+            if shouldround:
+                ans = ans._fix(context)
+            return ans
         if not self:
             exp = max(exp, other._exp - context.prec-1)
-            ans = other._rescale(exp, watchexp=0, context=context)
+            ans = other._rescale(exp, context.rounding)
             if shouldround:
                 ans = ans._fix(context)
             return ans
         if not other:
             exp = max(exp, self._exp - context.prec-1)
-            ans = self._rescale(exp, watchexp=0, context=context)
+            ans = self._rescale(exp, context.rounding)
             if shouldround:
                 ans = ans._fix(context)
             return ans
@@ -1042,10 +1025,10 @@
         if op1.sign != op2.sign:
             # Equal and opposite
             if op1.int == op2.int:
-                if exp < context.Etiny():
-                    exp = context.Etiny()
-                    context._raise_error(Clamped)
-                return Decimal((negativezero, (0,), exp))
+                ans = Decimal((negativezero, (0,), exp))
+                if shouldround:
+                    ans = ans._fix(context)
+                return ans
             if op1.int < op2.int:
                 op1, op2 = op2, op1
                 # OK, now abs(op1) > abs(op2)
@@ -1076,7 +1059,7 @@
     __radd__ = __add__
 
     def __sub__(self, other, context=None):
-        """Return self + (-other)"""
+        """Return self - other"""
         other = _convert_other(other)
         if other is NotImplemented:
             return other
@@ -1086,41 +1069,28 @@
             if ans:
                 return ans
 
-        # -Decimal(0) = Decimal(0), which we don't want since
-        # (-0 - 0 = -0 + (-0) = -0, but -0 + 0 = 0.)
-        # so we change the sign directly to a copy
-        tmp = Decimal(other)
-        tmp._sign = 1-tmp._sign
-
-        return self.__add__(tmp, context=context)
+        # self - other is computed as self + other.copy_negate()
+        return self.__add__(other.copy_negate(), context=context)
 
     def __rsub__(self, other, context=None):
-        """Return other + (-self)"""
+        """Return other - self"""
         other = _convert_other(other)
         if other is NotImplemented:
             return other
 
-        tmp = Decimal(self)
-        tmp._sign = 1 - tmp._sign
-        return other.__add__(tmp, context=context)
+        return other.__sub__(self, context=context)
 
-    def _increment(self, round=1, context=None):
+    def _increment(self):
         """Special case of add, adding 1eExponent
 
         Since it is common, (rounding, for example) this adds
         (sign)*one E self._exp to the number more efficiently than add.
 
+        Assumes that self is nonspecial.
+
         For example:
         Decimal('5.624e10')._increment() == Decimal('5.625e10')
         """
-        if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-
-            # Must be infinite, and incrementing makes no difference
-            return Decimal(self)
-
         L = list(self._int)
         L[-1] += 1
         spot = len(L)-1
@@ -1131,13 +1101,7 @@
                 break
             L[spot-1] += 1
             spot -= 1
-        ans = Decimal((self._sign, L, self._exp))
-
-        if context is None:
-            context = getcontext()
-        if round and context._rounding_decision == ALWAYS_ROUND:
-            ans = ans._fix(context)
-        return ans
+        return Decimal((self._sign, L, self._exp))
 
     def __mul__(self, other, context=None):
         """Return self * other.
@@ -1205,25 +1169,9 @@
 
     def __truediv__(self, other, context=None):
         """Return self / other."""
-        return self._divide(other, context=context)
-
-    def _divide(self, other, divmod = 0, context=None):
-        """Return a / b, to context.prec precision.
-
-        divmod:
-        0 => true division
-        1 => (a //b, a%b)
-        2 => a //b
-        3 => a%b
-
-        Actually, if divmod is 2 or 3 a tuple is returned, but errors for
-        computing the other value are not raised.
-        """
         other = _convert_other(other)
         if other is NotImplemented:
-            if divmod in (0, 1):
-                return NotImplemented
-            return (NotImplemented, NotImplemented)
+            return NotImplemented
 
         if context is None:
             context = getcontext()
@@ -1233,144 +1181,84 @@
         if self._is_special or other._is_special:
             ans = self._check_nans(other, context)
             if ans:
-                if divmod:
-                    return (ans, ans)
                 return ans
 
             if self._isinfinity() and other._isinfinity():
-                if divmod:
-                    return (context._raise_error(InvalidOperation,
-                                            '(+-)INF // (+-)INF'),
-                            context._raise_error(InvalidOperation,
-                                            '(+-)INF % (+-)INF'))
                 return context._raise_error(InvalidOperation, '(+-)INF/(+-)INF')
 
             if self._isinfinity():
-                if divmod == 1:
-                    return (Infsign[sign],
-                            context._raise_error(InvalidOperation, 'INF % x'))
-                elif divmod == 2:
-                    return (Infsign[sign], NaN)
-                elif divmod == 3:
-                    return (Infsign[sign],
-                            context._raise_error(InvalidOperation, 'INF % x'))
                 return Infsign[sign]
 
             if other._isinfinity():
-                if divmod:
-                    return (Decimal((sign, (0,), 0)), Decimal(self))
                 context._raise_error(Clamped, 'Division by infinity')
                 return Decimal((sign, (0,), context.Etiny()))
 
         # Special cases for zeroes
-        if not self and not other:
-            if divmod:
-                return context._raise_error(DivisionUndefined, '0 / 0', 1)
-            return context._raise_error(DivisionUndefined, '0 / 0')
-
-        if not self:
-            if divmod:
-                otherside = Decimal(self)
-                otherside._exp = min(self._exp, other._exp)
-                return (Decimal((sign, (0,), 0)),  otherside)
-            exp = self._exp - other._exp
-            if exp < context.Etiny():
-                exp = context.Etiny()
-                context._raise_error(Clamped, '0e-x / y')
-            if exp > context.Emax:
-                exp = context.Emax
-                context._raise_error(Clamped, '0e+x / y')
-            return Decimal( (sign, (0,), exp) )
-
         if not other:
-            if divmod:
-                return context._raise_error(DivisionByZero, 'divmod(x,0)',
-                                           sign, 1)
+            if not self:
+                return context._raise_error(DivisionUndefined, '0 / 0')
             return context._raise_error(DivisionByZero, 'x / 0', sign)
 
-        # OK, so neither = 0, INF or NaN
-        shouldround = context._rounding_decision == ALWAYS_ROUND
+        if not self:
+            exp = self._exp - other._exp
+            coeff = 0
+        else:
+            # OK, so neither = 0, INF or NaN
+            shift = len(other._int) - len(self._int) + context.prec + 1
+            exp = self._exp - other._exp - shift
+            op1 = _WorkRep(self)
+            op2 = _WorkRep(other)
+            if shift >= 0:
+                coeff, remainder = divmod(op1.int * 10**shift, op2.int)
+            else:
+                coeff, remainder = divmod(op1.int, op2.int * 10**-shift)
+            if remainder:
+                # result is not exact; adjust to ensure correct rounding
+                if coeff % 5 == 0:
+                    coeff += 1
+            else:
+                # result is exact; get as close to ideal exponent as possible
+                ideal_exp = self._exp - other._exp
+                while exp < ideal_exp and coeff % 10 == 0:
+                    coeff //= 10
+                    exp += 1
 
-        # If we're dividing into ints, and self < other, stop.
-        # self.__abs__(0) does not round.
-        if divmod and (self.__abs__(0, context) < other.__abs__(0, context)):
+        ans = Decimal((sign, list(map(int, str(coeff))), exp))
+        return ans._fix(context)
 
-            if divmod == 1 or divmod == 3:
-                exp = min(self._exp, other._exp)
-                ans2 = self._rescale(exp, context=context, watchexp=0)
-                if shouldround:
-                    ans2 = ans2._fix(context)
-                return (Decimal( (sign, (0,), 0) ),
-                        ans2)
+    def _divide(self, other, context):
+        """Return (self // other, self % other), to context.prec precision.
 
-            elif divmod == 2:
-                # Don't round the mod part, if we don't need it.
-                return (Decimal( (sign, (0,), 0) ), Decimal(self))
+        Assumes that neither self nor other is a NaN, that self is not
+        infinite and that other is nonzero.
+        """
+        sign = self._sign ^ other._sign
+        if other._isinfinity():
+            ideal_exp = self._exp
+        else:
+            ideal_exp = min(self._exp, other._exp)
 
-        op1 = _WorkRep(self)
-        op2 = _WorkRep(other)
-        op1, op2, adjust = _adjust_coefficients(op1, op2)
-        res = _WorkRep( (sign, 0, (op1.exp - op2.exp)) )
-        if divmod and res.exp > context.prec + 1:
-            return context._raise_error(DivisionImpossible)
+        expdiff = self.adjusted() - other.adjusted()
+        if not self or other._isinfinity() or expdiff <= -2:
+            return (Decimal((sign, (0,), 0)),
+                    self._rescale(ideal_exp, context.rounding))
+        if expdiff <= context.prec:
+            op1 = _WorkRep(self)
+            op2 = _WorkRep(other)
+            if op1.exp >= op2.exp:
+                op1.int *= 10**(op1.exp - op2.exp)
+            else:
+                op2.int *= 10**(op2.exp - op1.exp)
+            q, r = divmod(op1.int, op2.int)
+            if q < 10**context.prec:
+                return (Decimal((sign, list(map(int, str(q))), 0)),
+                        Decimal((self._sign, list(map(int, str(r))),
+                                ideal_exp)))
 
-        prec_limit = 10 ** context.prec
-        while 1:
-            while op2.int <= op1.int:
-                res.int += 1
-                op1.int -= op2.int
-            if res.exp == 0 and divmod:
-                if res.int >= prec_limit and shouldround:
-                    return context._raise_error(DivisionImpossible)
-                otherside = Decimal(op1)
-                frozen = context._ignore_all_flags()
-
-                exp = min(self._exp, other._exp)
-                otherside = otherside._rescale(exp, context=context, watchexp=0)
-                context._regard_flags(*frozen)
-                if shouldround:
-                    otherside = otherside._fix(context)
-                return (Decimal(res), otherside)
-
-            if op1.int == 0 and adjust >= 0 and not divmod:
-                break
-            if res.int >= prec_limit and shouldround:
-                if divmod:
-                    return context._raise_error(DivisionImpossible)
-                shouldround=1
-                # Really, the answer is a bit higher, so adding a one to
-                # the end will make sure the rounding is right.
-                if op1.int != 0:
-                    res.int *= 10
-                    res.int += 1
-                    res.exp -= 1
-
-                break
-            res.int *= 10
-            res.exp -= 1
-            adjust += 1
-            op1.int *= 10
-            op1.exp -= 1
-
-            if res.exp == 0 and divmod and op2.int > op1.int:
-                # Solves an error in precision.  Same as a previous block.
-
-                if res.int >= prec_limit and shouldround:
-                    return context._raise_error(DivisionImpossible)
-                otherside = Decimal(op1)
-                frozen = context._ignore_all_flags()
-
-                exp = min(self._exp, other._exp)
-                otherside = otherside._rescale(exp, context=context)
-
-                context._regard_flags(*frozen)
-
-                return (Decimal(res), otherside)
-
-        ans = Decimal(res)
-        if shouldround:
-            ans = ans._fix(context)
-        return ans
+        # Here the quotient is too large to be representable
+        ans = context._raise_error(DivisionImpossible,
+                                   'quotient too large in //, % or divmod')
+        return ans, ans
 
     def __rtruediv__(self, other, context=None):
         """Swaps self/other and returns __truediv__."""
@@ -1381,9 +1269,40 @@
 
     def __divmod__(self, other, context=None):
         """
-        (self // other, self % other)
+        Return (self // other, self % other)
         """
-        return self._divide(other, 1, context)
+        other = _convert_other(other)
+        if other is NotImplemented:
+            return other
+
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return (ans, ans)
+
+        sign = self._sign ^ other._sign
+        if self._isinfinity():
+            if other._isinfinity():
+                ans = context._raise_error(InvalidOperation, 'divmod(INF, INF)')
+                return ans, ans
+            else:
+                return (Infsign[sign],
+                        context._raise_error(InvalidOperation, 'INF % x'))
+
+        if not other:
+            if not self:
+                ans = context._raise_error(DivisionUndefined, 'divmod(0, 0)')
+                return ans, ans
+            else:
+                return (context._raise_error(DivisionByZero, 'x // 0', sign),
+                        context._raise_error(InvalidOperation, 'x % 0'))
+
+        quotient, remainder = self._divide(other, context)
+        if context._rounding_decision == ALWAYS_ROUND:
+            remainder = remainder._fix(context)
+        return quotient, remainder
 
     def __rdivmod__(self, other, context=None):
         """Swaps self/other and returns __divmod__."""
@@ -1400,15 +1319,25 @@
         if other is NotImplemented:
             return other
 
-        if self._is_special or other._is_special:
-            ans = self._check_nans(other, context)
-            if ans:
-                return ans
+        if context is None:
+            context = getcontext()
 
-        if self and not other:
-            return context._raise_error(InvalidOperation, 'x % 0')
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
 
-        return self._divide(other, 3, context)[1]
+        if self._isinfinity():
+            return context._raise_error(InvalidOperation, 'INF % x')
+        elif not other:
+            if self:
+                return context._raise_error(InvalidOperation, 'x % 0')
+            else:
+                return context._raise_error(DivisionUndefined, '0 % 0')
+
+        remainder = self._divide(other, context)[1]
+        if context._rounding_decision == ALWAYS_ROUND:
+            remainder = remainder._fix(context)
+        return remainder
 
     def __rmod__(self, other, context=None):
         """Swaps self/other and returns __mod__."""
@@ -1421,85 +1350,104 @@
         """
         Remainder nearest to 0-  abs(remainder-near) <= other/2
         """
+        if context is None:
+            context = getcontext()
+
+        other = _convert_other(other, raiseit=True)
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        # self == +/-infinity -> InvalidOperation
+        if self._isinfinity():
+            return context._raise_error(InvalidOperation,
+                                        'remainder_near(infinity, x)')
+
+        # other == 0 -> either InvalidOperation or DivisionUndefined
+        if not other:
+            if self:
+                return context._raise_error(InvalidOperation,
+                                            'remainder_near(x, 0)')
+            else:
+                return context._raise_error(DivisionUndefined,
+                                            'remainder_near(0, 0)')
+
+        # other = +/-infinity -> remainder = self
+        if other._isinfinity():
+            ans = Decimal(self)
+            return ans._fix(context)
+
+        # self = 0 -> remainder = self, with ideal exponent
+        ideal_exponent = min(self._exp, other._exp)
+        if not self:
+            ans = Decimal((self._sign, (0,), ideal_exponent))
+            return ans._fix(context)
+
+        # catch most cases of large or small quotient
+        expdiff = self.adjusted() - other.adjusted()
+        if expdiff >= context.prec + 1:
+            # expdiff >= prec+1 => abs(self/other) > 10**prec
+            return context._raise_error(DivisionImpossible)
+        if expdiff <= -2:
+            # expdiff <= -2 => abs(self/other) < 0.1
+            ans = self._rescale(ideal_exponent, context.rounding)
+            return ans._fix(context)
+
+        # adjust both arguments to have the same exponent, then divide
+        op1 = _WorkRep(self)
+        op2 = _WorkRep(other)
+        if op1.exp >= op2.exp:
+            op1.int *= 10**(op1.exp - op2.exp)
+        else:
+            op2.int *= 10**(op2.exp - op1.exp)
+        q, r = divmod(op1.int, op2.int)
+        # remainder is r*10**ideal_exponent; other is +/-op2.int *
+        # 10**ideal_exponent.   Apply correction to ensure that
+        # abs(remainder) <= abs(other)/2
+        if 2*r + (q&1) > op2.int:
+            r -= op2.int
+            q += 1
+
+        if q >= 10**context.prec:
+            return context._raise_error(DivisionImpossible)
+
+        # result has same sign as self unless r is negative
+        sign = self._sign
+        if r < 0:
+            sign = 1-sign
+            r = -r
+
+        ans = Decimal((sign, list(map(int, str(r))), ideal_exponent))
+        return ans._fix(context)
+
+    def __floordiv__(self, other, context=None):
+        """self // other"""
         other = _convert_other(other)
         if other is NotImplemented:
             return other
 
-        if self._is_special or other._is_special:
-            ans = self._check_nans(other, context)
-            if ans:
-                return ans
-        if self and not other:
-            return context._raise_error(InvalidOperation, 'x % 0')
-
         if context is None:
             context = getcontext()
-        # If DivisionImpossible causes an error, do not leave Rounded/Inexact
-        # ignored in the calling function.
-        context = context._shallow_copy()
-        flags = context._ignore_flags(Rounded, Inexact)
-        # Keep DivisionImpossible flags
-        (side, r) = self.__divmod__(other, context=context)
 
-        if r._isnan():
-            context._regard_flags(*flags)
-            return r
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
 
-        context = context._shallow_copy()
-        rounding = context._set_rounding_decision(NEVER_ROUND)
-
-        if other._sign:
-            comparison = other.__truediv__(Decimal(-2), context=context)
-        else:
-            comparison = other.__truediv__(Decimal(2), context=context)
-
-        context._set_rounding_decision(rounding)
-        context._regard_flags(*flags)
-
-        s1, s2 = r._sign, comparison._sign
-        r._sign, comparison._sign = 0, 0
-
-        if r < comparison:
-            r._sign, comparison._sign = s1, s2
-            # Get flags now
-            self.__divmod__(other, context=context)
-            return r._fix(context)
-        r._sign, comparison._sign = s1, s2
-
-        rounding = context._set_rounding_decision(NEVER_ROUND)
-
-        (side, r) = self.__divmod__(other, context=context)
-        context._set_rounding_decision(rounding)
-        if r._isnan():
-            return r
-
-        decrease = not side._iseven()
-        rounding = context._set_rounding_decision(NEVER_ROUND)
-        side = side.__abs__(context=context)
-        context._set_rounding_decision(rounding)
-
-        s1, s2 = r._sign, comparison._sign
-        r._sign, comparison._sign = 0, 0
-        if r > comparison or decrease and r == comparison:
-            r._sign, comparison._sign = s1, s2
-            context.prec += 1
-            numbsquant = len(side.__add__(Decimal(1), context=context)._int)
-            if numbsquant >= context.prec:
-                context.prec -= 1
-                return context._raise_error(DivisionImpossible)[1]
-            context.prec -= 1
-            if self._sign == other._sign:
-                r = r.__sub__(other, context=context)
+        if self._isinfinity():
+            if other._isinfinity():
+                return context._raise_error(InvalidOperation, 'INF // INF')
             else:
-                r = r.__add__(other, context=context)
-        else:
-            r._sign, comparison._sign = s1, s2
+                return Infsign[self._sign ^ other._sign]
 
-        return r._fix(context)
+        if not other:
+            if self:
+                return context._raise_error(DivisionByZero, 'x // 0',
+                                            self._sign ^ other._sign)
+            else:
+                return context._raise_error(DivisionUndefined, '0 // 0')
 
-    def __floordiv__(self, other, context=None):
-        """self // other"""
-        return self._divide(other, 2, context)[0]
+        return self._divide(other, context)[0]
 
     def __rfloordiv__(self, other, context=None):
         """Swaps self/other and returns __floordiv__."""
@@ -1519,22 +1467,27 @@
                 context = getcontext()
                 return context._raise_error(InvalidContext)
             elif self._isinfinity():
-                raise OverflowError("Cannot convert infinity to long")
+                raise OverflowError("Cannot convert infinity to int")
+        s = (-1)**self._sign
         if self._exp >= 0:
-            s = ''.join(map(str, self._int)) + '0'*self._exp
+            return s*int(''.join(map(str, self._int)))*10**self._exp
         else:
-            s = ''.join(map(str, self._int))[:self._exp]
-        if s == '':
-            s = '0'
-        sign = '-'*self._sign
-        return int(sign + s)
+            return s*int(''.join(map(str, self._int))[:self._exp] or '0')
 
-    def __long__(self):
-        """Converts to a long.
+    def _fix_nan(self, context):
+        """Decapitate the payload of a NaN to fit the context"""
+        payload = self._int
 
-        Equivalent to long(int(self))
-        """
-        return int(self.__int__())
+        # maximum length of payload is precision if _clamp=0,
+        # precision-1 if _clamp=1.
+        max_payload_len = context.prec - context._clamp
+        if len(payload) > max_payload_len:
+            pos = len(payload)-max_payload_len
+            while pos < len(payload) and payload[pos] == 0:
+                pos += 1
+            payload = payload[pos:]
+            return Decimal((self._sign, payload, self._exp))
+        return Decimal(self)
 
     def _fix(self, context):
         """Round if it is necessary to keep self within prec precision.
@@ -1545,303 +1498,649 @@
         self - Decimal instance
         context - context used.
         """
-        if self._is_special:
-            return self
+
         if context is None:
             context = getcontext()
-        prec = context.prec
-        ans = self._fixexponents(context)
-        if len(ans._int) > prec:
-            ans = ans._round(prec, context=context)
-            ans = ans._fixexponents(context)
-        return ans
-
-    def _fixexponents(self, context):
-        """Fix the exponents and return a copy with the exponent in bounds.
-        Only call if known to not be a special value.
-        """
-        folddown = context._clamp
-        Emin = context.Emin
-        ans = self
-        ans_adjusted = ans.adjusted()
-        if ans_adjusted < Emin:
-            Etiny = context.Etiny()
-            if ans._exp < Etiny:
-                if not ans:
-                    ans = Decimal(self)
-                    ans._exp = Etiny
-                    context._raise_error(Clamped)
-                    return ans
-                ans = ans._rescale(Etiny, context=context)
-                # It isn't zero, and exp < Emin => subnormal
-                context._raise_error(Subnormal)
-                if context.flags[Inexact]:
-                    context._raise_error(Underflow)
-            else:
-                if ans:
-                    # Only raise subnormal if non-zero.
-                    context._raise_error(Subnormal)
-        else:
-            Etop = context.Etop()
-            if folddown and ans._exp > Etop:
-                context._raise_error(Clamped)
-                ans = ans._rescale(Etop, context=context)
-            else:
-                Emax = context.Emax
-                if ans_adjusted > Emax:
-                    if not ans:
-                        ans = Decimal(self)
-                        ans._exp = Emax
-                        context._raise_error(Clamped)
-                        return ans
-                    context._raise_error(Inexact)
-                    context._raise_error(Rounded)
-                    c = context._raise_error(Overflow, 'above Emax', ans._sign)
-                    return c
-        return ans
-
-    def _round(self, prec=None, rounding=None, context=None):
-        """Returns a rounded version of self.
-
-        You can specify the precision or rounding method.  Otherwise, the
-        context determines it.
-        """
 
         if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-
-            if self._isinfinity():
+            if self._isnan():
+                # decapitate payload if necessary
+                return self._fix_nan(context)
+            else:
+                # self is +/-Infinity; return unaltered
                 return Decimal(self)
 
-        if context is None:
-            context = getcontext()
-
-        if rounding is None:
-            rounding = context.rounding
-        if prec is None:
-            prec = context.prec
-
+        # if self is zero then exponent should be between Etiny and
+        # Emax if _clamp==0, and between Etiny and Etop if _clamp==1.
+        Etiny = context.Etiny()
+        Etop = context.Etop()
         if not self:
-            if prec <= 0:
-                dig = (0,)
-                exp = len(self._int) - prec + self._exp
+            exp_max = [context.Emax, Etop][context._clamp]
+            new_exp = min(max(self._exp, Etiny), exp_max)
+            if new_exp != self._exp:
+                context._raise_error(Clamped)
+                return Decimal((self._sign, (0,), new_exp))
             else:
-                dig = (0,) * prec
-                exp = len(self._int) + self._exp - prec
-            ans = Decimal((self._sign, dig, exp))
+                return Decimal(self)
+
+        # exp_min is the smallest allowable exponent of the result,
+        # equal to max(self.adjusted()-context.prec+1, Etiny)
+        exp_min = len(self._int) + self._exp - context.prec
+        if exp_min > Etop:
+            # overflow: exp_min > Etop iff self.adjusted() > Emax
+            context._raise_error(Inexact)
             context._raise_error(Rounded)
-            return ans
+            return context._raise_error(Overflow, 'above Emax', self._sign)
+        self_is_subnormal = exp_min < Etiny
+        if self_is_subnormal:
+            context._raise_error(Subnormal)
+            exp_min = Etiny
 
-        if prec == 0:
-            temp = Decimal(self)
-            temp._int = (0,)+temp._int
-            prec = 1
-        elif prec < 0:
-            exp = self._exp + len(self._int) - prec - 1
-            temp = Decimal( (self._sign, (0, 1), exp))
-            prec = 1
-        else:
-            temp = Decimal(self)
-
-        numdigits = len(temp._int)
-        if prec == numdigits:
-            return temp
-
-        # See if we need to extend precision
-        expdiff = prec - numdigits
-        if expdiff > 0:
-            tmp = list(temp._int)
-            tmp.extend([0] * expdiff)
-            ans =  Decimal( (temp._sign, tmp, temp._exp - expdiff))
-            return ans
-
-        # OK, but maybe all the lost digits are 0.
-        lostdigits = self._int[expdiff:]
-        if lostdigits == (0,) * len(lostdigits):
-            ans = Decimal( (temp._sign, temp._int[:prec], temp._exp - expdiff))
-            # Rounded, but not Inexact
+        # round if self has too many digits
+        if self._exp < exp_min:
             context._raise_error(Rounded)
+            ans = self._rescale(exp_min, context.rounding)
+            if ans != self:
+                context._raise_error(Inexact)
+                if self_is_subnormal:
+                    context._raise_error(Underflow)
+                    if not ans:
+                        # raise Clamped on underflow to 0
+                        context._raise_error(Clamped)
+                elif len(ans._int) == context.prec+1:
+                    # we get here only if rescaling rounds the
+                    # cofficient up to exactly 10**context.prec
+                    if ans._exp < Etop:
+                        ans = Decimal((ans._sign, ans._int[:-1], ans._exp+1))
+                    else:
+                        # Inexact and Rounded have already been raised
+                        ans = context._raise_error(Overflow, 'above Emax',
+                                                   self._sign)
             return ans
 
-        # Okay, let's round and lose data
+        # fold down if _clamp == 1 and self has too few digits
+        if context._clamp == 1 and self._exp > Etop:
+            context._raise_error(Clamped)
+            self_padded = self._int + (0,)*(self._exp - Etop)
+            return Decimal((self._sign, self_padded, Etop))
 
-        this_function = getattr(temp, self._pick_rounding_function[rounding])
-        # Now we've got the rounding function
-
-        if prec != context.prec:
-            context = context._shallow_copy()
-            context.prec = prec
-        ans = this_function(prec, expdiff, context)
-        context._raise_error(Rounded)
-        context._raise_error(Inexact, 'Changed in rounding')
-
-        return ans
+        # here self was representable to begin with; return unchanged
+        return Decimal(self)
 
     _pick_rounding_function = {}
 
-    def _round_down(self, prec, expdiff, context):
+    # for each of the rounding functions below:
+    #   self is a finite, nonzero Decimal
+    #   prec is an integer satisfying 0 <= prec < len(self._int)
+    # the rounded result will have exponent self._exp + len(self._int) - prec;
+
+    def _round_down(self, prec):
         """Also known as round-towards-0, truncate."""
-        return Decimal( (self._sign, self._int[:prec], self._exp - expdiff) )
+        newexp = self._exp + len(self._int) - prec
+        return Decimal((self._sign, self._int[:prec] or (0,), newexp))
 
-    def _round_half_up(self, prec, expdiff, context, tmp = None):
-        """Rounds 5 up (away from 0)"""
-
-        if tmp is None:
-            tmp = Decimal( (self._sign,self._int[:prec], self._exp - expdiff))
-        if self._int[prec] >= 5:
-            tmp = tmp._increment(round=0, context=context)
-            if len(tmp._int) > prec:
-                return Decimal( (tmp._sign, tmp._int[:-1], tmp._exp + 1))
-        return tmp
-
-    def _round_half_even(self, prec, expdiff, context):
-        """Round 5 to even, rest to nearest."""
-
-        tmp = Decimal( (self._sign, self._int[:prec], self._exp - expdiff))
-        half = (self._int[prec] == 5)
-        if half:
-            for digit in self._int[prec+1:]:
-                if digit != 0:
-                    half = 0
-                    break
-        if half:
-            if self._int[prec-1] & 1 == 0:
-                return tmp
-        return self._round_half_up(prec, expdiff, context, tmp)
-
-    def _round_half_down(self, prec, expdiff, context):
-        """Round 5 down"""
-
-        tmp = Decimal( (self._sign, self._int[:prec], self._exp - expdiff))
-        half = (self._int[prec] == 5)
-        if half:
-            for digit in self._int[prec+1:]:
-                if digit != 0:
-                    half = 0
-                    break
-        if half:
-            return tmp
-        return self._round_half_up(prec, expdiff, context, tmp)
-
-    def _round_up(self, prec, expdiff, context):
+    def _round_up(self, prec):
         """Rounds away from 0."""
-        tmp = Decimal( (self._sign, self._int[:prec], self._exp - expdiff) )
+        newexp = self._exp + len(self._int) - prec
+        tmp = Decimal((self._sign, self._int[:prec] or (0,), newexp))
         for digit in self._int[prec:]:
             if digit != 0:
-                tmp = tmp._increment(round=1, context=context)
-                if len(tmp._int) > prec:
-                    return Decimal( (tmp._sign, tmp._int[:-1], tmp._exp + 1))
-                else:
-                    return tmp
+                return tmp._increment()
         return tmp
 
-    def _round_ceiling(self, prec, expdiff, context):
+    def _round_half_up(self, prec):
+        """Rounds 5 up (away from 0)"""
+        if self._int[prec] >= 5:
+            return self._round_up(prec)
+        else:
+            return self._round_down(prec)
+
+    def _round_half_down(self, prec):
+        """Round 5 down"""
+        if self._int[prec] == 5:
+            for digit in self._int[prec+1:]:
+                if digit != 0:
+                    break
+            else:
+                return self._round_down(prec)
+        return self._round_half_up(prec)
+
+    def _round_half_even(self, prec):
+        """Round 5 to even, rest to nearest."""
+        if prec and self._int[prec-1] & 1:
+            return self._round_half_up(prec)
+        else:
+            return self._round_half_down(prec)
+
+    def _round_ceiling(self, prec):
         """Rounds up (not away from 0 if negative.)"""
         if self._sign:
-            return self._round_down(prec, expdiff, context)
+            return self._round_down(prec)
         else:
-            return self._round_up(prec, expdiff, context)
+            return self._round_up(prec)
 
-    def _round_floor(self, prec, expdiff, context):
+    def _round_floor(self, prec):
         """Rounds down (not towards 0 if negative)"""
         if not self._sign:
-            return self._round_down(prec, expdiff, context)
+            return self._round_down(prec)
         else:
-            return self._round_up(prec, expdiff, context)
+            return self._round_up(prec)
 
-    def __pow__(self, n, modulo = None, context=None):
-        """Return self ** n (mod modulo)
+    def _round_05up(self, prec):
+        """Round down unless digit prec-1 is 0 or 5."""
+        if prec == 0 or self._int[prec-1] in (0, 5):
+            return self._round_up(prec)
+        else:
+            return self._round_down(prec)
 
-        If modulo is None (default), don't take it mod modulo.
+    def fma(self, other, third, context=None):
+        """Fused multiply-add.
+
+        Returns self*other+third with no rounding of the intermediate
+        product self*other.
+
+        self and other are multiplied together, with no rounding of
+        the result.  The third operand is then added to the result,
+        and a single final rounding is performed.
         """
-        n = _convert_other(n)
-        if n is NotImplemented:
-            return n
+
+        other = _convert_other(other, raiseit=True)
+        third = _convert_other(third, raiseit=True)
 
         if context is None:
             context = getcontext()
 
-        if self._is_special or n._is_special or n.adjusted() > 8:
-            # Because the spot << doesn't work with really big exponents
-            if n._isinfinity() or n.adjusted() > 8:
-                return context._raise_error(InvalidOperation, 'x ** INF')
+        # do self*other in fresh context with no traps and no rounding
+        mul_context = Context(traps=[], flags=[],
+                              _rounding_decision=NEVER_ROUND)
+        product = self.__mul__(other, mul_context)
 
-            ans = self._check_nans(n, context)
-            if ans:
-                return ans
+        if mul_context.flags[InvalidOperation]:
+            # reraise in current context
+            return context._raise_error(InvalidOperation,
+                                        'invalid multiplication in fma',
+                                        1, product)
 
-        if not n._isinteger():
-            return context._raise_error(InvalidOperation, 'x ** (non-integer)')
+        ans = product.__add__(third, context)
+        return ans
 
-        if not self and not n:
-            return context._raise_error(InvalidOperation, '0 ** 0')
+    def _power_modulo(self, other, modulo, context=None):
+        """Three argument version of __pow__"""
 
-        if not n:
-            return Decimal(1)
+        # if can't convert other and modulo to Decimal, raise
+        # TypeError; there's no point returning NotImplemented (no
+        # equivalent of __rpow__ for three argument pow)
+        other = _convert_other(other, raiseit=True)
+        modulo = _convert_other(modulo, raiseit=True)
 
-        if self == Decimal(1):
-            return Decimal(1)
+        if context is None:
+            context = getcontext()
 
-        sign = self._sign and not n._iseven()
-        n = int(n)
+        # deal with NaNs: if there are any sNaNs then first one wins,
+        # (i.e. behaviour for NaNs is identical to that of fma)
+        self_is_nan = self._isnan()
+        other_is_nan = other._isnan()
+        modulo_is_nan = modulo._isnan()
+        if self_is_nan or other_is_nan or modulo_is_nan:
+            if self_is_nan == 2:
+                return context._raise_error(InvalidOperation, 'sNaN',
+                                        1, self)
+            if other_is_nan == 2:
+                return context._raise_error(InvalidOperation, 'sNaN',
+                                        1, other)
+            if modulo_is_nan == 2:
+                return context._raise_error(InvalidOperation, 'sNaN',
+                                        1, modulo)
+            if self_is_nan:
+                return self._fix_nan(context)
+            if other_is_nan:
+                return other._fix_nan(context)
+            return modulo._fix_nan(context)
 
+        # check inputs: we apply same restrictions as Python's pow()
+        if not (self._isinteger() and
+                other._isinteger() and
+                modulo._isinteger()):
+            return context._raise_error(InvalidOperation,
+                                        'pow() 3rd argument not allowed '
+                                        'unless all arguments are integers')
+        if other < 0:
+            return context._raise_error(InvalidOperation,
+                                        'pow() 2nd argument cannot be '
+                                        'negative when 3rd argument specified')
+        if not modulo:
+            return context._raise_error(InvalidOperation,
+                                        'pow() 3rd argument cannot be 0')
+
+        # additional restriction for decimal: the modulus must be less
+        # than 10**prec in absolute value
+        if modulo.adjusted() >= context.prec:
+            return context._raise_error(InvalidOperation,
+                                        'insufficient precision: pow() 3rd '
+                                        'argument must not have more than '
+                                        'precision digits')
+
+        # define 0**0 == NaN, for consistency with two-argument pow
+        # (even though it hurts!)
+        if not other and not self:
+            return context._raise_error(InvalidOperation,
+                                        'at least one of pow() 1st argument '
+                                        'and 2nd argument must be nonzero ;'
+                                        '0**0 is not defined')
+
+        # compute sign of result
+        if other._iseven():
+            sign = 0
+        else:
+            sign = self._sign
+
+        # convert modulo to a Python integer, and self and other to
+        # Decimal integers (i.e. force their exponents to be >= 0)
+        modulo = abs(int(modulo))
+        base = _WorkRep(self.to_integral_value())
+        exponent = _WorkRep(other.to_integral_value())
+
+        # compute result using integer pow()
+        base = (base.int % modulo * pow(10, base.exp, modulo)) % modulo
+        for i in range(exponent.exp):
+            base = pow(base, 10, modulo)
+        base = pow(base, exponent.int, modulo)
+
+        return Decimal((sign, list(map(int, str(base))), 0))
+
+    def _power_exact(self, other, p):
+        """Attempt to compute self**other exactly.
+
+        Given Decimals self and other and an integer p, attempt to
+        compute an exact result for the power self**other, with p
+        digits of precision.  Return None if self**other is not
+        exactly representable in p digits.
+
+        Assumes that elimination of special cases has already been
+        performed: self and other must both be nonspecial; self must
+        be positive and not numerically equal to 1; other must be
+        nonzero.  For efficiency, other._exp should not be too large,
+        so that 10**abs(other._exp) is a feasible calculation."""
+
+        # In the comments below, we write x for the value of self and
+        # y for the value of other.  Write x = xc*10**xe and y =
+        # yc*10**ye.
+
+        # The main purpose of this method is to identify the *failure*
+        # of x**y to be exactly representable with as little effort as
+        # possible.  So we look for cheap and easy tests that
+        # eliminate the possibility of x**y being exact.  Only if all
+        # these tests are passed do we go on to actually compute x**y.
+
+        # Here's the main idea.  First normalize both x and y.  We
+        # express y as a rational m/n, with m and n relatively prime
+        # and n>0.  Then for x**y to be exactly representable (at
+        # *any* precision), xc must be the nth power of a positive
+        # integer and xe must be divisible by n.  If m is negative
+        # then additionally xc must be a power of either 2 or 5, hence
+        # a power of 2**n or 5**n.
+        #
+        # There's a limit to how small |y| can be: if y=m/n as above
+        # then:
+        #
+        #  (1) if xc != 1 then for the result to be representable we
+        #      need xc**(1/n) >= 2, and hence also xc**|y| >= 2.  So
+        #      if |y| <= 1/nbits(xc) then xc < 2**nbits(xc) <=
+        #      2**(1/|y|), hence xc**|y| < 2 and the result is not
+        #      representable.
+        #
+        #  (2) if xe != 0, |xe|*(1/n) >= 1, so |xe|*|y| >= 1.  Hence if
+        #      |y| < 1/|xe| then the result is not representable.
+        #
+        # Note that since x is not equal to 1, at least one of (1) and
+        # (2) must apply.  Now |y| < 1/nbits(xc) iff |yc|*nbits(xc) <
+        # 10**-ye iff len(str(|yc|*nbits(xc)) <= -ye.
+        #
+        # There's also a limit to how large y can be, at least if it's
+        # positive: the normalized result will have coefficient xc**y,
+        # so if it's representable then xc**y < 10**p, and y <
+        # p/log10(xc).  Hence if y*log10(xc) >= p then the result is
+        # not exactly representable.
+
+        # if len(str(abs(yc*xe)) <= -ye then abs(yc*xe) < 10**-ye,
+        # so |y| < 1/xe and the result is not representable.
+        # Similarly, len(str(abs(yc)*xc_bits)) <= -ye implies |y|
+        # < 1/nbits(xc).
+
+        x = _WorkRep(self)
+        xc, xe = x.int, x.exp
+        while xc % 10 == 0:
+            xc //= 10
+            xe += 1
+
+        y = _WorkRep(other)
+        yc, ye = y.int, y.exp
+        while yc % 10 == 0:
+            yc //= 10
+            ye += 1
+
+        # case where xc == 1: result is 10**(xe*y), with xe*y
+        # required to be an integer
+        if xc == 1:
+            if ye >= 0:
+                exponent = xe*yc*10**ye
+            else:
+                exponent, remainder = divmod(xe*yc, 10**-ye)
+                if remainder:
+                    return None
+            if y.sign == 1:
+                exponent = -exponent
+            # if other is a nonnegative integer, use ideal exponent
+            if other._isinteger() and other._sign == 0:
+                ideal_exponent = self._exp*int(other)
+                zeros = min(exponent-ideal_exponent, p-1)
+            else:
+                zeros = 0
+            return Decimal((0, (1,) + (0,)*zeros, exponent-zeros))
+
+        # case where y is negative: xc must be either a power
+        # of 2 or a power of 5.
+        if y.sign == 1:
+            last_digit = xc % 10
+            if last_digit in (2,4,6,8):
+                # quick test for power of 2
+                if xc & -xc != xc:
+                    return None
+                # now xc is a power of 2; e is its exponent
+                e = _nbits(xc)-1
+                # find e*y and xe*y; both must be integers
+                if ye >= 0:
+                    y_as_int = yc*10**ye
+                    e = e*y_as_int
+                    xe = xe*y_as_int
+                else:
+                    ten_pow = 10**-ye
+                    e, remainder = divmod(e*yc, ten_pow)
+                    if remainder:
+                        return None
+                    xe, remainder = divmod(xe*yc, ten_pow)
+                    if remainder:
+                        return None
+
+                if e*65 >= p*93: # 93/65 > log(10)/log(5)
+                    return None
+                xc = 5**e
+
+            elif last_digit == 5:
+                # e >= log_5(xc) if xc is a power of 5; we have
+                # equality all the way up to xc=5**2658
+                e = _nbits(xc)*28//65
+                xc, remainder = divmod(5**e, xc)
+                if remainder:
+                    return None
+                while xc % 5 == 0:
+                    xc //= 5
+                    e -= 1
+                if ye >= 0:
+                    y_as_integer = yc*10**ye
+                    e = e*y_as_integer
+                    xe = xe*y_as_integer
+                else:
+                    ten_pow = 10**-ye
+                    e, remainder = divmod(e*yc, ten_pow)
+                    if remainder:
+                        return None
+                    xe, remainder = divmod(xe*yc, ten_pow)
+                    if remainder:
+                        return None
+                if e*3 >= p*10: # 10/3 > log(10)/log(2)
+                    return None
+                xc = 2**e
+            else:
+                return None
+
+            if xc >= 10**p:
+                return None
+            xe = -e-xe
+            return Decimal((0, list(map(int, str(xc))), xe))
+
+        # now y is positive; find m and n such that y = m/n
+        if ye >= 0:
+            m, n = yc*10**ye, 1
+        else:
+            if xe != 0 and len(str(abs(yc*xe))) <= -ye:
+                return None
+            xc_bits = _nbits(xc)
+            if xc != 1 and len(str(abs(yc)*xc_bits)) <= -ye:
+                return None
+            m, n = yc, 10**(-ye)
+            while m % 2 == n % 2 == 0:
+                m //= 2
+                n //= 2
+            while m % 5 == n % 5 == 0:
+                m //= 5
+                n //= 5
+
+        # compute nth root of xc*10**xe
+        if n > 1:
+            # if 1 < xc < 2**n then xc isn't an nth power
+            if xc != 1 and xc_bits <= n:
+                return None
+
+            xe, rem = divmod(xe, n)
+            if rem != 0:
+                return None
+
+            # compute nth root of xc using Newton's method
+            a = 1 << -(-_nbits(xc)//n) # initial estimate
+            while True:
+                q, r = divmod(xc, a**(n-1))
+                if a <= q:
+                    break
+                else:
+                    a = (a*(n-1) + q)//n
+            if not (a == q and r == 0):
+                return None
+            xc = a
+
+        # now xc*10**xe is the nth root of the original xc*10**xe
+        # compute mth power of xc*10**xe
+
+        # if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
+        # 10**p and the result is not representable.
+        if xc > 1 and m > p*100//_log10_lb(xc):
+            return None
+        xc = xc**m
+        xe *= m
+        if xc > 10**p:
+            return None
+
+        # by this point the result *is* exactly representable
+        # adjust the exponent to get as close as possible to the ideal
+        # exponent, if necessary
+        str_xc = str(xc)
+        if other._isinteger() and other._sign == 0:
+            ideal_exponent = self._exp*int(other)
+            zeros = min(xe-ideal_exponent, p-len(str_xc))
+        else:
+            zeros = 0
+        return Decimal((0, list(map(int, str_xc))+[0,]*zeros, xe-zeros))
+
+    def __pow__(self, other, modulo=None, context=None):
+        """Return self ** other [ % modulo].
+
+        With two arguments, compute self**other.
+
+        With three arguments, compute (self**other) % modulo.  For the
+        three argument form, the following restrictions on the
+        arguments hold:
+
+         - all three arguments must be integral
+         - other must be nonnegative
+         - either self or other (or both) must be nonzero
+         - modulo must be nonzero and must have at most p digits,
+           where p is the context precision.
+
+        If any of these restrictions is violated the InvalidOperation
+        flag is raised.
+
+        The result of pow(self, other, modulo) is identical to the
+        result that would be obtained by computing (self**other) %
+        modulo with unbounded precision, but is computed more
+        efficiently.  It is always exact.
+        """
+
+        if modulo is not None:
+            return self._power_modulo(other, modulo, context)
+
+        other = _convert_other(other)
+        if other is NotImplemented:
+            return other
+
+        if context is None:
+            context = getcontext()
+
+        # either argument is a NaN => result is NaN
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
+        if not other:
+            if not self:
+                return context._raise_error(InvalidOperation, '0 ** 0')
+            else:
+                return Dec_p1
+
+        # result has sign 1 iff self._sign is 1 and other is an odd integer
+        result_sign = 0
+        if self._sign == 1:
+            if other._isinteger():
+                if not other._iseven():
+                    result_sign = 1
+            else:
+                # -ve**noninteger = NaN
+                # (-0)**noninteger = 0**noninteger
+                if self:
+                    return context._raise_error(InvalidOperation,
+                        'x ** y with x negative and y not an integer')
+            # negate self, without doing any unwanted rounding
+            self = Decimal((0, self._int, self._exp))
+
+        # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
+        if not self:
+            if other._sign == 0:
+                return Decimal((result_sign, (0,), 0))
+            else:
+                return Infsign[result_sign]
+
+        # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
         if self._isinfinity():
-            if modulo:
-                return context._raise_error(InvalidOperation, 'INF % x')
-            if n > 0:
-                return Infsign[sign]
-            return Decimal( (sign, (0,), 0) )
+            if other._sign == 0:
+                return Infsign[result_sign]
+            else:
+                return Decimal((result_sign, (0,), 0))
 
-        # With ludicrously large exponent, just raise an overflow
-        # and return inf.
-        if not modulo and n > 0 and \
-           (self._exp + len(self._int) - 1) * n > context.Emax and self:
+        # 1**other = 1, but the choice of exponent and the flags
+        # depend on the exponent of self, and on whether other is a
+        # positive integer, a negative integer, or neither
+        if self == Dec_p1:
+            if other._isinteger():
+                # exp = max(self._exp*max(int(other), 0),
+                # 1-context.prec) but evaluating int(other) directly
+                # is dangerous until we know other is small (other
+                # could be 1e999999999)
+                if other._sign == 1:
+                    multiplier = 0
+                elif other > context.prec:
+                    multiplier = context.prec
+                else:
+                    multiplier = int(other)
 
-            tmp = Decimal('inf')
-            tmp._sign = sign
-            context._raise_error(Rounded)
+                exp = self._exp * multiplier
+                if exp < 1-context.prec:
+                    exp = 1-context.prec
+                    context._raise_error(Rounded)
+            else:
+                context._raise_error(Inexact)
+                context._raise_error(Rounded)
+                exp = 1-context.prec
+
+            return Decimal((result_sign, (1,)+(0,)*-exp, exp))
+
+        # compute adjusted exponent of self
+        self_adj = self.adjusted()
+
+        # self ** infinity is infinity if self > 1, 0 if self < 1
+        # self ** -infinity is infinity if self < 1, 0 if self > 1
+        if other._isinfinity():
+            if (other._sign == 0) == (self_adj < 0):
+                return Decimal((result_sign, (0,), 0))
+            else:
+                return Infsign[result_sign]
+
+        # from here on, the result always goes through the call
+        # to _fix at the end of this function.
+        ans = None
+
+        # crude test to catch cases of extreme overflow/underflow.  If
+        # log10(self)*other >= 10**bound and bound >= len(str(Emax))
+        # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
+        # self**other >= 10**(Emax+1), so overflow occurs.  The test
+        # for underflow is similar.
+        bound = self._log10_exp_bound() + other.adjusted()
+        if (self_adj >= 0) == (other._sign == 0):
+            # self > 1 and other +ve, or self < 1 and other -ve
+            # possibility of overflow
+            if bound >= len(str(context.Emax)):
+                ans = Decimal((result_sign, (1,), context.Emax+1))
+        else:
+            # self > 1 and other -ve, or self < 1 and other +ve
+            # possibility of underflow to 0
+            Etiny = context.Etiny()
+            if bound >= len(str(-Etiny)):
+                ans = Decimal((result_sign, (1,), Etiny-1))
+
+        # try for an exact result with precision +1
+        if ans is None:
+            ans = self._power_exact(other, context.prec + 1)
+            if ans is not None and result_sign == 1:
+                ans = Decimal((1, ans._int, ans._exp))
+
+        # usual case: inexact result, x**y computed directly as exp(y*log(x))
+        if ans is None:
+            p = context.prec
+            x = _WorkRep(self)
+            xc, xe = x.int, x.exp
+            y = _WorkRep(other)
+            yc, ye = y.int, y.exp
+            if y.sign == 1:
+                yc = -yc
+
+            # compute correctly rounded result:  start with precision +3,
+            # then increase precision until result is unambiguously roundable
+            extra = 3
+            while True:
+                coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
+                if coeff % (5*10**(len(str(coeff))-p-1)):
+                    break
+                extra += 3
+
+            ans = Decimal((result_sign, list(map(int, str(coeff))), exp))
+
+        # the specification says that for non-integer other we need to
+        # raise Inexact, even when the result is actually exact.  In
+        # the same way, we need to raise Underflow here if the result
+        # is subnormal.  (The call to _fix will take care of raising
+        # Rounded and Subnormal, as usual.)
+        if not other._isinteger():
             context._raise_error(Inexact)
-            context._raise_error(Overflow, 'Big power', sign)
-            return tmp
+            # pad with zeros up to length context.prec+1 if necessary
+            if len(ans._int) <= context.prec:
+                expdiff = context.prec+1 - len(ans._int)
+                ans = Decimal((ans._sign, ans._int+(0,)*expdiff, ans._exp-expdiff))
+            if ans.adjusted() < context.Emin:
+                context._raise_error(Underflow)
 
-        elength = len(str(abs(n)))
-        firstprec = context.prec
-
-        if not modulo and firstprec + elength + 1 > DefaultContext.Emax:
-            return context._raise_error(Overflow, 'Too much precision.', sign)
-
-        mul = Decimal(self)
-        val = Decimal(1)
-        context = context._shallow_copy()
-        context.prec = firstprec + elength + 1
-        if n < 0:
-            # n is a long now, not Decimal instance
-            n = -n
-            mul = Decimal(1).__truediv__(mul, context=context)
-
-        spot = 1
-        while spot <= n:
-            spot <<= 1
-
-        spot >>= 1
-        # spot is the highest power of 2 less than n
-        while spot:
-            val = val.__mul__(val, context=context)
-            if val._isinfinity():
-                val = Infsign[sign]
-                break
-            if spot & n:
-                val = val.__mul__(mul, context=context)
-            if modulo is not None:
-                val = val.__mod__(modulo, context=context)
-            spot >>= 1
-        context.prec = firstprec
-
-        if context._rounding_decision == ALWAYS_ROUND:
-            return val._fix(context)
-        return val
+        # unlike exp, ln and log10, the power function respects the
+        # rounding mode; no need to use ROUND_HALF_EVEN here
+        ans = ans._fix(context)
+        return ans
 
     def __rpow__(self, other, context=None):
         """Swaps self/other and returns __pow__."""
@@ -1853,6 +2152,9 @@
     def normalize(self, context=None):
         """Normalize- strip trailing 0s, change anything equal to 0 to 0e0"""
 
+        if context is None:
+            context = getcontext()
+
         if self._is_special:
             ans = self._check_nans(context=context)
             if ans:
@@ -1864,19 +2166,26 @@
 
         if not dup:
             return Decimal( (dup._sign, (0,), 0) )
+        exp_max = [context.Emax, context.Etop()][context._clamp]
         end = len(dup._int)
         exp = dup._exp
-        while dup._int[end-1] == 0:
+        while dup._int[end-1] == 0 and exp < exp_max:
             exp += 1
             end -= 1
         return Decimal( (dup._sign, dup._int[:end], exp) )
 
-
-    def quantize(self, exp, rounding=None, context=None, watchexp=1):
+    def quantize(self, exp, rounding=None, context=None, watchexp=True):
         """Quantize self so its exponent is the same as that of exp.
 
         Similar to self._rescale(exp._exp) but with error checking.
         """
+        exp = _convert_other(exp, raiseit=True)
+
+        if context is None:
+            context = getcontext()
+        if rounding is None:
+            rounding = context.rounding
+
         if self._is_special or exp._is_special:
             ans = self._check_nans(exp, context)
             if ans:
@@ -1884,12 +2193,56 @@
 
             if exp._isinfinity() or self._isinfinity():
                 if exp._isinfinity() and self._isinfinity():
-                    return self  # if both are inf, it is OK
-                if context is None:
-                    context = getcontext()
+                    return Decimal(self)  # if both are inf, it is OK
                 return context._raise_error(InvalidOperation,
                                         'quantize with one INF')
-        return self._rescale(exp._exp, rounding, context, watchexp)
+
+        # if we're not watching exponents, do a simple rescale
+        if not watchexp:
+            ans = self._rescale(exp._exp, rounding)
+            # raise Inexact and Rounded where appropriate
+            if ans._exp > self._exp:
+                context._raise_error(Rounded)
+                if ans != self:
+                    context._raise_error(Inexact)
+            return ans
+
+        # exp._exp should be between Etiny and Emax
+        if not (context.Etiny() <= exp._exp <= context.Emax):
+            return context._raise_error(InvalidOperation,
+                   'target exponent out of bounds in quantize')
+
+        if not self:
+            ans = Decimal((self._sign, (0,), exp._exp))
+            return ans._fix(context)
+
+        self_adjusted = self.adjusted()
+        if self_adjusted > context.Emax:
+            return context._raise_error(InvalidOperation,
+                                        'exponent of quantize result too large for current context')
+        if self_adjusted - exp._exp + 1 > context.prec:
+            return context._raise_error(InvalidOperation,
+                                        'quantize result has too many digits for current context')
+
+        ans = self._rescale(exp._exp, rounding)
+        if ans.adjusted() > context.Emax:
+            return context._raise_error(InvalidOperation,
+                                        'exponent of quantize result too large for current context')
+        if len(ans._int) > context.prec:
+            return context._raise_error(InvalidOperation,
+                                        'quantize result has too many digits for current context')
+
+        # raise appropriate flags
+        if ans._exp > self._exp:
+            context._raise_error(Rounded)
+            if ans != self:
+                context._raise_error(Inexact)
+        if ans and ans.adjusted() < context.Emin:
+            context._raise_error(Subnormal)
+
+        # call to fix takes care of any necessary folddown
+        ans = ans._fix(context)
+        return ans
 
     def same_quantum(self, other):
         """Test whether self and other have the same exponent.
@@ -1903,83 +2256,85 @@
                 return self._isinfinity() and other._isinfinity() and True
         return self._exp == other._exp
 
-    def _rescale(self, exp, rounding=None, context=None, watchexp=1):
-        """Rescales so that the exponent is exp.
+    def _rescale(self, exp, rounding):
+        """Rescale self so that the exponent is exp, either by padding with zeros
+        or by truncating digits, using the given rounding mode.
+
+        Specials are returned without change.  This operation is
+        quiet: it raises no flags, and uses no information from the
+        context.
 
         exp = exp to scale to (an integer)
-        rounding = rounding version
-        watchexp: if set (default) an error is returned if exp is greater
-        than Emax or less than Etiny.
+        rounding = rounding mode
         """
-        if context is None:
-            context = getcontext()
-
         if self._is_special:
-            if self._isinfinity():
-                return context._raise_error(InvalidOperation, 'rescale with an INF')
-
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-
-        if watchexp and (context.Emax  < exp or context.Etiny() > exp):
-            return context._raise_error(InvalidOperation, 'rescale(a, INF)')
-
+            return Decimal(self)
         if not self:
-            ans = Decimal(self)
-            ans._int = (0,)
-            ans._exp = exp
-            return ans
+            return Decimal((self._sign, (0,), exp))
 
-        diff = self._exp - exp
-        digits = len(self._int) + diff
+        if self._exp >= exp:
+            # pad answer with zeros if necessary
+            return Decimal((self._sign, self._int + (0,)*(self._exp - exp), exp))
 
-        if watchexp and digits > context.prec:
-            return context._raise_error(InvalidOperation, 'Rescale > prec')
-
-        tmp = Decimal(self)
-        tmp._int = (0,) + tmp._int
-        digits += 1
-
+        # too many digits; round and lose data.  If self.adjusted() <
+        # exp-1, replace self by 10**(exp-1) before rounding
+        digits = len(self._int) + self._exp - exp
         if digits < 0:
-            tmp._exp = -digits + tmp._exp
-            tmp._int = (0,1)
-            digits = 1
-        tmp = tmp._round(digits, rounding, context=context)
+            self = Decimal((self._sign, (1,), exp-1))
+            digits = 0
+        this_function = getattr(self, self._pick_rounding_function[rounding])
+        return this_function(digits)
 
-        if tmp._int[0] == 0 and len(tmp._int) > 1:
-            tmp._int = tmp._int[1:]
-        tmp._exp = exp
+    def to_integral_exact(self, rounding=None, context=None):
+        """Rounds to a nearby integer.
 
-        tmp_adjusted = tmp.adjusted()
-        if tmp and tmp_adjusted < context.Emin:
-            context._raise_error(Subnormal)
-        elif tmp and tmp_adjusted > context.Emax:
-            return context._raise_error(InvalidOperation, 'rescale(a, INF)')
-        return tmp
+        If no rounding mode is specified, take the rounding mode from
+        the context.  This method raises the Rounded and Inexact flags
+        when appropriate.
 
-    def to_integral(self, rounding=None, context=None):
-        """Rounds to the nearest integer, without raising inexact, rounded."""
+        See also: to_integral_value, which does exactly the same as
+        this method except that it doesn't raise Inexact or Rounded.
+        """
         if self._is_special:
             ans = self._check_nans(context=context)
             if ans:
                 return ans
-            return self
+            return Decimal(self)
         if self._exp >= 0:
-            return self
+            return Decimal(self)
+        if not self:
+            return Decimal((self._sign, (0,), 0))
         if context is None:
             context = getcontext()
-        flags = context._ignore_flags(Rounded, Inexact)
-        ans = self._rescale(0, rounding, context=context)
-        context._regard_flags(flags)
+        if rounding is None:
+            rounding = context.rounding
+        context._raise_error(Rounded)
+        ans = self._rescale(0, rounding)
+        if ans != self:
+            context._raise_error(Inexact)
         return ans
 
-    def sqrt(self, context=None):
-        """Return the square root of self.
+    def to_integral_value(self, rounding=None, context=None):
+        """Rounds to the nearest integer, without raising inexact, rounded."""
+        if context is None:
+            context = getcontext()
+        if rounding is None:
+            rounding = context.rounding
+        if self._is_special:
+            ans = self._check_nans(context=context)
+            if ans:
+                return ans
+            return Decimal(self)
+        if self._exp >= 0:
+            return Decimal(self)
+        else:
+            return self._rescale(0, rounding)
 
-        Uses a converging algorithm (Xn+1 = 0.5*(Xn + self / Xn))
-        Should quadratically approach the right answer.
-        """
+    # the method name changed, but we provide also the old one, for compatibility
+    to_integral = to_integral_value
+
+    def sqrt(self, context=None):
+        """Return the square root of self."""
         if self._is_special:
             ans = self._check_nans(context=context)
             if ans:
@@ -1989,16 +2344,9 @@
                 return Decimal(self)
 
         if not self:
-            # exponent = self._exp / 2, using round_down.
-            # if self._exp < 0:
-            #    exp = (self._exp+1) // 2
-            # else:
-            exp = (self._exp) // 2
-            if self._sign == 1:
-                # sqrt(-0) = -0
-                return Decimal( (1, (0,), exp))
-            else:
-                return Decimal( (0, (0,), exp))
+            # exponent = self._exp // 2.  sqrt(-0) = -0
+            ans = Decimal((self._sign, (0,), self._exp // 2))
+            return ans._fix(context)
 
         if context is None:
             context = getcontext()
@@ -2006,104 +2354,94 @@
         if self._sign == 1:
             return context._raise_error(InvalidOperation, 'sqrt(-x), x > 0')
 
-        tmp = Decimal(self)
+        # At this point self represents a positive number.  Let p be
+        # the desired precision and express self in the form c*100**e
+        # with c a positive real number and e an integer, c and e
+        # being chosen so that 100**(p-1) <= c < 100**p.  Then the
+        # (exact) square root of self is sqrt(c)*10**e, and 10**(p-1)
+        # <= sqrt(c) < 10**p, so the closest representable Decimal at
+        # precision p is n*10**e where n = round_half_even(sqrt(c)),
+        # the closest integer to sqrt(c) with the even integer chosen
+        # in the case of a tie.
+        #
+        # To ensure correct rounding in all cases, we use the
+        # following trick: we compute the square root to an extra
+        # place (precision p+1 instead of precision p), rounding down.
+        # Then, if the result is inexact and its last digit is 0 or 5,
+        # we increase the last digit to 1 or 6 respectively; if it's
+        # exact we leave the last digit alone.  Now the final round to
+        # p places (or fewer in the case of underflow) will round
+        # correctly and raise the appropriate flags.
 
-        expadd = tmp._exp // 2
-        if tmp._exp & 1:
-            tmp._int += (0,)
-            tmp._exp = 0
+        # use an extra digit of precision
+        prec = context.prec+1
+
+        # write argument in the form c*100**e where e = self._exp//2
+        # is the 'ideal' exponent, to be used if the square root is
+        # exactly representable.  l is the number of 'digits' of c in
+        # base 100, so that 100**(l-1) <= c < 100**l.
+        op = _WorkRep(self)
+        e = op.exp >> 1
+        if op.exp & 1:
+            c = op.int * 10
+            l = (len(self._int) >> 1) + 1
         else:
-            tmp._exp = 0
+            c = op.int
+            l = len(self._int)+1 >> 1
 
-        context = context._shallow_copy()
-        flags = context._ignore_all_flags()
-        firstprec = context.prec
-        context.prec = 3
-        if tmp.adjusted() & 1 == 0:
-            ans = Decimal( (0, (8,1,9), tmp.adjusted()  - 2) )
-            ans = ans.__add__(tmp.__mul__(Decimal((0, (2,5,9), -2)),
-                                          context=context), context=context)
-            ans._exp -= 1 + tmp.adjusted() // 2
+        # rescale so that c has exactly prec base 100 'digits'
+        shift = prec-l
+        if shift >= 0:
+            c *= 100**shift
+            exact = True
         else:
-            ans = Decimal( (0, (2,5,9), tmp._exp + len(tmp._int)- 3) )
-            ans = ans.__add__(tmp.__mul__(Decimal((0, (8,1,9), -3)),
-                                          context=context), context=context)
-            ans._exp -= 1 + tmp.adjusted()  // 2
+            c, remainder = divmod(c, 100**-shift)
+            exact = not remainder
+        e -= shift
 
-        # ans is now a linear approximation.
-        Emax, Emin = context.Emax, context.Emin
-        context.Emax, context.Emin = DefaultContext.Emax, DefaultContext.Emin
-
-        half = Decimal('0.5')
-
-        maxp = firstprec + 2
-        rounding = context._set_rounding(ROUND_HALF_EVEN)
-        while 1:
-            context.prec = min(2*context.prec - 2, maxp)
-            ans = half.__mul__(ans.__add__(tmp.__truediv__(ans, context=context),
-                                           context=context), context=context)
-            if context.prec == maxp:
+        # find n = floor(sqrt(c)) using Newton's method
+        n = 10**prec
+        while True:
+            q = c//n
+            if n <= q:
                 break
+            else:
+                n = n + q >> 1
+        exact = exact and n*n == c
 
-        # Round to the answer's precision-- the only error can be 1 ulp.
-        context.prec = firstprec
-        prevexp = ans.adjusted()
-        ans = ans._round(context=context)
-
-        # Now, check if the other last digits are better.
-        context.prec = firstprec + 1
-        # In case we rounded up another digit and we should actually go lower.
-        if prevexp != ans.adjusted():
-            ans._int += (0,)
-            ans._exp -= 1
-
-
-        lower = ans.__sub__(Decimal((0, (5,), ans._exp-1)), context=context)
-        context._set_rounding(ROUND_UP)
-        if lower.__mul__(lower, context=context) > (tmp):
-            ans = ans.__sub__(Decimal((0, (1,), ans._exp)), context=context)
-
+        if exact:
+            # result is exact; rescale to use ideal exponent e
+            if shift >= 0:
+                # assert n % 10**shift == 0
+                n //= 10**shift
+            else:
+                n *= 10**-shift
+            e += shift
         else:
-            upper = ans.__add__(Decimal((0, (5,), ans._exp-1)),context=context)
-            context._set_rounding(ROUND_DOWN)
-            if upper.__mul__(upper, context=context) < tmp:
-                ans = ans.__add__(Decimal((0, (1,), ans._exp)),context=context)
+            # result is not exact; fix last digit as described above
+            if n % 5 == 0:
+                n += 1
 
-        ans._exp += expadd
+        ans = Decimal((0, list(map(int, str(n))), e))
 
-        context.prec = firstprec
-        context.rounding = rounding
+        # round, and fit to current context
+        context = context._shallow_copy()
+        rounding = context._set_rounding(ROUND_HALF_EVEN)
         ans = ans._fix(context)
+        context.rounding = rounding
 
-        rounding = context._set_rounding_decision(NEVER_ROUND)
-        if not ans.__mul__(ans, context=context) == self:
-            # Only rounded/inexact if here.
-            context._regard_flags(flags)
-            context._raise_error(Rounded)
-            context._raise_error(Inexact)
-        else:
-            # Exact answer, so let's set the exponent right.
-            # if self._exp < 0:
-            #    exp = (self._exp +1)// 2
-            # else:
-            exp = self._exp // 2
-            context.prec += ans._exp - exp
-            ans = ans._rescale(exp, context=context)
-            context.prec = firstprec
-            context._regard_flags(flags)
-        context.Emax, context.Emin = Emax, Emin
-
-        return ans._fix(context)
+        return ans
 
     def max(self, other, context=None):
         """Returns the larger value.
 
-        like max(self, other) except if one is not a number, returns
+        Like max(self, other) except if one is not a number, returns
         NaN (and signals if one is sNaN).  Also rounds.
         """
-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
+        other = _convert_other(other, raiseit=True)
+
+        if context is None:
+            context = getcontext()
 
         if self._is_special or other._is_special:
             # If one operand is a quiet NaN and the other is number, then the
@@ -2112,12 +2450,11 @@
             on = other._isnan()
             if sn or on:
                 if on == 1 and sn != 2:
-                    return self
+                    return self._fix_nan(context)
                 if sn == 1 and on != 2:
-                    return other
+                    return other._fix_nan(context)
                 return self._check_nans(other, context)
 
-        ans = self
         c = self.__cmp__(other)
         if c == 0:
             # If both operands are finite and equal in numerical value
@@ -2127,19 +2464,14 @@
             # positive sign and min returns the operand with the negative sign
             #
             # If the signs are the same then the exponent is used to select
-            # the result.
-            if self._sign != other._sign:
-                if self._sign:
-                    ans = other
-            elif self._exp < other._exp and not self._sign:
-                ans = other
-            elif self._exp > other._exp and self._sign:
-                ans = other
-        elif c == -1:
-            ans = other
+            # the result.  This is exactly the ordering used in compare_total.
+            c = self.compare_total(other)
 
-        if context is None:
-            context = getcontext()
+        if c == -1:
+            ans = other
+        else:
+            ans = self
+
         if context._rounding_decision == ALWAYS_ROUND:
             return ans._fix(context)
         return ans
@@ -2150,9 +2482,10 @@
         Like min(self, other) except if one is not a number, returns
         NaN (and signals if one is sNaN).  Also rounds.
         """
-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
+        other = _convert_other(other, raiseit=True)
+
+        if context is None:
+            context = getcontext()
 
         if self._is_special or other._is_special:
             # If one operand is a quiet NaN and the other is number, then the
@@ -2161,49 +2494,37 @@
             on = other._isnan()
             if sn or on:
                 if on == 1 and sn != 2:
-                    return self
+                    return self._fix_nan(context)
                 if sn == 1 and on != 2:
-                    return other
+                    return other._fix_nan(context)
                 return self._check_nans(other, context)
 
-        ans = self
         c = self.__cmp__(other)
         if c == 0:
-            # If both operands are finite and equal in numerical value
-            # then an ordering is applied:
-            #
-            # If the signs differ then max returns the operand with the
-            # positive sign and min returns the operand with the negative sign
-            #
-            # If the signs are the same then the exponent is used to select
-            # the result.
-            if self._sign != other._sign:
-                if other._sign:
-                    ans = other
-            elif self._exp > other._exp and not self._sign:
-                ans = other
-            elif self._exp < other._exp and self._sign:
-                ans = other
-        elif c == 1:
+            c = self.compare_total(other)
+
+        if c == -1:
+            ans = self
+        else:
             ans = other
 
-        if context is None:
-            context = getcontext()
         if context._rounding_decision == ALWAYS_ROUND:
             return ans._fix(context)
         return ans
 
     def _isinteger(self):
         """Returns whether self is an integer"""
+        if self._is_special:
+            return False
         if self._exp >= 0:
             return True
         rest = self._int[self._exp:]
         return rest == (0,)*len(rest)
 
     def _iseven(self):
-        """Returns 1 if self is even.  Assumes self is an integer."""
-        if self._exp > 0:
-            return 1
+        """Returns True if self is even.  Assumes self is an integer."""
+        if not self or self._exp > 0:
+            return True
         return self._int[-1+self._exp] & 1 == 0
 
     def adjusted(self):
@@ -2214,6 +2535,875 @@
         except TypeError:
             return 0
 
+    def canonical(self, context=None):
+        """Returns the same Decimal object.
+
+        As we do not have different encodings for the same number, the
+        received object already is in its canonical form.
+        """
+        return self
+
+    def compare_signal(self, other, context=None):
+        """Compares self to the other operand numerically.
+
+        It's pretty much like compare(), but all NaNs signal, with signaling
+        NaNs taking precedence over quiet NaNs.
+        """
+        if context is None:
+            context = getcontext()
+
+        self_is_nan = self._isnan()
+        other_is_nan = other._isnan()
+        if self_is_nan == 2:
+            return context._raise_error(InvalidOperation, 'sNaN',
+                                        1, self)
+        if other_is_nan == 2:
+            return context._raise_error(InvalidOperation, 'sNaN',
+                                        1, other)
+        if self_is_nan:
+            return context._raise_error(InvalidOperation, 'NaN in compare_signal',
+                                        1, self)
+        if other_is_nan:
+            return context._raise_error(InvalidOperation, 'NaN in compare_signal',
+                                        1, other)
+        return self.compare(other, context=context)
+
+    def compare_total(self, other):
+        """Compares self to other using the abstract representations.
+
+        This is not like the standard compare, which use their numerical
+        value. Note that a total ordering is defined for all possible abstract
+        representations.
+        """
+        # if one is negative and the other is positive, it's easy
+        if self._sign and not other._sign:
+            return Dec_n1
+        if not self._sign and other._sign:
+            return Dec_p1
+        sign = self._sign
+
+        # let's handle both NaN types
+        self_nan = self._isnan()
+        other_nan = other._isnan()
+        if self_nan or other_nan:
+            if self_nan == other_nan:
+                if self._int < other._int:
+                    if sign:
+                        return Dec_p1
+                    else:
+                        return Dec_n1
+                if self._int > other._int:
+                    if sign:
+                        return Dec_n1
+                    else:
+                        return Dec_p1
+                return Dec_0
+
+            if sign:
+                if self_nan == 1:
+                    return Dec_n1
+                if other_nan == 1:
+                    return Dec_p1
+                if self_nan == 2:
+                    return Dec_n1
+                if other_nan == 2:
+                    return Dec_p1
+            else:
+                if self_nan == 1:
+                    return Dec_p1
+                if other_nan == 1:
+                    return Dec_n1
+                if self_nan == 2:
+                    return Dec_p1
+                if other_nan == 2:
+                    return Dec_n1
+
+        if self < other:
+            return Dec_n1
+        if self > other:
+            return Dec_p1
+
+        if self._exp < other._exp:
+            if sign:
+                return Dec_p1
+            else:
+                return Dec_n1
+        if self._exp > other._exp:
+            if sign:
+                return Dec_n1
+            else:
+                return Dec_p1
+        return Dec_0
+
+
+    def compare_total_mag(self, other):
+        """Compares self to other using abstract repr., ignoring sign.
+
+        Like compare_total, but with operand's sign ignored and assumed to be 0.
+        """
+        s = self.copy_abs()
+        o = other.copy_abs()
+        return s.compare_total(o)
+
+    def copy_abs(self):
+        """Returns a copy with the sign set to 0. """
+        return Decimal((0, self._int, self._exp))
+
+    def copy_negate(self):
+        """Returns a copy with the sign inverted."""
+        if self._sign:
+            return Decimal((0, self._int, self._exp))
+        else:
+            return Decimal((1, self._int, self._exp))
+
+    def copy_sign(self, other):
+        """Returns self with the sign of other."""
+        return Decimal((other._sign, self._int, self._exp))
+
+    def exp(self, context=None):
+        """Returns e ** self."""
+
+        if context is None:
+            context = getcontext()
+
+        # exp(NaN) = NaN
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        # exp(-Infinity) = 0
+        if self._isinfinity() == -1:
+            return Dec_0
+
+        # exp(0) = 1
+        if not self:
+            return Dec_p1
+
+        # exp(Infinity) = Infinity
+        if self._isinfinity() == 1:
+            return Decimal(self)
+
+        # the result is now guaranteed to be inexact (the true
+        # mathematical result is transcendental). There's no need to
+        # raise Rounded and Inexact here---they'll always be raised as
+        # a result of the call to _fix.
+        p = context.prec
+        adj = self.adjusted()
+
+        # we only need to do any computation for quite a small range
+        # of adjusted exponents---for example, -29 <= adj <= 10 for
+        # the default context.  For smaller exponent the result is
+        # indistinguishable from 1 at the given precision, while for
+        # larger exponent the result either overflows or underflows.
+        if self._sign == 0 and adj > len(str((context.Emax+1)*3)):
+            # overflow
+            ans = Decimal((0, (1,), context.Emax+1))
+        elif self._sign == 1 and adj > len(str((-context.Etiny()+1)*3)):
+            # underflow to 0
+            ans = Decimal((0, (1,), context.Etiny()-1))
+        elif self._sign == 0 and adj < -p:
+            # p+1 digits; final round will raise correct flags
+            ans = Decimal((0, (1,) + (0,)*(p-1) + (1,), -p))
+        elif self._sign == 1 and adj < -p-1:
+            # p+1 digits; final round will raise correct flags
+            ans = Decimal((0, (9,)*(p+1), -p-1))
+        # general case
+        else:
+            op = _WorkRep(self)
+            c, e = op.int, op.exp
+            if op.sign == 1:
+                c = -c
+
+            # compute correctly rounded result: increase precision by
+            # 3 digits at a time until we get an unambiguously
+            # roundable result
+            extra = 3
+            while True:
+                coeff, exp = _dexp(c, e, p+extra)
+                if coeff % (5*10**(len(str(coeff))-p-1)):
+                    break
+                extra += 3
+
+            ans = Decimal((0, list(map(int, str(coeff))), exp))
+
+        # at this stage, ans should round correctly with *any*
+        # rounding mode, not just with ROUND_HALF_EVEN
+        context = context._shallow_copy()
+        rounding = context._set_rounding(ROUND_HALF_EVEN)
+        ans = ans._fix(context)
+        context.rounding = rounding
+
+        return ans
+
+    def is_canonical(self):
+        """Returns 1 if self is canonical; otherwise returns 0."""
+        return Dec_p1
+
+    def is_finite(self):
+        """Returns 1 if self is finite, otherwise returns 0.
+
+        For it to be finite, it must be neither infinite nor a NaN.
+        """
+        if self._is_special:
+            return Dec_0
+        else:
+            return Dec_p1
+
+    def is_infinite(self):
+        """Returns 1 if self is an Infinite, otherwise returns 0."""
+        if self._isinfinity():
+            return Dec_p1
+        else:
+            return Dec_0
+
+    def is_nan(self):
+        """Returns 1 if self is qNaN or sNaN, otherwise returns 0."""
+        if self._isnan():
+            return Dec_p1
+        else:
+            return Dec_0
+
+    def is_normal(self, context=None):
+        """Returns 1 if self is a normal number, otherwise returns 0."""
+        if self._is_special:
+            return Dec_0
+        if not self:
+            return Dec_0
+        if context is None:
+            context = getcontext()
+        if context.Emin <= self.adjusted() <= context.Emax:
+            return Dec_p1
+        else:
+            return Dec_0
+
+    def is_qnan(self):
+        """Returns 1 if self is a quiet NaN, otherwise returns 0."""
+        if self._isnan() == 1:
+            return Dec_p1
+        else:
+            return Dec_0
+
+    def is_signed(self):
+        """Returns 1 if self is negative, otherwise returns 0."""
+        return Decimal(self._sign)
+
+    def is_snan(self):
+        """Returns 1 if self is a signaling NaN, otherwise returns 0."""
+        if self._isnan() == 2:
+            return Dec_p1
+        else:
+            return Dec_0
+
+    def is_subnormal(self, context=None):
+        """Returns 1 if self is subnormal, otherwise returns 0."""
+        if self._is_special:
+            return Dec_0
+        if not self:
+            return Dec_0
+        if context is None:
+            context = getcontext()
+
+        r = self._exp + len(self._int)
+        if r <= context.Emin:
+            return Dec_p1
+        return Dec_0
+
+    def is_zero(self):
+        """Returns 1 if self is a zero, otherwise returns 0."""
+        if self:
+            return Dec_0
+        else:
+            return Dec_p1
+
+    def _ln_exp_bound(self):
+        """Compute a lower bound for the adjusted exponent of self.ln().
+        In other words, compute r such that self.ln() >= 10**r.  Assumes
+        that self is finite and positive and that self != 1.
+        """
+
+        # for 0.1 <= x <= 10 we use the inequalities 1-1/x <= ln(x) <= x-1
+        adj = self._exp + len(self._int) - 1
+        if adj >= 1:
+            # argument >= 10; we use 23/10 = 2.3 as a lower bound for ln(10)
+            return len(str(adj*23//10)) - 1
+        if adj <= -2:
+            # argument <= 0.1
+            return len(str((-1-adj)*23//10)) - 1
+        op = _WorkRep(self)
+        c, e = op.int, op.exp
+        if adj == 0:
+            # 1 < self < 10
+            num = str(c-10**-e)
+            den = str(c)
+            return len(num) - len(den) - (num < den)
+        # adj == -1, 0.1 <= self < 1
+        return e + len(str(10**-e - c)) - 1
+
+
+    def ln(self, context=None):
+        """Returns the natural (base e) logarithm of self."""
+
+        if context is None:
+            context = getcontext()
+
+        # ln(NaN) = NaN
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        # ln(0.0) == -Infinity
+        if not self:
+            return negInf
+
+        # ln(Infinity) = Infinity
+        if self._isinfinity() == 1:
+            return Inf
+
+        # ln(1.0) == 0.0
+        if self == Dec_p1:
+            return Dec_0
+
+        # ln(negative) raises InvalidOperation
+        if self._sign == 1:
+            return context._raise_error(InvalidOperation,
+                                        'ln of a negative value')
+
+        # result is irrational, so necessarily inexact
+        op = _WorkRep(self)
+        c, e = op.int, op.exp
+        p = context.prec
+
+        # correctly rounded result: repeatedly increase precision by 3
+        # until we get an unambiguously roundable result
+        places = p - self._ln_exp_bound() + 2 # at least p+3 places
+        while True:
+            coeff = _dlog(c, e, places)
+            # assert len(str(abs(coeff)))-p >= 1
+            if coeff % (5*10**(len(str(abs(coeff)))-p-1)):
+                break
+            places += 3
+        ans = Decimal((int(coeff<0), list(map(int, str(abs(coeff)))), -places))
+
+        context = context._shallow_copy()
+        rounding = context._set_rounding(ROUND_HALF_EVEN)
+        ans = ans._fix(context)
+        context.rounding = rounding
+        return ans
+
+    def _log10_exp_bound(self):
+        """Compute a lower bound for the adjusted exponent of self.log10().
+        In other words, find r such that self.log10() >= 10**r.
+        Assumes that self is finite and positive and that self != 1.
+        """
+
+        # For x >= 10 or x < 0.1 we only need a bound on the integer
+        # part of log10(self), and this comes directly from the
+        # exponent of x.  For 0.1 <= x <= 10 we use the inequalities
+        # 1-1/x <= log(x) <= x-1. If x > 1 we have |log10(x)| >
+        # (1-1/x)/2.31 > 0.  If x < 1 then |log10(x)| > (1-x)/2.31 > 0
+
+        adj = self._exp + len(self._int) - 1
+        if adj >= 1:
+            # self >= 10
+            return len(str(adj))-1
+        if adj <= -2:
+            # self < 0.1
+            return len(str(-1-adj))-1
+        op = _WorkRep(self)
+        c, e = op.int, op.exp
+        if adj == 0:
+            # 1 < self < 10
+            num = str(c-10**-e)
+            den = str(231*c)
+            return len(num) - len(den) - (num < den) + 2
+        # adj == -1, 0.1 <= self < 1
+        num = str(10**-e-c)
+        return len(num) + e - (num < "231") - 1
+
+    def log10(self, context=None):
+        """Returns the base 10 logarithm of self."""
+
+        if context is None:
+            context = getcontext()
+
+        # log10(NaN) = NaN
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        # log10(0.0) == -Infinity
+        if not self:
+            return negInf
+
+        # log10(Infinity) = Infinity
+        if self._isinfinity() == 1:
+            return Inf
+
+        # log10(negative or -Infinity) raises InvalidOperation
+        if self._sign == 1:
+            return context._raise_error(InvalidOperation,
+                                        'log10 of a negative value')
+
+        # log10(10**n) = n
+        if self._int[0] == 1 and self._int[1:] == (0,)*(len(self._int) - 1):
+            # answer may need rounding
+            ans = Decimal(self._exp + len(self._int) - 1)
+        else:
+            # result is irrational, so necessarily inexact
+            op = _WorkRep(self)
+            c, e = op.int, op.exp
+            p = context.prec
+
+            # correctly rounded result: repeatedly increase precision
+            # until result is unambiguously roundable
+            places = p-self._log10_exp_bound()+2
+            while True:
+                coeff = _dlog10(c, e, places)
+                # assert len(str(abs(coeff)))-p >= 1
+                if coeff % (5*10**(len(str(abs(coeff)))-p-1)):
+                    break
+                places += 3
+            ans = Decimal((int(coeff<0), list(map(int, str(abs(coeff)))),
+                          -places))
+
+        context = context._shallow_copy()
+        rounding = context._set_rounding(ROUND_HALF_EVEN)
+        ans = ans._fix(context)
+        context.rounding = rounding
+        return ans
+
+    def logb(self, context=None):
+        """ Returns the exponent of the magnitude of self's MSD.
+
+        The result is the integer which is the exponent of the magnitude
+        of the most significant digit of self (as though it were truncated
+        to a single digit while maintaining the value of that digit and
+        without limiting the resulting exponent).
+        """
+        # logb(NaN) = NaN
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        if context is None:
+            context = getcontext()
+
+        # logb(+/-Inf) = +Inf
+        if self._isinfinity():
+            return Inf
+
+        # logb(0) = -Inf, DivisionByZero
+        if not self:
+            return context._raise_error(DivisionByZero, 'logb(0)', 1)
+
+        # otherwise, simply return the adjusted exponent of self, as a
+        # Decimal.  Note that no attempt is made to fit the result
+        # into the current context.
+        return Decimal(self.adjusted())
+
+    def _islogical(self):
+        """Return True if self is a logical operand.
+
+        For being logical, it must be a finite numbers with a sign of 0,
+        an exponent of 0, and a coefficient whose digits must all be
+        either 0 or 1.
+        """
+        if self._sign != 0 or self._exp != 0:
+            return False
+        for dig in self._int:
+            if dig not in (0, 1):
+                return False
+        return True
+
+    def _fill_logical(self, context, opa, opb):
+        dif = context.prec - len(opa)
+        if dif > 0:
+            opa = (0,)*dif + opa
+        elif dif < 0:
+            opa = opa[-context.prec:]
+        dif = context.prec - len(opb)
+        if dif > 0:
+            opb = (0,)*dif + opb
+        elif dif < 0:
+            opb = opb[-context.prec:]
+        return opa, opb
+
+    def logical_and(self, other, context=None):
+        """Applies an 'and' operation between self and other's digits."""
+        if context is None:
+            context = getcontext()
+        if not self._islogical() or not other._islogical():
+            return context._raise_error(InvalidOperation)
+
+        # fill to context.prec
+        (opa, opb) = self._fill_logical(context, self._int, other._int)
+
+        # make the operation, and clean starting zeroes
+        result = [a&b for a,b in zip(opa,opb)]
+        for i,d in enumerate(result):
+            if d == 1:
+                break
+        result = tuple(result[i:])
+
+        # if empty, we must have at least a zero
+        if not result:
+            result = (0,)
+        return Decimal((0, result, 0))
+
+    def logical_invert(self, context=None):
+        """Invert all its digits."""
+        if context is None:
+            context = getcontext()
+        return self.logical_xor(Decimal((0,(1,)*context.prec,0)), context)
+
+    def logical_or(self, other, context=None):
+        """Applies an 'or' operation between self and other's digits."""
+        if context is None:
+            context = getcontext()
+        if not self._islogical() or not other._islogical():
+            return context._raise_error(InvalidOperation)
+
+        # fill to context.prec
+        (opa, opb) = self._fill_logical(context, self._int, other._int)
+
+        # make the operation, and clean starting zeroes
+        result = [a|b for a,b in zip(opa,opb)]
+        for i,d in enumerate(result):
+            if d == 1:
+                break
+        result = tuple(result[i:])
+
+        # if empty, we must have at least a zero
+        if not result:
+            result = (0,)
+        return Decimal((0, result, 0))
+
+    def logical_xor(self, other, context=None):
+        """Applies an 'xor' operation between self and other's digits."""
+        if context is None:
+            context = getcontext()
+        if not self._islogical() or not other._islogical():
+            return context._raise_error(InvalidOperation)
+
+        # fill to context.prec
+        (opa, opb) = self._fill_logical(context, self._int, other._int)
+
+        # make the operation, and clean starting zeroes
+        result = [a^b for a,b in zip(opa,opb)]
+        for i,d in enumerate(result):
+            if d == 1:
+                break
+        result = tuple(result[i:])
+
+        # if empty, we must have at least a zero
+        if not result:
+            result = (0,)
+        return Decimal((0, result, 0))
+
+    def max_mag(self, other, context=None):
+        """Compares the values numerically with their sign ignored."""
+        other = _convert_other(other, raiseit=True)
+
+        if context is None:
+            context = getcontext()
+
+        if self._is_special or other._is_special:
+            # If one operand is a quiet NaN and the other is number, then the
+            # number is always returned
+            sn = self._isnan()
+            on = other._isnan()
+            if sn or on:
+                if on == 1 and sn != 2:
+                    return self._fix_nan(context)
+                if sn == 1 and on != 2:
+                    return other._fix_nan(context)
+                return self._check_nans(other, context)
+
+        c = self.copy_abs().__cmp__(other.copy_abs())
+        if c == 0:
+            c = self.compare_total(other)
+
+        if c == -1:
+            ans = other
+        else:
+            ans = self
+
+        if context._rounding_decision == ALWAYS_ROUND:
+            return ans._fix(context)
+        return ans
+
+    def min_mag(self, other, context=None):
+        """Compares the values numerically with their sign ignored."""
+        other = _convert_other(other, raiseit=True)
+
+        if context is None:
+            context = getcontext()
+
+        if self._is_special or other._is_special:
+            # If one operand is a quiet NaN and the other is number, then the
+            # number is always returned
+            sn = self._isnan()
+            on = other._isnan()
+            if sn or on:
+                if on == 1 and sn != 2:
+                    return self._fix_nan(context)
+                if sn == 1 and on != 2:
+                    return other._fix_nan(context)
+                return self._check_nans(other, context)
+
+        c = self.copy_abs().__cmp__(other.copy_abs())
+        if c == 0:
+            c = self.compare_total(other)
+
+        if c == -1:
+            ans = self
+        else:
+            ans = other
+
+        if context._rounding_decision == ALWAYS_ROUND:
+            return ans._fix(context)
+        return ans
+
+    def next_minus(self, context=None):
+        """Returns the largest representable number smaller than itself."""
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        if self._isinfinity() == -1:
+            return negInf
+        if self._isinfinity() == 1:
+            return Decimal((0, (9,)*context.prec, context.Etop()))
+
+        context = context.copy()
+        context._set_rounding(ROUND_FLOOR)
+        context._ignore_all_flags()
+        new_self = self._fix(context)
+        if new_self != self:
+            return new_self
+        return self.__sub__(Decimal((0, (1,), context.Etiny()-1)), context)
+
+    def next_plus(self, context=None):
+        """Returns the smallest representable number larger than itself."""
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(context=context)
+        if ans:
+            return ans
+
+        if self._isinfinity() == 1:
+            return Inf
+        if self._isinfinity() == -1:
+            return Decimal((1, (9,)*context.prec, context.Etop()))
+
+        context = context.copy()
+        context._set_rounding(ROUND_CEILING)
+        context._ignore_all_flags()
+        new_self = self._fix(context)
+        if new_self != self:
+            return new_self
+        return self.__add__(Decimal((0, (1,), context.Etiny()-1)), context)
+
+    def next_toward(self, other, context=None):
+        """Returns the number closest to self, in the direction towards other.
+
+        The result is the closest representable number to self
+        (excluding self) that is in the direction towards other,
+        unless both have the same value.  If the two operands are
+        numerically equal, then the result is a copy of self with the
+        sign set to be the same as the sign of other.
+        """
+        other = _convert_other(other, raiseit=True)
+
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        comparison = self.__cmp__(other)
+        if comparison == 0:
+            return Decimal((other._sign, self._int, self._exp))
+
+        if comparison == -1:
+            ans = self.next_plus(context)
+        else: # comparison == 1
+            ans = self.next_minus(context)
+
+        # decide which flags to raise using value of ans
+        if ans._isinfinity():
+            context._raise_error(Overflow,
+                                 'Infinite result from next_toward',
+                                 ans._sign)
+            context._raise_error(Rounded)
+            context._raise_error(Inexact)
+        elif ans.adjusted() < context.Emin:
+            context._raise_error(Underflow)
+            context._raise_error(Subnormal)
+            context._raise_error(Rounded)
+            context._raise_error(Inexact)
+            # if precision == 1 then we don't raise Clamped for a
+            # result 0E-Etiny.
+            if not ans:
+                context._raise_error(Clamped)
+
+        return ans
+
+    def number_class(self, context=None):
+        """Returns an indication of the class of self.
+
+        The class is one of the following strings:
+          -sNaN
+          -NaN
+          -Infinity
+          -Normal
+          -Subnormal
+          -Zero
+          +Zero
+          +Subnormal
+          +Normal
+          +Infinity
+        """
+        if self.is_snan():
+            return "sNaN"
+        if self.is_qnan():
+            return "NaN"
+        inf = self._isinfinity()
+        if inf == 1:
+            return "+Infinity"
+        if inf == -1:
+            return "-Infinity"
+        if self.is_zero():
+            if self._sign:
+                return "-Zero"
+            else:
+                return "+Zero"
+        if context is None:
+            context = getcontext()
+        if self.is_subnormal(context=context):
+            if self._sign:
+                return "-Subnormal"
+            else:
+                return "+Subnormal"
+        # just a normal, regular, boring number, :)
+        if self._sign:
+            return "-Normal"
+        else:
+            return "+Normal"
+
+    def radix(self):
+        """Just returns 10, as this is Decimal, :)"""
+        return Decimal(10)
+
+    def rotate(self, other, context=None):
+        """Returns a rotated copy of self, value-of-other times."""
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        if other._exp != 0:
+            return context._raise_error(InvalidOperation)
+        if not (-context.prec <= int(other) <= context.prec):
+            return context._raise_error(InvalidOperation)
+
+        if self._isinfinity():
+            return Decimal(self)
+
+        # get values, pad if necessary
+        torot = int(other)
+        rotdig = self._int
+        topad = context.prec - len(rotdig)
+        if topad:
+            rotdig = ((0,)*topad) + rotdig
+
+        # let's rotate!
+        rotated = rotdig[torot:] + rotdig[:torot]
+
+        # clean starting zeroes
+        for i,d in enumerate(rotated):
+            if d != 0:
+                break
+        rotated = rotated[i:]
+
+        return Decimal((self._sign, rotated, self._exp))
+
+
+    def scaleb (self, other, context=None):
+        """Returns self operand after adding the second value to its exp."""
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        if other._exp != 0:
+            return context._raise_error(InvalidOperation)
+        liminf = -2 * (context.Emax + context.prec)
+        limsup =  2 * (context.Emax + context.prec)
+        if not (liminf <= int(other) <= limsup):
+            return context._raise_error(InvalidOperation)
+
+        if self._isinfinity():
+            return Decimal(self)
+
+        d = Decimal((self._sign, self._int, self._exp + int(other)))
+        d = d._fix(context)
+        return d
+
+    def shift(self, other, context=None):
+        """Returns a shifted copy of self, value-of-other times."""
+        if context is None:
+            context = getcontext()
+
+        ans = self._check_nans(other, context)
+        if ans:
+            return ans
+
+        if other._exp != 0:
+            return context._raise_error(InvalidOperation)
+        if not (-context.prec <= int(other) <= context.prec):
+            return context._raise_error(InvalidOperation)
+
+        if self._isinfinity():
+            return Decimal(self)
+
+        # get values, pad if necessary
+        torot = int(other)
+        if not torot:
+            return Decimal(self)
+        rotdig = self._int
+        topad = context.prec - len(rotdig)
+        if topad:
+            rotdig = ((0,)*topad) + rotdig
+
+        # let's shift!
+        if torot < 0:
+            rotated = rotdig[:torot]
+        else:
+            rotated = (rotdig + ((0,) * torot))
+            rotated = rotated[-context.prec:]
+
+        # clean starting zeroes
+        if rotated:
+            for i,d in enumerate(rotated):
+                if d != 0:
+                    break
+            rotated = rotated[i:]
+        else:
+            rotated = (0,)
+
+        return Decimal((self._sign, rotated, self._exp))
+
+
     # Support for pickling, copy, and deepcopy
     def __reduce__(self):
         return (self.__class__, (str(self),))
@@ -2426,6 +3616,9 @@
     def create_decimal(self, num='0'):
         """Creates a new Decimal instance but using self as context."""
         d = Decimal(num, context=self)
+        if d._isnan() and len(d._int) > self.prec - self._clamp:
+            return self._raise_error(ConversionSyntax,
+                                     "diagnostic info too long in NaN")
         return d._fix(self)
 
     # Methods
@@ -2460,6 +3653,17 @@
     def _apply(self, a):
         return str(a._fix(self))
 
+    def canonical(self, a):
+        """Returns the same Decimal object.
+
+        As we do not have different encodings for the same number, the
+        received object already is in its canonical form.
+
+        >>> ExtendedContext.canonical(Decimal('2.50'))
+        Decimal("2.50")
+        """
+        return a.canonical(context=self)
+
     def compare(self, a, b):
         """Compares values numerically.
 
@@ -2489,6 +3693,110 @@
         """
         return a.compare(b, context=self)
 
+    def compare_signal(self, a, b):
+        """Compares the values of the two operands numerically.
+
+        It's pretty much like compare(), but all NaNs signal, with signaling
+        NaNs taking precedence over quiet NaNs.
+
+        >>> c = ExtendedContext
+        >>> c.compare_signal(Decimal('2.1'), Decimal('3'))
+        Decimal("-1")
+        >>> c.compare_signal(Decimal('2.1'), Decimal('2.1'))
+        Decimal("0")
+        >>> c.flags[InvalidOperation] = 0
+        >>> print(c.flags[InvalidOperation])
+        0
+        >>> c.compare_signal(Decimal('NaN'), Decimal('2.1'))
+        Decimal("NaN")
+        >>> print(c.flags[InvalidOperation])
+        1
+        >>> c.flags[InvalidOperation] = 0
+        >>> print(c.flags[InvalidOperation])
+        0
+        >>> c.compare_signal(Decimal('sNaN'), Decimal('2.1'))
+        Decimal("NaN")
+        >>> print(c.flags[InvalidOperation])
+        1
+        """
+        return a.compare_signal(b, context=self)
+
+    def compare_total(self, a, b):
+        """Compares two operands using their abstract representation.
+
+        This is not like the standard compare, which use their numerical
+        value. Note that a total ordering is defined for all possible abstract
+        representations.
+
+        >>> ExtendedContext.compare_total(Decimal('12.73'), Decimal('127.9'))
+        Decimal("-1")
+        >>> ExtendedContext.compare_total(Decimal('-127'),  Decimal('12'))
+        Decimal("-1")
+        >>> ExtendedContext.compare_total(Decimal('12.30'), Decimal('12.3'))
+        Decimal("-1")
+        >>> ExtendedContext.compare_total(Decimal('12.30'), Decimal('12.30'))
+        Decimal("0")
+        >>> ExtendedContext.compare_total(Decimal('12.3'),  Decimal('12.300'))
+        Decimal("1")
+        >>> ExtendedContext.compare_total(Decimal('12.3'),  Decimal('NaN'))
+        Decimal("-1")
+        """
+        return a.compare_total(b)
+
+    def compare_total_mag(self, a, b):
+        """Compares two operands using their abstract representation ignoring sign.
+
+        Like compare_total, but with operand's sign ignored and assumed to be 0.
+        """
+        return a.compare_total_mag(b)
+
+    def copy_abs(self, a):
+        """Returns a copy of the operand with the sign set to 0.
+
+        >>> ExtendedContext.copy_abs(Decimal('2.1'))
+        Decimal("2.1")
+        >>> ExtendedContext.copy_abs(Decimal('-100'))
+        Decimal("100")
+        """
+        return a.copy_abs()
+
+    def copy_decimal(self, a):
+        """Returns a copy of the decimal objet.
+
+        >>> ExtendedContext.copy_decimal(Decimal('2.1'))
+        Decimal("2.1")
+        >>> ExtendedContext.copy_decimal(Decimal('-1.00'))
+        Decimal("-1.00")
+        """
+        return Decimal(a)
+
+    def copy_negate(self, a):
+        """Returns a copy of the operand with the sign inverted.
+
+        >>> ExtendedContext.copy_negate(Decimal('101.5'))
+        Decimal("-101.5")
+        >>> ExtendedContext.copy_negate(Decimal('-101.5'))
+        Decimal("101.5")
+        """
+        return a.copy_negate()
+
+    def copy_sign(self, a, b):
+        """Copies the second operand's sign to the first one.
+
+        In detail, it returns a copy of the first operand with the sign
+        equal to the sign of the second operand.
+
+        >>> ExtendedContext.copy_sign(Decimal( '1.50'), Decimal('7.33'))
+        Decimal("1.50")
+        >>> ExtendedContext.copy_sign(Decimal('-1.50'), Decimal('7.33'))
+        Decimal("1.50")
+        >>> ExtendedContext.copy_sign(Decimal( '1.50'), Decimal('-7.33'))
+        Decimal("-1.50")
+        >>> ExtendedContext.copy_sign(Decimal('-1.50'), Decimal('-7.33'))
+        Decimal("-1.50")
+        """
+        return a.copy_sign(b)
+
     def divide(self, a, b):
         """Decimal division in a specified context.
 
@@ -2530,6 +3838,316 @@
     def divmod(self, a, b):
         return a.__divmod__(b, context=self)
 
+    def exp(self, a):
+        """Returns e ** a.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.exp(Decimal('-Infinity'))
+        Decimal("0")
+        >>> c.exp(Decimal('-1'))
+        Decimal("0.367879441")
+        >>> c.exp(Decimal('0'))
+        Decimal("1")
+        >>> c.exp(Decimal('1'))
+        Decimal("2.71828183")
+        >>> c.exp(Decimal('0.693147181'))
+        Decimal("2.00000000")
+        >>> c.exp(Decimal('+Infinity'))
+        Decimal("Infinity")
+        """
+        return a.exp(context=self)
+
+    def fma(self, a, b, c):
+        """Returns a multiplied by b, plus c.
+
+        The first two operands are multiplied together, using multiply,
+        the third operand is then added to the result of that
+        multiplication, using add, all with only one final rounding.
+
+        >>> ExtendedContext.fma(Decimal('3'), Decimal('5'), Decimal('7'))
+        Decimal("22")
+        >>> ExtendedContext.fma(Decimal('3'), Decimal('-5'), Decimal('7'))
+        Decimal("-8")
+        >>> ExtendedContext.fma(Decimal('888565290'), Decimal('1557.96930'), Decimal('-86087.7578'))
+        Decimal("1.38435736E+12")
+        """
+        return a.fma(b, c, context=self)
+
+    def is_canonical(self, a):
+        """Returns 1 if the operand is canonical; otherwise returns 0.
+
+        >>> ExtendedContext.is_canonical(Decimal('2.50'))
+        Decimal("1")
+        """
+        return Dec_p1
+
+    def is_finite(self, a):
+        """Returns 1 if the operand is finite, otherwise returns 0.
+
+        For it to be finite, it must be neither infinite nor a NaN.
+
+        >>> ExtendedContext.is_finite(Decimal('2.50'))
+        Decimal("1")
+        >>> ExtendedContext.is_finite(Decimal('-0.3'))
+        Decimal("1")
+        >>> ExtendedContext.is_finite(Decimal('0'))
+        Decimal("1")
+        >>> ExtendedContext.is_finite(Decimal('Inf'))
+        Decimal("0")
+        >>> ExtendedContext.is_finite(Decimal('NaN'))
+        Decimal("0")
+        """
+        return a.is_finite()
+
+    def is_infinite(self, a):
+        """Returns 1 if the operand is an Infinite, otherwise returns 0.
+
+        >>> ExtendedContext.is_infinite(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_infinite(Decimal('-Inf'))
+        Decimal("1")
+        >>> ExtendedContext.is_infinite(Decimal('NaN'))
+        Decimal("0")
+        """
+        return a.is_infinite()
+
+    def is_nan(self, a):
+        """Returns 1 if the operand is qNaN or sNaN, otherwise returns 0.
+
+        >>> ExtendedContext.is_nan(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_nan(Decimal('NaN'))
+        Decimal("1")
+        >>> ExtendedContext.is_nan(Decimal('-sNaN'))
+        Decimal("1")
+        """
+        return a.is_nan()
+
+    def is_normal(self, a):
+        """Returns 1 if the operand is a normal number, otherwise returns 0.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.is_normal(Decimal('2.50'))
+        Decimal("1")
+        >>> c.is_normal(Decimal('0.1E-999'))
+        Decimal("0")
+        >>> c.is_normal(Decimal('0.00'))
+        Decimal("0")
+        >>> c.is_normal(Decimal('-Inf'))
+        Decimal("0")
+        >>> c.is_normal(Decimal('NaN'))
+        Decimal("0")
+        """
+        return a.is_normal(context=self)
+
+    def is_qnan(self, a):
+        """Returns 1 if the operand is a quiet NaN, otherwise returns 0.
+
+        >>> ExtendedContext.is_qnan(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_qnan(Decimal('NaN'))
+        Decimal("1")
+        >>> ExtendedContext.is_qnan(Decimal('sNaN'))
+        Decimal("0")
+        """
+        return a.is_qnan()
+
+    def is_signed(self, a):
+        """Returns 1 if the operand is negative, otherwise returns 0.
+
+        >>> ExtendedContext.is_signed(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_signed(Decimal('-12'))
+        Decimal("1")
+        >>> ExtendedContext.is_signed(Decimal('-0'))
+        Decimal("1")
+        """
+        return a.is_signed()
+
+    def is_snan(self, a):
+        """Returns 1 if the operand is a signaling NaN, otherwise returns 0.
+
+        >>> ExtendedContext.is_snan(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_snan(Decimal('NaN'))
+        Decimal("0")
+        >>> ExtendedContext.is_snan(Decimal('sNaN'))
+        Decimal("1")
+        """
+        return a.is_snan()
+
+    def is_subnormal(self, a):
+        """Returns 1 if the operand is subnormal, otherwise returns 0.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.is_subnormal(Decimal('2.50'))
+        Decimal("0")
+        >>> c.is_subnormal(Decimal('0.1E-999'))
+        Decimal("1")
+        >>> c.is_subnormal(Decimal('0.00'))
+        Decimal("0")
+        >>> c.is_subnormal(Decimal('-Inf'))
+        Decimal("0")
+        >>> c.is_subnormal(Decimal('NaN'))
+        Decimal("0")
+        """
+        return a.is_subnormal(context=self)
+
+    def is_zero(self, a):
+        """Returns 1 if the operand is a zero, otherwise returns 0.
+
+        >>> ExtendedContext.is_zero(Decimal('0'))
+        Decimal("1")
+        >>> ExtendedContext.is_zero(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.is_zero(Decimal('-0E+2'))
+        Decimal("1")
+        """
+        return a.is_zero()
+
+    def ln(self, a):
+        """Returns the natural (base e) logarithm of the operand.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.ln(Decimal('0'))
+        Decimal("-Infinity")
+        >>> c.ln(Decimal('1.000'))
+        Decimal("0")
+        >>> c.ln(Decimal('2.71828183'))
+        Decimal("1.00000000")
+        >>> c.ln(Decimal('10'))
+        Decimal("2.30258509")
+        >>> c.ln(Decimal('+Infinity'))
+        Decimal("Infinity")
+        """
+        return a.ln(context=self)
+
+    def log10(self, a):
+        """Returns the base 10 logarithm of the operand.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.log10(Decimal('0'))
+        Decimal("-Infinity")
+        >>> c.log10(Decimal('0.001'))
+        Decimal("-3")
+        >>> c.log10(Decimal('1.000'))
+        Decimal("0")
+        >>> c.log10(Decimal('2'))
+        Decimal("0.301029996")
+        >>> c.log10(Decimal('10'))
+        Decimal("1")
+        >>> c.log10(Decimal('70'))
+        Decimal("1.84509804")
+        >>> c.log10(Decimal('+Infinity'))
+        Decimal("Infinity")
+        """
+        return a.log10(context=self)
+
+    def logb(self, a):
+        """ Returns the exponent of the magnitude of the operand's MSD.
+
+        The result is the integer which is the exponent of the magnitude
+        of the most significant digit of the operand (as though the
+        operand were truncated to a single digit while maintaining the
+        value of that digit and without limiting the resulting exponent).
+
+        >>> ExtendedContext.logb(Decimal('250'))
+        Decimal("2")
+        >>> ExtendedContext.logb(Decimal('2.50'))
+        Decimal("0")
+        >>> ExtendedContext.logb(Decimal('0.03'))
+        Decimal("-2")
+        >>> ExtendedContext.logb(Decimal('0'))
+        Decimal("-Infinity")
+        """
+        return a.logb(context=self)
+
+    def logical_and(self, a, b):
+        """Applies the logical operation 'and' between each operand's digits.
+
+        The operands must be both logical numbers.
+
+        >>> ExtendedContext.logical_and(Decimal('0'), Decimal('0'))
+        Decimal("0")
+        >>> ExtendedContext.logical_and(Decimal('0'), Decimal('1'))
+        Decimal("0")
+        >>> ExtendedContext.logical_and(Decimal('1'), Decimal('0'))
+        Decimal("0")
+        >>> ExtendedContext.logical_and(Decimal('1'), Decimal('1'))
+        Decimal("1")
+        >>> ExtendedContext.logical_and(Decimal('1100'), Decimal('1010'))
+        Decimal("1000")
+        >>> ExtendedContext.logical_and(Decimal('1111'), Decimal('10'))
+        Decimal("10")
+        """
+        return a.logical_and(b, context=self)
+
+    def logical_invert(self, a):
+        """Invert all the digits in the operand.
+
+        The operand must be a logical number.
+
+        >>> ExtendedContext.logical_invert(Decimal('0'))
+        Decimal("111111111")
+        >>> ExtendedContext.logical_invert(Decimal('1'))
+        Decimal("111111110")
+        >>> ExtendedContext.logical_invert(Decimal('111111111'))
+        Decimal("0")
+        >>> ExtendedContext.logical_invert(Decimal('101010101'))
+        Decimal("10101010")
+        """
+        return a.logical_invert(context=self)
+
+    def logical_or(self, a, b):
+        """Applies the logical operation 'or' between each operand's digits.
+
+        The operands must be both logical numbers.
+
+        >>> ExtendedContext.logical_or(Decimal('0'), Decimal('0'))
+        Decimal("0")
+        >>> ExtendedContext.logical_or(Decimal('0'), Decimal('1'))
+        Decimal("1")
+        >>> ExtendedContext.logical_or(Decimal('1'), Decimal('0'))
+        Decimal("1")
+        >>> ExtendedContext.logical_or(Decimal('1'), Decimal('1'))
+        Decimal("1")
+        >>> ExtendedContext.logical_or(Decimal('1100'), Decimal('1010'))
+        Decimal("1110")
+        >>> ExtendedContext.logical_or(Decimal('1110'), Decimal('10'))
+        Decimal("1110")
+        """
+        return a.logical_or(b, context=self)
+
+    def logical_xor(self, a, b):
+        """Applies the logical operation 'xor' between each operand's digits.
+
+        The operands must be both logical numbers.
+
+        >>> ExtendedContext.logical_xor(Decimal('0'), Decimal('0'))
+        Decimal("0")
+        >>> ExtendedContext.logical_xor(Decimal('0'), Decimal('1'))
+        Decimal("1")
+        >>> ExtendedContext.logical_xor(Decimal('1'), Decimal('0'))
+        Decimal("1")
+        >>> ExtendedContext.logical_xor(Decimal('1'), Decimal('1'))
+        Decimal("0")
+        >>> ExtendedContext.logical_xor(Decimal('1100'), Decimal('1010'))
+        Decimal("110")
+        >>> ExtendedContext.logical_xor(Decimal('1111'), Decimal('10'))
+        Decimal("1101")
+        """
+        return a.logical_xor(b, context=self)
+
     def max(self, a,b):
         """max compares two values numerically and returns the maximum.
 
@@ -2550,6 +4168,10 @@
         """
         return a.max(b, context=self)
 
+    def max_mag(self, a, b):
+        """Compares the values numerically with their sign ignored."""
+        return a.max_mag(b, context=self)
+
     def min(self, a,b):
         """min compares two values numerically and returns the minimum.
 
@@ -2570,6 +4192,10 @@
         """
         return a.min(b, context=self)
 
+    def min_mag(self, a, b):
+        """Compares the values numerically with their sign ignored."""
+        return a.min_mag(b, context=self)
+
     def minus(self, a):
         """Minus corresponds to unary prefix minus in Python.
 
@@ -2605,6 +4231,68 @@
         """
         return a.__mul__(b, context=self)
 
+    def next_minus(self, a):
+        """Returns the largest representable number smaller than a.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> ExtendedContext.next_minus(Decimal('1'))
+        Decimal("0.999999999")
+        >>> c.next_minus(Decimal('1E-1007'))
+        Decimal("0E-1007")
+        >>> ExtendedContext.next_minus(Decimal('-1.00000003'))
+        Decimal("-1.00000004")
+        >>> c.next_minus(Decimal('Infinity'))
+        Decimal("9.99999999E+999")
+        """
+        return a.next_minus(context=self)
+
+    def next_plus(self, a):
+        """Returns the smallest representable number larger than a.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> ExtendedContext.next_plus(Decimal('1'))
+        Decimal("1.00000001")
+        >>> c.next_plus(Decimal('-1E-1007'))
+        Decimal("-0E-1007")
+        >>> ExtendedContext.next_plus(Decimal('-1.00000003'))
+        Decimal("-1.00000002")
+        >>> c.next_plus(Decimal('-Infinity'))
+        Decimal("-9.99999999E+999")
+        """
+        return a.next_plus(context=self)
+
+    def next_toward(self, a, b):
+        """Returns the number closest to a, in direction towards b.
+
+        The result is the closest representable number from the first
+        operand (but not the first operand) that is in the direction
+        towards the second operand, unless the operands have the same
+        value.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.next_toward(Decimal('1'), Decimal('2'))
+        Decimal("1.00000001")
+        >>> c.next_toward(Decimal('-1E-1007'), Decimal('1'))
+        Decimal("-0E-1007")
+        >>> c.next_toward(Decimal('-1.00000003'), Decimal('0'))
+        Decimal("-1.00000002")
+        >>> c.next_toward(Decimal('1'), Decimal('0'))
+        Decimal("0.999999999")
+        >>> c.next_toward(Decimal('1E-1007'), Decimal('-100'))
+        Decimal("0E-1007")
+        >>> c.next_toward(Decimal('-1.00000003'), Decimal('-10'))
+        Decimal("-1.00000004")
+        >>> c.next_toward(Decimal('0.00'), Decimal('-0.0000'))
+        Decimal("-0.00")
+        """
+        return a.next_toward(b, context=self)
+
     def normalize(self, a):
         """normalize reduces an operand to its simplest form.
 
@@ -2626,6 +4314,53 @@
         """
         return a.normalize(context=self)
 
+    def number_class(self, a):
+        """Returns an indication of the class of the operand.
+
+        The class is one of the following strings:
+          -sNaN
+          -NaN
+          -Infinity
+          -Normal
+          -Subnormal
+          -Zero
+          +Zero
+          +Subnormal
+          +Normal
+          +Infinity
+
+        >>> c = Context(ExtendedContext)
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.number_class(Decimal('Infinity'))
+        '+Infinity'
+        >>> c.number_class(Decimal('1E-10'))
+        '+Normal'
+        >>> c.number_class(Decimal('2.50'))
+        '+Normal'
+        >>> c.number_class(Decimal('0.1E-999'))
+        '+Subnormal'
+        >>> c.number_class(Decimal('0'))
+        '+Zero'
+        >>> c.number_class(Decimal('-0'))
+        '-Zero'
+        >>> c.number_class(Decimal('-0.1E-999'))
+        '-Subnormal'
+        >>> c.number_class(Decimal('-1E-10'))
+        '-Normal'
+        >>> c.number_class(Decimal('-2.50'))
+        '-Normal'
+        >>> c.number_class(Decimal('-Infinity'))
+        '-Infinity'
+        >>> c.number_class(Decimal('NaN'))
+        'NaN'
+        >>> c.number_class(Decimal('-NaN'))
+        'NaN'
+        >>> c.number_class(Decimal('sNaN'))
+        'sNaN'
+        """
+        return a.number_class(context=self)
+
     def plus(self, a):
         """Plus corresponds to unary prefix plus in Python.
 
@@ -2643,49 +4378,69 @@
     def power(self, a, b, modulo=None):
         """Raises a to the power of b, to modulo if given.
 
-        The right-hand operand must be a whole number whose integer part (after
-        any exponent has been applied) has no more than 9 digits and whose
-        fractional part (if any) is all zeros before any rounding.  The operand
-        may be positive, negative, or zero; if negative, the absolute value of
-        the power is used, and the left-hand operand is inverted (divided into
-        1) before use.
+        With two arguments, compute a**b.  If a is negative then b
+        must be integral.  The result will be inexact unless b is
+        integral and the result is finite and can be expressed exactly
+        in 'precision' digits.
 
-        If the increased precision needed for the intermediate calculations
-        exceeds the capabilities of the implementation then an Invalid
-        operation condition is raised.
+        With three arguments, compute (a**b) % modulo.  For the
+        three argument form, the following restrictions on the
+        arguments hold:
 
-        If, when raising to a negative power, an underflow occurs during the
-        division into 1, the operation is not halted at that point but
-        continues.
+         - all three arguments must be integral
+         - b must be nonnegative
+         - at least one of a or b must be nonzero
+         - modulo must be nonzero and have at most 'precision' digits
 
-        >>> ExtendedContext.power(Decimal('2'), Decimal('3'))
+        The result of pow(a, b, modulo) is identical to the result
+        that would be obtained by computing (a**b) % modulo with
+        unbounded precision, but is computed more efficiently.  It is
+        always exact.
+
+        >>> c = ExtendedContext.copy()
+        >>> c.Emin = -999
+        >>> c.Emax = 999
+        >>> c.power(Decimal('2'), Decimal('3'))
         Decimal("8")
-        >>> ExtendedContext.power(Decimal('2'), Decimal('-3'))
+        >>> c.power(Decimal('-2'), Decimal('3'))
+        Decimal("-8")
+        >>> c.power(Decimal('2'), Decimal('-3'))
         Decimal("0.125")
-        >>> ExtendedContext.power(Decimal('1.7'), Decimal('8'))
+        >>> c.power(Decimal('1.7'), Decimal('8'))
         Decimal("69.7575744")
-        >>> ExtendedContext.power(Decimal('Infinity'), Decimal('-2'))
+        >>> c.power(Decimal('10'), Decimal('0.301029996'))
+        Decimal("2.00000000")
+        >>> c.power(Decimal('Infinity'), Decimal('-1'))
         Decimal("0")
-        >>> ExtendedContext.power(Decimal('Infinity'), Decimal('-1'))
-        Decimal("0")
-        >>> ExtendedContext.power(Decimal('Infinity'), Decimal('0'))
+        >>> c.power(Decimal('Infinity'), Decimal('0'))
         Decimal("1")
-        >>> ExtendedContext.power(Decimal('Infinity'), Decimal('1'))
+        >>> c.power(Decimal('Infinity'), Decimal('1'))
         Decimal("Infinity")
-        >>> ExtendedContext.power(Decimal('Infinity'), Decimal('2'))
-        Decimal("Infinity")
-        >>> ExtendedContext.power(Decimal('-Infinity'), Decimal('-2'))
-        Decimal("0")
-        >>> ExtendedContext.power(Decimal('-Infinity'), Decimal('-1'))
+        >>> c.power(Decimal('-Infinity'), Decimal('-1'))
         Decimal("-0")
-        >>> ExtendedContext.power(Decimal('-Infinity'), Decimal('0'))
+        >>> c.power(Decimal('-Infinity'), Decimal('0'))
         Decimal("1")
-        >>> ExtendedContext.power(Decimal('-Infinity'), Decimal('1'))
+        >>> c.power(Decimal('-Infinity'), Decimal('1'))
         Decimal("-Infinity")
-        >>> ExtendedContext.power(Decimal('-Infinity'), Decimal('2'))
+        >>> c.power(Decimal('-Infinity'), Decimal('2'))
         Decimal("Infinity")
-        >>> ExtendedContext.power(Decimal('0'), Decimal('0'))
+        >>> c.power(Decimal('0'), Decimal('0'))
         Decimal("NaN")
+
+        >>> c.power(Decimal('3'), Decimal('7'), Decimal('16'))
+        Decimal("11")
+        >>> c.power(Decimal('-3'), Decimal('7'), Decimal('16'))
+        Decimal("-11")
+        >>> c.power(Decimal('-3'), Decimal('8'), Decimal('16'))
+        Decimal("1")
+        >>> c.power(Decimal('3'), Decimal('7'), Decimal('-16'))
+        Decimal("11")
+        >>> c.power(Decimal('23E12345'), Decimal('67E189'), Decimal('123456789'))
+        Decimal("11729830")
+        >>> c.power(Decimal('-0'), Decimal('17'), Decimal('1729'))
+        Decimal("-0")
+        >>> c.power(Decimal('-23'), Decimal('0'), Decimal('65537'))
+        Decimal("1")
         """
         return a.__pow__(b, modulo, context=self)
 
@@ -2740,6 +4495,14 @@
         """
         return a.quantize(b, context=self)
 
+    def radix(self):
+        """Just returns 10, as this is Decimal, :)
+
+        >>> ExtendedContext.radix()
+        Decimal("10")
+        """
+        return Decimal(10)
+
     def remainder(self, a, b):
         """Returns the remainder from integer division.
 
@@ -2794,6 +4557,28 @@
         """
         return a.remainder_near(b, context=self)
 
+    def rotate(self, a, b):
+        """Returns a rotated copy of a, b times.
+
+        The coefficient of the result is a rotated copy of the digits in
+        the coefficient of the first operand.  The number of places of
+        rotation is taken from the absolute value of the second operand,
+        with the rotation being to the left if the second operand is
+        positive or to the right otherwise.
+
+        >>> ExtendedContext.rotate(Decimal('34'), Decimal('8'))
+        Decimal("400000003")
+        >>> ExtendedContext.rotate(Decimal('12'), Decimal('9'))
+        Decimal("12")
+        >>> ExtendedContext.rotate(Decimal('123456789'), Decimal('-2'))
+        Decimal("891234567")
+        >>> ExtendedContext.rotate(Decimal('123456789'), Decimal('0'))
+        Decimal("123456789")
+        >>> ExtendedContext.rotate(Decimal('123456789'), Decimal('+2'))
+        Decimal("345678912")
+        """
+        return a.rotate(b, context=self)
+
     def same_quantum(self, a, b):
         """Returns True if the two operands have the same exponent.
 
@@ -2811,6 +4596,41 @@
         """
         return a.same_quantum(b)
 
+    def scaleb (self, a, b):
+        """Returns the first operand after adding the second value its exp.
+
+        >>> ExtendedContext.scaleb(Decimal('7.50'), Decimal('-2'))
+        Decimal("0.0750")
+        >>> ExtendedContext.scaleb(Decimal('7.50'), Decimal('0'))
+        Decimal("7.50")
+        >>> ExtendedContext.scaleb(Decimal('7.50'), Decimal('3'))
+        Decimal("7.50E+3")
+        """
+        return a.scaleb (b, context=self)
+
+    def shift(self, a, b):
+        """Returns a shifted copy of a, b times.
+
+        The coefficient of the result is a shifted copy of the digits
+        in the coefficient of the first operand.  The number of places
+        to shift is taken from the absolute value of the second operand,
+        with the shift being to the left if the second operand is
+        positive or to the right otherwise.  Digits shifted into the
+        coefficient are zeros.
+
+        >>> ExtendedContext.shift(Decimal('34'), Decimal('8'))
+        Decimal("400000000")
+        >>> ExtendedContext.shift(Decimal('12'), Decimal('9'))
+        Decimal("0")
+        >>> ExtendedContext.shift(Decimal('123456789'), Decimal('-2'))
+        Decimal("1234567")
+        >>> ExtendedContext.shift(Decimal('123456789'), Decimal('0'))
+        Decimal("123456789")
+        >>> ExtendedContext.shift(Decimal('123456789'), Decimal('+2'))
+        Decimal("345678900")
+        """
+        return a.shift(b, context=self)
+
     def sqrt(self, a):
         """Square root of a non-negative number to context precision.
 
@@ -2866,7 +4686,36 @@
         """
         return a.__str__(context=self)
 
-    def to_integral(self, a):
+    def to_integral_exact(self, a):
+        """Rounds to an integer.
+
+        When the operand has a negative exponent, the result is the same
+        as using the quantize() operation using the given operand as the
+        left-hand-operand, 1E+0 as the right-hand-operand, and the precision
+        of the operand as the precision setting; Inexact and Rounded flags
+        are allowed in this operation.  The rounding mode is taken from the
+        context.
+
+        >>> ExtendedContext.to_integral_exact(Decimal('2.1'))
+        Decimal("2")
+        >>> ExtendedContext.to_integral_exact(Decimal('100'))
+        Decimal("100")
+        >>> ExtendedContext.to_integral_exact(Decimal('100.0'))
+        Decimal("100")
+        >>> ExtendedContext.to_integral_exact(Decimal('101.5'))
+        Decimal("102")
+        >>> ExtendedContext.to_integral_exact(Decimal('-101.5'))
+        Decimal("-102")
+        >>> ExtendedContext.to_integral_exact(Decimal('10E+5'))
+        Decimal("1.0E+6")
+        >>> ExtendedContext.to_integral_exact(Decimal('7.89E+77'))
+        Decimal("7.89E+77")
+        >>> ExtendedContext.to_integral_exact(Decimal('-Inf'))
+        Decimal("-Infinity")
+        """
+        return a.to_integral_exact(context=self)
+
+    def to_integral_value(self, a):
         """Rounds to an integer.
 
         When the operand has a negative exponent, the result is the same
@@ -2875,29 +4724,32 @@
         of the operand as the precision setting, except that no flags will
         be set.  The rounding mode is taken from the context.
 
-        >>> ExtendedContext.to_integral(Decimal('2.1'))
+        >>> ExtendedContext.to_integral_value(Decimal('2.1'))
         Decimal("2")
-        >>> ExtendedContext.to_integral(Decimal('100'))
+        >>> ExtendedContext.to_integral_value(Decimal('100'))
         Decimal("100")
-        >>> ExtendedContext.to_integral(Decimal('100.0'))
+        >>> ExtendedContext.to_integral_value(Decimal('100.0'))
         Decimal("100")
-        >>> ExtendedContext.to_integral(Decimal('101.5'))
+        >>> ExtendedContext.to_integral_value(Decimal('101.5'))
         Decimal("102")
-        >>> ExtendedContext.to_integral(Decimal('-101.5'))
+        >>> ExtendedContext.to_integral_value(Decimal('-101.5'))
         Decimal("-102")
-        >>> ExtendedContext.to_integral(Decimal('10E+5'))
+        >>> ExtendedContext.to_integral_value(Decimal('10E+5'))
         Decimal("1.0E+6")
-        >>> ExtendedContext.to_integral(Decimal('7.89E+77'))
+        >>> ExtendedContext.to_integral_value(Decimal('7.89E+77'))
         Decimal("7.89E+77")
-        >>> ExtendedContext.to_integral(Decimal('-Inf'))
+        >>> ExtendedContext.to_integral_value(Decimal('-Inf'))
         Decimal("-Infinity")
         """
-        return a.to_integral(context=self)
+        return a.to_integral_value(context=self)
+
+    # the method name changed, but we provide also the old one, for compatibility
+    to_integral = to_integral_value
 
 class _WorkRep(object):
     __slots__ = ('sign','int','exp')
     # sign: 0 or 1
-    # int:  int or long
+    # int:  int
     # exp:  None, int, or string
 
     def __init__(self, value=None):
@@ -2930,66 +4782,338 @@
 
     Done during addition.
     """
-    # Yes, the exponent is a long, but the difference between exponents
-    # must be an int-- otherwise you'd get a big memory problem.
-    numdigits = int(op1.exp - op2.exp)
-    if numdigits < 0:
-        numdigits = -numdigits
+    if op1.exp < op2.exp:
         tmp = op2
         other = op1
     else:
         tmp = op1
         other = op2
 
-
-    if shouldround and numdigits > prec + 1:
-        # Big difference in exponents - check the adjusted exponents
+    # Let exp = min(tmp.exp - 1, tmp.adjusted() - precision - 1).
+    # Then adding 10**exp to tmp has the same effect (after rounding)
+    # as adding any positive quantity smaller than 10**exp; similarly
+    # for subtraction.  So if other is smaller than 10**exp we replace
+    # it with 10**exp.  This avoids tmp.exp - other.exp getting too large.
+    if shouldround:
         tmp_len = len(str(tmp.int))
         other_len = len(str(other.int))
-        if numdigits > (other_len + prec + 1 - tmp_len):
-            # If the difference in adjusted exps is > prec+1, we know
-            # other is insignificant, so might as well put a 1 after the
-            # precision (since this is only for addition).  Also stops
-            # use of massive longs.
-
-            extend = prec + 2 - tmp_len
-            if extend <= 0:
-                extend = 1
-            tmp.int *= 10 ** extend
-            tmp.exp -= extend
+        exp = tmp.exp + min(-1, tmp_len - prec - 2)
+        if other_len + other.exp - 1 < exp:
             other.int = 1
-            other.exp = tmp.exp
-            return op1, op2
+            other.exp = exp
 
-    tmp.int *= 10 ** numdigits
-    tmp.exp -= numdigits
+    tmp.int *= 10 ** (tmp.exp - other.exp)
+    tmp.exp = other.exp
     return op1, op2
 
-def _adjust_coefficients(op1, op2):
-    """Adjust op1, op2 so that op2.int * 10 > op1.int >= op2.int.
+##### Integer arithmetic functions used by ln, log10, exp and __pow__ #####
 
-    Returns the adjusted op1, op2 as well as the change in op1.exp-op2.exp.
-
-    Used on _WorkRep instances during division.
+# This function from Tim Peters was taken from here:
+# http://mail.python.org/pipermail/python-list/1999-July/007758.html
+# The correction being in the function definition is for speed, and
+# the whole function is not resolved with math.log because of avoiding
+# the use of floats.
+def _nbits(n, correction = {
+        '0': 4, '1': 3, '2': 2, '3': 2,
+        '4': 1, '5': 1, '6': 1, '7': 1,
+        '8': 0, '9': 0, 'a': 0, 'b': 0,
+        'c': 0, 'd': 0, 'e': 0, 'f': 0}):
+    """Number of bits in binary representation of the positive integer n,
+    or 0 if n == 0.
     """
-    adjust = 0
-    # If op1 is smaller, make it larger
-    while op2.int > op1.int:
-        op1.int *= 10
-        op1.exp -= 1
-        adjust += 1
+    if n < 0:
+        raise ValueError("The argument to _nbits should be nonnegative.")
+    hex_n = "%x" % n
+    return 4*len(hex_n) - correction[hex_n[0]]
 
-    # If op2 is too small, make it larger
-    while op1.int >= (10 * op2.int):
-        op2.int *= 10
-        op2.exp -= 1
-        adjust -= 1
+def _sqrt_nearest(n, a):
+    """Closest integer to the square root of the positive integer n.  a is
+    an initial approximation to the square root.  Any positive integer
+    will do for a, but the closer a is to the square root of n the
+    faster convergence will be.
 
-    return op1, op2, adjust
+    """
+    if n <= 0 or a <= 0:
+        raise ValueError("Both arguments to _sqrt_nearest should be positive.")
+
+    b=0
+    while a != b:
+        b, a = a, a--n//a>>1
+    return a
+
+def _rshift_nearest(x, shift):
+    """Given an integer x and a nonnegative integer shift, return closest
+    integer to x / 2**shift; use round-to-even in case of a tie.
+
+    """
+    b, q = 1 << shift, x >> shift
+    return q + (2*(x & (b-1)) + (q&1) > b)
+
+def _div_nearest(a, b):
+    """Closest integer to a/b, a and b positive integers; rounds to even
+    in the case of a tie.
+
+    """
+    q, r = divmod(a, b)
+    return q + (2*r + (q&1) > b)
+
+def _ilog(x, M, L = 8):
+    """Integer approximation to M*log(x/M), with absolute error boundable
+    in terms only of x/M.
+
+    Given positive integers x and M, return an integer approximation to
+    M * log(x/M).  For L = 8 and 0.1 <= x/M <= 10 the difference
+    between the approximation and the exact result is at most 22.  For
+    L = 8 and 1.0 <= x/M <= 10.0 the difference is at most 15.  In
+    both cases these are upper bounds on the error; it will usually be
+    much smaller."""
+
+    # The basic algorithm is the following: let log1p be the function
+    # log1p(x) = log(1+x).  Then log(x/M) = log1p((x-M)/M).  We use
+    # the reduction
+    #
+    #    log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
+    #
+    # repeatedly until the argument to log1p is small (< 2**-L in
+    # absolute value).  For small y we can use the Taylor series
+    # expansion
+    #
+    #    log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
+    #
+    # truncating at T such that y**T is small enough.  The whole
+    # computation is carried out in a form of fixed-point arithmetic,
+    # with a real number z being represented by an integer
+    # approximation to z*M.  To avoid loss of precision, the y below
+    # is actually an integer approximation to 2**R*y*M, where R is the
+    # number of reductions performed so far.
+
+    y = x-M
+    # argument reduction; R = number of reductions performed
+    R = 0
+    while (R <= L and abs(y) << L-R >= M or
+           R > L and abs(y) >> R-L >= M):
+        y = _div_nearest((M*y) << 1,
+                         M + _sqrt_nearest(M*(M+_rshift_nearest(y, R)), M))
+        R += 1
+
+    # Taylor series with T terms
+    T = -int(-10*len(str(M))//(3*L))
+    yshift = _rshift_nearest(y, R)
+    w = _div_nearest(M, T)
+    for k in range(T-1, 0, -1):
+        w = _div_nearest(M, k) - _div_nearest(yshift*w, M)
+
+    return _div_nearest(w*y, M)
+
+def _dlog10(c, e, p):
+    """Given integers c, e and p with c > 0, p >= 0, compute an integer
+    approximation to 10**p * log10(c*10**e), with an absolute error of
+    at most 1.  Assumes that c*10**e is not exactly 1."""
+
+    # increase precision by 2; compensate for this by dividing
+    # final result by 100
+    p += 2
+
+    # write c*10**e as d*10**f with either:
+    #   f >= 0 and 1 <= d <= 10, or
+    #   f <= 0 and 0.1 <= d <= 1.
+    # Thus for c*10**e close to 1, f = 0
+    l = len(str(c))
+    f = e+l - (e+l >= 1)
+
+    if p > 0:
+        M = 10**p
+        k = e+p-f
+        if k >= 0:
+            c *= 10**k
+        else:
+            c = _div_nearest(c, 10**-k)
+
+        log_d = _ilog(c, M) # error < 5 + 22 = 27
+        log_10 = _ilog(10*M, M) # error < 15
+        log_d = _div_nearest(log_d*M, log_10)
+        log_tenpower = f*M # exact
+    else:
+        log_d = 0  # error < 2.31
+        log_tenpower = div_nearest(f, 10**-p) # error < 0.5
+
+    return _div_nearest(log_tenpower+log_d, 100)
+
+def _dlog(c, e, p):
+    """Given integers c, e and p with c > 0, compute an integer
+    approximation to 10**p * log(c*10**e), with an absolute error of
+    at most 1.  Assumes that c*10**e is not exactly 1."""
+
+    # Increase precision by 2. The precision increase is compensated
+    # for at the end with a division by 100.
+    p += 2
+
+    # rewrite c*10**e as d*10**f with either f >= 0 and 1 <= d <= 10,
+    # or f <= 0 and 0.1 <= d <= 1.  Then we can compute 10**p * log(c*10**e)
+    # as 10**p * log(d) + 10**p*f * log(10).
+    l = len(str(c))
+    f = e+l - (e+l >= 1)
+
+    # compute approximation to 10**p*log(d), with error < 27
+    if p > 0:
+        k = e+p-f
+        if k >= 0:
+            c *= 10**k
+        else:
+            c = _div_nearest(c, 10**-k)  # error of <= 0.5 in c
+
+        # _ilog magnifies existing error in c by a factor of at most 10
+        log_d = _ilog(c, 10**p) # error < 5 + 22 = 27
+    else:
+        # p <= 0: just approximate the whole thing by 0; error < 2.31
+        log_d = 0
+
+    # compute approximation to 10**p*f*log(10), with error < 17
+    if f:
+        sign_f = [-1, 1][f > 0]
+        if p >= 0:
+            M = 10**p * abs(f)
+        else:
+            M = _div_nearest(abs(f), 10**-p) # M = 10**p*|f|, error <= 0.5
+
+        if M:
+            f_log_ten = sign_f*_ilog(10*M, M)   # M*log(10), error <= 1.2 + 15 < 17
+        else:
+            f_log_ten = 0
+    else:
+        f_log_ten = 0
+
+    # error in sum < 17+27 = 44; error after division < 0.44 + 0.5 < 1
+    return _div_nearest(f_log_ten + log_d, 100)
+
+def _iexp(x, M, L=8):
+    """Given integers x and M, M > 0, such that x/M is small in absolute
+    value, compute an integer approximation to M*exp(x/M).  For 0 <=
+    x/M <= 2.4, the absolute error in the result is bounded by 60 (and
+    is usually much smaller)."""
+
+    # Algorithm: to compute exp(z) for a real number z, first divide z
+    # by a suitable power R of 2 so that |z/2**R| < 2**-L.  Then
+    # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
+    # series
+    #
+    #     expm1(x) = x + x**2/2! + x**3/3! + ...
+    #
+    # Now use the identity
+    #
+    #     expm1(2x) = expm1(x)*(expm1(x)+2)
+    #
+    # R times to compute the sequence expm1(z/2**R),
+    # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).
+
+    # Find R such that x/2**R/M <= 2**-L
+    R = _nbits((x<<L)//M)
+
+    # Taylor series.  (2**L)**T > M
+    T = -int(-10*len(str(M))//(3*L))
+    y = _div_nearest(x, T)
+    Mshift = M<<R
+    for i in range(T-1, 0, -1):
+        y = _div_nearest(x*(Mshift + y), Mshift * i)
+
+    # Expansion
+    for k in range(R-1, -1, -1):
+        Mshift = M<<(k+2)
+        y = _div_nearest(y*(y+Mshift), Mshift)
+
+    return M+y
+
+def _dexp(c, e, p):
+    """Compute an approximation to exp(c*10**e), with p decimal places of
+    precision.
+
+    Returns d, f such that:
+
+      10**(p-1) <= d <= 10**p, and
+      (d-1)*10**f < exp(c*10**e) < (d+1)*10**f
+
+    In other words, d*10**f is an approximation to exp(c*10**e) with p
+    digits of precision, and with an error in d of at most 1.  This is
+    almost, but not quite, the same as the error being < 1ulp: when d
+    = 10**(p-1) the error could be up to 10 ulp."""
+
+    # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision
+    p += 2
+
+    # compute log10 with extra precision = adjusted exponent of c*10**e
+    extra = max(0, e + len(str(c)) - 1)
+    q = p + extra
+    log10 = _dlog(10, 0, q)  # error <= 1
+
+    # compute quotient c*10**e/(log10/10**q) = c*10**(e+q)/log10,
+    # rounding down
+    shift = e+q
+    if shift >= 0:
+        cshift = c*10**shift
+    else:
+        cshift = c//10**-shift
+    quot, rem = divmod(cshift, log10)
+
+    # reduce remainder back to original precision
+    rem = _div_nearest(rem, 10**extra)
+
+    # error in result of _iexp < 120;  error after division < 0.62
+    return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
+
+def _dpower(xc, xe, yc, ye, p):
+    """Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and
+    y = yc*10**ye, compute x**y.  Returns a pair of integers (c, e) such that:
+
+      10**(p-1) <= c <= 10**p, and
+      (c-1)*10**e < x**y < (c+1)*10**e
+
+    in other words, c*10**e is an approximation to x**y with p digits
+    of precision, and with an error in c of at most 1.  (This is
+    almost, but not quite, the same as the error being < 1ulp: when c
+    == 10**(p-1) we can only guarantee error < 10ulp.)
+
+    We assume that: x is positive and not equal to 1, and y is nonzero.
+    """
+
+    # Find b such that 10**(b-1) <= |y| <= 10**b
+    b = len(str(abs(yc))) + ye
+
+    # log(x) = lxc*10**(-p-b-1), to p+b+1 places after the decimal point
+    lxc = _dlog(xc, xe, p+b+1)
+
+    # compute product y*log(x) = yc*lxc*10**(-p-b-1+ye) = pc*10**(-p-1)
+    shift = ye-b
+    if shift >= 0:
+        pc = lxc*yc*10**shift
+    else:
+        pc = _div_nearest(lxc*yc, 10**-shift)
+
+    if pc == 0:
+        # we prefer a result that isn't exactly 1; this makes it
+        # easier to compute a correctly rounded result in __pow__
+        if ((len(str(xc)) + xe >= 1) == (yc > 0)): # if x**y > 1:
+            coeff, exp = 10**(p-1)+1, 1-p
+        else:
+            coeff, exp = 10**p-1, -p
+    else:
+        coeff, exp = _dexp(pc, -(p+1), p+1)
+        coeff = _div_nearest(coeff, 10)
+        exp += 1
+
+    return coeff, exp
+
+def _log10_lb(c, correction = {
+        '1': 100, '2': 70, '3': 53, '4': 40, '5': 31,
+        '6': 23, '7': 16, '8': 10, '9': 5}):
+    """Compute a lower bound for 100*log10(c) for a positive integer c."""
+    if c <= 0:
+        raise ValueError("The argument to _log10_lb should be nonnegative.")
+    str_c = str(c)
+    return 100*len(str_c) - correction[str_c[0]]
 
 ##### Helper Functions ####################################################
 
-def _convert_other(other):
+def _convert_other(other, raiseit=False):
     """Convert other to Decimal.
 
     Verifies that it's ok to use in an implicit construction.
@@ -2998,6 +5122,8 @@
         return other
     if isinstance(other, int):
         return Decimal(other)
+    if raiseit:
+        raise TypeError("Unable to convert %s to Decimal" % other)
     return NotImplemented
 
 _infinity_map = {
@@ -3085,12 +5211,16 @@
 # Reusable defaults
 Inf = Decimal('Inf')
 negInf = Decimal('-Inf')
+NaN = Decimal('NaN')
+Dec_0 = Decimal(0)
+Dec_p1 = Decimal(1)
+Dec_n1 = Decimal(-1)
+Dec_p2 = Decimal(2)
+Dec_n2 = Decimal(-2)
 
 # Infsign[sign] is infinity w/ that sign
 Infsign = (Inf, negInf)
 
-NaN = Decimal('NaN')
-
 
 ##### crud for parsing strings #############################################
 import re