author Stefan Israelsson Tampe Tue, 17 Apr 2018 14:02:56 +0000 (16:02 +0200) committer Stefan Israelsson Tampe Tue, 17 Apr 2018 14:02:56 +0000 (16:02 +0200)

index ca5e32e5ea2eb6cfcb3346d469f205700f0b9e8f..4fab02d11495c48fa3bc6ce22c565be43893b856 100644 (file)
@@ -11,7 +11,7 @@
(define-inlinable (xy v seed)
(modulo
(logxor seed
-           (+ (py-hash v)
+           (+ v
#x9e3779b9
(ash seed 6)
(ash seed -2)))
@@ -42,7 +42,7 @@
s))))

(define-method (py-hash (x <p>))
-  (aif it (ref x '__hash__)
-       (it)
-       (hash x complexity)))
+  (aif it (pk 'it (ref x '__hash__))
+       (pk 'hash (it))
+       (next-method)))

index 2cf1ce574698e37ba1c48372d3b4889666759030..a4c31c745bba2f9d358cb0af0241361b6311da14 100644 (file)
(begin
(set! xc (/ xc n))
(set! xe (+ xe 1))
-                    (lp)))))
+                    (lp))))))

(let* ((x  (_WorkRep self))
(xc (ref x 'int))

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

-        # now y is positive; find m and n such that y = m/n
-        if ye >= 0:
-            m, n = yc*10**ye, 1
-        else:
-            if xe != 0 and len(str(abs(yc*xe))) <= -ye:
-                return None
-            xc_bits = _nbits(xc)
-            if xc != 1 and len(str(abs(yc)*xc_bits)) <= -ye:
-                return None
-            m, n = yc, 10**(-ye)
-            while m % 2 == n % 2 == 0:
-                m //= 2
-                n //= 2
-            while m % 5 == n % 5 == 0:
-                m //= 5
-                n //= 5
-
-        # compute nth root of xc*10**xe
-        if n > 1:
-            # if 1 < xc < 2**n then xc isn't an nth power
-            if xc != 1 and xc_bits <= n:
-                return None
-
-            xe, rem = divmod(xe, n)
-            if rem != 0:
-                return None
-
-            # compute nth root of xc using Newton's method
-            a = 1 << -(-_nbits(xc)//n) # initial estimate
-            while True:
-                q, r = divmod(xc, a**(n-1))
-                if a <= q:
-                    break
-                else:
-                    a = (a*(n-1) + q)//n
-            if not (a == q and r == 0):
-                return None
-            xc = a
-
-        # now xc*10**xe is the nth root of the original xc*10**xe
-        # compute mth power of xc*10**xe
-
-        # if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
-        # 10**p and the result is not representable.
-        if xc > 1 and m > p*100//_log10_lb(xc):
-            return None
-        xc = xc**m
-        xe *= m
-        if xc > 10**p:
-            return None
-
-        # by this point the result *is* exactly representable
-        # adjust the exponent to get as close as possible to the ideal
-        # exponent, if necessary
-        str_xc = str(xc)
-        if other._isinteger() and other._sign == 0:
-            ideal_exponent = self._exp*int(other)
-            zeros = min(xe-ideal_exponent, p-len(str_xc))
-        else:
-            zeros = 0
-        return _dec_from_triple(0, str_xc+'0'*zeros, xe-zeros)
+        ;; 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) it
+            (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

-    def __pow__(self, other, modulo=None, context=None):
-        """Return self ** other [ % modulo].
+        ;; 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.

result that would be obtained by computing (self**other) %
modulo with unbounded precision, but is computed more
efficiently.  It is always exact.
-        """
-
-        if modulo is not None:
-            return self._power_modulo(other, modulo, context)
-
-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
-
-        if context is None:
-            context = getcontext()
-
-        # either argument is a NaN => result is NaN
-        ans = self._check_nans(other, context)
-        if ans:
-            return ans
-
-        # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
-        if not other:
-            if not self:
-                return context._raise_error(InvalidOperation, '0 ** 0')
-            else:
-                return _One
-
-        # result has sign 1 iff self._sign is 1 and other is an odd integer
-        result_sign = 0
-        if self._sign == 1:
-            if other._isinteger():
-                if not other._iseven():
-                    result_sign = 1
-            else:
-                # -ve**noninteger = NaN
-                # (-0)**noninteger = 0**noninteger
-                if self:
-                    return context._raise_error(InvalidOperation,
-                        'x ** y with x negative and y not an integer')
-            # negate self, without doing any unwanted rounding
-            self = self.copy_negate()
-
-        # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
-        if not self:
-            if other._sign == 0:
-                return _dec_from_triple(result_sign, '0', 0)
-            else:
-                return _SignedInfinity[result_sign]
-
-        # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
-        if self._isinfinity():
-            if other._sign == 0:
-                return _SignedInfinity[result_sign]
-            else:
-                return _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
-        if self == _One:
-            if other._isinteger():
-                # exp = max(self._exp*max(int(other), 0),
-                # 1-context.prec) but evaluating int(other) directly
-                # is dangerous until we know other is small (other
-                # could be 1e999999999)
-                if other._sign == 1:
-                    multiplier = 0
-                elif other > context.prec:
-                    multiplier = context.prec
-                else:
-                    multiplier = int(other)
-
-                exp = self._exp * multiplier
-                if exp < 1-context.prec:
-                    exp = 1-context.prec
-                    context._raise_error(Rounded)
-            else:
-                context._raise_error(Inexact)
-                context._raise_error(Rounded)
-                exp = 1-context.prec
-
-            return _dec_from_triple(result_sign, '1'+'0'*-exp, exp)
-
-        # compute adjusted exponent of self
+        "

-        # self ** infinity is infinity if self > 1, 0 if self < 1
-        # self ** -infinity is infinity if self < 1, 0 if self > 1
-        if other._isinfinity():
-            if (other._sign == 0) == (self_adj < 0):
-                return _dec_from_triple(result_sign, '0', 0)
-            else:
-                return _SignedInfinity[result_sign]
-
-        # from here on, the result always goes through the call
-        # to _fix at the end of this function.
-        ans = None
-        exact = False
-
-        # crude test to catch cases of extreme overflow/underflow.  If
-        # log10(self)*other >= 10**bound and bound >= len(str(Emax))
-        # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
-        # self**other >= 10**(Emax+1), so overflow occurs.  The test
-        # for underflow is similar.
-        bound = self._log10_exp_bound() + other.adjusted()
-        if (self_adj >= 0) == (other._sign == 0):
-            # self > 1 and other +ve, or self < 1 and other -ve
-            # possibility of overflow
-            if bound >= len(str(context.Emax)):
-                ans = _dec_from_triple(result_sign, '1', context.Emax+1)
-        else:
-            # self > 1 and other -ve, or self < 1 and other +ve
-            # possibility of underflow to 0
-            Etiny = context.Etiny()
-            if bound >= len(str(-Etiny)):
-                ans = _dec_from_triple(result_sign, '1', Etiny-1)
-
-        # try for an exact result with precision +1
-        if ans is None:
-            ans = self._power_exact(other, context.prec + 1)
-            if ans is not None:
-                if result_sign == 1:
-                    ans = _dec_from_triple(1, ans._int, ans._exp)
-                exact = True
-
-        # usual case: inexact result, x**y computed directly as exp(y*log(x))
-        if ans is None:
-            p = context.prec
-            x = _WorkRep(self)
-            xc, xe = x.int, x.exp
-            y = _WorkRep(other)
-            yc, ye = y.int, y.exp
-            if y.sign == 1:
-                yc = -yc
-
-            # compute correctly rounded result:  start with precision +3,
-            # then increase precision until result is unambiguously roundable
-            extra = 3
-            while True:
-                coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
-                if coeff % (5*10**(len(str(coeff))-p-1)):
-                    break
-                extra += 3
+       (twix
+        ((not (eq= modulo None)) it
+         ((ref self '_power_modulo) other modulo context))

-            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 exact and not other._isinteger():
-            # pad with zeros up to length context.prec+1 if necessary; this
-            # ensures that the Rounded signal will be raised.
-            if len(ans._int) <= context.prec:
-                expdiff = context.prec + 1 - len(ans._int)
-                ans = _dec_from_triple(ans._sign, ans._int+'0'*expdiff,
-                                       ans._exp-expdiff)
-
-            # create a copy of the current context, with cleared flags/traps
-            newcontext = context.copy()
-            newcontext.clear_flags()
-            for exception in _signals:
-                newcontext.traps[exception] = 0
-
-            # round in the new context
-            ans = ans._fix(newcontext)
-
-            # raise Inexact, and if necessary, Underflow
-            newcontext._raise_error(Inexact)
-            if newcontext.flags[Subnormal]:
-                newcontext._raise_error(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 newcontext.flags[Overflow]:
-                context._raise_error(Overflow, 'above Emax', ans._sign)
-            for exception in Underflow, Subnormal, Inexact, Rounded, Clamped:
-                if newcontext.flags[exception]:
-                    context._raise_error(exception)
+        ((norm-op other ) it it)
+        (let (get-context context))

-        else:
-          ans = ans._fix(context)
+        ;; either argument is a NaN => result is NaN
+        (bin-special self other context)

-        return ans
+        ;; 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
+        ((not (bool other))
+         (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)

-    def __rpow__(self, other, context=None):
-        """Swaps self/other and returns __pow__."""
-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
-        return other.__pow__(self, context=context)
+        ;; 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
+
+        ;; 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))
+        (if (eq? (= (ref other '_sign) 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 False)
+             (bound (+ ((ref self '_log10_exp_bound))
+
+         ;; 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 (strcoeff) 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)))

-    def normalize(self, context=None):
-        """Normalize- strip trailing 0s, change anything equal to 0 to 0e0"""
+        ans))))

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

-        if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-
-        dup = self._fix(context)
-        if dup._isinfinity():
-            return dup
-
-        if not dup:
-            return _dec_from_triple(dup._sign, '0', 0)
-        exp_max = [context.Emax, context.Etop()][context.clamp]
-        end = len(dup._int)
-        exp = dup._exp
-        while dup._int[end-1] == '0' and exp < exp_max:
-            exp += 1
-            end -= 1
-        return _dec_from_triple(dup._sign, dup._int[:end], exp)
-
-    def quantize(self, exp, rounding=None, context=None):
-        """Quantize self so its exponent is the same as that of exp.
+    (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)
+        (let ((dup ((ref self _fix) context))))
+        (((dup '_isinfinity)) it dup)
+        ((not (bool dup))
+         (_dec_from_triple (reg 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.
-        """
-        exp = _convert_other(exp, raiseit=True)
-
-        if context is None:
-            context = getcontext()
-        if rounding is None:
-            rounding = context.rounding
-
-        if self._is_special or exp._is_special:
-            ans = self._check_nans(exp, context)
-            if ans:
-                return ans
-
-            if exp._isinfinity() or self._isinfinity():
-                if exp._isinfinity() and self._isinfinity():
-                    return Decimal(self)  # if both are inf, it is OK
-                return context._raise_error(InvalidOperation,
-                                        'quantize with one INF')
-
-        # exp._exp should be between Etiny and Emax
-        if not (context.Etiny() <= exp._exp <= context.Emax):
-            return context._raise_error(InvalidOperation,
-                   'target exponent out of bounds in quantize')
-
-        if not self:
-            ans = _dec_from_triple(self._sign, '0', exp._exp)
-            return ans._fix(context)
-
-            return context._raise_error(InvalidOperation,
-                                        'exponent of quantize result too large for current context')
-        if self_adjusted - exp._exp + 1 > context.prec:
-            return context._raise_error(InvalidOperation,
-                                        'quantize result has too many digits for current context')
+        "
+       (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"))

-        ans = self._rescale(exp._exp, rounding)
-            return context._raise_error(InvalidOperation,
-                                        'exponent of quantize result too large for current context')
-        if len(ans._int) > context.prec:
-            return context._raise_error(InvalidOperation,
-                                        'quantize result has too many digits for current context')
+        ((not (bool self)) it
+        (let ((ans (_dec_from_triple (ref self '_sign) "0" _eexp)))
+          ((ref ans '_fix) context)))

-        # raise appropriate flags
-        if ans and ans.adjusted() < context.Emin:
-            context._raise_error(Subnormal)
-        if ans._exp > self._exp:
-            if ans != self:
-                context._raise_error(Inexact)
-            context._raise_error(Rounded)
+             (prec          (cx-prec context))))

-        # call to fix takes care of any necessary folddown, and
-        # signals Clamped if necessary
-        ans = ans._fix(context)
-        return ans
+       ((> 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"))

-    def same_quantum(self, other, context=None):
-        """Return True if self and other have the same exponent; otherwise
+       (let ((ans ((ref self '_rescale) _eexp rounding))))
+
+       (if (> ((ref ans 'adjusted)) Emax) it
+           ((cx-error context) InvalidOperation,
+            "exponent of quantize result too large for current context"))
+
+        ((> (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 (> (reg 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.
-        """
-        other = _convert_other(other, raiseit=True)
-        if self._is_special or other._is_special:
-            return (self.is_nan() and other.is_nan() or
-                    self.is_infinite() and other.is_infinite())
-        return self._exp == other._exp
+        "
+        (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))))

-    def _rescale(self, exp, rounding):
-        """Rescale self so that the exponent is exp, either by padding with zeros
+    (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

exp = exp to scale to (an integer)
rounding = rounding mode
-        """
-        if self._is_special:
-            return Decimal(self)
-        if not self:
-            return _dec_from_triple(self._sign, '0', exp)
-
-        if self._exp >= exp:
-            return _dec_from_triple(self._sign,
-                                        self._int + '0'*(self._exp - exp), exp)
-
-        # too many digits; round and lose data.  If self.adjusted() <
-        # exp-1, replace self by 10**(exp-1) before rounding
-        digits = len(self._int) + self._exp - exp
-        if digits < 0:
-            self = _dec_from_triple(self._sign, '1', exp-1)
-            digits = 0
-        this_function = self._pick_rounding_function[rounding]
-        changed = this_function(self, digits)
-        coeff = self._int[:digits] or '0'
-        if changed == 1:
-            coeff = str(int(coeff)+1)
-        return _dec_from_triple(self._sign, coeff, exp)
-
-    def _round(self, places, rounding):
-        """Round a nonzero, nonspecial Decimal to a fixed number of
+        "
+
+       (cond
+        ((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)
+            (_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.

-        """
-        if places <= 0:
-            raise ValueError("argument should be at least 1 in _round")
-        if self._is_special or not self:
-            return Decimal(self)
-        # 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.
-        return ans
-
-    def to_integral_exact(self, rounding=None, context=None):
-        """Rounds to a nearby integer.
+        "
+        (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.
+               (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

See also: to_integral_value, which does exactly the same as
this method except that it doesn't raise Inexact or Rounded.
-        """
-        if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-            return Decimal(self)
-        if self._exp >= 0:
-            return Decimal(self)
-        if not self:
-            return _dec_from_triple(self._sign, '0', 0)
-        if context is None:
-            context = getcontext()
-        if rounding is None:
-            rounding = context.rounding
-        ans = self._rescale(0, rounding)
-        if ans != self:
-            context._raise_error(Inexact)
-        context._raise_error(Rounded)
-        return ans
+        "
+       (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 (boool 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)

-    def to_integral_value(self, rounding=None, context=None):
-        """Rounds to the nearest integer, without raising inexact, rounded."""
-        if context is None:
-            context = getcontext()
-        if rounding is None:
-            rounding = context.rounding
-        if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
-            return Decimal(self)
-        if self._exp >= 0:
-            return Decimal(self)
-        else:
-            return self._rescale(0, rounding)
+          ans)))))

-    # the method name changed, but we provide also the old one, for compatibility
-    to_integral = to_integral_value
+    (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))))))

-    def sqrt(self, context=None):
-        """Return the square root of self."""
-        if context is None:
-            context = getcontext()
+    ;; the method name changed, but we provide also the old one,
+    ;; for compatibility
+    (define to_integral to_integral_value)

-        if self._is_special:
-            ans = self._check_nans(context=context)
-            if ans:
-                return ans
+    (define sqrt
+      (lam (self (= context None))
+       "Return the square root of self."
+       (twix
+        (let (get-context context))

-            if self._isinfinity() and self._sign == 0:
-                return Decimal(self)
+        ((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)))

-        if not self:
-            # exponent = self._exp // 2.  sqrt(-0) = -0
-            ans = _dec_from_triple(self._sign, '0', self._exp // 2)
-            return ans._fix(context)
+             #f) it it)

-        if self._sign == 1:
-            return context._raise_error(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
-        prec = context.prec+1
-
-        # write argument in the form c*100**e where e = self._exp//2
-        # is the 'ideal' exponent, to be used if the square root is
-        # exactly representable.  l is the number of 'digits' of c in
-        # base 100, so that 100**(l-1) <= c < 100**l.
-        op = _WorkRep(self)
-        e = op.exp >> 1
-        if op.exp & 1:
-            c = op.int * 10
-            l = (len(self._int) >> 1) + 1
-        else:
-            c = op.int
-            l = len(self._int)+1 >> 1
-
-        # rescale so that c has exactly prec base 100 'digits'
-        shift = prec-l
-        if shift >= 0:
-            c *= 100**shift
-            exact = True
-        else:
-            c, remainder = divmod(c, 100**-shift)
-            exact = not remainder
-        e -= shift
+        ((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)
+         ((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))))

-        # find n = floor(sqrt(c)) using Newton's method
-        n = 10**prec
-        while True:
-            q = c//n
-            if n <= q:
-                break
-            else:
-                n = n + q >> 1
-        exact = exact and n*n == c
-
-        if exact:
-            # result is exact; rescale to use ideal exponent e
-            if shift >= 0:
-                # assert n % 10**shift == 0
-                n //= 10**shift
-            else:
-                n *= 10**-shift
-            e += shift
-        else:
-            # result is not exact; fix last digit as described above
-            if n % 5 == 0:
-                n += 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))))

ans = _dec_from_triple(0, str(n), e)

return ans._fix(context)

-    def min(self, other, context=None):
+   def min(self, other, context=None):
"""Returns the smaller value.

Like min(self, other) except if one is not a number, returns
@@ -5731,3 +5849,746 @@ def _normalize(op1, op2, prec = 0):
return op1, op2

+##### Integer arithmetic functions used by ln, log10, exp and __pow__ #####
+
+_nbits = int.bit_length
+
+def _decimal_lshift_exact(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
+
+    """
+    if n == 0:
+        return 0
+    elif e >= 0:
+        return n * 10**e
+    else:
+        # val_n = largest power of 10 dividing n.
+        str_n = str(abs(n))
+        val_n = len(str_n) - len(str_n.rstrip('0'))
+        return None if val_n < -e else n // 10**-e
+
+def _sqrt_nearest(n, a):
+    """Closest integer to the square root of the positive integer n.  a is
+    an initial approximation to the square root.  Any positive integer
+    will do for a, but the closer a is to the square root of n the
+    faster convergence will be.
+
+    """
+    if n <= 0 or a <= 0:
+        raise ValueError("Both arguments to _sqrt_nearest should be positive.")
+
+    b=0
+    while a != b:
+        b, a = a, a--n//a>>1
+    return a
+
+def _rshift_nearest(x, shift):
+    """Given an integer x and a nonnegative integer shift, return closest
+    integer to x / 2**shift; use round-to-even in case of a tie.
+
+    """
+    b, q = 1 << shift, x >> shift
+    return q + (2*(x & (b-1)) + (q&1) > b)
+
+def _div_nearest(a, b):
+    """Closest integer to a/b, a and b positive integers; rounds to even
+    in the case of a tie.
+
+    """
+    q, r = divmod(a, b)
+    return q + (2*r + (q&1) > b)
+
+def _ilog(x, M, L = 8):
+    """Integer approximation to M*log(x/M), with absolute error boundable
+    in terms only of x/M.
+
+    Given positive integers x and M, return an integer approximation to
+    M * log(x/M).  For L = 8 and 0.1 <= x/M <= 10 the difference
+    between the approximation and the exact result is at most 22.  For
+    L = 8 and 1.0 <= x/M <= 10.0 the difference is at most 15.  In
+    both cases these are upper bounds on the error; it will usually be
+    much smaller."""
+
+    # The basic algorithm is the following: let log1p be the function
+    # log1p(x) = log(1+x).  Then log(x/M) = log1p((x-M)/M).  We use
+    # the reduction
+    #
+    #    log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
+    #
+    # repeatedly until the argument to log1p is small (< 2**-L in
+    # absolute value).  For small y we can use the Taylor series
+    # expansion
+    #
+    #    log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
+    #
+    # truncating at T such that y**T is small enough.  The whole
+    # computation is carried out in a form of fixed-point arithmetic,
+    # with a real number z being represented by an integer
+    # approximation to z*M.  To avoid loss of precision, the y below
+    # is actually an integer approximation to 2**R*y*M, where R is the
+    # number of reductions performed so far.
+
+    y = x-M
+    # argument reduction; R = number of reductions performed
+    R = 0
+    while (R <= L and abs(y) << L-R >= M or
+           R > L and abs(y) >> R-L >= M):
+        y = _div_nearest((M*y) << 1,
+                         M + _sqrt_nearest(M*(M+_rshift_nearest(y, R)), M))
+        R += 1
+
+    # Taylor series with T terms
+    T = -int(-10*len(str(M))//(3*L))
+    yshift = _rshift_nearest(y, R)
+    w = _div_nearest(M, T)
+    for k in range(T-1, 0, -1):
+        w = _div_nearest(M, k) - _div_nearest(yshift*w, M)
+
+    return _div_nearest(w*y, M)
+
+def _dlog10(c, e, p):
+    """Given integers c, e and p with c > 0, p >= 0, compute an integer
+    approximation to 10**p * log10(c*10**e), with an absolute error of
+    at most 1.  Assumes that c*10**e is not exactly 1."""
+
+    # increase precision by 2; compensate for this by dividing
+    # final result by 100
+    p += 2
+
+    # write c*10**e as d*10**f with either:
+    #   f >= 0 and 1 <= d <= 10, or
+    #   f <= 0 and 0.1 <= d <= 1.
+    # Thus for c*10**e close to 1, f = 0
+    l = len(str(c))
+    f = e+l - (e+l >= 1)
+
+    if p > 0:
+        M = 10**p
+        k = e+p-f
+        if k >= 0:
+            c *= 10**k
+        else:
+            c = _div_nearest(c, 10**-k)
+
+        log_d = _ilog(c, M) # error < 5 + 22 = 27
+        log_10 = _log10_digits(p) # error < 1
+        log_d = _div_nearest(log_d*M, log_10)
+        log_tenpower = f*M # exact
+    else:
+        log_d = 0  # error < 2.31
+        log_tenpower = _div_nearest(f, 10**-p) # error < 0.5
+
+    return _div_nearest(log_tenpower+log_d, 100)
+
+def _dlog(c, e, p):
+    """Given integers c, e and p with c > 0, compute an integer
+    approximation to 10**p * log(c*10**e), with an absolute error of
+    at most 1.  Assumes that c*10**e is not exactly 1."""
+
+    # Increase precision by 2. The precision increase is compensated
+    # for at the end with a division by 100.
+    p += 2
+
+    # rewrite c*10**e as d*10**f with either f >= 0 and 1 <= d <= 10,
+    # or f <= 0 and 0.1 <= d <= 1.  Then we can compute 10**p * log(c*10**e)
+    # as 10**p * log(d) + 10**p*f * log(10).
+    l = len(str(c))
+    f = e+l - (e+l >= 1)
+
+    # compute approximation to 10**p*log(d), with error < 27
+    if p > 0:
+        k = e+p-f
+        if k >= 0:
+            c *= 10**k
+        else:
+            c = _div_nearest(c, 10**-k)  # error of <= 0.5 in c
+
+        # _ilog magnifies existing error in c by a factor of at most 10
+        log_d = _ilog(c, 10**p) # error < 5 + 22 = 27
+    else:
+        # p <= 0: just approximate the whole thing by 0; error < 2.31
+        log_d = 0
+
+    # compute approximation to f*10**p*log(10), with error < 11.
+    if f:
+        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
+            f_log_ten = _div_nearest(f*_log10_digits(p+extra), 10**extra)
+        else:
+            f_log_ten = 0
+    else:
+        f_log_ten = 0
+
+    # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1
+    return _div_nearest(f_log_ten + log_d, 100)
+
+class _Log10Memoize(object):
+    """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__."""
+    def __init__(self):
+        self.digits = "23025850929940456840179914546843642076011014886"
+
+    def getdigits(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(self.digits):
+            # compute p+3, p+6, p+9, ... digits; continue until at
+            # least one of the extra digits is nonzero
+            extra = 3
+            while True:
+                # compute p+extra digits, correct to within 1ulp
+                M = 10**(p+extra+2)
+                digits = str(_div_nearest(_ilog(10*M, M), 100))
+                if digits[-extra:] != '0'*extra:
+                    break
+                extra += 3
+            # keep all reliable digits so far; remove trailing zeros
+            # and next nonzero digit
+            self.digits = digits.rstrip('0')[:-1]
+        return int(self.digits[:p+1])
+
+_log10_digits = _Log10Memoize().getdigits
+
+def _iexp(x, M, L=8):
+    """Given integers x and M, M > 0, such that x/M is small in absolute
+    value, compute an integer approximation to M*exp(x/M).  For 0 <=
+    x/M <= 2.4, the absolute error in the result is bounded by 60 (and
+    is usually much smaller)."""
+
+    # Algorithm: to compute exp(z) for a real number z, first divide z
+    # by a suitable power R of 2 so that |z/2**R| < 2**-L.  Then
+    # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
+    # series
+    #
+    #     expm1(x) = x + x**2/2! + x**3/3! + ...
+    #
+    # Now use the identity
+    #
+    #     expm1(2x) = expm1(x)*(expm1(x)+2)
+    #
+    # R times to compute the sequence expm1(z/2**R),
+    # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).
+
+    # Find R such that x/2**R/M <= 2**-L
+    R = _nbits((x<<L)//M)
+
+    # Taylor series.  (2**L)**T > M
+    T = -int(-10*len(str(M))//(3*L))
+    y = _div_nearest(x, T)
+    Mshift = M<<R
+    for i in range(T-1, 0, -1):
+        y = _div_nearest(x*(Mshift + y), Mshift * i)
+
+    # Expansion
+    for k in range(R-1, -1, -1):
+        Mshift = M<<(k+2)
+        y = _div_nearest(y*(y+Mshift), Mshift)
+
+    return M+y
+
+def _dexp(c, e, p):
+    """Compute an approximation to exp(c*10**e), with p decimal places of
+    precision.
+
+    Returns 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
+    p += 2
+
+    # compute log(10) with extra precision = adjusted exponent of c*10**e
+    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
+    shift = e+q
+    if shift >= 0:
+        cshift = c*10**shift
+    else:
+        cshift = c//10**-shift
+    quot, rem = divmod(cshift, _log10_digits(q))
+
+    # reduce remainder back to original precision
+    rem = _div_nearest(rem, 10**extra)
+
+    # error in result of _iexp < 120;  error after division < 0.62
+    return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
+
+def _dpower(xc, xe, yc, ye, p):
+    """Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and
+    y = yc*10**ye, compute x**y.  Returns a pair of integers (c, e) such that:
+
+      10**(p-1) <= c <= 10**p, and
+      (c-1)*10**e < x**y < (c+1)*10**e
+
+    in other words, c*10**e is an approximation to x**y with p digits
+    of precision, and with an error in c of at most 1.  (This is
+    almost, but not quite, the same as the error being < 1ulp: when c
+    == 10**(p-1) we can only guarantee error < 10ulp.)
+
+    We assume that: x is positive and not equal to 1, and y is nonzero.
+    """
+
+    # Find b such that 10**(b-1) <= |y| <= 10**b
+    b = len(str(abs(yc))) + ye
+
+    # log(x) = lxc*10**(-p-b-1), to p+b+1 places after the decimal point
+    lxc = _dlog(xc, xe, p+b+1)
+
+    # compute product y*log(x) = yc*lxc*10**(-p-b-1+ye) = pc*10**(-p-1)
+    shift = ye-b
+    if shift >= 0:
+        pc = lxc*yc*10**shift
+    else:
+        pc = _div_nearest(lxc*yc, 10**-shift)
+
+    if pc == 0:
+        # we prefer a result that isn't exactly 1; this makes it
+        # easier to compute a correctly rounded result in __pow__
+        if ((len(str(xc)) + xe >= 1) == (yc > 0)): # if x**y > 1:
+            coeff, exp = 10**(p-1)+1, 1-p
+        else:
+            coeff, exp = 10**p-1, -p
+    else:
+        coeff, exp = _dexp(pc, -(p+1), p+1)
+        coeff = _div_nearest(coeff, 10)
+        exp += 1
+
+    return coeff, exp
+
+def _log10_lb(c, correction = {
+        '1': 100, '2': 70, '3': 53, '4': 40, '5': 31,
+        '6': 23, '7': 16, '8': 10, '9': 5}):
+    """Compute a lower bound for 100*log10(c) for a positive integer c."""
+    if c <= 0:
+        raise ValueError("The argument to _log10_lb should be nonnegative.")
+    str_c = str(c)
+    return 100*len(str_c) - correction[str_c]
+
+##### Helper Functions ####################################################
+
+def _convert_other(other, raiseit=False, allow_float=False):
+    """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).
+
+    """
+    if isinstance(other, Decimal):
+        return other
+    if isinstance(other, int):
+        return Decimal(other)
+    if allow_float and isinstance(other, float):
+        return Decimal.from_float(other)
+
+    if raiseit:
+        raise TypeError("Unable to convert %s to Decimal" % other)
+    return NotImplemented
+
+def _convert_for_comparison(self, other, equality_op=False):
+    """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".
+
+    """
+    if isinstance(other, Decimal):
+        return 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.
+    if isinstance(other, _numbers.Rational):
+        if not self._is_special:
+            self = _dec_from_triple(self._sign,
+                                    str(int(self._int) * other.denominator),
+                                    self._exp)
+        return self, Decimal(other.numerator)
+
+    # Comparisons with float and complex types.  == and != comparisons
+    # with complex numbers should succeed, returning either True or False
+    # as appropriate.  Other comparisons return NotImplemented.
+    if equality_op and isinstance(other, _numbers.Complex) and other.imag == 0:
+        other = other.real
+    if isinstance(other, float):
+        context = getcontext()
+        if equality_op:
+            context.flags[FloatOperation] = 1
+        else:
+            context._raise_error(FloatOperation,
+                "strict semantics for mixing floats and Decimals are enabled")
+        return self, Decimal.from_float(other)
+    return NotImplemented, NotImplemented
+
+
+##### Setup Specific Contexts ############################################
+
+# The default context prototype used by Context()
+# Is mutable, so that new contexts can have different default values
+
+DefaultContext = Context(
+        prec=28, rounding=ROUND_HALF_EVEN,
+        traps=[DivisionByZero, Overflow, InvalidOperation],
+        flags=[],
+        Emax=999999,
+        Emin=-999999,
+        capitals=1,
+        clamp=0
+)
+
+# 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.
+
+BasicContext = Context(
+        prec=9, rounding=ROUND_HALF_UP,
+        traps=[DivisionByZero, Overflow, InvalidOperation, Clamped, Underflow],
+        flags=[],
+)
+
+ExtendedContext = Context(
+        prec=9, rounding=ROUND_HALF_EVEN,
+        traps=[],
+        flags=[],
+)
+
+
+##### crud for parsing strings #############################################
+#
+# Regular expression used for parsing numeric strings.  Additional
+#
+# 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.
+
+import re
+_parser = re.compile(r"""        # 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
+""", re.VERBOSE | re.IGNORECASE).match
+
+_all_zeros = re.compile('0*\$').match
+_exact_half = re.compile('50*\$').match
+
+##### 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][#][minimumwidth][,][.precision][type]
+
+_parse_format_specifier_regex = re.compile(r"""\A
+(?:
+   (?P<fill>.)?
+   (?P<align>[<>=^])
+)?
+(?P<sign>[-+ ])?
+(?P<alt>\#)?
+(?P<minimumwidth>(?!0)\d+)?
+(?P<thousands_sep>,)?
+(?:\.(?P<precision>0|(?!0)\d+))?
+(?P<type>[eEfFgGn%])?
+\Z
+""", re.VERBOSE|re.DOTALL)
+
+del re
+
+# 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.
+try:
+    import locale as _locale
+except ImportError:
+    pass
+
+def _parse_format_specifier(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
+      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
+
+    """
+    m = _parse_format_specifier_regex.match(format_spec)
+    if m is None:
+        raise ValueError("Invalid format specifier: " + format_spec)
+
+    # get the dictionary
+    format_dict = m.groupdict()
+
+    # is requested, the fill and align fields should be absent.
+    fill = format_dict['fill']
+    align = format_dict['align']
+        if fill is not None:
+            raise ValueError("Fill character conflicts with '0'"
+                             " in format specifier: " + format_spec)
+        if align is not None:
+            raise ValueError("Alignment conflicts with '0' in "
+                             "format specifier: " + format_spec)
+    format_dict['fill'] = fill or ' '
+    # 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.
+    format_dict['align'] = align or '>'
+
+    # default sign handling: '-' for negative, '' for positive
+    if format_dict['sign'] is None:
+        format_dict['sign'] = '-'
+
+    # minimumwidth defaults to 0; precision remains None if not given
+    format_dict['minimumwidth'] = int(format_dict['minimumwidth'] or '0')
+    if format_dict['precision'] is not None:
+        format_dict['precision'] = int(format_dict['precision'])
+
+    # 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 format_dict['precision'] == 0:
+        if format_dict['type'] is None or format_dict['type'] in 'gGn':
+            format_dict['precision'] = 1
+
+    # determine thousands separator, grouping, and decimal separator, and
+    # add appropriate entries to format_dict
+    if format_dict['type'] == 'n':
+        # apart from separators, 'n' behaves just like 'g'
+        format_dict['type'] = 'g'
+        if _localeconv is None:
+            _localeconv = _locale.localeconv()
+        if format_dict['thousands_sep'] is not None:
+            raise ValueError("Explicit thousands separator conflicts with "
+                             "'n' type in format specifier: " + format_spec)
+        format_dict['thousands_sep'] = _localeconv['thousands_sep']
+        format_dict['grouping'] = _localeconv['grouping']
+        format_dict['decimal_point'] = _localeconv['decimal_point']
+    else:
+        if format_dict['thousands_sep'] is None:
+            format_dict['thousands_sep'] = ''
+        format_dict['grouping'] = [3, 0]
+        format_dict['decimal_point'] = '.'
+
+    return format_dict
+
+def _format_align(sign, body, spec):
+    """Given an unpadded, non-aligned numeric string 'body' and sign
+    format specifier dictionary 'spec' (as produced by
+    parse_format_specifier).
+
+    """
+    # how much extra space do we have to play with?
+    minimumwidth = spec['minimumwidth']
+    fill = spec['fill']
+    padding = fill*(minimumwidth - len(sign) - len(body))
+
+    align = spec['align']
+    if align == '<':
+        result = sign + body + padding
+    elif align == '>':
+        result = padding + sign + body
+    elif align == '=':
+        result = sign + padding + body
+    elif align == '^':
+    else:
+        raise ValueError('Unrecognised alignment field')
+
+    return result
+
+def _group_lengths(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 + 
+    #   (3) list of positive integers + [locale.CHAR_MAX], or
+
+    from itertools import chain, repeat
+    if not grouping:
+        return []
+    elif grouping[-1] == 0 and len(grouping) >= 2:
+        return chain(grouping[:-1], repeat(grouping[-2]))
+    elif grouping[-1] == _locale.CHAR_MAX:
+        return grouping[:-1]
+    else:
+        raise ValueError('unrecognised format for grouping')
+
+def _insert_thousands_sep(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.
+
+    """
+
+    sep = spec['thousands_sep']
+    grouping = spec['grouping']
+
+    groups = []
+    for l in _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
+        l = min(max(len(digits), min_width, 1), l)
+        groups.append('0'*(l - len(digits)) + digits[-l:])
+        digits = digits[:-l]
+        min_width -= l
+        if not digits and min_width <= 0:
+            break
+        min_width -= len(sep)
+    else:
+        l = max(len(digits), min_width, 1)
+        groups.append('0'*(l - len(digits)) + digits[-l:])
+    return sep.join(reversed(groups))
+
+def _format_sign(is_negative, spec):
+    """Determine sign character."""
+
+    if is_negative:
+        return '-'
+    elif spec['sign'] in ' +':
+        return spec['sign']
+    else:
+        return ''
+
+def _format_number(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
+      fill and align if necessary
+    """
+
+    sign = _format_sign(is_negative, spec)
+
+    if fracpart or spec['alt']:
+        fracpart = spec['decimal_point'] + fracpart
+
+    if exp != 0 or spec['type'] in 'eE':
+        echar = {'E': 'E', 'e': 'e', 'G': 'E', 'g': 'e'}[spec['type']]
+        fracpart += "{0}{1:+}".format(echar, exp)
+    if spec['type'] == '%':
+        fracpart += '%'
+
+        min_width = spec['minimumwidth'] - len(fracpart) - len(sign)
+    else:
+        min_width = 0
+    intpart = _insert_thousands_sep(intpart, spec, min_width)
+
+    return _format_align(sign, intpart+fracpart, spec)
+
+
+##### Useful Constants (internal use only) ################################
+
+# Reusable defaults
+_Infinity = Decimal('Inf')
+_NegativeInfinity = Decimal('-Inf')
+_NaN = Decimal('NaN')
+_Zero = Decimal(0)
+_One = Decimal(1)
+_NegativeOne = Decimal(-1)
+
+# _SignedInfinity[sign] is infinity w/ that sign
+_SignedInfinity = (_Infinity, _NegativeInfinity)
+
+# Constants related to the hash implementation;  hash(x) is based
+# on the reduction of x modulo _PyHASH_MODULUS
+_PyHASH_MODULUS = sys.hash_info.modulus
+# hash values to use for positive and negative infinities, and nans
+_PyHASH_INF = sys.hash_info.inf
+_PyHASH_NAN = sys.hash_info.nan
+
+# _PyHASH_10INV is the inverse of 10 modulo the prime _PyHASH_MODULUS
+_PyHASH_10INV = pow(10, _PyHASH_MODULUS - 2, _PyHASH_MODULUS)
+del sys
index 479c035336de2a291d8916914e8efef6b84e8772..37ed5f081a6811a023272967fda14490af889a86 100644 (file)
@@ -408,12 +408,9 @@ explicitly tell it to not update etc.
#f)))
(if (or (not f) (eq? f not-implemented))
(gox xx (mrefx xx key l))
-          (catch #t
-            (lambda ()
-              (f xx key))
-            (lambda x
-              (gox xx (mrefx xx key l))))))))
-
+         (catch #t
+                (lambda () (gox xx (f xx key)))
+                (lambda z  (if (pair? l) (car l) #f)))))))

(define-syntax-rule (mref x key l)
(let ((xx x))
@@ -921,18 +918,19 @@ explicitly tell it to not update etc.
#'(let ()
(define name
(letruc ((dname (make-up dval)) ...)
-                    (make-p-class 'name doc
-                                   parents
-                                   (lambda (dict)
-                                     (pylist-set! dict 'dname dname)
-                                     ...
-                                     (values)))))
-
-           (begin
-             (module-define! (current-module) 'ddname (ref name 'dname))
-             (name-object ddname))
-           ...
-
+                     (let ((ret
+                            (make-p-class 'name doc
+                                          parents
+                                          (lambda (dict)
+                                            (pylist-set! dict 'dname dname)
+                                            ...
+                                            (values)))))
+                       (begin
+                         (module-define! (current-module) 'ddname dname)
+                         (name-object ddname))
+                       ...
+                       ret)))
+
(module-define! (current-module) 'nname (ref name '__goops__))
(name-object nname)
(name-object name)
@@ -970,25 +968,21 @@ explicitly tell it to not update etc.
(define name
(letruc ((dname (make-up dval)) ...)
body
-                      (make-p-class 'name doc
-                                    parents
-                                    (lambda (dict)
-                                      (pylist-set! dict 'dname dname)
-                                      ...
-                                      (values)))))
-
-            (pk 1)
-           (begin
-              (pk 'ddname 'dname)
-             (module-define! (current-module) 'ddname (ref name 'dname))
-              (pk '*)
-             (name-object ddname))
-           ...
-            (pk 2)
+                      (let ((ret
+                            (make-p-class 'name doc
+                                          parents
+                                          (lambda (dict)
+                                            (pylist-set! dict 'dname dname)
+                                            ...
+                                            (values)))))
+                       (begin
+                         (module-define! (current-module) 'ddname dname)
+                         (name-object ddname))
+                       ...
+                       ret)))
(module-define! (current-module) 'nname (ref name '__goops__))
(name-object nname)
(name-object name)
-            (pk 3)
name))))))

(define-syntax mk-p-class-noname