(define-module (language python module decimal)
  #:use-module ((language python module collections) #:select (namedtuple))
  #:use-module ((language python module itertools) #:select (chain repeat))
  #:use-module ((language python module sys) #:select (maxsize hash_info))
  #:use-module (language python module)
  #:use-module ((language python module python) #:select
		(isinstance str float int tuple classmethod pow property
			    complex range reversed (format . str-format)))
  #:use-module (language python list)
  #:use-module (language python number)
  #:use-module (language python string)
  #:use-module (language python for)
  #:use-module (language python try)
  #:use-module (language python hash)
  #:use-module (language python dict)
  #:use-module (language python def)
  #:use-module (language python exceptions)
  #:use-module (language python bool)
  #:use-module (language python module)
  #:use-module (oop pf-objects)
  #:use-module (language python module re)
  #:use-module (ice-9 control)
  #:use-module ((ice-9 match) #:select ((match . ice:match)))
  #:export
  (        ;; Two major classes
   Decimal Context

	   ;; Named tuple representation
	   DecimalTuple

	   ;; Contexts
	   DefaultContext BasicContext ExtendedContext
		    
	   ;; Exceptions
	   DecimalException Clamped InvalidOperation DivisionByZero
	   Inexact Rounded Subnormal Overflow Underflow
	   FloatOperation

	   ;; Exceptional conditions that trigger InvalidOperation
	   DivisionImpossible InvalidContext ConversionSyntax DivisionUndefined

	   ;; 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_05UP

	   ;; Functions for manipulating contexts
	   setcontext getcontext localcontext

	   ;; Limits for the C version for compatibility
	   MAX_PREC  MAX_EMAX MIN_EMIN MIN_ETINY))

(define-syntax-rule (aif it p . l) (let ((it p)) (if it . l)))

#|
This is the copyright information of the file ported over to scheme
# Copyright (c) 2004 Python Software Foundation.
# All rights reserved.

# Written by Eric Price <eprice at tjhsst.edu>
#    and Facundo Batista <facundo at taniquetil.com.ar>
#    and Raymond Hettinger <python at rcn.com>
#    and Aahz <aahz at pobox.com>
#    and Tim Peters

# This module should be kept in sync with the latest updates of the
# IBM specification as it evolves.  Those updates will be treated
# as bug fixes (deviation from the spec is a compatibility, usability
# bug) and will be backported.  At this point the spec is stabilizing
# and the updates are becoming fewer, smaller, and less significant.
|#

(define guile:modulo (@ (guile) modulo))

(define __name__   "decimal")
(define __xname__  __name__)
(define __version__ "1.70")
;; Highest version of the spec this complies with
;; See http://speleotrove.com/decimal/

(define DecimalTuple (namedtuple "DecimalTuple" "sign,digits,exponent"))

;; Rounding
(define ROUND_DOWN      'ROUND_DOWN)
(define ROUND_HALF_UP   'ROUND_HALF_UP)
(define ROUND_HALF_EVEN 'ROUND_HALF_EVEN)
(define ROUND_CEILING   'ROUND_CEILING)
(define ROUND_FLOOR     'ROUND_FLOOR)
(define ROUND_UP        'ROUND_UP)
(define ROUND_HALF_DOWN 'ROUND_HALF_DOWN)
(define ROUND_05UP      'ROUND_05UP)

;; Compatibility with the C version
(define MAX_PREC 425000000)
(define MAX_EMAX 425000000)
(define MIN_EMIN -425000000)

(if (= maxsize  (- (ash 1 63) 1))
    (begin
      (set! MAX_PREC 999999999999999999)
      (set! MAX_EMAX 999999999999999999)
      (set! MIN_EMIN -999999999999999999)))

(define MIN_ETINY  (- MIN_EMIN (- MAX_PREC 1)))

;; Context
(define-inlinable (cx-prec     x) (rawref x 'prec))
(define-inlinable (cx-emax     x) (rawref x 'Emax))
(define-inlinable (cx-emin     x) (rawref x 'Emin))
(define-inlinable (cx-etiny    x) ((ref x 'Etiny)))
(define-inlinable (cx-etop     x) ((ref x 'Etop)))
(define-inlinable (cx-copy     x) ((ref x 'copy)))
(define-inlinable (cx-clear_flags x) ((ref x 'clear_flags)))
(define-inlinable (cx-raise    x) (ref x '_raise_error))
(define-inlinable (cx-error    x) (ref x '_raise_error))
(define-inlinable (cx-capitals x) (rawref x 'capitals))
(define-inlinable (cx-cap      x) (rawref x 'capitals))
(define-inlinable (cx-rounding x) (rawref x 'rounding))
(define-inlinable (cx-clamp    x) (rawref x 'clamp))
(define-inlinable (cx-traps    x) (rawref x 'traps))
(define-inlinable (cx-flags    x) (rawref x 'flags))

;; Errors

(define-python-class DecimalException (ArithmeticError)
    "Base exception class.

    Used exceptions derive from this.
    If an exception derives from another exception besides this (such as
    Underflow (Inexact, Rounded, Subnormal) that indicates that it is only
    called if the others are present.  This isn't actually used for
    anything, though.

    handle  -- Called when context._raise_error is called and the
               trap_enabler is not set.  First argument is self, second is the
               context.  More arguments can be given, those being after
               the explanation in _raise_error (For example,
               context._raise_error(NewError, '(-x)!', self._sign) would
               call NewError().handle(context, self._sign).)

    To define a new exception, it should be sufficient to have it derive
    from DecimalException.
    "
    
    (define handle
      (lambda (self context . args)
        (values))))


(define-python-class Clamped (DecimalException)
    "Exponent of a 0 changed to fit bounds.

    This occurs and signals clamped if the exponent of a result has been
    altered in order to fit the constraints of a specific concrete
    representation.  This may occur when the exponent of a zero result would
    be outside the bounds of a representation, or when a large normal
    number would have an encoded exponent that cannot be represented.  In
    this latter case, the exponent is reduced to fit and the corresponding
    number of zero digits are appended to the coefficient ('fold-down').
    ")

(define-python-class InvalidOperation (DecimalException)
    "An invalid operation was performed.

    Various bad things cause this:

    Something creates a signaling NaN
    -INF + INF
    0 * (+-)INF
    (+-)INF / (+-)INF
    x % 0
    (+-)INF % x
    x._rescale( non-integer )
    sqrt(-x) , x > 0
    0 ** 0
    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.
    "
    (define handle
      (lambda (self context . args)
        (if (bool args)
	    (let ((ans  (_dec_from_triple
			 (ref (car args) '_sign)
			 (ref (car args) '_int)
			 "n" #t)))
	      ((ref ans '_fix_nan) context))
	    _NaN))))

(define-python-class ConversionSyntax (InvalidOperation)
    "Trying to convert badly formed string.

    This occurs and signals invalid-operation if a string is being
    converted to a number and it does not conform to the numeric string
    syntax.  The result is [0,qNaN].
    "
    (define handle
      (lambda x _NaN)))

(define-python-class DivisionByZero (DecimalException ZeroDivisionError)
    "Division by 0.

    This occurs and signals division-by-zero if division of a finite number
    by zero was attempted (during a divide-integer or divide operation, or a
    power operation with negative right-hand operand), and the dividend was
    not zero.

    The result of the operation is [sign,inf], where sign is the exclusive
    or of the signs of the operands for divide, or is 1 for an odd power of
    -0, for power.
    "

    (define handle
      (lambda (self context sign . args)
        (pylist-ref _SignedInfinity sign))))

(define-python-class DivisionImpossible (InvalidOperation)
    "Cannot perform the division adequately.

    This occurs and signals invalid-operation if the integer result of a
    divide-integer or remainder operation had too many digits (would be
    longer than precision).  The result is [0,qNaN].
    "

    (define handle
      (lambda x _NaN)))

(define-python-class DivisionUndefined (InvalidOperation ZeroDivisionError)
    "Undefined result of division.

    This occurs and signals invalid-operation if division by zero was
    attempted (during a divide-integer, divide, or remainder operation), and
    the dividend is also zero.  The result is [0,qNaN].
    "

    (define handle
      (lambda x _NaN)))

(define-python-class Inexact (DecimalException)
    "Had to round, losing information.

    This occurs and signals inexact whenever the result of an operation is
    not exact (that is, it needed to be rounded and any discarded digits
    were non-zero), or if an overflow or underflow condition occurs.  The
    result in all cases is unchanged.

    The inexact signal may be tested (or trapped) to determine if a given
    operation (or sequence of operations) was inexact.
    ")

(define-python-class InvalidContext (InvalidOperation)
    "Invalid context.  Unknown rounding, for example.

    This occurs and signals invalid-operation if an invalid context was
    detected during an operation.  This can occur if contexts are not checked
    on creation and either the precision exceeds the capability of the
    underlying concrete representation or an unknown or unsupported rounding
    was specified.  These aspects of the context need only be checked when
    the values are required to be used.  The result is [0,qNaN].
    "

    (define handle
      (lambda x _NaN)))

(define-python-class Rounded (DecimalException)
    "Number got rounded (not  necessarily changed during rounding).

    This occurs and signals rounded whenever the result of an operation is
    rounded (that is, some zero or non-zero digits were discarded from the
    coefficient), or if an overflow or underflow condition occurs.  The
    result in all cases is unchanged.

    The rounded signal may be tested (or trapped) to determine if a given
    operation (or sequence of operations) caused a loss of precision.
    ")

(define-python-class Subnormal (DecimalException)
    "Exponent < Emin before rounding.

    This occurs and signals subnormal whenever the result of a conversion or
    operation is subnormal (that is, its adjusted exponent is less than
    Emin, before any rounding).  The result in all cases is unchanged.

    The subnormal signal may be tested (or trapped) to determine if a given
    or operation (or sequence of operations) yielded a subnormal result.
    ")

(define-python-class Overflow (Inexact Rounded)
    "Numerical overflow.

    This occurs and signals overflow if the adjusted exponent of a result
    (from a conversion or from an operation that is not an attempt to divide
    by zero), after rounding, would be greater than the largest value that
    can be handled by the implementation (the value Emax).

    The result depends on the rounding mode:

    For round-half-up and round-half-even (and for round-half-down and
    round-up, if implemented), the result of the operation is [sign,inf],
    where sign is the sign of the intermediate result.  For round-down, the
    result is the largest finite number that can be represented in the
    current precision, with the sign of the intermediate result.  For
    round-ceiling, the result is the same as for round-down if the sign of
    the intermediate result is 1, or is [0,inf] otherwise.  For round-floor,
    the result is the same as for round-down if the sign of the intermediate
    result is 0, or is [1,inf] otherwise.  In all cases, Inexact and Rounded
    will also be raised.
    "

    (define handle
      (let ((l (list ROUND_HALF_UP ROUND_HALF_EVEN
		     ROUND_HALF_DOWN ROUND_UP)))
	(lambda (self context sign . args)
	  (let/ec ret
	    (if (memq (ref context 'rounding) l)
		(ret (pylist-ref _SignedInfinity sign)))
            
	    (if (= sign 0)
		(if (eq? (ref context 'rounding) ROUND_CEILING)
		    (ret (pylist-ref _SignedInfinity sign))
		    (ret (_dec_from_triple
			  sign
			  (* "9" (cx-prec context))
			  (+ (- (cx-emax context) (cx-prec context)) 1)))))

	    (if (= sign 1)
		(if (eq? (ref context 'rounding) ROUND_FLOOR)
		    (ret (pylist-ref _SignedInfinity sign))
		    (ret (_dec_from_triple
			  sign
			  (* "9" (cx-prec context))
			  (+ (- (cx-emax context) (cx-prec context)) 1))))))))))


(define-python-class Underflow (Inexact Rounded Subnormal)
    "Numerical underflow with result rounded to 0.

    This occurs and signals underflow if a result is inexact and the
    adjusted exponent of the result would be smaller (more negative) than
    the smallest value that can be handled by the implementation (the value
    Emin).  That is, the result is both inexact and subnormal.

    The result after an underflow will be a subnormal number rounded, if
    necessary, so that its exponent is not less than Etiny.  This may result
    in 0 with the sign of the intermediate result and an exponent of Etiny.

    In all cases, Inexact, Rounded, and Subnormal will also be raised.
    ")

(define-python-class FloatOperation (DecimalException TypeError)
    "Enable stricter semantics for mixing floats and Decimals.

    If the signal is not trapped (default), mixing floats and Decimals is
    permitted in the Decimal() constructor, context.create_decimal() and
    all comparison operators. Both conversion and comparisons are exact.
    Any occurrence of a mixed operation is silently recorded by setting
    FloatOperation in the context flags.  Explicit conversions with
    Decimal.from_float() or context.create_decimal_from_float() do not
    set the flag.

    Otherwise (the signal is trapped), only equality comparisons and explicit
    conversions are silent. All other mixed operations raise FloatOperation.
    ")

;; List of public traps and flags
(define _signals
  (list Clamped DivisionByZero Inexact Overflow Rounded
	Underflow InvalidOperation Subnormal FloatOperation))

;; Map conditions (per the spec) to signals
(define _condition_map
  `((,ConversionSyntax   . ,InvalidOperation)
    (,DivisionImpossible . ,InvalidOperation)
    (,DivisionUndefined  . ,InvalidOperation)
    (,InvalidContext     . ,InvalidOperation)))

;; Valid rounding modes
(define _rounding_modes
  (list ROUND_DOWN ROUND_HALF_UP ROUND_HALF_EVEN ROUND_CEILING
	ROUND_FLOOR ROUND_UP ROUND_HALF_DOWN ROUND_05UP))

;; ##### Context Functions ##################################################

;; The getcontext() and setcontext() function manage access to a thread-local
;; current context.
(define *context* (make-fluid #f))
(define (getcontext)
  (fluid-ref *context*))
(define (setcontext context)
  (fluid-set! *context* context))

;; ##### Decimal class #######################################################

;; Do not subclass Decimal from numbers.Real and do not register it as such
;; (because Decimals are not interoperable with floats).  See the notes in
;; numbers.py for more detail.

(define _dec_from_triple
  (lam (sign coefficient exponent (= special #f))
    "Create a decimal instance directly, without any validation,
    normalization (e.g. removal of leading zeros) or argument
    conversion.

    This function is for *internal use only*.
    "
    (Decimal sign coefficient exponent special)))

(define (get-parsed-sign m)
  (if (equal? ((ref m 'group) "sign") "-")
      1
      0))

(define (get-parsed-int m)
  ((ref m 'group) "int"))

(define (get-parsed-frac m)
  ((ref m 'group) "frac"))

(define (get-parsed-exp m)
  ((ref m 'group) "exp"))

(define (get-parsed-diag m)
  ((ref m 'group) "diag"))

(define (get-parsed-sig m)
  ((ref m 'group) "signal"))

(def (_mk self __init__ (= value "0") (= context None))
       "Create a decimal point instance.

        >>> Decimal('3.14')              # string input
        Decimal('3.14')
        >>> Decimal((0, (3, 1, 4), -2))  # tuple (sign, digit_tuple, exponent)
        Decimal('3.14')
        >>> Decimal(314)                 # int
        Decimal('314')
        >>> Decimal(Decimal(314))        # another decimal instance
        Decimal('314')
        >>> Decimal('  3.14  \\n')        # leading and trailing whitespace okay
        Decimal('3.14')
        "

        ;; Note that the coefficient, self._int, is actually stored as
        ;; a string rather than as a tuple of digits.  This speeds up
        ;; the "digits to integer" and "integer to digits" conversions
        ;; that are used in almost every arithmetic operation on
        ;; Decimals.  This is an internal detail: the as_tuple function
        ;; and the Decimal constructor still deal with tuples of
        ;; digits.

        ;; From a string
        ;; REs insist on real strings, so we can too.
        (cond
	 ((isinstance value str)
	    (let ((m (_parser (scm-str str))))
	      (if (not m)
                (let ((context (if (eq? context None)
				   (getcontext)
				   context)))
		  ((cx-raise context)
		   ConversionSyntax
		   (+ "Invalid literal for Decimal: " value))))

	      (let ((sign     (get-parsed-sign m))
		    (intpart  (get-parsed-int  m))
		    (fracpart (get-parsed-frac m))
		    (exp      (get-parsed-exp  m))
		    (diag     (get-parsed-diag m))
		    (signal   (get-parsed-sig  m)))
		
		(set self 'sign sign)
		
		(if (not (eq? intpart None))
		    (begin
		      (set self '_int (str (int (+ intpart fracpart))))
		      (set self '_exp (- exp (len fracpart)))
		      (set self '_is_special #f))
		    (begin
		      (if (not (eq? diag None))
			  (begin
			    ;; NaN
			    (set self '_int
				 (py-lstrip (str (int (if (bool diag)
							  diag
							  "0")))
					    "0"))
			    (if signal
				(set self '_exp "N")
				(set self '_exp "n")))
			  (begin
			    ;; infinity
			    (set self '_int "0")
			    (set self '_exp "F")))
		      (set self '_is_special #t))))))
	 
	;; From an integer	   
        ((isinstance value int)
	 (if (>= value 0)
	     (set self '_sign 0)
	     (set self '_sign 1))
	 (set self '_exp 0)
         (set self '_int (str (abs value)))
	 (set self '_is_special #f))
	
        ;; From another decimal
        ((isinstance value Decimal)
            (set self '_exp        (ref value '_exp       ))
            (set self '_sign       (ref value '_sign      ))
            (set self '_int        (ref value '_int       ))
            (set self '_is_special (ref value '_is_special)))

        ;; From an internal working value
        ((isinstance value _WorkRep)
	 (set self '_exp        (int (ref value '_exp)))
	 (set self '_sign       (ref value '_sign))
	 (set self '_int        (str (ref value 'int)))
	 (set self '_is_special #f))

        ;; tuple/list conversion (possibly from as_tuple())
        ((isinstance value (list list tuple))
            (if (not (= (len value) 3))
                (raise (ValueError
			(+ "Invalid tuple size in creation of Decimal "
			   "from list or tuple.  The list or tuple "
			   "should have exactly three elements."))))
	    ;; # process sign.  The isinstance test rejects floats
	    (let ((v0 (pylist-ref value 0))
		  (v1 (pylist-ref value 1))
		  (v2 (pylist-ref value 2)))
	      (if (not (and (isinstance v0 int)
			    (or (= v0 0) (= v0 1))))
		  (raise (ValueError
			  (+ "Invalid sign.  The first value in the tuple "
			     "should be an integer; either 0 for a "
			     "positive number or 1 for a negative number."))))
	      (set self '_sign v0)
            (if (eq? v2 'F)
		(begin
		  (set self '_int "0")
		  (set self '_exp v2)
		  (set self 'is_special #t))
		(let ((digits (py-list)))
		  ;; process and validate the digits in value[1]
		  (for ((digit : v1)) ()
		       (if (and (isinstance digit int)
				(<= 0 digit)
				(<= digit 9))
			   ;; skip leading zeros
			   (if (or (bool digits) (> digit 0))
			       (pylist-append! digits digit))
			   (raise (ValueError
				   (+ "The second value in the tuple must "
				      "be composed of integers in the range "
				      "0 through 9.")))))
		  
		  (cond
		   ((or (eq? v2 'n) (eq? v2 'N))
		    (begin
		      ;; NaN: digits form the diagnostic
		      (set self '_int  (py-join "" (map str digits)))
		      (set self '_exp  v2)
			(set self '_is_special #t)))
		   ((isinstance v2  int)
		    ;; finite number: digits give the coefficient
		    (set self '_int (py-join "" (map str digits)))
		    (set self '_exp v2)
		    (set self '_is_special #f))
		   (else
		    (raise (ValueError
			    (+ "The third value in the tuple must "
			       "be an integer, or one of the "
			       "strings 'F', 'n', 'N'.")))))))))

	((isinstance value float)
	 (let ((context (if (eq? context None)
			    (getcontext)
			    context)))
	   ((cx-error context)
	    FloatOperation
	    (+ "strict semantics for mixing floats and Decimals are "
	       "enabled"))
	   
	   (__init__ self ((ref Decimal 'from_float) value))))

	(else
	 (raise (TypeError
		 (format #f "Cannot convert ~a to Decimal" value))))))

(define-inlinable (divmod x y)
  (values (quotient x y) (modulo x y)))

(define-syntax twix
  (syntax-rules (when let let* if)
    ((_ a) a)
    ((_ (let () a ...) . l)
     (begin a ... (twix . l)))

    ((_ (let ((a aa) ...) b ...) . l)
     (let ((a aa) ...) b ... (twix . l)))
    ((_ (let (a ...)) . l)
     (a ... (twix . l)))
    ((_ (let* (a ...) b ...) . l)
     (let* (a ...) b ... (twix . l)))
    ((_ (if . u) . l)
     (begin (if . u) (twix . l)))
    ((_ (when . u) . l)
     (begin (when . u) (twix . l)))
    ((_ (a it code ...) . l)
     (aif it a (begin code ...) (twix . l)))))

(define-syntax-rule (norm-op self op)
  (begin
    (set! op ((ref self '_convert_other) op))
    (if (eq? op NotImplemented)
	op
	#f)))

(define-syntax-rule (get-context context code)
  (let ((context (if (eq? context None)
		     (getcontext)
		     context)))
    code))

(define-syntax-rule (un-special self context)
  (if ((ref self '_is_special))
      (let ((ans ((ref self '_check_nans) #:context context)))
	(if (bool ans)
	    ans
	    #f))
      #f))

(define-syntax-rule (bin-special o1 o2 context)
  (if (or (ref o1  '_is_special)
	  (ref o2  '_is_special))
      (or (un-special o1 context) (un-special o2 context))))

(define-syntax-rule (add-special self other context)
  (or (bin-special self other context)
      (if ((ref self '_isinfinity))
	  ;; If both INF, same sign =>
	  ;; same as both, opposite => error.
	  (if (and (not (= (ref self '_sign) (ref other '_sign)))
		   ((ref other '_isinfinity)))
	      ((cx-error context) InvalidOperation "-INF + INF")
	      (Decimal self))	      
	  (if ((ref other '_isinfinity))
	      (Decimal other) ; Can't both be infinity here
	      #f))))

(define-syntax-rule (mul-special self other context resultsign)
  (if (or (ref self '_is_special) (ref other '_is_special))
      (twix
       ((bin-special self other context) it it)
	       
       ((if ((ref self '_isinfinity))
	    (if (not (bool other))
		((cx-error context) InvalidOperation "(+-)INF * 0")
		(pylist-ref _SignedInfinity resultsign))
	    #f) it it)

       (if ((ref other '_isinfinity))
	   (if (not (bool self))
	       ((cx-error context) InvalidOperation "(+-)INF * 0")
	       (pylist-ref _SignedInfinity resultsign))
	   #f))
      #f))

(define-syntax-rule (div-special self other context sign)
  (if (or (ref self '_is_special) (ref other '_is_special))
      (twix
       ((bin-special self other context) it it)
       
       ((and ((ref self '_isinfinity)) ((ref other '_isinfinity))) it
	((cx-error context) InvalidOperation "(+-)INF/(+-)INF"))

       (((ref self '_isinfinity)) it
	(pylist-ref _SignedInfinity sign))

       (((ref other '_isinfinity)) it
	((cx-error context) Clamped "Division by infinity")
	(_dec_from_triple sign  "0" (cx-etiny context))))))
        

(define-python-class Decimal (object)
    "Floating point class for decimal arithmetic."
    
    
    ;; Generally, the value of the Decimal instance is given by
    ;;  (-1)**_sign * _int * 10**_exp
    ;; Special values are signified by _is_special == True

    (define __init__
      (case-lambda
       ((self sign coefficient exponent special)
	(set self '_sign       sign)
	(set self '_int        coefficient)
	(set self '_exp        exponent)
	(set self '_is_special special))
       
       ((self)
	(_mk self __init__))
       ((self a)
	(_mk self __init__ a))
       ((self a b)
	(_mk self __init__ a b))))

    (define from_float
      (classmethod
       (lambda (cls f)
        "Converts a float to a decimal number, exactly.

        Note that Decimal.from_float(0.1) is not the same as Decimal('0.1').
        Since 0.1 is not exactly representable in binary floating point, the
        value is stored as the nearest representable value which is
        0x1.999999999999ap-4.  The exact equivalent of the value in decimal
        is 0.1000000000000000055511151231257827021181583404541015625.

        >>> Decimal.from_float(0.1)
        Decimal('0.1000000000000000055511151231257827021181583404541015625')
        >>> Decimal.from_float(float('nan'))
        Decimal('NaN')
        >>> Decimal.from_float(float('inf'))
        Decimal('Infinity')
        >>> Decimal.from_float(-float('inf'))
        Decimal('-Infinity')
        >>> Decimal.from_float(-0.0)
        Decimal('-0')

        "
	(define (frexp x)
	  (if (< x 0) (set! x (- x)))
	  
	  (let lp ((l (string->list (format #f "~e" x))) (r1 '()))
	    (ice:match l
	       ((#\. . l)
		(let lp ((l l) (r2 '()))
		  (ice:match l
		    ((#\E . l)
		     (let* ((n (length r1))
			    (m (list->string (append (reverse r1) (reverse r2))))
			    (e (+ (- n 1) (string->number (list->string l)))))
		       (cons m e)))
		    ((x . l)
		     (lp l (cons x r2))))))
	       
	       ((x . l)
		(lp l (cons x r1))))))
		
		    
        (cond
	 ((isinstance f  int)                ; handle integer inputs
	  (cls f))
	 ((not (isinstance f float))
	  (raise (TypeError "argument must be int or float.")))
	 ((or (inf? f) (nan? f))
	  (cls (cond
		((nan? f)          "")
		((eq? f (inf))     "")
		((eq? f (- (inf))) ""))))
	 (else
	  (let* ((sign (if (>= f 0) 0 1))
		 (me   (frexp f))
		 (m    (car me))
		 (e    (cdr me))
		 (res  (_dec_from_triple sign m e)))
	    (if (eq? cls Decimal)
		res
		(cls res))))))))
    
    (define _isnan
      (lambda (self)
        "Returns whether the number is not actually one.

        0 if a number
        1 if NaN
        2 if sNaN
        "
        (if (ref self '_is_special)
            (let ((exp (ref self '_exp)))
	      (cond
	       ((eq? exp 'n) 1)
	       ((eq? exp 'N) 2)
	       (else         0)))
	    0)))
    
    (define _isinfinity 
      (lambda (self)
        "Returns whether the number is infinite

        0 if finite or not a number
        1 if +INF
        -1 if -INF
        "
        (if (eq? (ref self '_exp) 'F)
            (if (eq? (ref self '_sign) 1)
                -1
		1)
	    0)))
    
    (define _check_nans
      (lam (self (= other None) (= context None))
        "Returns whether the number is not actually one.

        if self, other are sNaN, signal
        if self, other are NaN return nan
        return 0

        Done before operations.
        "

        (let ((self_is_nan ((ref self '_isnan)))
	      (other_is_nan
	       (if (eq? other None)
		   #f
		   ((ref other '_isnan)))))
	  
	  (if (or self_is_nan other_is_nan)
	      (let ((context (if (eq? context None)
				 (getcontext)
				 context)))
		(cond
		 ((eq? self_is_nan 2)
		  ((cx-error context) InvalidOperation "sNaN" self))
		 ((eq? other_is_nan 2)
		  ((cx-error context) InvalidOperation "sNaN" other))
		 (self_is_nan
		  ((ref self '_fix_nan) context))
		 (else
		  ((ref other '_fix_nan) context))))
	      0))))
		 

    (define _compare_check_nans
      (lambda (self other context)
        "Version of _check_nans used for the signaling comparisons
        compare_signal, __le__, __lt__, __ge__, __gt__.

        Signal InvalidOperation if either self or other is a (quiet
        or signaling) NaN.  Signaling NaNs take precedence over quiet
        NaNs.

        Return 0 if neither operand is a NaN.

        "
        (let ((context (if (eq? context None)
			   (getcontext)
			   context)))

	  (if (or (ref self  '_is_special)
		  (ref other '_is_special))
	      (cond
	       (((ref self 'is_snan))
		((cx-error context)
		 InvalidOperation
		 "comparison involving sNaN" self))

 	       (((ref other 'is_snan))
		((cx-error context)
		 InvalidOperation
		 "comparison involving sNaN" other))
		
	       (((ref self 'is_qnan))
		((cx-error context)
		 InvalidOperation
		 "comparison involving NaN" self))

 	       (((ref other 'is_qnan))
		((cx-error context)
		 InvalidOperation
		 "comparison involving NaN" other))
	       
	       (else 0))
	      0))))

    (define __bool__
      (lambda (self)
        "Return True if self is nonzero; otherwise return False.

        NaNs and infinities are considered nonzero.
        "
        (or (ref self '_is_special) (not (equal? (ref self '_int) "0")))))

    (define _cmp
      (lambda (self other)
        "Compare the two non-NaN decimal instances self and other.

        Returns -1 if self < other, 0 if self == other and 1
        if self > other.  This routine is for internal use only."

	(let ((self_sign  (ref self  '_sign))
	      (other_sign (ref other '_sign)))
	  (cond
	   ((or (ref self '_is_special) (ref other '_is_special))
	    (let ((self_inf  ((ref self  '_isinfinity)))
		  (other_inf ((ref other '_isinfinity))))
	      (cond
	       ((eq? self_inf other_inf) 0)
	       ((< self_inf other_inf)  -1)
	       (else                     1))))
	    
	   ;; check for zeros;  Decimal('0') == Decimal('-0')
	   ((not (bool self))
	    (if (not (bool other))
		0
		(let ((s (ref other '_sign)))
		  (if (= s 0)
		      -1
		      1))))
	   ((not (bool other))
	    (let ((s (ref self '_sign)))
	      (if (= s 0)
		  1
		  -1)))
	  
	   ((< other_sign self_sign)
	    -1)
	   ((< self_sign other_sign)
	    1)

	   (else
	    (let ((self_adjusted  ((ref self  'adjusted)))
		  (other_adjusted ((ref other 'adjusted)))
		  (self_exp       (ref self  '_exp))
		  (other_exp      (ref other '_exp)))
	      (cond
	       ((= self_adjusted other_adjusted)
		(let ((self_padded  (+ (ref self '_int)
				       (* "0" (- self_exp  other_exp))))
		      (other_padded (+ (ref other '_int)
				       (* "0" (- other_exp  self_exp)))))
		  (cond
		   ((equal? self_padded other_padded)
		    0)
		   ((< self_padded other_padded)
		    (if (= self_sign 0)
			-1
			1))
		   (else
		    (if (= self_sign 0)
			1
			-1)))))
	       ((> self_adjusted other_adjusted)
		(if (= self_sign 0)
		    1
		    -1))
	       (else
		(if (= self_sign 0)
		    -1
		    1)))))))))
    
    ;; Note: The Decimal standard doesn't cover rich comparisons for
    ;; Decimals.  In particular, the specification is silent on the
    ;; subject of what should happen for a comparison involving a NaN.
    ;; We take the following approach:
    ;;
    ;;   == comparisons involving a quiet NaN always return False
    ;;   != comparisons involving a quiet NaN always return True
    ;;   == or != comparisons involving a signaling NaN signal
    ;;      InvalidOperation, and return False or True as above if the
    ;;      InvalidOperation is not trapped.
    ;;   <, >, <= and >= comparisons involving a (quiet or signaling)
    ;;      NaN signal InvalidOperation, and return False if the
    ;;      InvalidOperation is not trapped.
    ;;
    ;; This behavior is designed to conform as closely as possible to
    ;; that specified by IEEE 754.
    
    (define __eq__
      (lam (self other (= context None))
	   (let* ((so    (_convert_for_comparison self other #:equality_op #t))
		  (self  (car  so))
		  (other (cadr so)))
	     
	     (cond
	      ((eq? other NotImplemented)
	       other)
	      ((bool ((ref self '_check_nans) other context))
	       #f)
	      (else (= ((ref self '_cmp) other) 0))))))

    (define _xlt 
      (lambda (<)
	(lam (self other (= context None))
	     (let* ((so (_convert_for_comparison self other #:equality_op #t))
		    (self  (car  so))
		    (other (cadr so)))
	     
	       (cond
		((eq? other NotImplemented)
		 other)
		((bool ((ref self '_compare_check_nans) other context))
		 #f)
		(else (< ((ref self '_cmp) other) 0)))))))
      
    (define __lt__ (lambda x (apply (_xlt <  ) x)))
    (define __le__ (lambda x (apply (_xlt <= ) x)))
    (define __gt__ (lambda x (apply (_xlt >  ) x)))
    (define __ge__ (lambda x (apply (_xlt >= ) x)))

    (define compare 
      (lam (self other (= context None))
        "Compare self to other.  Return a decimal value:

        a or b is a NaN ==> Decimal('NaN')
        a < b           ==> Decimal('-1')
        a == b          ==> Decimal('0')
        a > b           ==> Decimal('1')
        "
        (let ((other (_convert_other other #:raiseit #t)))
	  ;; Compare(NaN, NaN) = NaN
	  (if (or (ref self '_is_special)
		  (and (bool other)
		       (ref other '_is_special)))
	      (aif it ((ref self '_check_nans) other context)
		   it
		   (Decimal ((ref self '_cmp) other)))))))

    (define __hash__
      (lambda (self)
	"x.__hash__() <==> hash(x)"

        ;; In order to make sure that the hash of a Decimal instance
        ;; agrees with the hash of a numerically equal integer, float
        ;; or Fraction, we follow the rules for numeric hashes outlined
        ;; in the documentation.  (See library docs, 'Built-in Types').
        (cond
	 ((ref self '_is_special)
	  (cond
	   (((ref self 'is_snan))
	    (raise (TypeError "Cannot hash a signaling NaN value.")))
	   (((ref self 'is_snan))
	    (hash (nan) pyhash-N))
	   ((= 1 (ref self '_sign))
	    (hash (- (inf)) pyhash-N))
	   (else
	    (hash (inf) pyhash-N))))
	    
	 (else
	  (let* ((exp (ref self '_exp))
		 (exp_hash
		  (if (>= exp 0)
		      (pow 10            exp     pyhash-N)
		      (pow _PyHASH_10INV (- exp) pyhash-N)))
		 
		 (hash_
		  (modulo (* (int (ref self '_int)) exp_hash)
			   pyhash-N))

		 (ans
		  (if (>= self 0) hash_ (- hash_))))
	    (if (= ans -1) -2 ans))))))
    
    (define as_tuple
      (lambda (self)
        "Represents the number as a triple tuple.

        To show the internals exactly as they are.
        "
	(DecimalTuple (ref self '_sign)
		      (tuple (map int (ref self '_int)))
		      (ref self '_exp))))

    (define as_integer_ratio
      (lambda (self)
        "Express a finite Decimal instance in the form n / d.

        Returns a pair (n, d) of integers.  When called on an infinity
        or NaN, raises OverflowError or ValueError respectively.

        >>> Decimal('3.14').as_integer_ratio()
        (157, 50)
        >>> Decimal('-123e5').as_integer_ratio()
        (-12300000, 1)
        >>> Decimal('0.00').as_integer_ratio()
        (0, 1)
        "
        (if (ref self '_is_special)
            (if ((ref self 'is_nan))
                (raise (ValueError
			"cannot convert NaN to integer ratio"))
                (raise (OverflowError
			"cannot convert Infinity to integer ratio"))))

        (if (not (bool self))
            (values 0 1)
	    (let* ((s (ref self '_sign))
		   (n (int (ref self '_int)))
		   (e (ref self '_exp))
		   (x
		    (* n (if (> exp 0)
			     (expt 10 exp)
			     (/ 1 (expt 10 (- expt)))))))
	      (values (numerator x)
		      (denominator x))))))
    
    (define __repr__
      (lambda (self)
        "Represents the number as an instance of Decimal."
        ;# Invariant:  eval(repr(d)) == d
        (format #f "Decimal('~a')" (str self))))

    (define __str__
      (lam (self  (= eng #f) (= context None))
        "Return string representation of the number in scientific notation.

        Captures all of the information in the underlying representation.
        "
	(let* ((sign         (if (= (ref self '_sign) 0) "" "-"))
	       (exp          (ref self '_exp))
	       (i            (ref self '_int))
	       (leftdigits   (+ exp (len i)))
	       (dotplace     #f)
	       (intpart      #f)
	       (fracpart     #f)
	       (exppart      #f))
	  
	  (cond
	   ((ref self '_is_special)
            (cond
	     ((eq? (ref self '_exp) 'F)
	      (+ sign "Infinity"))
	     ((eq? (ref self '_exp) 'n)
	      (+ sign  "NaN" (ref self '_int)))
	     (else ; self._exp == 'N'
	      (+ sign  "sNaN" (ref self '_int)))))
	   (else
	    ;; dotplace is number of digits of self._int to the left of the
	    ;; decimal point in the mantissa of the output string (that is,
	    ;; after adjusting the exponent)
	    (cond
	     ((and (<= exp 0) (> leftdigits  -6))
	      ;; no exponent required
	      (set! dotplace leftdigits))
	     
	     ((not eng)
	      ;; usual scientific notation: 1 digit on left of the point
	      (set! dotplace 1))
	     
	     ((equal? i "0")
	      ;; engineering notation, zero
	      (set! dotplace (- (modulo (+ leftdigits 1) 3) 1)))
	     (else
	      ;; engineering notation, nonzero
	      (set! dotplace (- (modulo (+ leftdigits 1) 3) 1))))

	    (cond
	     ((<= dotplace 0)
	      (set! intpart "0")
	      (set! fracpart (+ "." + (* "0" (- dotplace)) + i)))
	     ((>= dotplace (len i))
	      (set! intpart (+ i (* "0" (- dotplace (len i)))))
	      (set! fracpart ""))
	     (else
	      (set! intpart (pylist-slice i None dotplace None))
	      (set! fracpart (+ '.' (pylist-slice i dotplace None None)))))

	    
	    (cond
	     ((= leftdigits dotplace)
	      (set! exp ""))
	     (else
	      (let ((context (if (eq? context None)
				 (getcontext)
				 context)))
		(set! exp
		      (+ (pylist-ref '("e" "E") (cx-capitals context))
			 (format #f "~@d" (- leftdigits dotplace)))))))
	    (+ sign intpart fracpart exp))))))
    
    (define to_eng_string
      (lam (self (= context None))
        "Convert to a string, using engineering notation if an exponent is needed.
        Engineering notation has an exponent which is a multiple of 3.  This
        can leave up to 3 digits to the left of the decimal place and may
        require the addition of either one or two trailing zeros.
        "
        ((ref self '__str__) #:eng #t #:contect context)))

    (define __neg__
      (lam (self (= context None))
        "Returns a copy with the sign switched.

        Rounds, if it has reason.
        "
	(twix 
	  ((un-special self context) it it)
	  (let* ((context (if (eq? context None)
			      (getcontext)
			      context))
		 (ans     (if (and (not (bool self))
				   (not (eq? (cx-rounding context)
					     ROUND_FLOOR)))
			      ;; -Decimal('0') is Decimal('0'),
			      ;; not Decimal('-0'), except
			      ;; in ROUND_FLOOR rounding mode.
			      ((ref self 'copy_abs))
			      ((ref self 'copy_negate)))))
	    
	    ((ref ans '_fix) context)))))

    (define __pos__
      (lam (self (= context None))
        "Returns a copy, unless it is a sNaN.

        Rounds the number (if more than precision digits)
        "
	(twix
	 ((un-special self context) it it)
	 
	 (let* ((context (if (eq? context None)
			     (getcontext)
			     context))
		(ans     (if (and (not (bool self))
				  (not (eq? (cx-rounding context)
					    ROUND_FLOOR)))
			     ;; -Decimal('0') is Decimal('0'),
			     ;; not Decimal('-0'), except
			     ;; in ROUND_FLOOR rounding mode.
			     ((ref self 'copy_abs))
			     (Decimal self))))

	   ((ref ans '_fix) context)))))

    (define __abs__
      (lam (self  (= round #t) (= context None))
	"Returns the absolute value of self.

        If the keyword argument 'round' is false, do not round.  The
        expression self.__abs__(round=False) is equivalent to
        self.copy_abs().
        "
	(twix
	 ((not (bool round)) it
	  ((ref self 'copy_abs)))

	 ((un-special self context) it it)

	 (if (= (ref self '_sign) 1)
	     ((ref self '__neg__) #:context context)
	     ((ref self '__pos__) #:context context)))))
    

    (define __add__
      (lam (self other (= context None))
        "Returns self + other.

        -INF + INF (or the reverse) cause InvalidOperation errors.
        "
	(twix
	 ((norm-op self other) it it)
	 
	 (let (get-context context))
	 
	 ((add-special self other context) it it)
	 		 
	 (let* ((negativezero 0)
		(self_sign    (ref self  '_sign))
		(other_sign   (ref other '_sign))
		(self_exp     (ref self  '_sign))
		(other_exp    (ref other '_sign))
		(prec         (cx-prec context))
		(exp          (min self_exp other_exp))		    
		(sign         #f)
		(ans          #f))
	  
	   (if (and (eq? (cx-rounding context) ROUND_FLOOR)
		    (not (= self_sign other_sign)))
	       ;; If the answer is 0, the sign should be negative,
	       ;; in this case.
	       (set! negativezero 1)))

	 ((if (and (not (bool self)) (not (bool other)))
	      (begin
		(set! sign (min self_sign other_sign))
		(if (= negativezero 1)
		    (set! sign 1))
		(set! ans (_dec_from_triple sign "0" exp))
		(set! ans ((ref ans '_fix) context))
		ans)
	      #f) it it)

	 ((if (not (bool self))
	      (begin
		(set! exp (max exp (- other_exp prec 1)))
		(set! ans ((ref other '_rescale) exp
			   (cx-rounding context)))
		(set! ans ((ref ans '_fix) context))
		ans)
	      #f) it it)

	 ((if (not (bool other))
	      (begin
		(set! exp (max exp (- self_exp prec 1)))
		(set! ans ((ref self '_rescale) exp
			   (cx-rounding context)))
		(set! ans ((ref ans '_fix) context))
		ans)
	      #f) it it)
		
		
	 (let* ((op1    (_WorkRep self))
		(op2    (_WorkRep other))
		(ab     (_normalize op1 op2 prec))
		(op1_i  (car  ab))
		(op2_i  (cadr ab))
		(result (_WorkRep))))

	 ((cond
	   ((not (= (ref op1 'sign) (ref op2 'sign)))
	    ;; Equal and opposite
	    (twix
	     ((= op1_i op2_i) it
	      (set! ans (_dec_from_triple negativezero "0" exp))
	      (set! ans ((ref ans '_fix) context))
	      ans)
		    
	     (begin
	       (if (< op1_i op2_i)
		   (let ((t op1))
		     (set! op1 op2)
		     (set! op2 t)))
		    
	       (if (= (ref op1 'sign) 1)			
		   (let ((t (ref op1 'sign)))
		     (set result 'sign 1)
		     (set op1 'sign (ref op2 'sign))
		     (set op2 'sign t))
		   (set result 'sign 0))
	       #f)))
	   ((= (ref op1 'sign) 1)
	    (set result 'sign 1)
	    #f)

	   (else
	    (set result 'sign 0)
	    #f)) it it)

	 (begin
	   (if (= (ref op2 'sign) 0)
	       (set result 'int (+ (ref op1 'int) (ref op2 'int)))
	       (set result 'int (- (ref op1 'int) (ref op2 'int))))

	   (set result 'exp (ref op1 'exp))
	   (set! ans (Decimal result))
	   ((ref ans '_fix) context)))))

    (define __radd__ (lambda x (apply __add__ x)))

    (define __sub__
      (lam (self other (= context None))
	"Return self - other"
	(twix
	 ((norm-op self other)            it it)
	 ((bin-special self other context) it it)
	 ((ref self '__add__)
	  ((ref other 'copy_negate)) #:context context))))

    (define  __rsub__
      (lam (self other (= context None))
        "Return other - self"
    	(twix
	 ((norm-op self other) it it)
	 ((ref 'other '__sub__) self  #:context context))))

    (define __mul__
      (lam (self other (= context None))
        "Return self * other.

        (+-) INF * 0 (or its reverse) raise InvalidOperation.
        "
	(twix
	 ((norm-op self other) it it)	 
	 (let (get-context context))
	 
	 (let ((resultsign (logxor (ref self  '_sign)
				   (ref other '_sign)))))

	 ((mul-special self other context resultsign) it it)

	 (let ((resultexp (+ (ref self '_exp) (ref other '_exp)))))

	 ;; Special case for multiplying by zero
	 ((or (not (bool self)) (not (bool other))) it
	  (let ((ans (_dec_from_triple resultsign "0" resultexp)))
	    ((ref ans '_fix) context)))
	 
	 ;; Special case for multiplying by power of 10
	 ((equal? (ref self '_int) "1") it
	  (let ((ans (_dec_from_triple resultsign (ref other '_int) resultexp)))
	    ((ref ans '_fix) context)))

	 ((equal? (ref other '_int) "1") it
	  (let ((ans (_dec_from_triple resultsign (ref self '_int) resultexp)))
	    ((ref ans '_fix) context)))

	 (let* ((op1 (_WorkRep self))
		(op2 (_WorkRep other))
		(ans (_dec_from_triple resultsign
				       (str (* (ref op1 'int) (ref op2 'int)))
				       resultexp)))
	   ((ref ans '_fix) context)))))
  
    (define __rmul__ __mul__)

    (define __truediv__
      (lam (self other (= context None))
	"Return self / other."
	(twix
	 ((norm-op self other) it it)	 
	 (let (get-context context))
   
	 (let ((sign (logxor (ref self  '_sign)
			     (ref other '_sign)))))

	 ((div-special self other context sign) it it)	 

	 ;; Special cases for zeroes
        ((if (not (bool other))
	     (if (not (bool self))
                 ((cx-error context) DivisionUndefined "0 / 0")
		 ((cx-error context) DivisionByZero    "x / 0" sign))
	     #f) it it)

	(let ((exp    #f)
	      (coeff  #f)
	      (prec   (cx-prec context))
	      (nself  (len (ref self '_int)))
	      (nother (len (ref other '_int))))
	  (if (not (bool self))
	      (begin
		(set! exp   (- (ref self '_exp) (ref other '_exp)))
		(set! coeff 0))
	      ;; OK, so neither = 0, INF or NaN
	      (let ((shift (+ nother (- nself) prec 1))
		    (op1   (_WorkRep self))
		    (op2   (_WorkRep other)))
		(set! exp (- (ref self '_exp) (ref other '_exp) shift))
		(call-with-values
		    (lambda ()
		      (if (>= shift 0)
			  (divmod (* (ref op1 'int) (expt 10 shift))
				  (ref op2 'int))
			  (divmod (ref op1 'int)
				  (* (ref op2 'int) (expt 10 shift)))))
		  (lambda (coeff- remainder)
		    (set! coeff
			  (if (not (= remainder 0))
			      ;; result is not exact adjust to ensure
			      ;; correct rounding
			      (if (= (modulo coeff- 5) 0)
				  (+ coeff- 1)
				  coeff)
			      (let ((ideal_exp (- (ref self '_exp)
						  (ref other '_exp))))
				(let lp ((coeff- coeff-) (exp- exp))
				  (if (and (< exp- ideal_exp)
					   (= (modulo coeff 10) 0))
				      (lp (/ coeff 10) (+ exp- 1))
				      (begin
					(set exp exp-)
					coeff))))))))))
			      
			    
	  (let ((ans (_dec_from_triple sign (str coeff) exp)))
	    ((ref ans '_fix) context))))))

    (define _divide
      (lambda (self other context)	
        "Return (self // other, self % other), to context.prec precision.

        Assumes that neither self nor other is a NaN, that self is not
        infinite and that other is nonzero.
        "
	(apply values
	(twix
	 (let (let ((sign
		     (logxor (ref self  '_sign)
			     (ref other '_sign)))
		    (ideal_exp
		     (if ((ref other '_isinfinity))
			 (ref self '_exp)
			 (min (ref self 'exp) (ref other '_exp))))
		    (expdiff
		     (- ((ref self 'adjusted)) ((ref other 'adjusted)))))))
	 
	 ((or (not (bool self))
	      ((ref other '_isinfinity))
	      (<= expdiff -1)) it
	   (list (_dec_from_triple sign "0" 0)
		 ((ref self '_rescale) ideal_exp (cx-rounding context))))
	 
	 ((if (<= expdiff  (cx-prec context))
	      (let ((op1 (_WorkRep self))
		    (op2 (_WorkRep other)))
		(if (>= (ref op1 'exp) (ref op2 'exp))
		    (set op1 'int (* (ref op1 'int)
				     (expt 10 (- (ref op1 'exp)
						 (ref op2 'exp)))))
		    (set op2 'int (* (ref op2 'int)
				     (expt 10 (- (ref op2 'exp)
						 (ref op1 'exp))))))
		(call-with-values (lambda () (divmod (ref op1 'int)
						     (ref op2 'int)))
		  (lambda (q r)		
		    (if (< q (expt 10 (cx-prec context)))
			(list (_dec_from_triple sign (str q) 0)
			      (_dec_from_triple (ref self '_sign)
						(str r)
						ideal_exp))
			#f))))
	      #f) it it)
	 (begin
	   ;; Here the quotient is too large to be representable
	   (let ((ans ((cx-raise context) DivisionImpossible
		       "quotient too large in //, % or divmod")))
	     (list ans ans)))))))

    (define __rtruediv__
      (lam (self other (= context None))
	"Swaps self/other and returns __truediv__."
	(twix
	 ((norm-op self other) it it)	 
	 ((ref other '__truediv__) self #:context context))))
    
    (define __divmod__
      (lam (self other (= context None))
        "
        Return (self // other, self % other)
        "
	(apply values
	(twix
	 ((norm-op self other) it it)	 

	 (let (get-context context))
	 
	 ((add-special self other context) it it)
	 
	 (((ref self '_check_nans) other context) it
	  (list it it))

	 (let (let ((sign
		     (logxor (ref self  '_sign)
			     (ref other '_sign))))))

	 (((ref self '_isinfinity)) it
	  (if ((ref other '_isinfinity))
	      (let ((ans ((cx-error context) InvalidOperation
			  "divmod(INF, INF)")))
		(list ans ans))
	      (list (list-ref _SignedInfinity sign)
		    ((cx-raise context) InvalidOperation "INF % x"))))

	 ((not (bool other)) it
	  (if (not (bool self))
	      (let ((ans ((cx-error context) DivisionUndefined
			  "divmod(0, 0)")))
		(list ans ans))
	      (list ((cx-error context) DivisionByZero "x // 0" sign)
		    ((cx-error context) InvalidOperation "x % 0"))))
	     
	 (call-with-values (lambda () ((ref self '_divide) other context))
	   (lambda (quotient remainder)
	     (let ((remainder ((ref remainder '_fix) context)))
	       (list quotient remainder))))))))

    (define __rdivmod__
      (lam (self other (= context None))
	"Swaps self/other and returns __divmod__."
    	(twix
	 ((norm-op self other) it it)	 
	 ((ref other '__divmod__) self #:context context))))

    (define __mod__
      (lam (self other (= context None))
        "
        self % other
        "
	(twix
	 ((norm-op self other) it it)	 

	 (let (get-context context))

	 ((bin-special self other context) it it)

	 (((ref self '_isinfinity)) it
	  ((cx-error context) InvalidOperation "INF % x"))
	 
	 ((not (bool other)) it
          (if (bool self)
	      ((cx-error context) InvalidOperation  "x % 0")          
	      ((cx-error context) DivisionUndefined "0 % 0")))

	 (let* ((remainder ((ref self '_divide) other context)))
	   ((ref remainder '_fix) context)))))

    (define __rmod__
      (lam (self other (= context None))
	"Swaps self/other and returns __mod__."
    	(twix
	 ((norm-op self other) it it)	 
	 ((ref other '__mod__) self #:context context))))

    (define remainder_near
      (lam (self other (= context None))
        "
        Remainder nearest to 0-  abs(remainder-near) <= other/2
        "
	(twix
	 ((norm-op self other) it it)

	 (let (get-context context))

	 ((bin-special self other context) it it)

	 ;; self == +/-infinity -> InvalidOperation
	 (((ref self '_isinfinity)) it
	  ((cx-error context) InvalidOperation "remainder_near(infinity, x)"))

	 ;; other == 0 -> either InvalidOperation or DivisionUndefined
	 ((not (bool other)) it
	  (if (not (bool self))
	      ((cx-error context) InvalidOperation  "remainder_near(x, 0)")
	      ((cx-error context) DivisionUndefined "remainder_near(0, 0)")))

	 ;; other = +/-infinity -> remainder = self
	 (((ref other '_isinfinity)) it
	  (let ((ans (Decimal self)))
            ((ref ans '_fix) context)))

	 ;; self = 0 -> remainder = self, with ideal exponent
	 (let (let ((ideal_exponent (min (ref self '_exp) (ref other '_exp))))))

	 ((not (bool self)) it
	  (let ((ans (_dec_from_triple (ref self '_sign) "0" ideal_exponent)))
            ((ref ans '_fix) context)))

	 ;; catch most cases of large or small quotient
	 (let (let ((expdiff
		     (- ((ref self 'adjusted)) ((ref other 'adjusted)))))))
	 
        ((>= expdiff (+ (cx-prec context) 1)) it
	 ;; expdiff >= prec+1 => abs(self/other) > 10**prec
	 ((cx-error context) DivisionImpossible))
	
        ((<= expdiff -2) it
	 ;; expdiff <= -2 => abs(self/other) < 0.1
	 (let ((ans ((ref self '_rescale)
		     ideal_exponent (cx-rounding context))))
	   ((ref ans '_fix) context)))

        (let ((op1  (_WorkRep self))
	      (op2  (_WorkRep other)))

	  ;; adjust both arguments to have the same exponent, then divide
	  (if (>= (ref op1 'exp) (ref op2 'exp))
	      (set op1 'int (* (ref op1 'int)
			       (expt 10 (- (ref op1 'exp) (ref op2 'exp)))))
	      (set op2 'int (* (ref op2 'int)
			       (expt 10 (- (ref op2 'exp) (ref op1 'exp))))))
	  
	  (call-with-values (lambda () (divmod (ref op1 'int) (ref op2 'int)))
	    (lambda (q r)

	      ;; 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)  + (logand q 1)) (ref op2 'int))
		  (set! r (- r (ref op2 'int)))
		  (set! q (+ q 1)))

	      (if (>= q (expt 10 (cx-prec context)))
		  ((cx-error context) DivisionImpossible)
		  (let ((sign (ref self '_sign)))
		    (if (< r 0)
			(set! sign (- 1 sign))
			(set! r    (- r)))
		    (let ((ans (_dec_from_triple sign (str r) ideal_exponent)))
		      ((ref ans '_fix) context))))))))))

    (define __floordiv__
      (lam (self other (= context None))
        "self // other"
	(twix
	 ((norm-op self other) it it)

	 (let (get-context context))

	 ((bin-special self other context) it it)

	 (((ref self '_isinfinity)) it
	  (if ((ref other '_isinfinity))
               ((cx-error context) InvalidOperation "INF // INF")
	       (pylist-ref _SignedInfinity (logxor (ref self  '_sign)
						   (ref other '_sign)))))

	 ((not (bool other)) it
	  (if (bool self)
	      ((cx-error context) DivisionByZero    "x // 0"
	       (logxor (ref self  '_sign) (ref other '_sign)))
	      ((cx-error context) DivisionUndefined "0 // 0")))

	 ((ref self '_divide) other context))))

    (define __rfloordiv__
      (lam (self other (= context None))
	"Swaps self/other and returns __floordiv__."
	(twix
	 ((norm-op self other) it it)	 
	 ((ref other '__floordiv__) self #:context context))))

    (define __float__
      (lambda (self)
        "Float representation."
        (if ((ref self '_isnan))
            (if ((ref self 'is_snan))
                (raise (ValueError "Cannot convert signaling NaN to float"))
		(if (= (ref self '_sign))
		    (- (nan))
		    (nan)))
	    (if ((ref self '_isspecial))
		(if (= (ref self '_sign))
		    (- (inf))
		    (inf))
		(float (str self))))))
    
    (define __int__
      (lambda (self)
        "Converts self to an int, truncating if necessary."
	(if ((ref self '_isnan))
	    (raise (ValueError "Cannot convert NaN to integer"))
	    (if ((ref self '_isspecial))
		(raise (OverflowError "Cannot convert infinity to integer"))
		(let ((s (if (= (ref self '_sign) 1) -1 1)))
		  (if (>= (ref self '_exp) 0)
		      (* s (int (ref self '_int)) (expt 10 (ref self '_exp)))
		      (* s (int (or (bool (pylist-slice (ref self '_int)
						    None (ref self '_exp) None))
				    "0")))))))))

    (define __trunc__ __int__)

    (define real
      (property (lambda (self) self)))
    
    (define imag
      (property
       (lambda (self)
	 (Decimal 0))))

    (define conjugate
      (lambda (self) self))
    
    (define __complex__
      (lambda (self)
        (complex (float self))))

    (define _fix_nan
      (lambda (self context)
        "Decapitate the payload of a NaN to fit the context"
        (let ((payload (ref self '_int))
	      
	      ;; maximum length of payload is precision if clamp=0,
	      ;; precision-1 if clamp=1.
	      (max_payload_len
	       (- (ref context 'prec)
		  (ref context 'clamp))))
	  
	  (if (> (len payload)  max_payload_len)
	      (let ((payload (py-lstrip
			      (pylist-slice payload
					    (- (len payload) max_payload_len)
					    None None)  "0")))
		(_dec_from_triple (ref self '_sign) payload (ref self '_exp)
				  #t))
	      (Decimal self)))))

    (define _fix 
      (lambda (self context)
        "Round if it is necessary to keep self within prec precision.

        Rounds and fixes the exponent.  Does not raise on a sNaN.

        Arguments:
        self - Decimal instance
        context - context used.
        "

	(twix
	 (((ref self '_is_special)) it
	  (if ((ref self '_isnan))
	      ;; decapitate payload if necessary
	      ((ref self '_fix_nan) context)

	      ;; self is +/-Infinity; return unaltered
	      (Decimal self)))

	 ;; if self is zero then exponent should be between Etiny and
	 ;; Emax if clamp==0, and between Etiny and Etop if clamp==1.
	 (let ((Etiny (cx-etiny context))
	       (Etop  (cx-etop  context))))
	 
	 ((not (bool self)) it
	  (let* ((exp_max (if (= (cx-clamp context) 0)
			      (cx-emax context)
			      Etop))
		 (new_exp (min (max (ref self '_exp) Etiny) exp_max)))
            (if (not (= new_exp (ref self '_exp)))
		(begin
		  ((cx-error context) Clamped)
		  (_dec_from_triple (ref self '_sign) "0" new_exp))
		(Decimal self))))

	 ;; exp_min is the smallest allowable exponent of the result,
	 ;; equal to max(self.adjusted()-context.prec+1, Etiny)
	 (let ((exp_min (+ (len (ref self '_int))
			   (ref self '_exp)
			   (- (cx-prec context))))))
	 ((> exp_min Etop) it
	  ;; overflow: exp_min > Etop iff self.adjusted() > Emax
	  (let ((ans ((cx-error context) Overflow "above Emax"
		      (ref self '_sign))))
            ((cx-error context) Inexact)
	    ((cx-error context) Rounded)
	    ans))

	 (let* ((self_is_subnormal (< exp_min Etiny))
		(exp_min           (if self_is_subnormal Etiny exp_min))))

	 ;; round if self has too many digits
	 ((< (ref self '_exp) exp_min) it
	  (let ((digits (+ (len (ref self '_int))
			   (ref self '_exp)
			   (- exp_min))))
            (if (< digits 0)
                (set! self (_dec_from_triple (ref self '_sign)
					     "1" (- exp_min 1)))
                (set! digits 0))
	    
            (let* ((ans #f)
		   (rounding_method (pylist-ref
				     (ref self '_pick_rounding_function)
				     (cx-rounding context)))
		   (changed (rounding_method self digits))
		   (coeff   (or (bool (pylist-slice (ref self '_int)
						    None digits None)) "0")))
	      (if (> changed  0)
		  (begin
		    (set! coeff (str (+ (int coeff) 1)))
		    (if (> (len coeff) (cx-prec context))
			(begin
			  (set! coeff (pylist-slice coeff None -1 None))
			  (set! exp_min (+ exp_min  1))))))

	      ;; check whether the rounding pushed the exponent out of range
	      (if (> exp_min  Etop)
		  (set! ans
			((cx-error context) Overflow "above Emax"
			 (ref self '_sign)))
		  (set! ans (_dec_from_triple (ref self '_sign) coeff exp_min)))

	      ;; raise the appropriate signals, taking care to respect
	      ;; the precedence described in the specification
	      (if (and changed self_is_subnormal)
		  ((cx-error context) Underflow))
	      (if self_is_subnormal
		  ((cx-error context) Subnormal))
	      (if changed
		  ((cx-error context) Inexact))

	      ((cx-error context) Rounded)

	      (if (not (bool ans))
		  ;; raise Clamped on underflow to 0
		  ((cx-error context) Clamped))

	      ans)))
	 (begin
	   (if self_is_subnormal
	       ((cx-error context) Subnormal))


	   ;; fold down if clamp == 1 and self has too few digits
	   (if (and (= (cx-clamp context) 1) (> (ref self '_exp) Etop))
	       (begin
		 ((cx-error context) Clamped)
		 (let ((self_padded  (+ (ref self '_int)
					(* "0"
					   (- (ref self '_exp) Etop)))))
		   (_dec_from_triple (ref self '_sign) self_padded Etop)))
	       
	       ;; here self was representable to begin with; return unchanged
	       (Decimal self))))))
    
    

    ;; for each of the rounding functions below:
    ;;   self is a finite, nonzero Decimal
    ;;   prec is an integer satisfying 0 <= prec < len(self._int)
    ;;
    ;; each function returns either -1, 0, or 1, as follows:
    ;;   1 indicates that self should be rounded up (away from zero)
    ;;   0 indicates that self should be truncated, and that all the
    ;;    digits to be truncated are zeros (so the value is unchanged)
    ;;  -1 indicates that there are nonzero digits to be truncated

    (define _round_down
      (lambda (self prec)
        "Also known as round-towards-0, truncate."
        (if (_all_zeros (ref self '_int) prec)
            0
	    -1)))

    (define _round_up
      (lambda (self prec)
        "Rounds away from 0."
        (- (_round_down self prec))))

    (define _round_half_up
      (lambda (self prec)
        "Rounds 5 up (away from 0)"
        (cond
	 ((in (pylist-ref (ref self '_int) prec) "56789")
	  1)
	 ((_all_zeros (ref self '_int) prec)
	  0)
	 (else -1))))
    
    (define _round_half_down
      (lambda (self prec)
        "Round 5 down"
        (if (_exact_half (ref self '_int) prec)
            -1
            (_round_half_up self prec))))

    (define _round_half_even
      (lambda (self prec)
        "Round 5 to even, rest to nearest."
        (if (and (_exact_half (ref self '_int) prec)
		 (or (= prec 0)
		     (in (pylist-ref (ref self '_int) (- prec 1)) "02468")))
            -1)
       (_round_half_up self prec)))

    (define _round_ceiling
      (lambda (self prec)
        "Rounds up (not away from 0 if negative.)"
        (if (= (ref self '_sign) 1)
            (_round_down self prec)
	    (- (_round_down self prec)))))

    (define _round_floor
      (lambda (self prec)
        "Rounds down (not towards 0 if negative)"
        (if (= (ref self '_sign) 1)
	    (- (_round_down self prec))
	    (_round_down self prec))))
    
    (define _round_05up
      (lambda (self prec)
        "Round down unless digit prec-1 is 0 or 5."
        (if (and prec (not (in (pylist-ref (ref self '_int) (- prec 1) "05"))))
	    (_round_down self prec)
	    (- (_round_down self prec)))))
    
    (define _pick_rounding_function
      (dict `((,ROUND_DOWN      . ,_round_down   )
	      (,ROUND_UP        . ,_round_up     )
	      (,ROUND_HALF_UP   . ,_round_half_up)
	      (,ROUND_HALF_DOWN . ,_round_half_down)
	      (,ROUND_HALF_EVEN . ,_round_half_even)
	      (,ROUND_CEILING   . ,_round_ceiling)
	      (,ROUND_FLOOR     . ,_round_floor)
	      (,ROUND_05UP      . ,_round_05up))))

    (define __round__
      (lam (self  (= n None))
        "Round self to the nearest integer, or to a given precision.

        If only one argument is supplied, round a finite Decimal
        instance self to the nearest integer.  If self is infinite or
        a NaN then a Python exception is raised.  If self is finite
        and lies exactly halfway between two integers then it is
        rounded to the integer with even last digit.

        >>> round(Decimal('123.456'))
        123
        >>> round(Decimal('-456.789'))
        -457
        >>> round(Decimal('-3.0'))
        -3
        >>> round(Decimal('2.5'))
        2
        >>> round(Decimal('3.5'))
        4
        >>> round(Decimal('Inf'))
        Traceback (most recent call last):
          ...
        OverflowError: cannot round an infinity
        >>> round(Decimal('NaN'))
        Traceback (most recent call last):
          ...
        ValueError: cannot round a NaN

        If a second argument n is supplied, self is rounded to n
        decimal places using the rounding mode for the current
        context.

        For an integer n, round(self, -n) is exactly equivalent to
        self.quantize(Decimal('1En')).

        >>> round(Decimal('123.456'), 0)
        Decimal('123')
        >>> round(Decimal('123.456'), 2)
        Decimal('123.46')
        >>> round(Decimal('123.456'), -2)
        Decimal('1E+2')
        >>> round(Decimal('-Infinity'), 37)
        Decimal('NaN')
        >>> round(Decimal('sNaN123'), 0)
        Decimal('NaN123')

        "
        (if (not (eq? n None))
            ;; two-argument form: use the equivalent quantize call
            (if (not (isinstance n int))
                (raise (TypeError
			"Second argument to round should be integral"))
		(let ((exp (_dec_from_triple 0 "1" (- n))))
		  ((ref self 'quantize) exp)))

	    ;; one-argument form§x
	    (if (ref self '_is_special)
		(if ((ref self 'is_nan))
		    (raise (ValueError    "cannot round a NaN"))
		    (raise (OverflowError "cannot round an infinity")))
		(int ((ref self '_rescale) 0 ROUND_HALF_EVEN))))))

    (define __floor__
      (lambda (self)
        "Return the floor of self, as an integer.

        For a finite Decimal instance self, return the greatest
        integer n such that n <= self.  If self is infinite or a NaN
        then a Python exception is raised.

        "
	(if (ref self '_is_special)
		(if ((ref self 'is_nan))
		    (raise (ValueError    "cannot round a NaN"))
		    (raise (OverflowError "cannot round an infinity")))
		(int ((ref self '_rescale) 0 ROUND_FLOOR)))))

    (define __ceil__
      (lambda (self)
        """Return the ceiling of self, as an integer.

        For a finite Decimal instance self, return the least integer n
        such that n >= self.  If self is infinite or a NaN then a
        Python exception is raised.

        """
	(if (ref self '_is_special)
	    (if ((ref self 'is_nan))
		(raise (ValueError    "cannot round a NaN"))
		(raise (OverflowError "cannot round an infinity")))
	    (int ((ref self '_rescale) 0 ROUND_CEILING)))))

    (define fma
      (lam (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.
        "
	(twix
	 (let ((other (_convert_other other #:raiseit #t))
	       (third (_convert_other third #:raiseit #t))
	       (fin   (lambda (product)
			((ref product '__add__) third context)))))
	 ;; compute product; raise InvalidOperation if either operand is
	 ;; a signaling NaN or if the product is zero times infinity.	 
        ((if (or (ref self '_is_special) (ref other '_is_special))
	     (twix
	      (let (get-context context))
	      ((equal? (ref self '_exp) "N")  it
	       ((cx-error context) InvalidOperation "sNaN" self))
	      ((equal? (ref other '_exp) "N") it
	       ((cx-error context) InvalidOperation "sNaN" other))
	      ((equal? (ref self '_exp) "n")  it
	       (fin self))
	      ((equal? (ref other '_exp) "n") it
	       (fin other))
	      ((equal? (ref self '_exp) "F")  it
	       (if (not (bool other))
		   ((cx-error context) InvalidOperation "INF * 0 in fma")
		   (pylist-ref _SignedInfinity
			       (logxor (ref self  '_sign)
				       (ref other '_sign)))))
	      ((equal? (ref other '_exp) "F")  it
	       (if (not (bool self))
		   ((cx-error context) InvalidOperation "0 * INF in fma")
		   (pylist-ref _SignedInfinity
			       (logxor (ref self  '_sign)
				       (ref other '_sign)))))
	      #f)) it it)
	
	(fin
	 (_dec_from_triple (logxor (ref self '_sign) (ref other '_sign))
			   (str (* (int (ref self '_int))
				   (int (ref other '_int))))
			   (+ (ref self '_exp) (ref other '_exp)))))))
    
    (define _power_modulo
      (lam (self other modulo (= context None))
        "Three argument version of __pow__"
	(twix
	 ((norm-op self other ) it it)
	 ((norm-op self modulo) it it)
	 (let (get-context context))

	 ;; deal with NaNs: if there are any sNaNs then first one wins,
	 ;; (i.e. behaviour for NaNs is identical to that of fma)
	 (let ((self_is_nan   (ref self   '_isnan))
	       (other_is_nan  (ref other  '_isnan))
	       (modulo_is_nan (ref modulo '_isnan))))
	 
	 ((or (bool self_is_nan) (bool other_is_nan) (bool modulo_is_nan)) it
	  (cond
	   ((= self_is_nan 2)
	    ((cx-error context) InvalidOperation "sNaN"  self))
	   ((= other_is_nan 2)
	    ((cx-error context) InvalidOperation "sNaN"  other))
	   ((modulo_is_nan 2)
	    ((cx-error context) InvalidOperation "sNaN"  modulo))
	   ((bool self_is_nan)
	    (_fix_nan self context))
	   ((bool other_is_nan)
	    (_fix_nan other context))
	   (else
	    (_fix_nan modulo context))))
         
	 ;;check inputs: we apply same restrictions as Python's pow()
	 ((not (and ((ref self   '_isinteger))
		    ((ref other  '_isinteger))
		    ((ref modulo '_isinteger)))) it
	  ((cx-error context) InvalidOperation
	   (+ "pow() 3rd argument not allowed "
	      "unless all arguments are integers")))
	 
	 ((< other 0) it
	  ((cx-error context) InvalidOperation
	   (+ "pow() 2nd argument cannot be "
	      "negative when 3rd argument specified")))

	 ((not (bool modulo)) it
	     ((cx-error context) InvalidOperation
	      "pow() 3rd argument cannot be 0"))

	 ;; additional restriction for decimal: the modulus must be less
	 ;; than 10**prec in absolute value
	 ((>= ((ref modulo 'adjusted)) (cx-prec context)) it
	  ((cx-error context) 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!)
	 ((and (not (bool other)) (not (bool self))) it
	  ((cx-error context) InvalidOperation
	   (+ "at least one of pow() 1st argument "
	      "and 2nd argument must be nonzero ;"
	      "0**0 is not defined")))

	 ;; compute sign of result
	 (let ((sign     (if ((ref other '_iseven))
			     0
			     (ref self '_sign)))
	       (base     (_WorkRep ((ref self  'to_integral_value))))
	       (exponent (_WorkRep ((ref other 'to_integral_value)))))


	   ;; convert modulo to a Python integer, and self and other to
	   ;; Decimal integers (i.e. force their exponents to be >= 0)
	   (set! modulo (abs (int modulo)))

	   ;; compute result using integer pow()
	   (set! base (guile:modulo
		       (* (guile:modulo (ref base 'int) modulo)
			  (modulo-expt 10 (ref base 'exp) modulo))
		       modulo))
	   
	   (let lp ((i (ref exponent 'exp)))
	     (if (> i 0)
		 (begin
		   (set! base (modulo-expt base 10 modulo))
		   (lp (- i 1)))))
	   
	   (set! base (modulo-expt base (ref exponent 'int) modulo))
	   
	   (_dec_from_triple sign (str base) 0)))))

    (define _power_exact 
      (lambda (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 abs(y) = yc*10**ye, with xc
        ;; and yc positive integers not divisible by 10.

        ;; 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.  Express y as a rational number 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 y 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).

	(twix
	 (let ()
	   (define-syntax-rule (clean xc xe n +)
	     (let lp ()
	       (if (= (modulo xc n) 0)
		   (begin
		     (set! xc (/ xc n))
		     (set! xe (+ xe 1))
		     (lp))))))

	 (let* ((x  (_WorkRep self))
		(xc (ref x 'int))
		(xe (ref x 'exp)))
	   (clean xc xe 10 +))
	 
	 (let* ((y  (_WorkRep other))
		(yc (ref y 'int))
		(ye (ref y 'exp)))
	   (clean yc ye 10 +))
	 
        ;; case where xc == 1: result is 10**(xe*y), with xe*y
        ;; required to be an integer
        ((= xc 1) it
	 (set! xe (* xe yc))

	 ;; result is now 10**(xe * 10**ye);  xe * 10**ye must be integral
	 (clean xe ye 10 +)
	 
	 (if (< ye 0)
	     None
	     (let ((exponent (* xe (expt 10 ye)))
		   (zeros    #f))
	       (if (= (ref y 'sign) 1)
		   (set! exponent (- exponent)))
	       
	       ;;if other is a nonnegative integer, use ideal exponent
	       (if (and ((ref other '_isinteger)) (= (ref other '_sign) 0))
		   (begin
		     (let ((ideal_exponent (* (ref self '_exp) (int other))))
		       (set! zeros (min (- exponent ideal_exponent) (- p 1)))))
		   (set! zeros 0))

	       (_dec_from_triple 0 (+ "1" (* "0" zeros)) (- exponent zeros)))))

        ;; case where y is negative: xc must be either a power
        ;; of 2 or a power of 5.
        ((= (ref y 'sign) 1) it
	 (let ((last_digit (modulo xc 10)))
	   (twix
	    (let ((e #f)))
            ((cond
	      ((= (modulo last_digit 2) 0)
	       ;; quick test for power of 2
	       (twix
		((not (= (logand xc (- xc)) xc)) it
		 None)
		;; now xc is a power of 2; e is its exponent
		(let () (set! e (- (_nbits xc) 1)))
		
		;; We now have:
		;;
		;;   x = 2**e * 10**xe, e > 0, and y < 0.
		;;
		;; The exact result is:
		;;
		;;   x**y = 5**(-e*y) * 10**(e*y + xe*y)
		;;
		;; provided that both e*y and xe*y are integers.
		;; Note that if
		;; 5**(-e*y) >= 10**p, then the result can't be expressed
		;; exactly with p digits of precision.
		;;
		;; Using the above, we can guard against large values of ye.
		;; 93/65 is an upper bound for log(10)/log(5), so if
		;;
		;;   ye >= len(str(93*p//65))
		;;
		;; then
		;;
		;;   -e*y >= -y >= 10**ye > 93*p/65 > p*log(10)/log(5),
		;;
		;; so 5**(-e*y) >= 10**p, and the coefficient of the result
		;; can't be expressed in p digits.
		
		;; emax >= largest e such that 5**e < 10**p.
		(let ((emax (quotient (* p 93) 65))))

		((>= ye (len (str emax))) it
		 None)

		;; Find -e*y and -xe*y; both must be integers
		(let ()
		  (set! e  (_decimal_lshift_exact (* e  yc) ye))
		  (set! xe (_decimal_lshift_exact (* xe yc) ye)))
		  
		((or (eq? e None) (eq? xe None)) it
		 None)

		((> e emax) it
		 None)

		(begin
		  (set! xc (expt 5 e))
		  #f)))
	    
	      ((= last_digit 5)
	       (twix
		;; e >= log_5(xc) if xc is a power of 5; we have
		;; equality all the way up to xc=5**2658
		(let* ((e         (quotient (* (_nbits xc) 28) 65))
		       (q         (expt 5 e))
		       (xz        xc)
		       (xc        (quotient q xz))
		       (remainder (modulo   q xz))))
		
		((not (= remainder 0)) it
		 None)
		
		(let () (clean xc e 5 -))
                
		;; Guard against large values of ye, using the same logic as in
		;; the 'xc is a power of 2' branch.  10/3 is an upper bound for
		;; log(10)/log(2).
		(let ((emax (quotient (* p 10) 3))))

		((>= ye (len (str emax))) it
		 None)

		(let ()
		  (set! e  (_decimal_lshift_exact (* e  yc) ye))
		  (set! xe (_decimal_lshift_exact (* xe yc) ye)))
		
		((or (eq? e None) (eq? xe None)) it
		 None)

		((> e emax) it
		 None)

		(begin
		  (set! xc (expt 2 e))
		  #f)))

	      (else
	       None)) it it)

	    ((>= xc (expt 10 p)) it it)

	    (begin
	      (set! xe (+ (- e) (- xe)))
	      (_dec_from_triple 0 (str xc) xe)))))

        ;; now y is positive; find m and n such that y = m/n
	(let ((m #f) (n #f) (xc_bits (_nbits xc))))
	((if (>= ye 0)
	     (begin
	       (set! m (* yc (expt 10 ye)))
	       (set! n 1)
	       #f)
	     (twix
	      ((and (not (= xe 0)) (<= (len (str (abs (* yc xe)))) (- ye))) it
	       None)
	   
	      ((and (not (= xc 1))
		    (<= (len (str (* (abs yc) xc_bits))) (- ye))) it
		None)

	      (begin
		(set! m yc)
		(set! n (expt 10 (- ye)))

		(let lp()
		  (if (and (= (modulo m 2) 0) (= (modulo n 2) 0))
		      (begin
			(set! m (quotient m 2))
			(set! n (quotient n 2)))))
		(let lp()
		  (if (and (= (modulo m 5) 0) (= (modulo n 5) 0))
		      (begin
			(set! m (quotient m 5))
			(set! n (quotient n 5)))))
		#f))) it it)
	
        ;; compute nth root of xc*10**xe
        ((if (> n 1)
	     (begin
	       (twix
		;; if 1 < xc < 2**n then xc isn't an nth power
		((and (not (= xc 1)) (<= xc_bits n)) it
		 None)

		((not (= (modulo xe n) 0)) it
		 None)

		(begin
		  (let ((a (ash 1 (- (quotient (- xc_bits) n)))))
		    (set! xe (quotient xe n))
		
		    ;; compute nth root of xc using Newton's method	       
		    (let lp ()
		      (let* ((x (expt a (- n 1)))
			     (q (quotient xc x))
			     (r (modulo   xc x)))
			(if (<= a q)
			    (if (not (and (= a q) (= r 0)))
				None
				(begin
				  (set! xc a)
				  #f))
			    (begin
			      (set! a (quotient (+ (* a (- n 1)) q) n))
			      (lp)))))))))
	     #f) it it)

        ;; 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.
        ((and (> xc 1) (> m (quotient (* p 100) (_log10_lb xc)))) it
	 None)

	(let ()
	  (set! xc (expt xc m))
	  (set! xe (xe * m)))
	
        ((> xc  (expt 10 p)) it
	 None)

	(begin
	  ;; by this point the result *is* exactly representable
	  ;; adjust the exponent to get as close as possible to the ideal
	  ;; exponent, if necessary
	  (let* ((str_xc (str xc))
		 (zeros (if (and ((ref other '_isinteger))
				 (= (ref other '_sign) 0))
			    (let ((ideal_exponent
				   (* (ref self '_exp) (int other))))
			      (min (- xe ideal_exponent)
				   (- p  (len str_xc))))
			    0)))
	    (_dec_from_triple 0  (+ str_xc (* '0' zeros)) (- xe zeros)))))))

    (define __pow__
      (lam (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.
        "

	(twix
	 ((not (eq? modulo None)) it
	  ((ref self '_power_modulo) other modulo context))

	 ((norm-op self other ) it it)
	 (let (get-context context))

	 ;; either argument is a NaN => result is NaN
	 ((bin-special self other context) it it)

	 ;; 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
	 ((not (bool other)) it
	  (if (not (bool self))
	      ((cx-error context) InvalidOperation "0 ** 0")
	      _One))
	 
	 ;; result has sign 1 iff self._sign is 1 and other is an odd integer
	 (let ((result_sign  0)))

	 ((if (= (ref self '_sign) 1)
	      (twix
		((if ((ref other '_isinteger))
		     (if (not ((ref other '_iseven)))
			 (begin
			   (set! result_sign 1)
			   #f)
			 #f)
		     ;; -ve**noninteger = NaN
		     ;; (-0)**noninteger = 0**noninteger
		     (if (bool self)
			 ((cx-error context) InvalidOperation 
			  "x ** y with x negative and y not an integer")
			 #f)) it it)
		(begin
		  ;; negate self, without doing any unwanted rounding
		  (set! self ((ref self 'copy_negate)))
		  #f))	      
	      #f) it it)

        ;; 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
        ((not (bool self)) it
	 (if (= (ref other '_sign) 0)
	     (_dec_from_triple result_sign "0" 0)
	     (pylist-ref _SignedInfinity result_sign)))

        ;; Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
        (((self '_isinfinity)) it
	 (if (= (ref other '_sign) 0)
	     (pylist-ref _SignedInfinity result_sign)
	     (_dec_from_triple result_sign "0" 0)))

        ;; 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
	(let ((prec (cx-prec context))))
	
        ((equal? self _One) it
	 (let ((exp  #f))	      
	   (if ((ref 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)
	       (let ((multiplier
		      (cond
		       ((= (ref other '_sign) 1)
			0)
		       ((> other prec)
			prec)
		       (else
			(int other)))))

		 (set! exp (* (ref self '_exp) multiplier))
		 (if (< exp (- 1 prec))
		     (begin
		       (set! exp (- 1 prec))
		       ((cx-error context) Rounded))))
	       (begin
		 ((cx-error context) Inexact)
		 ((cx-error context) Rounded)
		 (set! exp (- 1 prec))))

	   (_dec_from_triple result_sign  (+ "1" (* "0" (- exp)) exp))))

        ;; compute adjusted exponent of self
        (let ((self_adj ((ref self 'adjusted)))))

        ;; self ** infinity is infinity if self > 1, 0 if self < 1
        ;; self ** -infinity is infinity if self < 1, 0 if self > 1
        (((ref other '_isinfinity)) it
	 (if (eq? (= (ref other '_sign) 0)
		  (< self_adj           0))
	     (_dec_from_triple result_sign "0" 0)
	     (pylist-ref _SignedInfinity result_sign)))

        ;; from here on, the result always goes through the call
        ;; to _fix at the end of this function.
        (let ((ans   None)
	      (exact #f)
	      (bound (+ ((ref self '_log10_exp_bound))
			((ref other 'adjusted)))))

	  ;; 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.
        
        (if (eq? (>= self_adj 0) (= (ref other '_sign) 0))
            ;; self > 1 and other +ve, or self < 1 and other -ve
            ;; possibility of overflow
            (if (>= bound (len (str (cx-emax context))))
                (set! ans
		      (_dec_from_triple result_sign "1"
					(+ (cx-emax context) 1))))
	    
            ;; self > 1 and other -ve, or self < 1 and other +ve
            ;; possibility of underflow to 0
            (let ((Etiny (cx-etiny context)))
	      (if (>= bound (len (str (- Etiny))))
		  (set! ans (_dec_from_triple result_sign "1" (- Etiny 1))))))

        ;; try for an exact result with precision +1
        (when (eq? ans None)
	  (set! ans ((ref self '_power_exact) other (+ prec 1)))
	  (when (not (eq? ans None))	      
	    (if (= result_sign 1)
		(set! ans (_dec_from_triple 1  (ref ans '_int)
					    (ref ans '_exp))))
	    (set! exact #t)))

        ;; usual case: inexact result, x**y computed directly as exp(y*log(x))
        (when (eq? ans None)
	  (let* ((p  prec)
		     
		 (x  (_WorkRep self))
		 (xc (ref x 'int))
		 (xe (ref x 'exp))

		 (y  (_WorkRep other))
		 (yc (ref y 'int))
		 (ye (ref y 'exp)))

            (if (= (ref y 'sign) 1)
                (set! yc (- yc)))

            ;; compute correctly rounded result:  start with precision +3,
            ;; then increase precision until result is unambiguously roundable
	    (call-with-values
		(lambda ()
		  (let lp ((extra 3))
		    (call-with-values
			(lambda () (_dpower xc xe yc ye (+ p extra)))
		      (lambda (coeff exp)		  
			(if (modulo coeff
				    (* 5 (expt 10 (- (len (str coeff)) p 1))))
			    (values coeff exp)
			    (lp (+ extra 3)))))))
	      (lambda (coeff exp)	    
		(set! ans (_dec_from_triple result_sign (str coeff) exp))))))

        ;; unlike exp, ln and log10, the power function respects the
        ;; rounding mode; no need to switch to ROUND_HALF_EVEN here

        ;; There's a difficulty here when 'other' is not an integer and
        ;; the result is exact.  In this case, the specification
	;; requires that the Inexact flag be raised (in spite of
	;; exactness), but since the result is exact _fix won't do this
        ;; for us.  (Correspondingly, the Underflow signal should also
        ;; be raised for subnormal results.)  We can't directly raise
        ;; these signals either before or after calling _fix, since
	;; that would violate the precedence for signals.  So we wrap
        ;; the ._fix call in a temporary context, and reraise
        ;; afterwards.
        (if (and exact (not ((ref other '_isinteger))))
	    (begin
	      ;; pad with zeros up to length context.prec+1 if necessary; this
	      ;; ensures that the Rounded signal will be raised.
	      (if (<= (len (ref ans '_int)) prec)
		  (let ((expdiff (+ prec 1 (- (len (ref ans '_int))))))
		    (set! ans (_dec_from_triple (ref ans '_sign)
						(+ (ref ans '_int)
						   (* "0" expdiff))
						(- (ref ans '_exp) expdiff)))))

	      ;; create a copy of the current context, with cleared flags/traps
	      (let ((newcontext (cx-copy context)))
		(cx-clear_flags newcontext)

		(for ((exception : _signals)) ()
		     (pylist-set! (cx-traps newcontext) exception 0)
		     (values))

		;; round in the new context
		(set! ans ((ref ans '_fix) newcontext))

		;; raise Inexact, and if necessary, Underflow
		((cx-error newcontext) Inexact)
		(if (bool (pylist-ref (cx-flags newcontext) Subnormal))
		    ((cx-error newcontext) Underflow))

		;; propagate signals to the original context; _fix could
		;; have raised any of Overflow, Underflow, Subnormal,
		;; Inexact, Rounded, Clamped.  Overflow needs the correct
		;; arguments.  Note that the order of the exceptions is
		;; important here.
		(if (bool (pylist-ref (cx-flags newcontext) Overflow))
		    ((cx-error newcontext)
		     Overflow "above Emax" (ref ans '_sign)))
	      
		(for ((exception : (list Underflow Subnormal
					 Inexact   Rounded Clamped))) ()
		   (if (bool (pylist-ref (cx-flags newcontext) exception))
		       ((cx-error newcontext) exception)
		       (values))))
	    
	      (set! ans ((ref ans '_fix) context)))

	    ans)))))

    (define __rpow__
      (lam (self other (= context None))
	   "Swaps self/other and returns __pow__."
	   (twix
	    ((norm-op self other) it it)
	    ((ref 'other '__pow__) self  #:context context))))

    (define normalize
      (lam (self (= context None))
        "Normalize- strip trailing 0s, change anything equal to 0 to 0e0"
	
	(twix
	 (let (get-context context))
	 ((un-special self context) it it)
	 (let ((dup ((ref self _fix) context))))
	 (((dup '_isinfinity)) it dup)
	 ((not (bool dup)) it
	  (_dec_from_triple (ref dup '_sign) "0" 0))

	 (let* ((_int    (ref dup '_int))
		(exp_max (let ((i (cx-clamp context)))
			   (if (= i 0)
			       (cx-emax context)
			       (cx-etop context)))))
	   (let lp ((end (len _int)) (exp (ref dup '_exp)))
	     (if (and (equal? (pylist-ref _int (- end 1))
			      "0")
		      (< exp exp_max))
		 (lp (- end 1) (+ exp 1))
		 (_dec_from_triple
		  (ref dup '_sign)
		  (pylist-slice _int None end None)
		  exp)))))))

    (define quantize
      (lam (self exp (= rounding None) (= context None))
        "Quantize self so its exponent is the same as that of exp.

        Similar to self._rescale(exp._exp) but with error checking.
        "
	(twix
	 (let* ((exp      (_convert_other exp #:raiseit #t))
		(context  (if (eq? context  None) (getcontext) context))
		(rounding (if (eq? rounding None)
			      (cx-rounding context)
			      rounding))))
	  
        ((if (or ((self '_is_special)) ((exp '_is_special)))
	     (let ((ans ((ref self '_check_nans) exp context)))
	       (cond
		((bool ans)
		   ans)
		((or ((ref exp '_isinfinity)) ((ref self '_isinfinity)))
		 (if (and ((ref exp '_isinfinity))
			  ((ref self '_isinfinity)))
		     (Decimal self))  ; if both are inf, it is OK
		 ((cx-error context) InvalidOperation "quantize with one INF"))
		(else
		 #f)))
	     #f) it it)
	
        ;; exp._exp should be between Etiny and Emax
	(let ((_eexp (ref exp '_exp))
	      (Emax  (cx-emax context))))
	
        ((not (and (<= (cx-etiny context) _eexp) (<= _eexp Emax))) it
	 ((cx-error context) InvalidOperation
	  "target exponent out of bounds in quantize"))

        ((not (bool self)) it
	 (let ((ans (_dec_from_triple (ref self '_sign) "0" _eexp)))
	   ((ref ans '_fix) context)))

        (let ((self_adjusted ((ref self 'adjusted)))
	      (prec          (cx-prec context))))

	((> self_adjusted (cx-emax context)) it
	 ((cx-error context) InvalidOperation
	  "exponent of quantize result too large for current context"))
	
        ((> (+ self_adjusted (- _eexp)  1) prec) it
	 ((cx-error context) InvalidOperation
	  "quantize result has too many digits for current context"))

	(let ((ans ((ref self '_rescale) _eexp rounding))))
	 
	(if (> ((ref ans 'adjusted)) Emax)
	    ((cx-error context) InvalidOperation
	     "exponent of quantize result too large for current context"))
	
        (if (> (len (ref ans '_int)) prec)
	 ((cx-error context) InvalidOperation
	  "quantize result has too many digits for current context"))


	(begin
	  ;; raise appropriate flags
	  (if (and (bool ans) (< ((ref ans 'adjusted)) (cx-emin context)))
	      ((cx-error context) Subnormal))

	  (when (> (ref ans '_exp) (ref self '_exp))
	    (if (not (equal? ans self))
		((cx-error context) Inexact))
	    ((cx-error context) Rounded))
	            
	  ;; call to fix takes care of any necessary folddown, and
	  ;; signals Clamped if necessary
	  ((ref ans '_fix) context)))))
    
    (define same_quantum
      (lam (self other (= context None))
        "Return True if self and other have the same exponent; otherwise
        return False.

        If either operand is a special value, the following rules are used:
           * return True if both operands are infinities
           * return True if both operands are NaNs
           * otherwise, return False.
        "
        (let ((other (_convert_other other #:raiseit #t)))
	  (if (or (ref self '_is_special) (ref other '_is_special))
	      (or (and ((ref self 'is_nan))      ((ref other 'is_nan)))
		  (and ((ref self 'is_infinite)) ((ref other 'is_infinite))))))
	
        (= (ref self '_exp) (ref other '_exp))))

    (define _rescale
      (lam (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 mode
        "

	(twix
	 ((ref self '_is_special) it
	  (Decimal self))
	 
	 ((not (bool self)) it
	  (_dec_from_triple (ref self '_sign) "0" exp))

	 (let ((_exp  (ref self '_exp))
	       (_sign (ref self '_sign))
	       (_int  (ref self '_int))))

	 ((>= _exp exp) it
            ;; pad answer with zeros if necessary
            (_dec_from_triple _sign (+ _int (* "0" (- _exp exp))) exp))

	 
        ;; too many digits; round and lose data.  If self.adjusted() <
        ;; exp-1, replace self by 10**(exp-1) before rounding
        (let ((digits (+ (len _int) _exp (- exp))))
	  (if (< digits 0)
	      (set! self (_dec_from_triple _sign "1" (- exp 1)))
	      (set! digits 0))

	  (let* ((this_function (pylist-ref (ref self '_pick_rounding_function)
					    rounding))
		 (changed       (this_function self digits))
		 (coeff         (or (bool
				     (pylist-slice _int None digits None))
				    "0")))
	    (if (= changed 1)
		(set! coeff (str (+ (int coeff) 1))))

	    (_dec_from_triple _sign coeff exp))))))

    (define _round
      (lam (self places rounding)
        "Round a nonzero, nonspecial Decimal to a fixed number of
        significant figures, using the given rounding mode.

        Infinities, NaNs and zeros are returned unaltered.

        This operation is quiet: it raises no flags, and uses no
        information from the context.

        "
        (cond
	 ((<= places 0)
	  (raise (ValueError "argument should be at least 1 in _round")))
	 ((or (ref self '_is_special) (not (bool self)))
	  (Decimal self))
	 (else
	  (let ((ans ((ref self '_rescale) (+ ((self 'adjusted)) 1 (- places))
		      rounding)))
	    ;; it can happen that the rescale alters the adjusted exponent;
	    ;; for example when rounding 99.97 to 3 significant figures.
	    ;; When this happens we end up with an extra 0 at the end of
	    ;; the number; a second rescale fixes this.
	    (if (not (= ((ref ans 'adjusted)) ((ref self 'adjusted))))
		(set! ans ((ref ans '_rescale) (+ ((ans 'adjusted)) 1
						  (- places))
			   rounding)))
	    ans)))))

    (define to_integral_exact
      (lam (self (= rounding None) (= context None))
        "Rounds to a nearby integer.

        If no rounding mode is specified, take the rounding mode from
        the context.  This method raises the Rounded and Inexact flags
        when appropriate.

        See also: to_integral_value, which does exactly the same as
        this method except that it doesn't raise Inexact or Rounded.
        "
	(cond
        ((ref self '_is_special)
	 (let ((ans ((ref self '_check_nans) #:context context)))
	   (if (bool ans)
               ans
	       (Decimal self))))
	        
        ((>= (ref self '_exp) 0)
	 (Decimal self))
        ((not (bool self))	 
	 (_dec_from_triple (ref self '_sign) "0" 0))
	(else
	 (let* ((context  (if (eq? context  None) (getcontext) context))
		(rounding (if (eq? rounding None)
			      (cx-rounding context)
			      rounding))
		(ans      ((ref self '_rescale) 0 rounding)))
	   
	   (if (not (equal? ans self))
	       ((cx-error context) Inexact))
	   
	   ((cx-error context) Rounded)

	   ans)))))

    (define to_integral_value
      (lam (self (= rounding None) (= context None))
	"Rounds to the nearest integer, without raising inexact, rounded."
	(let* ((context  (if (eq? context  None) (getcontext) context))
	       (rounding (if (eq? rounding None)
			     (cx-rounding context)
			     rounding)))
	  
	  (cond
	   ((ref self '_is_special)
	    (let ((ans ((ref self '_check_nans) #:context context)))
	      (if (bool ans)
		  ans
		  (Decimal self))))
	    
	   ((>= (ref self '_exp) 0)
	    (Decimal self))
	    
	   (else
            ((ref self '_rescale) 0 rounding))))))

    ;; the method name changed, but we provide also the old one,
    ;; for compatibility
    (define to_integral to_integral_value)

    (define sqrt
      (lam (self (= context None))
	"Return the square root of self."
	(twix
	 (let (get-context context))	 

	 ((if (ref self '_is_special)
	      (let ((ans ((ref self '_check_nans) #:context context)))
		(if (bool ans)
		    ans
		    (if (and ((self '_isinfinity)) (= (ref self '_sign) 0))
			(Decimal self)
			#f)))

	      #f) it it)

	 ((not (bool self)) it
	  ;; exponent = self._exp // 2.  sqrt(-0) = -0
	  (let ((ans (_dec_from_triple (ref self '_sign) "0"
				       (quotient (ref self '_exp) 2))))
	    ((ref ans '_fix) context)))

	 ((= (ref self '_sign) 1) it
	  ((cx-error context) InvalidOperation "sqrt(-x), x > 0"))

        ;; 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.

        ;; use an extra digit of precision
        (let* ((prec  (+ (cx-prec context) 1))
	       (op    (_WorkRep self))
	       (e     (ash (ref op 'exp) -1))
	       (c     #f)
	       (l     #f)
	       (shift #f)
	       (exact #f))

	  ;; 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.
	  (if (= (logand (ref op 'exp)  1) 1)
	      (begin
		(set! c (* (ref op 'int) 10))
		(set! l (+ (ash (len (ref self '_int)) -1)  1)))
	      (begin
		(set! c (ref op 'int))
		(set! l (ash (+ (len (ref self '_int)) 1) -1))))

	  ;; rescale so that c has exactly prec base 100 'digits'
	  (set! shift (- prec l))
	  (if (>= shift 0)
	      (begin
		(set! c (* c (expt 100 shift)))
		(set! exact #t))
	      (let ((x (expt 100 (- shift))))
		(set! c (quotient c x))
		(let ((remainder (modulo c x)))
		  (set! exact (= remainder 0)))))

	  (set! e (- e shift))

	  ;; find n = floor(sqrt(c)) using Newton's method
	  (let ((n (let lp ((n  (expt 10 prec)))
		     (let ((q (quotient c n)))
		       (if (<= n q)
			   n
			   (lp (ash (+ n q) -1)))))))

	    (set! exact (and exact (= (* n n) c)))

	    (if exact
		;; result is exact; rescale to use ideal exponent e
		(begin
		  (if (>= shift 0)
		      ;; assert n % 10**shift == 0
		      (set! n (quotient n (expt 10 shift)))         
		      (set! n (*        n (expt 10 (- shift)))))
		  (set! e (+ e shift)))
		;; result is not exact; fix last digit as described above
		(if (= (modulo n 5) 0)
		    (set! n (+ n 1))))

	    (let* ((ans (_dec_from_triple 0 (str n) e))
		   ;; round, and fit to current context
		   (context  ((ref context '_shallow_copy)))
		   (rounding ((ref context '_set_rounding) ROUND_HALF_EVEN))
		   (ans      ((ref ans '_fix) context)))
	      (set context 'rounding rounding)
	      ans))))))
    
    (define max
      (lam (self other (= context None))
        "Returns the larger value.

        Like max(self, other) except if one is not a number, returns
        NaN (and signals if one is sNaN).  Also rounds.
        "
	(twix
	 (let ((other (_convert_other other #:raiseit #t))))
	 (let (get-context context))
	 ((if (or (ref self '_is_special) (ref other '_is_special))
	      (begin
		;; If one operand is a quiet NaN and the other is number, then
		;; the number is always returned
		(let ((sn ((ref self  '_isnan)))
		      (on ((ref other '_isnan))))
		  (if (or (bool sn) (bool on))
		      (if (and (= on 1) (= sn 0))
			  ((ref self '_fix) context)
			  (if (and (= sn 1) (= on 0))
			      ((ref other '_fix) context)
			      ((ref self '_check_nans) other context)))
		      #f)))
	      #f) it it)

	 (let ((c ((ref 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.  This is exactly the ordering used in
	       ;; compare_total.
	       (set! c ((ref self 'compare_total) other)))

	   (let ((ans (if (= c -1)
			  other
			  self)))

	     ((ref ans '_fix) context))))))

    (define min
      (lam (self other (= context None))
        "Returns the larger value.

        Like max(self, other) except if one is not a number, returns
        NaN (and signals if one is sNaN).  Also rounds.
        "
	(twix
	 (let ((other (_convert_other other #:raiseit #t))))
	 (let (get-context context))
	 ((if (or (ref self '_is_special) (ref other '_is_special))
	      (begin
		;; If one operand is a quiet NaN and the other is number, then
		;; the number is always returned
		(let ((sn ((ref self  '_isnan)))
		      (on ((ref other '_isnan))))
		  (if (or (bool sn) (bool on))
		      (if (and (= on 1) (= sn 0))
			  ((ref self '_fix) context)
			  (if (and (= sn 1) (= on 0))
			      ((ref other '_fix) context)
			      ((ref self '_check_nans) other context)))
		      #f)))
	      #f) it it)

	 (let ((c ((ref 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.  This is exactly the ordering used in
	       ;; compare_total.
	       (set! c ((ref self 'compare_total) other)))

	   (let ((ans (if (= c -1)
			  self
			  other)))

	     ((ref ans '_fix) context))))))

    (define _isinteger
      (lambda (self)
        "Returns whether self is an integer"
        (cond
	 ((bool (ref self '_is_special))
          #f)
	 ((>= (ref self '_exp) 0)
	  #t)
	 (else
	  (let ((rest (pylist-ref (ref self '_int) (ref self '_exp))))
	    (equal? rest (* "0" (len rest))))))))

    (define _iseven
      (lambda (self)
        "Returns True if self is even.  Assumes self is an integer."
        (if (or (not (bool self)) (> (ref self '_exp) 0))
            #t
	    (in (pylist-ref (ref self '_int) (+ -1 (ref self '_exp)))
		"02468"))))

    (define adjusted
      (lambda (self)
        "Return the adjusted exponent of self"
        (try
          (lambda () (+ (ref self '_exp) + (len (ref self '_int)) -1))
	  ;; If NaN or Infinity, self._exp is string
	  (#:except TypeError (lambda z 0)))))

    (define canonical
      (lambda (self)
        "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.
        "
	self))

    (define compare_signal
      (lam (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.
        "
        (let* ((other (_convert_other other #:raiseit #t))
	       (ans   ((ref self '_compare_check_nans) other context)))
	  (if (bool ans)
	      ans
	      ((ref self 'compare) other #:context context)))))

    (define compare_total
      (lam (self other (= context None))
        "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.
        "

	(twix
	 (let ((other (_convert_other other #:raiseit #t))))

	 ;; if one is negative and the other is positive, it's easy
        ((and (bool (ref self '_sign)) (not (bool other '_sign))) it
	 _NegativeOne)

	((and (not (bool (ref self '_sign))) (bool other '_sign)) it
	 _One)
	
        (let ((sign      (ref self '_sign))
	      ;; let's handle both NaN types
	      (self_nan  ((ref self  '_isnan)))
	      (other_nan ((ref other '_isnan)))))
	
        ((if (or (bool self_nan) (bool other_nan))
	     (if (= self_nan other_nan)
                ;; compare payloads as though they're integers
		 (let ((self_key  (list (len (ref self '_int))
					(ref self '_int)))
		       (other_key (list (len (ref other '_int))
					(ref other '_int))))
		   (cond
		    ((< self_key other_key)
		     (if (bool sign)
			 _One
			 _NegativeOne))
		    ((> self_key other_key)
		     (if sign
			 _NegativeOne
			 _One))
		    (else
		     _Zero)))

		 (if (bool sign)
		     (cond
		      ((= self_nan 1)
			 _NegativeOne)
		      ((= other_nan 1)
		       _One)
		      ((= self_nan 2)
		       _NegativeOne)
		      ((= other_nan 2)
		       _One)
		      (else
		       #f))		     
		     (cond
		      ((= self_nan 1)
		       _One)
		      ((= other_nan 1)
		       _NegativeOne)
		      ((= self_nan 2)
		       _One)
		      ((= other_nan 2)
		       _NegativeOne)
		      (else #f))))
	     #f) it it)

        ((< self other) it
	 _NegativeOne)

	((> self > other) it
	 _One)

        ((< (ref self '_exp) (ref other '_exp)) it
	 (if (bool sign)
	     _One
	     _NegativeOne))
	
        ((> (ref self '_exp) (ref other '_exp)) it
	 (if (bool sign)
	     _NegativeOne
	     _One))
	
	_Zero)))


    (define compare_total_mag
      (lam (self other (= context None))
        "Compares self to other using abstract repr., ignoring sign.

        Like compare_total, but with operand's sign ignored and assumed to be 0.
        "
        (let* ((other (_convert_other other #:raiseit #t))
	       (s     ((ref self  'copy_abs)))
	       (o     ((ref other 'copy_abs))))
	  ((ref s 'compare_total) o))))

    (define copy_abs
      (lambda (self)
        "Returns a copy with the sign set to 0. "
	(_dec_from_triple 0 (ref self '_int)
			  (ref self '_exp) (ref self '_is_special))))

    (define copy_negate
      (lambda (self)
        "Returns a copy with the sign inverted."
        (if (bool (ref self '_sign))
	    (_dec_from_triple 0 (ref self '_int)
			  (ref self '_exp) (ref self '_is_special))
	    (_dec_from_triple 1 (ref self '_int)
			      (ref self '_exp) (ref self '_is_special)))))
    
    (define copy_sign
      (lam (self other (= context None))
        "Returns self with the sign of other."
        (let ((other (_convert_other other #:raiseit #t)))
	  (_dec_from_triple (ref other 'sign) (ref self '_int)
			    (ref self '_exp) (ref self '_is_special)))))

    (define exp
      (lam (self (= context None))
        "Returns e ** self."

	(twix
	 (let (get-context context))

	 ;; exp(NaN) = NaN
	 (let ((ans ((ref self '_check_nans) #:context context))))
	     
	 (ans it it)

	 ;; exp(-Infinity) = 0
	 ((= ((ref self '_isinfinity)) -1) it
	  _Zero)

	 ;; exp(0) = 1
	 ((not (bool self)) it
	  _One)

	 ;; exp(Infinity) = Infinity
	 ((= ((ref self '_isinfinity)) 1) it
	  (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.
	 (let ((p   (cx-prec context))
	       (adj ((ref 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.
	 (let* ((sign  (ref self '_sign))
		(emax  (cx-emax context))
		(etiny (cx-etiny context))
		(ans
		 (cond
		  ((and (= sign 0)
			(> adj (len (str (* (+ emax 1) 3)))))
		   ;; overflow
		   (_dec_from_triple 0  "1" (+ emax 1)))
		  
		  ((and (= sign 1)
			(> adj (len (str (* (- 1 etiny) 3)))))
		   ;; underflow to 0
		   (_dec_from_triple 0 "1" (- etiny 1)))
		  ((and (= sign  0) (< adj (- p)))
		   ;; p+1 digits; final round will raise correct flags
		   (_dec_from_triple 0 (+ "1" (* "0" (- p 1)) "1") (- p)))
		  ((and (= sign 1) (< adj (- (+ p 1))))
		   ;; p+1 digits; final round will raise correct flags
		   (_dec_from_triple 0  (* "9" (+ p 1)) (- (+ p 1))))
		  (else
		   ;; general case
		   (let* ((op (_WorkRep self))
			  (c  (ref op 'int))
			  (e  (ref op 'exp)))
		     (if (= (ref op 'sign) 1)
			 (set! c (- c)))
		    
		     ;; compute correctly rounded result: increase precision by
		     ;; 3 digits at a time until we get an unambiguously
		     ;; roundable result
		    
		     (let lp ((extra 3))
		       (call-with-values (lambda () (_dexp c e (+ p extra)))
			 (lambda (coeff exp)
			   (if (not (= (modulo
					coeff
					(* 5 (expt 10
						   (- (len (str coeff)) p 1))))
				       0))
			      
			       (_dec_from_triple 0 (str coeff) exp)
			       (lp (+ extra 3))))))))))))

	 ;; at this stage, ans should round correctly with *any*
	 ;; rounding mode, not just with ROUND_HALF_EVEN
	
	 (let* ((context  ((ref context '_shallow_copy)))
		(rounding ((ref context '_set_rounding) ROUND_HALF_EVEN))
		(ans      ((ref ans '_fix) context)))
	  
	   (set context 'rounding rounding)
	   ans))))
    
    (define is_canonical
      (lambda (self)
        "Return True if self is canonical; otherwise return False.

        Currently, the encoding of a Decimal instance is always
        canonical, so this method returns True for any Decimal.
        "
        #t))

    (define is_finite
      (lambda (self)
        "Return True if self is finite; otherwise return False.

        A Decimal instance is considered finite if it is neither
        infinite nor a NaN.
        "
	(not (ref self '_is_special))))

    (define is_infinite
      (lambda (self)
        "Return True if self is infinite; otherwise return False."
        (equal? (ref self '_exp) "F")))

    (define is_nan
      (lambda (self)
        "Return True if self is a qNaN or sNaN; otherwise return False."
	(let ((e (ref self '_exp)))
	  (or (equal? e "n") (equal? e "N")))))

    (define is_normal
      (lam (self (= context None))
        "Return True if self is a normal number; otherwise return False."
        (if (or (bool (ref self '_is_special)) (not (bool self)))
            #f
	    (let ((context (if (eq? context None)
			       (getcontext)
			       context)))
	      (<= (cx-emin context) ((ref self 'adjusted)))))))

    (define is_qnan
      (lambda (self)
        "Return True if self is a quiet NaN; otherwise return False."
       (equal? (ref self '_exp) "n")))

    (define is_signed
      (lambda (self)
        "Return True if self is negative; otherwise return False."
       (= (ref self '_sign) 1 )))

    (define is_snan
      (lambda (self)
        "Return True if self is a signaling NaN; otherwise return False."
	(equal? (ref self '_exp) "N")))

    (define is_subnormal
      (lam (self (= context None))
	   "Return True if self is subnormal; otherwise return False."
	   (if (or (bool (ref self '_is_special)) (not (bool self)))
	       #f
	       (let ((context (if (eq? context None)
				  (getcontext)
				  context)))
		 (> (cx-emin context) ((ref self 'adjusted)))))))

    (define is_zero
      (lambda (self)
        "Return True if self is a zero; otherwise return False."
        (and (not (bool (ref self '_is_special)))
	     (equal? (ref self '_int) "0"))))

    (define _ln_exp_bound
      (lambda (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
        (let ((adj (+ (ref self '_exp) (len (ref self '_int)) (- 1))))
	  (cond
	   ((>= adj 1)
            ;; argument >= 10; we use 23/10 = 2.3 as a lower bound for ln(10)
            (- (len (str (floor-quotient (* adj 23) 10))) 1))
	   ((<= adj -2)
            ;; argument <= 0.1
            (- (len (str (floor-quotient (* (- (+ 1 adj)) 23) 10))) 1))
	   (else
	    (let* ((op (_WorkRep self))
		   (c  (ref op 'int))
		   (e  (ref op 'exp)))
	      (if (= adj 0)
		  ;; 1 < self < 10
		  (let ((num (str (- c (expt 10 (- e)))))
			(den (str c)))
		    (- (len num) (len den) - (< num den)))
		  ;; adj == -1, 0.1 <= self < 1
		  (+ e (len (str (- (expt 10 (- e)) c))) (- 1)))))))))


    (define ln
      (lam (self (= context None))
        "Returns the natural (base e) logarithm of self."

	(twix
	 (let (get-context context))

	 ;; ln(NaN) = NaN
	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 ;; ln(0.0) == -Infinity
         ((not (bool self)) it
	  _NegativeInfinity)

	 ;; ln(Infinity) = Infinity
	 ((= ((ref self '_isinfinity)) 1) it
          _Infinity)

	 ;; ln(1.0) == 0.0
	 ((equal? self _One) it
	     _Zero)

	 ;; ln(negative) raises InvalidOperation
	 (if (= (ref self '_sign) 1)
	  ((cx-error context) InvalidOperation
	   "ln of a negative value"))

	 ;; result is irrational, so necessarily inexact
	 (let* ((op (_WorkRep self))
		(c  (ref op 'int))
		(e  (ref op 'exp))
		(p  (cx-prec context))))
	   

	 ;; correctly rounded result: repeatedly increase precision by 3
	 ;; until we get an unambiguously roundable result
	 (let ((places (+ p
			  (- ((ref self '_ln_exp_bound)))
			  2)) ;; at least p+3 places
	       (ans #f))
	    (let lp ((places places))
	      (let ((coeff (_dlog c  e places)))
		;; assert len(str(abs(coeff)))-p >= 1
		(if (not (= (modulo coeff
				    (* 5 (expt 10 (- (len (str (abs coeff)))
						     p 1))))
			    0))
		    (set! ans (_dec_from_triple (int (< coeff 0))
						(str (abs coeff))
						(- places)))
						    
		    (lp (+ places 3)))))

	    (let* ((context  ((ref context '_shallow_copy)))
		   (rounding ((ref context '_set_rounding) ROUND_HALF_EVEN))
		   (ans      ((ref ans '_fix) context)))
	      (set context 'rounding rounding)
	      ans)))))

    (define _log10_exp_bound
      (lambda (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

        (let ((adj (+ (ref self '_exp) (len (ref self '_int)) (- 1))))
	  (cond
	   ((>= adj 1)
	    ;; self >= 10
	    (- (len (str adj)) 1))
	   ((<= adj -2)
	    ;;# self < 0.1
	    (- (len (str (- (+ 1 adj)))) 1))
	   (else
	    (let* ((op (_WorkRep self))
		   (c  (ref op 'int))
		   (e  (ref op 'exp)))
	      (if(= adj 0)
		 ;; 1 < self < 10
		 (let ((num (str (- c (expt 10 (- e)))))
		       (den (str (* 231 c))))
		   (+ (len num) (- (len den)) (- (< num den)) 2))	       
		 ;; adj == -1, 0.1 <= self < 1
		 (let ((num (str (- (expt 10 (- e)) c))))
		   (+ (len num) e (- (< num "231")) (- 1))))))))))

    (define log10
      (lam (self (= context None))
        "Returns the base 10 logarithm of self."
	(twix
	 (let (get-context context))

	 ;; log(NaN) = NaN
	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 ;; log10(0.0) == -Infinity
         ((not (bool self)) it
	  _NegativeInfinity)

	 ;; log10(Infinity) = Infinity
	 ((= ((ref self '_isinfinity)) 1) it
          _Infinity)

	 ;; log10(1.0) == 0.0
	 ((equal? self _One) it
	     _Zero)

	 ;; ln(negative) raises InvalidOperation
	 (if (= (ref self '_sign) 1)
	  ((cx-error context) InvalidOperation
	   "log10 of a negative value"))

	 (let ((ans #f)))
	 
	 ;; log10(10**n) = n
	 (begin
	   (if (and (equal? (string-ref (ref self '_int) 0) #\1)
		    (equal? (pylist-slice (ref self '_int) 1 None None)
			    (* "0" (- (len (ref self '_int)) 1))))
	       ;;answer may need rounding
	       (set! ans (Decimal (+ (ref self '_exp) (len (ref self '_int)) (- 1))))
	       ;; result is irrational, so necessarily inexact
	       (let* ((op (_WorkRep self))
		      (c  (ref op 'int))
		      (e  (ref op 'exp))
		      (p  (cx-prec context)))

		 ;; correctly rounded result: repeatedly increase precision
		 ;; until result is unambiguously roundable
		 (let lp ((places (+ p (- ((ref self '_log10_exp_bound))) 2)))
		   (let ((coeff (_dlog10 c e places)))
		     ;; assert len(str(abs(coeff)))-p >= 1
		     (if (not (= (modulo coeff
					 (* 5 (expt 10
						    (- (len (str (abs coeff)))
						       p 1))))
				 0))
			 (set! ans
			       (_dec_from_triple (int (< coeff 0))
						 (str (abs coeff))
						 (- places)))
			 (lp (+ places 3)))))))

	   (let* ((context  ((ref context '_shallow_copy)))
		  (rounding ((ref context '_set_rounding) ROUND_HALF_EVEN))
		  (ans      ((ref ans '_fix) context)))
	     (set context 'rounding rounding)
	     ans)))))

    (define logb
      (lam (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).
        "
	(twix
	 (let (get-context context))

	 ;; logb(NaN) = NaN
	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 ;; logb(+/-Inf) = +Inf
	 (((ref self '_isinfinity)) it
	  _Infinity)

	 ;; logb(0) = -Inf, DivisionByZero
	 ((not (bool self)) it
	  ((cx-error context) 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.
	 (let ((ans (Decimal ((ref self 'adjusted)))))
	   ((ref ans '_fix) context)))))

    (define _islogical
      (lambda (self)
        "Return True if self is a logical operand.

        For being logical, it must be a finite number with a sign of 0,
        an exponent of 0, and a coefficient whose digits must all be
        either 0 or 1.
        "
        (if (or (not (= (ref self '_sign) 0))
		(not (= (ref self '_exp)  0)))
            #f
	    (for ((dig : (ref self '_int))) ()
		 (if (not (or (equal? dig "0") (equal? dig "1")))
		     (break #f))
		 #:final #t))))
    
    (define _fill_logical
      (lambda (self context opa opb)
	(define (o opa dif)
	  (cond
	   ((> dif 0)
	    (* "0" dif) opa)
	   ((< dif 0)
	    (pylist-slice opa (- (cx-prec context) None None)))
	   (else
	    opa)))
	
        (let* ((dif (- (cx-prec context) (len opa)))
	       (opa (o opa dif))
	       (dif (- (cx-prec context) (len opb)))
	       (opb (o opb dif)))
	  (values opa opb))))

    (define logical_*
      (lambda (logand)
      (lam (self other (= context None))
	"Applies an 'and' operation between self and other's digits."

	(twix
	 (let (get-context context))	 
	 (let ((other (_convert_other other #:raiseit #t))))

	 (if (or (not ((ref self '_islogical)))  (not ((ref other '_islogical))))
	  ((cx-error context) InvalidOperation))

	 ;; fill to context.prec
	 (call-with-values
	     (lambda ()
	       ((ref self '_fill_logical)
		context (ref self '_int) (ref other '_int)))
	   (lambda (opa opb)
	     ;; make the operation and clean starting zeroes
	     (_dec_from_triple
	      0
	      (for ((a : opa) (b : opb)) ((l '()) (f #t))
		   (let ((i (logand (int a) (int b))))
		     (if (and f (= i 0))
			 (values l                #t)
			 (values (cons (str i) l) #f)))
		   #:final
		   (if (null? l)
		       "0"
		       (list->string (reverse l))))
	      0)))))))

    (define logical_and (lambda x (apply (logical_* logand) x)))
    (define logical_or  (lambda x (apply (logical_* logior) x)))
    (define logical_xor (lambda x (apply (logical_* logxor) x)))
    
    (define logical_invert
      (lam (self (= context None))
        "Invert all its digits."
        (let ((context (if (eq? context None)
			   (getcontext)
			   context)))
	  (logical_xor self
		       (_dec_from_triple 0 (* "1" (cx-prec context)) 0)
		       context))))

    (define x_mag
      (lambda (nott)
      (lam (self other (= context None))
        "Compares the values numerically with their sign ignored."
	(twix
	 (let ((other (_convert_other other  #:raiseit #t))))
	 (let (get-context context))	 

	 ((if (or (bool (ref self '_is_special)) (bool (other '_is_special)))
	      ;; If one operand is a quiet NaN and the other is number, then the
	      ;; number is always returned
	      (let ((sn ((ref self  '_isnan)))
		    (on ((ref other '_isnan))))
		(if (or (bool sn) (bool on))
		    (cond
		     ((and (= on 1) (= sn 0))
		      ((ref self '_fix) context))
		     ((and (= on 0) (= sn 1))
		      ((ref other '_fix) context))
		     (else
		      ((ref self '_check_nans) other context)))
		    #f))
	      #f) it it)

	 (let* ((s   ((ref self 'copy_abs)))
		(o   ((ref other 'copy_abs)))
		(c   ((ref s '_cmp) o))
		(c   (if (= c 0)
			 ((self 'compare_total) other)
			 c))
		(ans (if (nott (= c -1)) other self)))
		 
	   ((ref ans '_fix) context))))))

    (define max_mag (lambda y (apply (x_mag (lambda (x) x)) y)))
    (define min_mag (lambda y (apply (x_mag not) y)))
      
    (define next_minus
      (lam (self (= context None))
	"Returns the largest representable number smaller than itself."

	(twix
	 (let (get-context context))	 
	 
	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 ((= ((ref self '_isinfinity)) -1) IT
	  _NegativeInfinity)

	 ((= ((ref self '_isinfinity)) 1) IT
	  (_dec_from_triple 0  (* '9' (cx-prec context)) (cx-etop context)))

	
	 (let* ((context  ((ref context 'copy)))
		(rounding ((ref context '_set_rounding) ROUND_FLOOR)))
	   ((context '_ignore_all_flags))
	   (let ((new_self ((ref self '_fix) context)))
	     (if (not (equal? self new_self))
		 new_self
		 ((ref self '__sub__)
		  (_dec_from_triple 0 "1" (- (cx-etiny context) 1))
		  context)))))))
    
    (define next_plus
      (lam (self (= context None))
	"Returns the largest representable number smaller than itself."

	(twix
	 (let (get-context context))	 
	 
	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 ((= ((ref self '_isinfinity)) 1) it
	  _Infinity)

	 ((= ((ref self '_isinfinity)) -1) it
	  (_dec_from_triple 1  (* '9' (cx-prec context)) (cx-etop context)))

	
	 (let* ((context  ((ref context 'copy)))
		(rounding ((ref context '_set_rounding) ROUND_CEILING)))
	   ((context '_ignore_all_flags))
	   (let ((new_self ((ref self '_fix) context)))
	     (if (not (equal? self new_self))
		 new_self
		 ((ref self '__add__)
		  (_dec_from_triple 0 "1" (- (cx-etiny context) 1))
		  context)))))))
    
    (define next_toward
      (lam (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.
        "

	(twix
	 (let ((other (_convert_other other  #:raiseit #t))))
	 (let (get-context context))	 

	 (let ((ans ((ref self '_check_nans) #:context context))))
	 (ans it it)

	 (let ((comparison ((ref self '_cmp) other))))
	 
	 ((= comparison 0) it
	  ((ref self 'copy_sign) other))

	 (let ((ans (if (= comparison -1)
			((ref self 'next_plus)  context)
			;; comparison == 1
			((ref self 'next_minus) context))))

	   ;; decide which flags to raise using value of ans
	   (cond
	    (((ref ans '_isinfinity))
	     ((cx-error context) Overflow "Infinite result from next_toward"
	      (ref ans '_sign))
	     ((cx-error context) Inexact)
	     ((cx-error context) Rounded))
	      
	    ((< ((ref ans 'adjusted)) (cx-emin context))
	     ((cx-error context) Underflow)
	     ((cx-error context) Subnormal)
	     ((cx-error context) Inexact)
	     ((cx-error context) Rounded)
	     ;; if precision == 1 then we don't raise Clamped for a
	     ;; result 0E-Etiny.
	     (if (not (bool ans))
		 ((cx-error context) Clamped)))
	    (else #f))

	  ans))))

    (define number_class
      (lam (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
        "
	(twix
	 (((ref self 'is_snan)) it
	  "sNaN")
	 (((ref self 'is_qnan)) it
	  "NaN")
	 (let ((inf ((ref self '_isinfinity)))))
	 ((= inf 1)  it
	  "+Infinity")
	 ((= inf -1) it
	  "-Infinity")
        (((ref self 'is_zero)) it
	 (if (bool (ref self '_sign))
	     "-Zero"
	     "+Zero"))

	(let (get-context context))	 

        (((ref self 'is_subnormal) #:context context) it
	 (if (bool (ref self '_sign))
	     "-Subnormal"
	     "+Subnormal"))

	;; just a normal, regular, boring number, :)
        (if (bool (ref self '_sign))
            "-Normal"
	    "+Normal"))))

    (define radix
      (lambda (self)
        "Just returns 10, as this is Decimal"
        (Decimal 10)))

    (define rotate
      (lam (self other (= context None))
	"Returns a rotated copy of self, value-of-other times."
	(twix
	 (let (get-context context))	 	 
	 (let ((other (_convert_other other  #:raiseit #t))))

	 (let ((ans ((ref other '_check_nans) #:context context))))
	 (ans it it)
	 
	 ((not (= (ref other '_exp) 0)) it
	  ((cx-error context) InvalidOperation))

	 (let ((o (int other))
	       (p (cx-prec context))))
	 
	 ((not (and (<= (- p) o) (<= o p))) it	  
	  ((cx-error context) InvalidOperation))

        (((ref self '_isinfinity)) it
	 (Decimal self))

        ;; get values, pad if necessary
        (let* ((torot   (int other))
	       (rotdig  (ref self '_int))
	       (topad   (- p (len rotdig))))
	  (cond
	   ((> topad 0)
	    (set! rotdig (+ (* "0" topad) + rotdig)))
	   ((< topad 0)
            (set! rotdig (pylist-slice rotdig (- topad) None None)))
	   (else #f))

	  (let ((rotated (+ (pylist-slice rotdig torot None  None)
			    (pylist-slice rotdig None  torot None))))
	    (_dec_from_triple (ref self '_sign)
			      (or (bool (py-lstrip rotated "0")) "0")
			      (ref self '_exp)))))))

    (define scaleb
      (lam (self other (= context None))
	"Returns self operand after adding the second value to its exp."
	(twix
	 (let (get-context context))	 	 
	 (let ((other (_convert_other other  #:raiseit #t))))

	 (let ((ans ((ref other '_check_nans) #:context context))))
	 (ans it it)

	 ((not (= (ref other '_exp))) it
	  ((cx-error context) InvalidOperation))
	 
	 (let ((liminf  (* -2  (+ (cx-emax context) (cx-prec context))))
	       (limsup  (*  2  (+ (cx-emax context) (cx-prec context))))))

	 ((not (let ((o (int other)))
		 (and (<= liminf o)
		      (<= o limsup)))) it
	  ((cx-error context) InvalidOperation))

	 (((ref self '_isinfinity)) it
	  (Decimal self))

	 (let* ((d (_dec_from_triple (ref self '_sign)
				     (ref self '_int)
				     (+ (ref self '_exp) (int other))))
		(d ((ref d '_fix) context)))
	   d))))

    (define shift
      (lam (self other (= context None))
	"Returns a rotated copy of self, value-of-other times."
	(twix
	 (let (get-context context))	 	 
	 (let ((other (_convert_other other  #:raiseit #t))))

	 (let ((ans ((ref other '_check_nans) #:context context))))
	 (ans it it)
	 
	 ((not (= (ref other '_exp) 0)) it
	  ((cx-error context) InvalidOperation))

	 (let ((o (int other))
	       (p (cx-prec context))))
	 
	 ((not (and (<= (- p) o) (<= o p))) it	  
	  ((cx-error context) InvalidOperation))

        (((ref self '_isinfinity)) it
	 (Decimal self))

        ;; get values, pad if necessary
        (let* ((torot   (int other))
	       (rotdig  (ref self '_int))
	       (topad   (- p (len rotdig))))

	  (cond
	   ((> topad 0)
	    (set! rotdig (+ (* "0" topad) + rotdig)))
	   ((< topad 0)
            (set! rotdig (pylist-slice rotdig (- topad) None None)))
	   (else #f))

	  
	  ;; let's shift!
	  (let ((shifted (if (< torot 0)
            (pylist-slice rotdig None torot None)
            (pylist-slice (+ rotdig (* "0" torot))
            (- p) None None))))

	    (_dec_from_triple (ref self '_sign)
			      (or (bool (py-lstrip shifted "0")) "0")
			      (ref self '_exp)))))))

    ;; Support for pickling, copy, and deepcopy
    ;; def __reduce__(self):
    ;;    return (self.__class__, (str(self),))

    ;; def __copy__(self):
    ;;     if type(self) is Decimal:
    ;;        return self     # I'm immutable; therefore I am my own clone
    ;;    return self.__class__(str(self))

    ;; def __deepcopy__(self, memo):
    ;;    if type(self) is Decimal:
    ;;        return self     # My components are also immutable
    ;;    return self.__class__(str(self))

    ;; PEP 3101 support.  the _localeconv keyword argument should be
    ;; considered private: it's provided for ease of testing only.    
    (define __format__
      (lam (self specifier (= context None) (= _localeconv None))
        "Format a Decimal instance according to the given specifier.

        The specifier should be a standard format specifier, with the
        form described in PEP 3101.  Formatting types 'e', 'E', 'f',
        'F', 'g', 'G', 'n' and '%' are supported.  If the formatting
        type is omitted it defaults to 'g' or 'G', depending on the
        value of context.capitals.
        "

        ;; Note: PEP 3101 says that if the type is not present then
        ;; there should be at least one digit after the decimal point.
        ;; We take the liberty of ignoring this requirement for
        ;; Decimal---it's presumably there to make sure that
        ;; format(float, '') behaves similarly to str(float).
        (if (eq? context None)
            (set! context (getcontext)))

	(twix
	 (let ((spec
		(_parse_format_specifier specifier
					 #:_localeconv _localeconv))))

	 (let ((type (pylist-ref spec "type"))))
	 ;; special values don't care about the type or precision
	 ((bool (ref self '_is_special)) it
	  (let ((sign (_format_sign (ref self '_sign) spec))
		(body (str ((ref self 'copy_abs)))))
	    (if (equal? type "%")
		(set! body (+ body "%")))
	    (_format_align sign body spec)))

	 ;; a type of None defaults to 'g' or 'G', depending on context
	 (if (eq? type None)
	     (pylist-set! spec "type"
			  (if (= (cx-cap context) 0) "g" "G")))

	 (let ((type (pylist-ref spec "type"))))
	 ;; if type is '%', adjust exponent of self accordingly
	 (if (equal? type "%")
	     (set! self
		   (_dec_from_triple (ref self '_sign)
				     (ref self '_int)
				     (+ (ref self '_exp) 2))))

	 ;; round if necessary, taking rounding mode from the context
	 (let ((rounding  (cx-rounding context))
	       (precision (pylist-ref spec "precision")))
	   (if (not (eq? precision None))
	       (cond
		((in type "eE")
		 (set! self ((ref self '_round) (+ precision 1) rounding)))
		((in type "fF%")
		 (set! self ((ref self '_rescale) (- precision) rounding)))
		((and (in type "gG") (> (len (ref self '_int)) precision))
		 (set! self  ((ref self '_round) precision rounding)))
		(else #t))))
	  
	 ;; special case: zeros with a positive exponent can't be
	 ;; represented in fixed point; rescale them to 0e0.
	 (if (and (not (bool self)) (> (ref self '_exp) 0) (in type "fF%"))
	     (set! self ((ref self '_rescale) 0 rounding)))

	 ;; figure out placement of the decimal point
	 (let* ((leftdigits (+ (ref self '_exp) (len (ref self '_int))))
		(dotplace
		 (cond
		  ((in type "eE")
		   (if (and (not (bool self)) (not (eq? precision None)))
		       (- 1 precision)
		       1))
		 ((in type "fF%")
		  leftdigits)
		 ((in type "gG")
		  (if (and (<= (ref self '_exp) 0) (> leftdigits -6))
		      leftdigits
		      1))
		 (else
		  1))))
		 
	   ;; find digits before and after decimal point, and get exponent
	   (call-with-values
	       (lambda ()
		 (cond
		  ((< dotplace 0)
		   (values '0'
			   (+ (* "0" (- dotplace)) (ref self '_int))))
		  ((> dotplace (len (ref self '_int)))
		   (values (+ (ref self '_int) (* "0" (- dotplace
							 (len (ref self '_int)))))
			   ""))
		  (else
		   (values
		    (or (bool (pylist-slice (ref self '_int) None dotplace None))
			"0")
		    (pylist-slice (ref self '_int) dotplace None None)))))
	     (lambda (intpart fracpart)
	       (let ((exp (- leftdigits dotplace)))
		 ;; done with the decimal-specific stuff;  hand over the rest
		 ;; of the formatting to the _format_number function
		 (_format_number (ref self '_sign) intpart fracpart exp spec)))))))))

(define _dec_from_triple
  (lam (sign  coefficient  exponent (= special #f))
    "Create a decimal instance directly, without any validation
    normalization (e.g. removal of leading zeros) or argument
    conversion.

    This function is for *internal use only*.
    "

    (let ((self ((ref object '__new__) Decimal)))
      (set self '_sign       sign)
      (set self '_int        coefficient)      
      (set self '_exp        exponent)
      (set self '_is_special special)

      self)))

;; Register Decimal as a kind of Number (an abstract base class).
;; However, do not register it as Real (because Decimals are not
;; interoperable with floats).
;; _numbers.Number.register(Decimal)

;; ##### Context class #######################################################

(define-python-class _ContextManager (object)
    "Context manager class to support localcontext().

      Sets a copy of the supplied context in __enter__() and restores
      the previous decimal context in __exit__()
    "
    (define __init__
      (lambda (self new_context)
        (set self 'new_context ((ref new_context 'copy)))))
    
    (define __enter__
      (lambda (self)
        (set self 'saved_context (getcontext))
        (setcontext (ref self 'new_context))
        (ref self 'new_context)))
    
    (define __exit__
      (lambda (self t v tb)
        (setcontext (ref self 'saved_context)))))

(define DefaultContext #f)

(define-syntax-rule (setq s q m)
  (if (eq? q None)
      (if (bool m)
	  (set s 'q  (ref m 'q)))
      (set s 'q q)))

(define-python-class Context (object)
    "Contains the context for a Decimal instance.

    Contains:
    prec - precision (for use in rounding, division square roots..)
    rounding - rounding type (how you round)
    traps - If traps[exception] = 1, then the exception is
                    raised when it is caused.  Otherwise, a value is
                    substituted in.
    flags  - When an exception is caused, flags[exception] is set.
             (Whether or not the trap_enabler is set)
             Should be reset by user of Decimal instance.
    Emin -   Minimum exponent
    Emax -   Maximum exponent
    capitals -      If 1, 1*10^1 is printed as 1E+1.
                    If 0, printed as 1e1
    clamp -  If 1, change exponents if too high (Default 0)
    "

    (define __init__
      (lam (self (= prec  None) (= rounding None) (= Emin None)
		 (= Emax  None) (= capitals None) (= clamp None)
		 (= flags None) (= traps None)    (= _ignored_flags None))
        ;; Set defaults; for everything except flags and _ignored_flags,
        ;; inherit from DefaultContext.
        (let ((dc DefaultContext))
	  (setq self prec     dc)
	  (setq self rounding dc)
	  (setq self Emin     dc)
	  (setq self Emax     dc)
	  (setq self capitals dc)
	  (setq self clamp    dc)
	  
	  (set self '_ignored_flags
	       (if (eq? _ignored_flags None)
		   (py-list)
		   _ignored_flags))
	  (set self 'traps
	       (cond
		((eq? traps None)
		 ((ref (ref dc traps) 'copy)))
		((not (isinstance traps dict))
		 (dict (for ((s : (+ _signals traps))) ((l '()))
			    (cons (list s (int (in s traps))) l)
			    #:final (reverse l))))
		(else traps)))

	  (set self 'flags
	       (cond
		((eq? flags None)
		 ((ref dict 'fromkeys) _signals 0))
		((not (isinstance flags dict))
		 (dict (for ((s : (+ _signals flags))) ((l '()))
			    (cons (list s (int (in s flags))) l)
			    #:final (reverse l))))
		(else flags))))))

    (define _set_integer_check
      (lambda (self name value vmin vmax)
        (if (not (isinstance value int))
            (raise (TypeError (format #f "~a must be an integer" name))))
	
        (cond
	 ((equal? vmin "-inf")
	  (if (> value vmax)
	      (raise (ValueError (format #f "~a must be in [~a, ~a]. got: ~a"
					 name vmin vmax value)))))
	 ((equal? vmax "inf")
	  (if (< value vmin)
	      (raise (ValueError (format #f "~a must be in [~a, ~a]. got: ~a"
					 name vmin vmax value)))))
	 (else
	  (if (or (< value vmin) (> value vmax))
	      (raise (ValueError (format #f "~a must be in [~a, ~a]. got ~a"
					 name vmin vmax value))))))
        (rawset self (string->symbol name) value)))

    (define _set_signal_dict
      (lambda (self name d)
        (if (not (isinstance d dict))
            (raise (TypeError (format #f "~a must be a signal dict" d))))
	
        (for ((key : d)) ()
	     (if (not (in key _signals))
		 (raise (KeyError (format #f "~a is not a valid signal dict"
					  d)))))
        (for ((key : _signals)) ()
	     (if (not (in key d))
		 (raise (KeyError (format #f "~a is not a valid signal dict"
					  d)))))
        (rawset self (string->symbol name) d)))

    (define __setattr__
      (lambda (self name value)
        (cond
	 ((equal? name "prec")
          ((ref self '_set_integer_check) name value 1 "inf"))
	 ((equal? name "Emin")
	  ((ref self '_set_integer_check) name value "-inf" 0))
	 ((equal? name "Emax")
	  ((ref self '_set_integer_check) name value 0 "inf"))
	 ((equal? name "capitals")
	  ((ref self '_set_integer_check) name value 0 1))
	 ((equal? name "clamp")
	  ((ref self '_set_integer_check) name value 0 1))
	 ((equal? name "rounding")
	  (if (not (member value _rounding_modes))
	      ;; raise TypeError even for strings to have consistency
	      ;; among various implementations.
	      (raise (TypeError (format #f "~a: invalid rounding mode" value))))
	  (rawset self (string->symbol name) value))
	 ((or (equal? name "flags") (equal? name "traps"))
	  ((ref self '_set_signal_dict) name value))
	 ((equal? name "_ignored_flags")
	  (rawset self (string->symbol name) value))
	 (else
	  (raise (AttributeError
		  (format #f
			  "'decimal.Context' object has no attribute '~a'"
			  name)))))))

    (define __delattr__
      (lambda (self name)
        (raise (AttributeError (format #f "~a cannot be deleted" name)))))

    ;;# Support for pickling, copy, and deepcopy
    ;;def __reduce__(self):
    ;;    flags = [sig for sig, v in self.flags.items() if v]
    ;;    traps = [sig for sig, v in self.traps.items() if v]
    ;;    return (self.__class__,
    ;;            (self.prec, self.rounding, self.Emin, self.Emax,
    ;;             self.capitals, self.clamp, flags, traps))

    (define __repr__
      (lambda (self)
        "Show the current context."
	(format #f "Context(prec=~a, rounding=~a, Emin=~a, Emax=~a capitals=~a clamp=~a, flags=~a, traps=~a)"
		(ref self 'prec)
		(ref self 'rounding)
		(ref self 'Emin)
		(ref self 'Emax)
		(ref self 'capitals)
		(ref self 'clamp)
		(for ((k v : (ref self 'flags))) ((l '()))
		     (cons k l)
		     #:final (reverse l))
		(for ((k v : (ref self 'traps))) ((l '()))
		     (cons k l)
		     #:final (reverse l)))))
    
    (define clear_flags
      (lambda (self)
        "Reset all flags to zero"
        (for ((flag : (ref self 'flags))) ()
	     (pylist-set! (ref self 'flags) flag 0))))

    (define clear_traps
      (lambda (self)
        "Reset all traps to zero"
        (for ((flag : (ref self 'traps))) ()
	     (pylist-set! (ref self 'traps) flag 0))))

    (define _shallow_copy
      (lambda (self)
        "Returns a shallow copy from self."
        (Context (ref self 'prec)
		 (ref self 'rounding)
		 (ref self 'Emin)
		 (ref self 'Emax)
		 (ref self 'capitals)
		 (ref self 'clamp)
		 (ref self 'flags)
		 (ref self 'traps)
		 (ref self '_ignored_flags))))
    
    (define copy
      (lambda (self)
        "Returns a deep copy from self."
	(Context (ref self 'prec)
		 (ref self 'rounding)
		 (ref self 'Emin)
		 (ref self 'Emax)
		 (ref self 'capitals)
		 (ref self 'clamp)
		 ((ref (ref self 'flags) 'copy))
		 ((ref (ref self 'traps) 'copy))
		 (ref self '_ignored_flags))))

    (define __copy__ copy)

    (define _raise_error
      (lam (self condition (= explanation None) (* args))
        "Handles an error

        If the flag is in _ignored_flags, returns the default response.
        Otherwise, it sets the flag, then, if the corresponding
        trap_enabler is set, it reraises the exception.  Otherwise, it returns
        the default value after setting the flag.
        "
        (let ((error ((ref _condition_map 'get) condition condition)))
	  (if (in error (ref self '_ignored_flags))
	      ;; Don't touch the flag
	      (py-apply (ref (error) 'handle) self (* args))
	      (begin
		(pylist-set! (ref self 'flags) error 1)
		(if (not (bool (pylist-ref (ref self 'traps) error)))
		    ;; The errors define how to handle themselves.
		    (py-apply (ref (condition) 'handle) self (* args))
		    
		    ;; Errors should only be risked on copies of the context
		    ;; self._ignored_flags = []
		    (raise (error explanation))))))))

    (define _ignore_all_flags
      (lambda (self)
        "Ignore all flags, if they are raised"
        (py-apply (ref self '_ignore_flags) (* _signals))))

    (define _ignore_flags
      (lambda (self . flags)
        "Ignore the flags, if they are raised"
        ;; Do not mutate-- This way, copies of a context leave the original
        ;; alone.
        (set self '_ignored_flags (+ (ref self '_ignored_flags) (py-list flags)))
        (py-list flags)))

    (define _regard_flags
      (lambda (self . flags)
        "Stop ignoring the flags, if they are raised"
        (let ((flags
	       (if (and (pair? flags) (isinstance (car flags) (list tuple py-list)))
		   (car flags)
		   flags)))
	  (for ((flag : flags)) ()
	       ((ref (ref self '_ignored_flags) 'remove) flag)))))

    ;; We inherit object.__hash__, so we must deny this explicitly
    (define __hash__ None)

    (define Etiny
      (lambda (self)
        "Returns Etiny (= Emin - prec + 1)"
        (int (+ (- (ref self 'Emin) (ref self 'prec)) 1))))

    (define Etop
      (lambda (self)
        "Returns maximum exponent (= Emax - prec + 1)"
        (int (+ (- (ref self 'Emax) (ref self 'prec)) 1))))
    
    (define _set_rounding
      (lambda (self type)
        "Sets the rounding type.

        Sets the rounding type, and returns the current (previous)
        rounding type.  Often used like:

        context = context.copy()
        # so you don't change the calling context
        # if an error occurs in the middle.
        rounding = context._set_rounding(ROUND_UP)
        val = self.__sub__(other, context=context)
        context._set_rounding(rounding)

        This will make it round up for that operation.
        "
        (let ((rounding (ref self 'rounding)))
	  (set self 'rounding type)
	  rounding)))

    (define create_decimal
      (lam (self (= num "0"))
        "Creates a new Decimal instance but using self as context.

        This method implements the to-number operation of the
        IBM Decimal specification."

        (if (or (and (isinstance num str) (not (equal? num ((ref num 'strip)))))
		(in "_" num))
            ((ref self '_raise_error) ConversionSyntax
	     "trailing or leading whitespace and "
	     "underscores are not permitted.")
	    (let ((d (Decimal num #:context self)))
	      (if (and ((ref d '_isnan)) (> (len (ref d '_int)) (- (ref self 'prec)
								   (ref self 'clamp))))
		  ((ref self '_raise_error) ConversionSyntax
		   "diagnostic info too long in NaN")
		  ((ref d '_fix) self))))))

    (define create_decimal_from_float
      (lambda (self f)
        "Creates a new Decimal instance from a float but rounding using self
        as the context.

        >>> context = Context(prec=5, rounding=ROUND_DOWN)
        >>> context.create_decimal_from_float(3.1415926535897932)
        Decimal('3.1415')
        >>> context = Context(prec=5, traps=[Inexact])
        >>> context.create_decimal_from_float(3.1415926535897932)
        Traceback (most recent call last):
            ...
        decimal.Inexact: None

        "
        (let ((d ((ref Decimal 'from_float) f)))       ; An exact conversion
	  ((ref d '_fix) self))))                      ; Apply the context rounding

    ;; Methods
    (define abs
      (lambda (self a)
        "Returns the absolute value of the operand.

        If the operand is negative, the result is the same as using the minus
        operation on the operand.  Otherwise, the result is the same as using
        the plus operation on the operand.

        >>> ExtendedContext.abs(Decimal('2.1'))
        Decimal('2.1')
        >>> ExtendedContext.abs(Decimal('-100'))
        Decimal('100')
        >>> ExtendedContext.abs(Decimal('101.5'))
        Decimal('101.5')
        >>> ExtendedContext.abs(Decimal('-101.5'))
        Decimal('101.5')
        >>> ExtendedContext.abs(-1)
        Decimal('1')
        "
        (let ((a (_convert_other a #:raiseit #t)))
	  ((ref a '__abs__) #:context self))))

    (define add
      (lambda (self a b)
        "Return the sum of the two operands.

        >>> ExtendedContext.add(Decimal('12'), Decimal('7.00'))
        Decimal('19.00')
        >>> ExtendedContext.add(Decimal('1E+2'), Decimal('1.01E+4'))
        Decimal('1.02E+4')
        >>> ExtendedContext.add(1, Decimal(2))
        Decimal('3')
        >>> ExtendedContext.add(Decimal(8), 5)
        Decimal('13')
        >>> ExtendedContext.add(5, 5)
        Decimal('10')
        "
        (let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__add__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #f "Unable to convert ~a to Decimal" b)))
	      r))))

    (define _apply
      (lambda (self a)
        (str ((ref a '_fix) self))))

    (define canonical
      (lambda (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')
        "
        (if (not (isinstance a Decimal))
            (raise (TypeError "canonical requires a Decimal as an argument.")))

	((ref a 'canonical))))

    (define compare
      (lambda (self a b)
        "Compares values numerically.

        If the signs of the operands differ, a value representing each operand
        ('-1' if the operand is less than zero, '0' if the operand is zero or
        negative zero, or '1' if the operand is greater than zero) is used in
        place of that operand for the comparison instead of the actual
        operand.

        The comparison is then effected by subtracting the second operand from
        the first and then returning a value according to the result of the
        subtraction: '-1' if the result is less than zero, '0' if the result is
        zero or negative zero, or '1' if the result is greater than zero.

        >>> ExtendedContext.compare(Decimal('2.1'), Decimal('3'))
        Decimal('-1')
        >>> ExtendedContext.compare(Decimal('2.1'), Decimal('2.1'))
        Decimal('0')
        >>> ExtendedContext.compare(Decimal('2.1'), Decimal('2.10'))
        Decimal('0')
        >>> ExtendedContext.compare(Decimal('3'), Decimal('2.1'))
        Decimal('1')
        >>> ExtendedContext.compare(Decimal('2.1'), Decimal('-3'))
        Decimal('1')
        >>> ExtendedContext.compare(Decimal('-3'), Decimal('2.1'))
        Decimal('-1')
        >>> ExtendedContext.compare(1, 2)
        Decimal('-1')
        >>> ExtendedContext.compare(Decimal(1), 2)
        Decimal('-1')
        >>> ExtendedContext.compare(1, Decimal(2))
        Decimal('-1')
        "
        (let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'compare) b #:context self))))

    (define compare_signal
      (lambda (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
        >>> c.compare_signal(-1, 2)
        Decimal('-1')
        >>> c.compare_signal(Decimal(-1), 2)
        Decimal('-1')
        >>> c.compare_signal(-1, Decimal(2))
        Decimal('-1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'compare_signal) b #:context self))))

    (define compare_total
      (lambda (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')
        >>> ExtendedContext.compare_total(1, 2)
        Decimal('-1')
        >>> ExtendedContext.compare_total(Decimal(1), 2)
        Decimal('-1')
        >>> ExtendedContext.compare_total(1, Decimal(2))
        Decimal('-1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'compare_total) b #:context self))))

    (define compare_total_mag
      (lambda (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.
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'compare_total_mag) b #:context self))))

    (define copy_abs
      (lambda (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')
        >>> ExtendedContext.copy_abs(-1)
        Decimal('1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'copy_abs)))))
    
    (define copy_decimal
      (lambda (self a)
        "Returns a copy of the decimal object.

        >>> ExtendedContext.copy_decimal(Decimal('2.1'))
        Decimal('2.1')
        >>> ExtendedContext.copy_decimal(Decimal('-1.00'))
        Decimal('-1.00')
        >>> ExtendedContext.copy_decimal(1)
        Decimal('1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  (Decimal a))))
    
    (define copy_negate
      (lambda (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')
        >>> ExtendedContext.copy_negate(1)
        Decimal('-1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'copy_negate)))))
    
    (define copy_sign
      (lambda (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')
        >>> ExtendedContext.copy_sign(1, -2)
        Decimal('-1')
        >>> ExtendedContext.copy_sign(Decimal(1), -2)
        Decimal('-1')
        >>> ExtendedContext.copy_sign(1, Decimal(-2))
        Decimal('-1')
        "
	(let ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'copy_sign) b))))

    (define divide
      (lambda (self a b)
        "Decimal division in a specified context.

        >>> ExtendedContext.divide(Decimal('1'), Decimal('3'))
        Decimal('0.333333333')
        >>> ExtendedContext.divide(Decimal('2'), Decimal('3'))
        Decimal('0.666666667')
        >>> ExtendedContext.divide(Decimal('5'), Decimal('2'))
        Decimal('2.5')
        >>> ExtendedContext.divide(Decimal('1'), Decimal('10'))
        Decimal('0.1')
        >>> ExtendedContext.divide(Decimal('12'), Decimal('12'))
        Decimal('1')
        >>> ExtendedContext.divide(Decimal('8.00'), Decimal('2'))
        Decimal('4.00')
        >>> ExtendedContext.divide(Decimal('2.400'), Decimal('2.0'))
        Decimal('1.20')
        >>> ExtendedContext.divide(Decimal('1000'), Decimal('100'))
        Decimal('10')
        >>> ExtendedContext.divide(Decimal('1000'), Decimal('1'))
        Decimal('1000')
        >>> ExtendedContext.divide(Decimal('2.40E+6'), Decimal('2'))
        Decimal('1.20E+6')
        >>> ExtendedContext.divide(5, 5)
        Decimal('1')
        >>> ExtendedContext.divide(Decimal(5), 5)
        Decimal('1')
        >>> ExtendedContext.divide(5, Decimal(5))
        Decimal('1')
        "
        (let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__truediv__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))
    
    (define divide_int
      (lambda (self a b)
        "Divides two numbers and returns the integer part of the result.

        >>> ExtendedContext.divide_int(Decimal('2'), Decimal('3'))
        Decimal('0')
        >>> ExtendedContext.divide_int(Decimal('10'), Decimal('3'))
        Decimal('3')
        >>> ExtendedContext.divide_int(Decimal('1'), Decimal('0.3'))
        Decimal('3')
        >>> ExtendedContext.divide_int(10, 3)
        Decimal('3')
        >>> ExtendedContext.divide_int(Decimal(10), 3)
        Decimal('3')
        >>> ExtendedContext.divide_int(10, Decimal(3))
        Decimal('3')
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__floordiv__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))
    
    (define divmod
      (lambda (self a b)
        "Return (a // b, a % b).

        >>> ExtendedContext.divmod(Decimal(8), Decimal(3))
        (Decimal('2'), Decimal('2'))
        >>> ExtendedContext.divmod(Decimal(8), Decimal(4))
        (Decimal('2'), Decimal('0'))
        >>> ExtendedContext.divmod(8, 4)
        (Decimal('2'), Decimal('0'))
        >>> ExtendedContext.divmod(Decimal(8), 4)
        (Decimal('2'), Decimal('0'))
        >>> ExtendedContext.divmod(8, Decimal(4))
        (Decimal('2'), Decimal('0'))
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__divmod__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))

    (define exp
      (lambda (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')
        >>> c.exp(10)
        Decimal('22026.4658')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'exp) #:context self))))
    
    (define fma
      (lambda (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')
        >>> ExtendedContext.fma(1, 3, 4)
        Decimal('7')
        >>> ExtendedContext.fma(1, Decimal(3), 4)
        Decimal('7')
        >>> ExtendedContext.fma(1, 3, Decimal(4))
        Decimal('7')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'fma) b c #:context self))))

    (define is_canonical
      (lambda (self a)
        "Return True if the operand is canonical; otherwise return False.

        Currently, the encoding of a Decimal instance is always
        canonical, so this method returns True for any Decimal.

        >>> ExtendedContext.is_canonical(Decimal('2.50'))
        True
        "
        (if (not (isinstance a Decimal))
            (raise (TypeError "is_canonical requires a Decimal as an argument."))
	    ((ref a 'is_canonical)))))

    (define is_finite
      (lambda (self a)
        "Return True if the operand is finite; otherwise return False.

        A Decimal instance is considered finite if it is neither
        infinite nor a NaN.

        >>> ExtendedContext.is_finite(Decimal('2.50'))
        True
        >>> ExtendedContext.is_finite(Decimal('-0.3'))
        True
        >>> ExtendedContext.is_finite(Decimal('0'))
        True
        >>> ExtendedContext.is_finite(Decimal('Inf'))
        False
        >>> ExtendedContext.is_finite(Decimal('NaN'))
        False
        >>> ExtendedContext.is_finite(1)
        True
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'finite)))))

    (define is_infinite
      (lambda (self a)
        "Return True if the operand is infinite; otherwise return False.

        >>> ExtendedContext.is_infinite(Decimal('2.50'))
        False
        >>> ExtendedContext.is_infinite(Decimal('-Inf'))
        True
        >>> ExtendedContext.is_infinite(Decimal('NaN'))
        False
        >>> ExtendedContext.is_infinite(1)
        False
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_infinite)))))

    (define is_nan
      (lambda (self a)
        "Return True if the operand is a qNaN or sNaN;
        otherwise return False.

        >>> ExtendedContext.is_nan(Decimal('2.50'))
        False
        >>> ExtendedContext.is_nan(Decimal('NaN'))
        True
        >>> ExtendedContext.is_nan(Decimal('-sNaN'))
        True
        >>> ExtendedContext.is_nan(1)
        False
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_nan)))))

    (define is_normal
      (lambda (self a)
        "Return True if the operand is a normal number;
        otherwise return False.

        >>> c = ExtendedContext.copy()
        >>> c.Emin = -999
        >>> c.Emax = 999
        >>> c.is_normal(Decimal('2.50'))
        True
        >>> c.is_normal(Decimal('0.1E-999'))
        False
        >>> c.is_normal(Decimal('0.00'))
        False
        >>> c.is_normal(Decimal('-Inf'))
        False
        >>> c.is_normal(Decimal('NaN'))
        False
        >>> c.is_normal(1)
        True
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_normal) #:context self))))

    (define is_qnan
      (lambda (self a)
        "Return True if the operand is a quiet NaN; otherwise return False.

        >>> ExtendedContext.is_qnan(Decimal('2.50'))
        False
        >>> ExtendedContext.is_qnan(Decimal('NaN'))
        True
        >>> ExtendedContext.is_qnan(Decimal('sNaN'))
        False
        >>> ExtendedContext.is_qnan(1)
        False
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_qnan)))))

    (define is_signed
      (lambda (self a)
        "Return True if the operand is negative; otherwise return False.

        >>> ExtendedContext.is_signed(Decimal('2.50'))
        False
        >>> ExtendedContext.is_signed(Decimal('-12'))
        True
        >>> ExtendedContext.is_signed(Decimal('-0'))
        True
        >>> ExtendedContext.is_signed(8)
        False
        >>> ExtendedContext.is_signed(-8)
        True
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_signed)))))

        
    (define is_snan
      (lambda (self a)
        "Return True if the operand is a signaling NaN;
        otherwise return False.

        >>> ExtendedContext.is_snan(Decimal('2.50'))
        False
        >>> ExtendedContext.is_snan(Decimal('NaN'))
        False
        >>> ExtendedContext.is_snan(Decimal('sNaN'))
        True
        >>> ExtendedContext.is_snan(1)
        False
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_snan)))))


    (define is_subnormal
      (lambda (self a)
	"Return True if the operand is subnormal; otherwise return False.

        >>> c = ExtendedContext.copy()
        >>> c.Emin = -999
        >>> c.Emax = 999
        >>> c.is_subnormal(Decimal('2.50'))
        False
        >>> c.is_subnormal(Decimal('0.1E-999'))
        True
        >>> c.is_subnormal(Decimal('0.00'))
        False
        >>> c.is_subnormal(Decimal('-Inf'))
        False
        >>> c.is_subnormal(Decimal('NaN'))
        False
        >>> c.is_subnormal(1)
        False
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_subnormal) #:context self))))

    (define is_zero
      (lambda (self a)
        "Return True if the operand is a zero; otherwise return False.

        >>> ExtendedContext.is_zero(Decimal('0'))
        True
        >>> ExtendedContext.is_zero(Decimal('2.50'))
        False
        >>> ExtendedContext.is_zero(Decimal('-0E+2'))
        True
        >>> ExtendedContext.is_zero(1)
        False
        >>> ExtendedContext.is_zero(0)
        True
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'is_zero)))))

    (define ln
      (lambda (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')
        >>> c.ln(1)
        Decimal('0')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'ln) #:context self))))

    (define log10
      (lambda (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')
        >>> c.log10(0)
        Decimal('-Infinity')
        >>> c.log10(1)
        Decimal('0')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'log10) #:context self))))

    (define logb
      (lambda (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')
        >>> ExtendedContext.logb(1)
        Decimal('0')
        >>> ExtendedContext.logb(10)
        Decimal('1')
        >>> ExtendedContext.logb(100)
        Decimal('2')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'logb) #:context self))))

    (define logical_and
      (lambda (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')
        >>> ExtendedContext.logical_and(110, 1101)
        Decimal('100')
        >>> ExtendedContext.logical_and(Decimal(110), 1101)
        Decimal('100')
        >>> ExtendedContext.logical_and(110, Decimal(1101))
        Decimal('100')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'logical_and) b #:context self))))

    (define logical_invert
      (lambda (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')
        >>> ExtendedContext.logical_invert(1101)
        Decimal('111110010')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'logical_invert) #:context self))))

    (define logical_or
      (lambda (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')
        >>> ExtendedContext.logical_or(110, 1101)
        Decimal('1111')
        >>> ExtendedContext.logical_or(Decimal(110), 1101)
        Decimal('1111')
        >>> ExtendedContext.logical_or(110, Decimal(1101))
        Decimal('1111')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'logical_or) b #:context self))))
	
    (define logical_xor
      (lambda (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')
        >>> ExtendedContext.logical_xor(110, 1101)
        Decimal('1011')
        >>> ExtendedContext.logical_xor(Decimal(110), 1101)
        Decimal('1011')
        >>> ExtendedContext.logical_xor(110, Decimal(1101))
        Decimal('1011')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'logical_xor) b #:context self))))
	
    (define max
      (lambda (self a b)
        "max compares two values numerically and returns the maximum.

        If either operand is a NaN then the general rules apply.
        Otherwise, the operands are compared as though by the compare
        operation.  If they are numerically equal then the left-hand operand
        is chosen as the result.  Otherwise the maximum (closer to positive
        infinity) of the two operands is chosen as the result.

        >>> ExtendedContext.max(Decimal('3'), Decimal('2'))
        Decimal('3')
        >>> ExtendedContext.max(Decimal('-10'), Decimal('3'))
        Decimal('3')
        >>> ExtendedContext.max(Decimal('1.0'), Decimal('1'))
        Decimal('1')
        >>> ExtendedContext.max(Decimal('7'), Decimal('NaN'))
        Decimal('7')
        >>> ExtendedContext.max(1, 2)
        Decimal('2')
        >>> ExtendedContext.max(Decimal(1), 2)
        Decimal('2')
        >>> ExtendedContext.max(1, Decimal(2))
        Decimal('2')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'max) b #:context self))))
	
    (define max_mag
      (lambda (self a b)
        "Compares the values numerically with their sign ignored.

        >>> ExtendedContext.max_mag(Decimal('7'), Decimal('NaN'))
        Decimal('7')
        >>> ExtendedContext.max_mag(Decimal('7'), Decimal('-10'))
        Decimal('-10')
        >>> ExtendedContext.max_mag(1, -2)
        Decimal('-2')
        >>> ExtendedContext.max_mag(Decimal(1), -2)
        Decimal('-2')
        >>> ExtendedContext.max_mag(1, Decimal(-2))
        Decimal('-2')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'max_mag) b #:context self))))
	
    (define min
      (lambda (self a b)
        "min compares two values numerically and returns the minimum.

        If either operand is a NaN then the general rules apply.
        Otherwise, the operands are compared as though by the compare
        operation.  If they are numerically equal then the left-hand operand
        is chosen as the result.  Otherwise the minimum (closer to negative
        infinity) of the two operands is chosen as the result.

        >>> ExtendedContext.min(Decimal('3'), Decimal('2'))
        Decimal('2')
        >>> ExtendedContext.min(Decimal('-10'), Decimal('3'))
        Decimal('-10')
        >>> ExtendedContext.min(Decimal('1.0'), Decimal('1'))
        Decimal('1.0')
        >>> ExtendedContext.min(Decimal('7'), Decimal('NaN'))
        Decimal('7')
        >>> ExtendedContext.min(1, 2)
        Decimal('1')
        >>> ExtendedContext.min(Decimal(1), 2)
        Decimal('1')
        >>> ExtendedContext.min(1, Decimal(29))
        Decimal('1')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'min) b #:context self))))

    (define min_mag
      (lambda (self a b)
        "Compares the values numerically with their sign ignored.

        >>> ExtendedContext.min_mag(Decimal('3'), Decimal('-2'))
        Decimal('-2')
        >>> ExtendedContext.min_mag(Decimal('-3'), Decimal('NaN'))
        Decimal('-3')
        >>> ExtendedContext.min_mag(1, -2)
        Decimal('1')
        >>> ExtendedContext.min_mag(Decimal(1), -2)
        Decimal('1')
        >>> ExtendedContext.min_mag(1, Decimal(-2))
        Decimal('1')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'min_mag) b #:context self))))

    (define minus
      (lambda (self a)
        "Minus corresponds to unary prefix minus in Python.

        The operation is evaluated using the same rules as subtract; the
        operation minus(a) is calculated as subtract('0', a) where the '0'
        has the same exponent as the operand.

        >>> ExtendedContext.minus(Decimal('1.3'))
        Decimal('-1.3')
        >>> ExtendedContext.minus(Decimal('-1.3'))
        Decimal('1.3')
        >>> ExtendedContext.minus(1)
        Decimal('-1')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a '__neg__) #:context self))))

    (define multiply
      (lambda (self a b)
        "multiply multiplies two operands.

        If either operand is a special value then the general rules apply.
        Otherwise, the operands are multiplied together
        ('long multiplication'), resulting in a number which may be as long as
        the sum of the lengths of the two operands.

        >>> ExtendedContext.multiply(Decimal('1.20'), Decimal('3'))
        Decimal('3.60')
        >>> ExtendedContext.multiply(Decimal('7'), Decimal('3'))
        Decimal('21')
        >>> ExtendedContext.multiply(Decimal('0.9'), Decimal('0.8'))
        Decimal('0.72')
        >>> ExtendedContext.multiply(Decimal('0.9'), Decimal('-0'))
        Decimal('-0.0')
        >>> ExtendedContext.multiply(Decimal('654321'), Decimal('654321'))
        Decimal('4.28135971E+11')
        >>> ExtendedContext.multiply(7, 7)
        Decimal('49')
        >>> ExtendedContext.multiply(Decimal(7), 7)
        Decimal('49')
        >>> ExtendedContext.multiply(7, Decimal(7))
        Decimal('49')
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__mul__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))

    (define next_minus
      (lambda (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')
        >>> c.next_minus(1)
        Decimal('0.999999999')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'next_minus) #:context self))))

    (define next_plus
      (lambda (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')
        >>> c.next_plus(1)
        Decimal('1.00000001')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'next_plus) #:context self))))
    
    (define next_toward
      (lambda (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')
        >>> c.next_toward(0, 1)
        Decimal('1E-1007')
        >>> c.next_toward(Decimal(0), 1)
        Decimal('1E-1007')
        >>> c.next_toward(0, Decimal(1))
        Decimal('1E-1007')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'next_toward) b #:context self))))

    (define normalize
      (lambda (self a)
        "normalize reduces an operand to its simplest form.

        Essentially a plus operation with all trailing zeros removed from the
        result.

        >>> ExtendedContext.normalize(Decimal('2.1'))
        Decimal('2.1')
        >>> ExtendedContext.normalize(Decimal('-2.0'))
        Decimal('-2')
        >>> ExtendedContext.normalize(Decimal('1.200'))
        Decimal('1.2')
        >>> ExtendedContext.normalize(Decimal('-120'))
        Decimal('-1.2E+2')
        >>> ExtendedContext.normalize(Decimal('120.00'))
        Decimal('1.2E+2')
        >>> ExtendedContext.normalize(Decimal('0.00'))
        Decimal('0')
        >>> ExtendedContext.normalize(6)
        Decimal('6')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'normalize) #:context self))))

    (define number_class
      (lambda (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 = ExtendedContext.copy()
        >>> 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'
        >>> c.number_class(123)
        '+Normal'
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'number_class) #:context self))))

    (define plus
      (lambda (self a)
        "Plus corresponds to unary prefix plus in Python.

        The operation is evaluated using the same rules as add; the
        operation plus(a) is calculated as add('0', a) where the '0'
        has the same exponent as the operand.

        >>> ExtendedContext.plus(Decimal('1.3'))
        Decimal('1.3')
        >>> ExtendedContext.plus(Decimal('-1.3'))
        Decimal('-1.3')
        >>> ExtendedContext.plus(-1)
        Decimal('-1')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a '__pos__) #:context self))))

    (define power
      (lam (self a b (= modulo None))
        "Raises a to the power of b, to modulo if given.

        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.

        With three arguments, compute (a**b) % modulo.  For the
        three argument form, the following restrictions on the
        arguments hold:

         - 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

        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')
        >>> c.power(Decimal('-2'), Decimal('3'))
        Decimal('-8')
        >>> c.power(Decimal('2'), Decimal('-3'))
        Decimal('0.125')
        >>> c.power(Decimal('1.7'), Decimal('8'))
        Decimal('69.7575744')
        >>> c.power(Decimal('10'), Decimal('0.301029996'))
        Decimal('2.00000000')
        >>> c.power(Decimal('Infinity'), Decimal('-1'))
        Decimal('0')
        >>> c.power(Decimal('Infinity'), Decimal('0'))
        Decimal('1')
        >>> c.power(Decimal('Infinity'), Decimal('1'))
        Decimal('Infinity')
        >>> c.power(Decimal('-Infinity'), Decimal('-1'))
        Decimal('-0')
        >>> c.power(Decimal('-Infinity'), Decimal('0'))
        Decimal('1')
        >>> c.power(Decimal('-Infinity'), Decimal('1'))
        Decimal('-Infinity')
        >>> c.power(Decimal('-Infinity'), Decimal('2'))
        Decimal('Infinity')
        >>> 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')
        >>> ExtendedContext.power(7, 7)
        Decimal('823543')
        >>> ExtendedContext.power(Decimal(7), 7)
        Decimal('823543')
        >>> ExtendedContext.power(7, Decimal(7), 2)
        Decimal('1')
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__pow__) b modulo #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))

    (define quantize
      (lambda (self a b)
        "Returns a value equal to 'a' (rounded), having the exponent of 'b'.

        The coefficient of the result is derived from that of the left-hand
        operand.  It may be rounded using the current rounding setting (if the
        exponent is being increased), multiplied by a positive power of ten (if
        the exponent is being decreased), or is unchanged (if the exponent is
        already equal to that of the right-hand operand).

        Unlike other operations, if the length of the coefficient after the
        quantize operation would be greater than precision then an Invalid
        operation condition is raised.  This guarantees that, unless there is
        an error condition the exponent of the result of a quantize is always
        equal to that of the right-hand operand.

        Also unlike other operations, quantize will never raise Underflow, even
        if the result is subnormal and inexact.

        >>> ExtendedContext.quantize(Decimal('2.17'), Decimal('0.001'))
        Decimal('2.170')
        >>> ExtendedContext.quantize(Decimal('2.17'), Decimal('0.01'))
        Decimal('2.17')
        >>> ExtendedContext.quantize(Decimal('2.17'), Decimal('0.1'))
        Decimal('2.2')
        >>> ExtendedContext.quantize(Decimal('2.17'), Decimal('1e+0'))
        Decimal('2')
        >>> ExtendedContext.quantize(Decimal('2.17'), Decimal('1e+1'))
        Decimal('0E+1')
        >>> ExtendedContext.quantize(Decimal('-Inf'), Decimal('Infinity'))
        Decimal('-Infinity')
        >>> ExtendedContext.quantize(Decimal('2'), Decimal('Infinity'))
        Decimal('NaN')
        >>> ExtendedContext.quantize(Decimal('-0.1'), Decimal('1'))
        Decimal('-0')
        >>> ExtendedContext.quantize(Decimal('-0'), Decimal('1e+5'))
        Decimal('-0E+5')
        >>> ExtendedContext.quantize(Decimal('+35236450.6'), Decimal('1e-2'))
        Decimal('NaN')
        >>> ExtendedContext.quantize(Decimal('-35236450.6'), Decimal('1e-2'))
        Decimal('NaN')
        >>> ExtendedContext.quantize(Decimal('217'), Decimal('1e-1'))
        Decimal('217.0')
        >>> ExtendedContext.quantize(Decimal('217'), Decimal('1e-0'))
        Decimal('217')
        >>> ExtendedContext.quantize(Decimal('217'), Decimal('1e+1'))
        Decimal('2.2E+2')
        >>> ExtendedContext.quantize(Decimal('217'), Decimal('1e+2'))
        Decimal('2E+2')
        >>> ExtendedContext.quantize(1, 2)
        Decimal('1')
        >>> ExtendedContext.quantize(Decimal(1), 2)
        Decimal('1')
        >>> ExtendedContext.quantize(1, Decimal(2))
        Decimal('1')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'quantize) b #:context self))))

    (define radix
      (lambda (self)
        "Just returns 10, as this is Decimal, :)

        >>> ExtendedContext.radix()
        Decimal('10')
        "
        (Decimal 10)))

    (define remainder
      (lambda (self a b)
        "Returns the remainder from integer division.

        The result is the residue of the dividend after the operation of
        calculating integer division as described for divide-integer, rounded
        to precision digits if necessary.  The sign of the result, if
        non-zero, is the same as that of the original dividend.

        This operation will fail under the same conditions as integer division
        (that is, if integer division on the same two operands would fail, the
        remainder cannot be calculated).

        >>> ExtendedContext.remainder(Decimal('2.1'), Decimal('3'))
        Decimal('2.1')
        >>> ExtendedContext.remainder(Decimal('10'), Decimal('3'))
        Decimal('1')
        >>> ExtendedContext.remainder(Decimal('-10'), Decimal('3'))
        Decimal('-1')
        >>> ExtendedContext.remainder(Decimal('10.2'), Decimal('1'))
        Decimal('0.2')
        >>> ExtendedContext.remainder(Decimal('10'), Decimal('0.3'))
        Decimal('0.1')
        >>> ExtendedContext.remainder(Decimal('3.6'), Decimal('1.3'))
        Decimal('1.0')
        >>> ExtendedContext.remainder(22, 6)
        Decimal('4')
        >>> ExtendedContext.remainder(Decimal(22), 6)
        Decimal('4')
        >>> ExtendedContext.remainder(22, Decimal(6))
        Decimal('4')
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__mod__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))

    (define remainder_near
      (lambda (self a b)
        "Returns to be 'a - b * n', where n is the integer nearest the exact
        value of 'x / b' (if two integers are equally near then the even one
        is chosen).  If the result is equal to 0 then its sign will be the
        sign of a.

        This operation will fail under the same conditions as integer division
        (that is, if integer division on the same two operands would fail, the
        remainder cannot be calculated).

        >>> ExtendedContext.remainder_near(Decimal('2.1'), Decimal('3'))
        Decimal('-0.9')
        >>> ExtendedContext.remainder_near(Decimal('10'), Decimal('6'))
        Decimal('-2')
        >>> ExtendedContext.remainder_near(Decimal('10'), Decimal('3'))
        Decimal('1')
        >>> ExtendedContext.remainder_near(Decimal('-10'), Decimal('3'))
        Decimal('-1')
        >>> ExtendedContext.remainder_near(Decimal('10.2'), Decimal('1'))
        Decimal('0.2')
        >>> ExtendedContext.remainder_near(Decimal('10'), Decimal('0.3'))
        Decimal('0.1')
        >>> ExtendedContext.remainder_near(Decimal('3.6'), Decimal('1.3'))
        Decimal('-0.3')
        >>> ExtendedContext.remainder_near(3, 11)
        Decimal('3')
        >>> ExtendedContext.remainder_near(Decimal(3), 11)
        Decimal('3')
        >>> ExtendedContext.remainder_near(3, Decimal(11))
        Decimal('3')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'remainder_near) b #:context self))))

    (define rotate
      (lambda (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')
        >>> ExtendedContext.rotate(1333333, 1)
        Decimal('13333330')
        >>> ExtendedContext.rotate(Decimal(1333333), 1)
        Decimal('13333330')
        >>> ExtendedContext.rotate(1333333, Decimal(1))
        Decimal('13333330')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'rotate) b #:context self))))

    (define same_quantum
      (lambda (self a b)
        "Returns True if the two operands have the same exponent.

        The result is never affected by either the sign or the coefficient of
        either operand.

        >>> ExtendedContext.same_quantum(Decimal('2.17'), Decimal('0.001'))
        False
        >>> ExtendedContext.same_quantum(Decimal('2.17'), Decimal('0.01'))
        True
        >>> ExtendedContext.same_quantum(Decimal('2.17'), Decimal('1'))
        False
        >>> ExtendedContext.same_quantum(Decimal('Inf'), Decimal('-Inf'))
        True
        >>> ExtendedContext.same_quantum(10000, -1)
        True
        >>> ExtendedContext.same_quantum(Decimal(10000), -1)
        True
        >>> ExtendedContext.same_quantum(10000, Decimal(-1))
        True
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'same_quantum) b #:context self))))

    (define scaleb
      (lambda (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')
        >>> ExtendedContext.scaleb(1, 4)
        Decimal('1E+4')
        >>> ExtendedContext.scaleb(Decimal(1), 4)
        Decimal('1E+4')
        >>> ExtendedContext.scaleb(1, Decimal(4))
        Decimal('1E+4')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'scaleb) b #:context self))))

    (define shift
      (lambda (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')
        >>> ExtendedContext.shift(88888888, 2)
        Decimal('888888800')
        >>> ExtendedContext.shift(Decimal(88888888), 2)
        Decimal('888888800')
        >>> ExtendedContext.shift(88888888, Decimal(2))
        Decimal('888888800')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'shift) b #:context self))))

    (define sqrt
      (lambda (self a)
        "Square root of a non-negative number to context precision.

        If the result must be inexact, it is rounded using the round-half-even
        algorithm.

        >>> ExtendedContext.sqrt(Decimal('0'))
        Decimal('0')
        >>> ExtendedContext.sqrt(Decimal('-0'))
        Decimal('-0')
        >>> ExtendedContext.sqrt(Decimal('0.39'))
        Decimal('0.624499800')
        >>> ExtendedContext.sqrt(Decimal('100'))
        Decimal('10')
        >>> ExtendedContext.sqrt(Decimal('1'))
        Decimal('1')
        >>> ExtendedContext.sqrt(Decimal('1.0'))
        Decimal('1.0')
        >>> ExtendedContext.sqrt(Decimal('1.00'))
        Decimal('1.0')
        >>> ExtendedContext.sqrt(Decimal('7'))
        Decimal('2.64575131')
        >>> ExtendedContext.sqrt(Decimal('10'))
        Decimal('3.16227766')
        >>> ExtendedContext.sqrt(2)
        Decimal('1.41421356')
        >>> ExtendedContext.prec
        9
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'sqrt) #:context self))))

    (define subtract
      (lambda (self a b)
        "Return the difference between the two operands.

        >>> ExtendedContext.subtract(Decimal('1.3'), Decimal('1.07'))
        Decimal('0.23')
        >>> ExtendedContext.subtract(Decimal('1.3'), Decimal('1.30'))
        Decimal('0.00')
        >>> ExtendedContext.subtract(Decimal('1.3'), Decimal('2.07'))
        Decimal('-0.77')
        >>> ExtendedContext.subtract(8, 5)
        Decimal('3')
        >>> ExtendedContext.subtract(Decimal(8), 5)
        Decimal('3')
        >>> ExtendedContext.subtract(8, Decimal(5))
        Decimal('3')
        "
	(let* ((a (_convert_other a #:raiseit #t))
	       (r ((ref a '__sub__) b #:context self)))
	  (if (equal? r NotImplemented)
	      (raise (TypeError (format #t "Unable to convert ~a to Decimal" b)))
	      r))))

    (define to_eng_string
      (lambda (self a)
        "Convert to a string, using engineering notation if an exponent is needed.

        Engineering notation has an exponent which is a multiple of 3.  This
        can leave up to 3 digits to the left of the decimal place and may
        require the addition of either one or two trailing zeros.

        The operation is not affected by the context.

        >>> ExtendedContext.to_eng_string(Decimal('123E+1'))
        '1.23E+3'
        >>> ExtendedContext.to_eng_string(Decimal('123E+3'))
        '123E+3'
        >>> ExtendedContext.to_eng_string(Decimal('123E-10'))
        '12.3E-9'
        >>> ExtendedContext.to_eng_string(Decimal('-123E-12'))
        '-123E-12'
        >>> ExtendedContext.to_eng_string(Decimal('7E-7'))
        '700E-9'
        >>> ExtendedContext.to_eng_string(Decimal('7E+1'))
        '70'
        >>> ExtendedContext.to_eng_string(Decimal('0E+1'))
        '0.00E+3'

        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'to_eng_string) #:context self))))

    (define to_sci_string
      (lambda (self a)
        "Converts a number to a string, using scientific notation.

        The operation is not affected by the context.
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a '__str__) #:context self))))

    (define to_integral_exact
      (lambda (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')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'to_integral_exact) #:context self))))
  
    (define to_integral_value
      (lambda (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, except that no flags will
        be set.  The rounding mode is taken from the context.

        >>> ExtendedContext.to_integral_value(Decimal('2.1'))
        Decimal('2')
        >>> ExtendedContext.to_integral_value(Decimal('100'))
        Decimal('100')
        >>> ExtendedContext.to_integral_value(Decimal('100.0'))
        Decimal('100')
        >>> ExtendedContext.to_integral_value(Decimal('101.5'))
        Decimal('102')
        >>> ExtendedContext.to_integral_value(Decimal('-101.5'))
        Decimal('-102')
        >>> ExtendedContext.to_integral_value(Decimal('10E+5'))
        Decimal('1.0E+6')
        >>> ExtendedContext.to_integral_value(Decimal('7.89E+77'))
        Decimal('7.89E+77')
        >>> ExtendedContext.to_integral_value(Decimal('-Inf'))
        Decimal('-Infinity')
        "
	(let* ((a (_convert_other a #:raiseit #t)))
	  ((ref a 'to_integral_value) #:context self))))

    ;; the method name changed, but we provide also the old one, for compatibility
    (define to_integral to_integral_value))

(define-python-class _WorkRep ()
  (define __init__
    (lam (self (= value None))
	 (cond
	  ((eq? value None)
	   (set self 'sign None)
	   (set self 'int  0)
	   (set self 'exp  None))
	  ((isinstance value Decimal)
	   (set self 'sign  (ref value '_sign))
	   (set self 'int   (int (ref value '_int)))
	   (set self 'exp   (ref value '_exp)))
	  (else
	   ;; assert isinstance(value, tuple)
	   (set self 'sign (pylist-ref value 0))
	   (set self 'int  (pylist-ref value 1))
	   (set self 'exp  (pylist-ref value 2))))))
  
  (define __repr__
    (lambda (self)
        (format #f "(~a, ~a, ~a)" (ref self 'sign) (ref self 'int) (ref self 'exp))))

  (define __str__  __repr__))


(define _normalize
  (lam (op1 op2 (= prec 0))
    "Normalizes op1, op2 to have the same exp and length of coefficient.

    Done during addition.
    "
    (call-with-values
	(lambda ()
	  (if (< (ref op1 'exp) (ref op2 'exp))
	      (values op2 op1)
	      (values op1 op2)))
      (lambda (tmp other)
	;; 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.
	(let* ((tmp_len   (len (str (ref tmp   'int))))
	       (other_len (len (str (ref other 'int))))
	       (exp       (+ (ref tmp 'exp)
			     (min -1  (- tmp_len prec 2)))))
	  (when (< (+ other_len (ref other 'exp) -1) exp)
	    (set other 'int 1)
	    (set other 'exp exp))

	  (set tmp 'int (* (ref tmp 'int) (expt 10 (- (ref tmp 'exp) (ref other 'exp)))))
	  (set tmp 'exp (ref other 'exp))
	  (values op1 op2))))))


;;##### Integer arithmetic functions used by ln, log10, exp and __pow__ #####

(define _nbits (ref int 'bit_length))

(define _decimal_lshift_exact
  (lambda (n e)
    " Given integers n and e, return n * 10**e if it's an integer, else None.

    The computation is designed to avoid computing large powers of 10
    unnecessarily.

    >>> _decimal_lshift_exact(3, 4)
    30000
    >>> _decimal_lshift_exact(300, -999999999)  # returns None

    "
    (cond
     ((= n 0)
      0)
     ((>= e 0)
      (* n (expt 10 e)))
     (else
      ;; val_n = largest power of 10 dividing n.
      (let* ((str_n (str (abs n)))
	     (val_n (- (len str_n) (len ((ref str_n 'rstrip) "0")))))
	(if (< val_n (- e))
	    None
	    (floor-quotient n (expt 10 (- e)))))))))

(define _sqrt_nearest
  (lambda (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.

    "
    (if (or (<= n 0) (<= a 0))
        (raise (ValueError "Both arguments to _sqrt_nearest should be positive.")))

    (let lp ((b 0) (a a))
      (if (not (= a b))
	  (lp a (ash (- a (floor-quotient (- n) a)) -1))
	  a))))

(define _rshift_nearest
  (lambda (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.

    "
    (let ((b (ash 1 shift))
	  (q (ash x (- shift))))
      (+ q (if (> (+ (* 2 (logand x (- b 1))) (logand  q 1)) b) 1 0)))))

(define _div_nearest
  (lambda (a b)
    "Closest integer to a/b, a and b positive integers; rounds to even
    in the case of a tie.

    "
    (call-with-values (lambda () (py-divmod a b))
      (lambda (q r)
	(+ q  (if (> (+ (* 2 r) (logand q 1)) b) 1 0))))))

(define _ilog
  (lam (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.

    ;; argument reduction; R = number of reductions performed
    (call-with-values
	(lambda ()
	  (let lp ((y (- x M)) (R 0))
	    (if (>= (ash (abs y) (- L R)) M)
		(values
		 (_div_nearest (ash (* M y) 1)
			       (+ M (_sqrt_nearest (* M (+ M (_rshift_nearest y R))) M)))
		 (+ R 1))
		(values y R))))
      (lambda (y R)
	;; Taylor series with T terms
	(let* ((T      (- (int (* -10 (floor-quotient (len (str M)) (* 3 L))))))
	       (yshift (_rshift_nearest y R))
	       (w      (_div_nearest M T)))
	  (for ((k : (range (- T 1) 0 -1))) ((w w))
	       (- (_div_nearest M k) (_div_nearest (* yshift w) M))
	       #:final
	       (_div_nearest (* w y) M)))))))

(define _dlog10
  (lambda (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
    (set! p (+ 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
    (call-with-values
	(lambda ()
	  (let* ((l (len (str c)))
		 (f (- (+ e l) (if (>= (+ e l) 1) 1 0))))	    
	    (if (> p 0)
		(let* ((M      (expt 10 p))
		       (k      (- (+ e p) f))
		       (c      (if (>= k 0)
				   (* c (expt 10 k))
				   (_div_nearest c  (expt 10 (- k)))))
		       (log_d  (_ilog c M))        ;; error < 5 + 22 = 27
		       (log_10 (_log10_digits p))) ;; error < 1
		  (values (_div_nearest (* log_d M) log_10)
			  (* f M))) ;; exact
		(values 0  ;; error < 2.31
			(_div_nearest f  (expt 10 (- p))))))) ;; error < 0.5
      (lambda (log_d log_tenpower)
	(_div_nearest (+ log_tenpower log_d) 100)))))

(define _dlog
  (lambda (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.
    (set! p (+ 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).
    (let* ((l (len (str c)))
	   (f (- (+ e l) (if (>= (+ e l) 1) 1 0))))

      ;; compute approximation to 10**p*log(d), with error < 27
      (call-with-values
	  (lambda ()
	    (if (> p 0)
		(let* ((k (- (+ e p) f))
		       (c (if (>= k 0)
			      (* c (expt 10 k))
			      (_div_nearest c (expt 10 (- k))))))  ; error of <= 0.5 in c

		  ;; _ilog magnifies existing error in c by a factor of at most 10
		  (_ilog c (expt 10 p))) ; error < 5 + 22 = 27
		;; p <= 0: just approximate the whole thing by 0; error < 2.31
		0))
	(lambda (log_d)
	  (call-with-values
	      (lambda ()
		;; compute approximation to f*10**p*log(10), with error < 11.
		(if (not (= f 0))
		    (let ((extra (- (len (str (abs f))) 1)))
		      (if (>= (+ p extra) 0)
			  ;; error in f * _log10_digits(p+extra) < |f| * 1 = |f|
			  ;; after division error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
			  (_div_nearest (* f (_log10_digits (+ p extra))) (expt 10 extra))
			  0))
		    0))
	    (lambda (f_log_ten)
	      ;; error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1
	      (_div_nearest (+ f_log_ten log_d) 100))))))))

(define-python-class _Log10Memoize ()
    "Class to compute, store, and allow retrieval of, digits of the
    constant log(10) = 2.302585....  This constant is needed by
    Decimal.ln, Decimal.log10, Decimal.exp and Decimal.__pow__."
    
    (define __init__
      (lambda (self)
        (set self 'digits "23025850929940456840179914546843642076011014886")))

    (define getdigits
      (lambda (self p)
        "Given an integer p >= 0, return floor(10**p)*log(10).

        For example, self.getdigits(3) returns 2302.
        "
        ;; digits are stored as a string, for quick conversion to
        ;; integer in the case that we've already computed enough
        ;; digits; the stored digits should always be correct
        ;; (truncated, not rounded to nearest).
        (if (< p 0) (raise (ValueError "p should be nonnegative")))

        (if (>= p (len (ref self 'digits)))
            ;; compute p+3, p+6, p+9, ... digits; continue until at
            ;; least one of the extra digits is nonzero
	    (begin
	      (let lp ((extra 3))
		;; compute p+extra digits, correct to within 1ulp
		(let* ((M       (expt 10 (+ p extra 2)))
		       (digits  (str (_div_nearest (_ilog (* 10 M) M) 100))))
		  (if (not (equal? (pylist-slice (ref self 'digits) (- extra) None None)
				   (* '0' extra)))
		      #t
		      (lp (+ extra  3))))))
	    ;; keep all reliable digits so far; remove trailing zeros
	    ;; and next nonzero digit
	    (set self 'digits (pylist-slice (py-rstrip (ref self 'digits)) "0") None -1 None))
	
        (int (pylist-slice (ref self 'digits) None (+ p 1) None)))))

(define _log10_digits (ref (_Log10Memoize) 'getdigits))

(define _iexp
  (lam (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
    (let ((R (_nbits (floor-quotient (ash x L) M))))
      ;; Taylor series.  (2**L)**T > M
      (let* ((T  (- (int (floor-quotient (* -10 (len (str M))) (* 3 L)))))
	     (y1 (let ((Mshift (ash M R)))
		   (for ((i : (range (- T 1) 0 -1))) ((y (_div_nearest x T)))
			(_div_nearest (* x (+ Mshift y)) (* Mshift i))
			#:final y)))
				   

	     ;; Expansion
	     (y2 (for ((k : (range (- R 1) -1 -1))) ((y y1))
		      (let ((Mshift (ash M (+ k 2))))
			(_div_nearest (* y (+ y Mshift)) Mshift))
		      #:final y)))

	(+ M y2)))))

(define _dexp
  (lambda (c e p)
    "Compute an approximation to exp(c*10**e), with p decimal places of
    precision.

    Returns integers 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
    (set! p (+ p 2))

    ;; compute log(10) with extra precision = adjusted exponent of c*10**e
    (let* ((extra (max 0 (+ e (len (str c)) -1)))
	   (q     (+ p extra)))

      ;; compute quotient c*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q),
      ;; rounding down
      (let* ((shift  (+ e q))
	     (cshift (if (>= shift 0)
			 (* c (expt 10 shift))
			 (floor-quotient c (expt 10 (- shift))))))
	(call-with-values
	    (lambda ()
	      (divmod cshift (_log10_digits q)))
	  (lambda (quot rem)
	    ;; reduce remainder back to original precision
	    (set! rem (_div_nearest rem (expt 10 extra)))

	    ;; error in result of _iexp < 120;  error after division < 0.62
	    (values (_div_nearest (_iexp rem (expt 10 p)) 1000)
		    (+ quot (- p) 3))))))))

(define _dpower
  (lambda (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.
    "

    (let*
	;; 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))
	 (ps    (if (>= shift 0)
		    (* lxc yc (expt 10 shift))
		    (_div_nearest (* lxc yc) (expt 10 (- shift))))))

      (if (= ps 0)
	  ;; we prefer a result that isn't exactly 1; this makes it
	  ;; easier to compute a correctly rounded result in __pow__
	  (if (eq? (>= (+ (len (str xc)) xe) 1)
		   (>  yc 0)) ; if x**y > 1:
	      (values (+ (expt 10 (- p 1)) 1) (- 1 p))
	      (values (- (expt 10 p) 1) (- p)))
	  (call-with-values
	      (lambda ()
		(_dexp ps (- (+ p 1)) (+ p 1)))
	    (lambda (coeff exp)
	      (values (_div_nearest coeff 10)
		      (+ exp 1))))))))
(pk 3)
(define _corr (dict '(("1" . 100) ("2" . 70) ("3" . 53) ("4" . 40) ("5" . 31)
		      ("6" . 23 ) ("7" . 16) ("8" . 10) ("9" . 5))))
(define _log10_lb
  (lam (c (= correction _corr))
       "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.")))
    (let ((str_c (str c)))
      (- (* 100 (len str_c) (pylist-ref correction (pylist-ref str_c 0)))))))

;;#### Helper Functions ####################################################

(define _convert_other
  (lam (other (= raiseit #f) (= allow_float #f))
    "Convert other to Decimal.

    Verifies that it's ok to use in an implicit construction.
    If allow_float is true, allow conversion from float;  this
    is used in the comparison methods (__eq__ and friends).

    "
    (cond
     ((isinstance other  Decimal)
      other)
     ((isinstance other int)
      (Decimal other))
     ((and allow_float (isinstance other float))
      ((ref Decimal 'from_float) other))
     (raiseit
      (raise (TypeError (format #f "Unable to convert  ~a to Decimal" other))))
     (else
      NotImplemented))))

(define _convert_for_comparison
  (lam (self other (= equality_op #f))
    "Given a Decimal instance self and a Python object other, return
    a pair (s, o) of Decimal instances such that 's op o' is
    equivalent to 'self op other' for any of the 6 comparison
    operators 'op'.

    "
    (cond
     ((isinstance other Decimal)
      (values self other))

     ;; Comparison with a Rational instance (also includes integers):
     ;; self op n/d <=> self*d op n (for n and d integers, d positive).
     ;; A NaN or infinity can be left unchanged without affecting the
     ;; comparison result.
     ((isinstance other int)
      (if (not (bool (ref self '_is_special)))
          (values
	   (_dec_from_triple (ref self '_sign)
			     (* (str int (ref self '_int)) (ref other 'denominator))
			     (ref self '_exp))
	   (Decimal (ref other 'numerator)))
	  (values NotImplemented NotImplemented)))
	  
     ;; Comparisons with float and complex types.  == and != comparisons
     ;; with complex numbers should succeed, returning either True or False
     ;; as appropriate.  Other comparisons return NotImplemented.
     (else
      (let ((other (if (and equality_op
			    (isinstance other complex)
			    (= (ref other 'imag) 0))
		       (ref other 'real)
		       other)))
	(if (isinstance other float)
	    (let ((context (getcontext)))
	      (if equality_op
		  (pylist-set! (ref context 'flags) FloatOperation 1)
		  ((cx-error context) FloatOperation
		   "strict semantics for mixing floats and Decimals are enabled"))
	      (values self ((ref Decimal 'from_float) other)))
	    (values NotImplemented NotImplemented)))))))


;;##### Setup Specific Contexts ############################################

;; The default context prototype used by Context()
;; Is mutable, so that new contexts can have different default values

(define DefaultContext
  (Context
   #:prec     28
   #:rounding ROUND_HALF_EVEN
   #:traps    (list DivisionByZero Overflow InvalidOperation)
   #:flags    '()
   #:Emax     999999
   #:Emin    -999999
   #:capitals 1
   #:clamp    0))
(pk 3 1)
;; Pre-made alternate contexts offered by the specification
;; Don't change these; the user should be able to select these
;; contexts and be able to reproduce results from other implementations
;; of the spec.

(define BasicContext
  (Context
   #:prec     9
   #:rounding ROUND_HALF_UP
   #:traps    (list DivisionByZero Overflow InvalidOperation Clamped Underflow)
   #:flags    '()))
(pk 3 2)
(define ExtendedContext
  (Context
   #:prec     9
   #:rounding ROUND_HALF_EVEN
   #:traps    '()
   #:flags    '()))

;;##### crud for parsing strings #############################################
;;#
;;# Regular expression used for parsing numeric strings.  Additional
;;# comments:
;;#
;;# 1. Uncomment the two '\s*' lines to allow leading and/or trailing
;;# whitespace.  But note that the specification disallows whitespace in
;;# a numeric string.
;;#
;;# 2. For finite numbers (not infinities and NaNs) the body of the
;;# number between the optional sign and the optional exponent must have
;;# at least one decimal digit, possibly after the decimal point.  The
;;# lookahead expression '(?=\d|\.\d)' checks this.
(pk 5)
(define _parser
  (ref
   (compile "        # A numeric string consists of:
#    \\s*
    (?P<sign>[-+])?              # an optional sign, followed by either...
    (
        (?=\\d|\\.\\d)           # ...a number (with at least one digit)
        (?P<int>\\d*)            # having a (possibly empty) integer part
        (\\.(?P<frac>\\d*))?      # followed by an optional fractional part
        (E(?P<exp>[-+]?\\d+))?   # followed by an optional exponent, or...
    |
        Inf(inity)?              # ...an infinity, or...
    |
        (?P<signal>s)?           # ...an (optionally signaling)
        NaN                      # NaN
        (?P<diag>\\d*)           # with (possibly empty) diagnostic info.
    )
#    \\s*
    \\Z
" (logior VERBOSE IGNORECASE))
   'match))
(pk 6)
(define _all_zeros  (ref (compile "0*$" ) 'match))
(pk 7)
(define _exact_half (ref (compile "50*$") 'match))
(pk 8)
;;##### PEP3101 support functions ##############################################
;;# The functions in this section have little to do with the Decimal
;;# class, and could potentially be reused or adapted for other pure
;;# Python numeric classes that want to implement __format__
;;#
;;# A format specifier for Decimal looks like:
;;#
;;#   [[fill]align][sign][#][0][minimumwidth][,][.precision][type]

(define _parse_format_specifier_regex
  (compile "\\A
(?:
   (?P<fill>.)?
   (?P<align>[<>=^])
)?
(?P<sign>[-+ ])?
(?P<alt>\\#)?
(?P<zeropad>0)?
(?P<minimumwidth>(?!0)\\d+)?
(?P<thousands_sep>,)?
(?:\\.(?P<precision>0|(?!0)\\d+))?
(?P<type>[eEfFgGn%])?
\\Z
" (logior VERBOSE DOTALL)))
(pk 9)
;; The locale module is only needed for the 'n' format specifier.  The
;; rest of the PEP 3101 code functions quite happily without it, so we
;; don't care too much if locale isn't present.
(define _locale (Module "locale"))
(pk 10)
(define _parse_format_specifier
  (lam (format_spec (= _localeconv None))
    "Parse and validate a format specifier.

    Turns a standard numeric format specifier into a dict, with the
    following entries:

      fill: fill character to pad field to minimum width
      align: alignment type, either '<', '>', '=' or '^'
      sign: either '+', '-' or ' '
      minimumwidth: nonnegative integer giving minimum width
      zeropad: boolean, indicating whether to pad with zeros
      thousands_sep: string to use as thousands separator, or ''
      grouping: grouping for thousands separators, in format
        used by localeconv
      decimal_point: string to use for decimal point
      precision: nonnegative integer giving precision or None
      type: one of the characters 'eEfFgG%', or None

    "
    (let* ((m (let ((m ((ref _parse_format_specifier_regex 'match) format_spec)))
		(if (eq? m None)
		    (raise (ValueError (+ "Invalid format specifier: " format_spec))))
		m))

	   ;; get the dictionary
	   (format_dict ((ref m 'groupdict)))

	   ;; zeropad; defaults for fill and alignment.  If zero padding
	   ;; is requested, the fill and align fields should be absent.
	   (fill  (pylist-ref format_dict "fill"))
	   (minw  (pylist-ref format_dict "minimumwidth"))
	   (sign  (pylist-ref format_dict "sign"))
	   (prec  (pylist-ref format_dict "precition"))
	   (sepM  (pylist-ref format_dict "thousands_sep"))
	   (type  (pylist-ref format_dict "type"))
	   (align (pylist-ref format_dict "align")))

      (pylist-set! format_dict "zeropad" (not (eq? (pylist-ref format_dict "zeropad") None)))
      (when (pylist-ref format_dict "zeropad")
	 (if (not (eq? fill  None))
	     (raise (ValueError (+ "Fill character conflicts with '0'"
				   " in format specifier: " format_spec))))
	 (if (not (eq? align  None))
	     (raise (ValueError (+ "Alignment conflicts with '0' in "
				   "format specifier: " format_spec)))))

      (pylist-set! format_dict "fill" (or (bool fill) " "))
      
      ;; PEP 3101 originally specified that the default alignment should
      ;; be left;  it was later agreed that right-aligned makes more sense
      ;; for numeric types.  See http://bugs.python.org/issue6857.

      (pylist-set! format_dict "align" (or (bool align) ">"))

      ;; default sign handling: '-' for negative, '' for positive
      (pylist-set! format_dict "sign" (or (bool sign) "-"))
      
      ;; minimumwidth defaults to 0; precision remains None if not given
      (pylist-set! format_dict "minimumwidth" (int (if (eq? minw None) "0" minw)))
      (if (not (eq? prec None))
	  (pylist-set! format_dict "precision" (let ((w (int prec))) (set! prec w) w)))

      ;; if format type is 'g' or 'G' then a precision of 0 makes little
      ;; sense; convert it to 1.  Same if format type is unspecified.
      (if (equal? prec 0)
	  (if (or (eq? type None) (in type "gGn"))
	      (pylist-set! format_dict "precision" 1)))

      ;; determine thousands separator, grouping, and decimal separator, and
      ;; add appropriate entries to format_dict
      (if (equal? type "n")
	  (begin
	    ;; apart from separators, 'n' behaves just like 'g'
	    (pylist-set! format_dict "type" "g")
	    (if (eq? _localeconv None)
		(set! _localeconv ((ref _locale 'localeconv))))
	    (if (not (eq? sepM None))
		(raise (ValueError (+ "Explicit thousands separator conflicts with "
				      "'n' type in format specifier: " format_spec))))
            (pylist-set! format_dict "thousands_sep" (pylist-ref _localeconv "thousands_sep"))
	    (pylist-set! format_dict "grouping"      (pylist-ref _localeconv "grouping"))
	    (pylist-set! format_dict "decimal_point" (pylist-ref _localeconv "decimal_point")))
	  (begin
	    (if (eq? sepM None)
		(pylist-set! format_dict "thousands_sep" ""))
	    (pylist-set! format_dict "grouping"      (list 3 0))
	    (pylist-set! format_dict "decimal_point" ".")))

      format_dict)))

(define _format_align
  (lambda (sign body spec)
    "Given an unpadded, non-aligned numeric string 'body' and sign
    string 'sign', add padding and alignment conforming to the given
    format specifier dictionary 'spec' (as produced by
    parse_format_specifier).

    "
    ;; how much extra space do we have to play with?
    (let* ((minimumwidth (pylist-ref spec "minimumwidth"))
	   (fill         (pylist-ref spec "fill"))
	   (padding      (* fill (- minimumwidth (len sign) (len body))))
	   (align        (pylist-ref spec "align")))
      (cond
       ((equal? align "<")
        (+ sign body padding))
       ((equal? align ">")
        (+ padding sign body))
       ((equal? align "=")
        (+ sign padding body))
       ((equal? align "^")
        (let* ((half (floor-quotient (len padding) 2))
	       (pad1  (pylist-slice padding None half None))
	       (pad2  (pylist-slice padding half None None)))
	  (+ pad1 sign  body pad2)))
       (else
        (raise (ValueError "Unrecognised alignment field")))))))

(define _group_lengths
  (lambda (grouping)
    "Convert a localeconv-style grouping into a (possibly infinite)
    iterable of integers representing group lengths.

    "
    ;; The result from localeconv()['grouping'], and the input to this
    ;; function should be a list of integers in one of the
    ;; following three forms:
    ;;
    ;;   (1) an empty list, or
    ;;   (2) nonempty list of positive integers + [0]
    ;;   (3) list of positive integers + [locale.CHAR_MAX], or

    (cond
     ((not (bool grouping))
      '())
     ((and (= (pylist-ref grouping -1) 0) (>= (len grouping) 2))
      (chain (pylist-slice grouping None -1 None) (repeat (pylist-ref grouping -2))))
     ((= (pylist-ref grouping -1) (ref _locale 'CHAR_MAX))
      (pylist-slice grouping None -1 None))
     (else
      (raise (ValueError "unrecognised format for grouping"))))))

(define _insert_thousands_sep
  (lam (digits spec (= min_width 1))
    "Insert thousands separators into a digit string.

    spec is a dictionary whose keys should include 'thousands_sep' and
    'grouping'; typically it's the result of parsing the format
    specifier using _parse_format_specifier.

    The min_width keyword argument gives the minimum length of the
    result, which will be padded on the left with zeros if necessary.

    If necessary, the zero padding adds an extra '0' on the left to
    avoid a leading thousands separator.  For example, inserting
    commas every three digits in '123456', with min_width=8, gives
    '0,123,456', even though that has length 9.

    "

    (let ((sep       (pylist-ref spec "thousands_sep"))
	  (grouping  (pylist-ref spec "grouping"))
	  (groups    (pylist)))

      (for ((l : (_group_lengths grouping))) ()
	   (if (<= l 0)
	       (raise (ValueError "group length should be positive")))
	   ;; max(..., 1) forces at least 1 digit to the left of a separator
	   (let ((l (min (max (len digits) min_width 1) l)))
	     ((ref groups 'append) (+ (* '0' (- l (len digits)))
				      (pylist-slice digits (- l) None None)))
	     (set! digits (pylist-slice digits None (- l) None))
	     (set! min_width (- min_width l))
	     (if (and (= 0 digits) (<= min_width  0))
		 (break))
	     (set! min_width (- min_width (len sep))))
	   #:final
	   (let ((l (max (len digits) min_width 1)))
	     ((ref groups 'append) (+ (* "0" (- l (len digits)))
				      (pylist-slice digits (- l) None None)))))
      ((ref sep 'join) (reversed groups)))))

(define _format_sign
  (lam (is_negative spec)
    "Determine sign character."

    (cond
     ((bool is_negative)
      "-")
     ((in (pylist-ref spec "sign") " +")
      (pylist-ref spec "sign"))
     (else
      ""))))
(pk 11)
(define typed (dict '(("E" .  "E") ("e" . "e") ("G" . "E") ("g" . "e"))))
(pk 12)
(define _format_number
  (lambda (is_negative intpart fracpart exp spec)
    "Format a number, given the following data:

    is_negative: true if the number is negative, else false
    intpart: string of digits that must appear before the decimal point
    fracpart: string of digits that must come after the point
    exp: exponent, as an integer
    spec: dictionary resulting from parsing the format specifier

    This function uses the information in spec to:
      insert separators (decimal separator and thousands separators)
      format the sign
      format the exponent
      add trailing '%' for the '%' type
      zero-pad if necessary
      fill and align if necessary
    "

    (let ((sign  (_format_sign is_negative  spec)))

      (if (or (bool fracpart) (bool (pylist-ref spec "alt")))
	  (set! fracpart (+ (pylist-ref spec "decimal_point") fracpart)))

      (if (or (not (= exp 0)) (in (pylist-ref spec "type") "eEgG"))
	  (let ((echar (pylist-ref typed (pylist-ref spec "type")))) 
	    (set! fracpart (+ fracpart (str-format "{0}{1:+}" echar exp)))))
    
      (if (equal? (pylist-ref spec "type") "%")
	  (set! fracpart (+ fracpart "%")))

      (let* ((min_width
	      (if (bool (pylist-ref spec "zeropad"))
		  (- (pylist-ref spec "minimumwidth") (len fracpart) (len sign))
		  0))
	     (intpart (_insert_thousands_sep intpart spec min_width)))


	(_format_align sign  (+ intpart fracpart) spec)))))


;;##### Useful Constants (internal use only) ################################

;; Reusable defaults
(pk 13 1)
(define _Infinity         (Decimal "Inf"))
(pk 13 2)
(define _NegativeInfinity (Decimal "-Inf"))
(pk 13 3)
(define _NaN              (Decimal "NaN"))
(pk 13 4)
(define _Zero             (Decimal 0))
(pk 13 5)
(define _One              (Decimal 1))
(pk 13 6)
(define _NegativeOne      (Decimal -1))
(pk 13 7)
;; _SignedInfinity[sign] is infinity w/ that sign
(define _SignedInfinity  (list _Infinity _NegativeInfinity))
(pk 13 8)
;; Constants related to the hash implementation;  hash(x) is based
;; on the reduction of x modulo _PyHASH_MODULUS
(define _PyHASH_MODULUS (ref hash_info 'modulus))
(pk 13 9)
;; hash values to use for positive and negative infinities, and nans
(define _PyHASH_INF (ref hash_info 'inf))
(pk 13 10)
(define _PyHASH_NAN (ref hash_info 'nan))
(pk 13 11)
;; _PyHASH_10INV is the inverse of 10 modulo the prime _PyHASH_MODULUS
(define _PyHASH_10INV (pow 10 (- _PyHASH_MODULUS 2) _PyHASH_MODULUS))
(pk 13 12)