author Stefan Israelsson Tampe Mon, 16 Apr 2018 13:55:30 +0000 (15:55 +0200) committer Stefan Israelsson Tampe Mon, 16 Apr 2018 13:55:30 +0000 (15:55 +0200)

@@ -2,6 +2,8 @@
#:use-module ((language python module collections) #:select (namedtuple))
#:export ())

+(define guile:modulo (@ (guile) moduolo))
+
(define __name__   "decimal")
(define __xname__  __name__)
(define __version__ "1.70")
(lambda (self context)
"Decapitate the payload of a NaN to fit the context"
-
+
;; maximum length of payload is precision if clamp=0,
;; precision-1 if clamp=1.
(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
-
-    def _round_down(self, prec):
-        """Also known as round-towards-0, truncate."""
-        if _all_zeros(self._int, prec):
-            return 0
-        else:
-            return -1
-
-    def _round_up(self, prec):
-        """Rounds away from 0."""
-        return -self._round_down(prec)
-
-    def _round_half_up(self, prec):
-        """Rounds 5 up (away from 0)"""
-        if self._int[prec] in '56789':
-            return 1
-        elif _all_zeros(self._int, prec):
-            return 0
-        else:
-            return -1
-
-    def _round_half_down(self, prec):
-        """Round 5 down"""
-        if _exact_half(self._int, prec):
-            return -1
-        else:
-            return self._round_half_up(prec)
-
-    def _round_half_even(self, prec):
-        """Round 5 to even, rest to nearest."""
-        if _exact_half(self._int, prec) and \
-                (prec == 0 or self._int[prec-1] in '02468'):
-            return -1
-        else:
-            return self._round_half_up(prec)
-
-    def _round_ceiling(self, prec):
-        """Rounds up (not away from 0 if negative.)"""
-        if self._sign:
-            return self._round_down(prec)
-        else:
-            return -self._round_down(prec)
-
-    def _round_floor(self, prec):
-        """Rounds down (not towards 0 if negative)"""
-        if not self._sign:
-            return self._round_down(prec)
-        else:
-            return -self._round_down(prec)
-
-    def _round_05up(self, prec):
-        """Round down unless digit prec-1 is 0 or 5."""
-        if prec and self._int[prec-1] not in '05':
-            return self._round_down(prec)
-        else:
-            return -self._round_down(prec)
-
-    _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,
-    )
-
-    def __round__(self, n=None):
-        """Round self to the nearest integer, or to a given precision.
+    ;; 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
>>> round(Decimal('sNaN123'), 0)
Decimal('NaN123')

-        """
-        if n is not None:
-            # two-argument form: use the equivalent quantize call
-            if not isinstance(n, int):
-                raise TypeError('Second argument to round should be integral')
-            exp = _dec_from_triple(0, '1', -n)
-            return self.quantize(exp)
-
-        # one-argument form
-        if self._is_special:
-            if self.is_nan():
-                raise ValueError("cannot round a NaN")
-            else:
-                raise OverflowError("cannot round an infinity")
-        return int(self._rescale(0, ROUND_HALF_EVEN))
-
-    def __floor__(self):
-        """Return the floor of self, as an integer.
+        "
+        (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
+           (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 self._is_special:
-            if self.is_nan():
-                raise ValueError("cannot round a NaN")
-            else:
-                raise OverflowError("cannot round an infinity")
-        return int(self._rescale(0, ROUND_FLOOR))
+        "
+       (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)))))

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

For a finite Decimal instance self, return the least integer n
Python exception is raised.

"""
-        if self._is_special:
-            if self.is_nan():
-                raise ValueError("cannot round a NaN")
-            else:
-                raise OverflowError("cannot round an infinity")
-        return int(self._rescale(0, ROUND_CEILING))
+       (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)))))

-    def fma(self, other, third, context=None):
+    (define fma
+      (lam (self other third (= context None))

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))
+             ((equals? (ref self '_exp) "N")  it
+              ((cx-error context) InvalidOperation "sNaN" self))
+             ((equals? (ref other '_exp) "N") it
+              ((cx-error context) InvalidOperation "sNaN" other))
+             ((equals? (ref self '_exp) "n")  it
+              (fin self))
+             ((equals? (ref other '_exp) "n") it
+              (fin other))
+             ((equals? (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)))))
+             ((equals? (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 other ) it it)
+        ((norm-op modulo) it it)
+        (let (get-context context))

-        other = _convert_other(other, raiseit=True)
-        third = _convert_other(third, raiseit=True)
+        ;; 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)))
+         ((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)))))

-        # compute product; raise InvalidOperation if either operand is
-        # a signaling NaN or if the product is zero times infinity.
-        if self._is_special or other._is_special:
-            if context is None:
-                context = getcontext()
-            if self._exp == 'N':
-                return context._raise_error(InvalidOperation, 'sNaN', self)
-            if other._exp == 'N':
-                return context._raise_error(InvalidOperation, 'sNaN', other)
-            if self._exp == 'n':
-                product = self
-            elif other._exp == 'n':
-                product = other
-            elif self._exp == 'F':
-                if not other:
-                    return context._raise_error(InvalidOperation,
-                                                'INF * 0 in fma')
-                product = _SignedInfinity[self._sign ^ other._sign]
-            elif other._exp == 'F':
-                if not self:
-                    return context._raise_error(InvalidOperation,
-                                                '0 * INF in fma')
-                product = _SignedInfinity[self._sign ^ other._sign]
-        else:
-            product = _dec_from_triple(self._sign ^ other._sign,
-                                       str(int(self._int) * int(other._int)),
-                                       self._exp + other._exp)
+    (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.

-    def _power_modulo(self, other, modulo, context=None):
-        """Three argument version of __pow__"""
+        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).

-        other = _convert_other(other)
-        if other is NotImplemented:
-            return other
-        modulo = _convert_other(modulo)
-        if modulo is NotImplemented:
-            return modulo
+       (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))

-        if context is None:
-            context = getcontext()
+        ;; 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
+            ((cond
+             ((= (modulo last_digit 2) 0)
+              ;; quick test for power of 2
+              (twix
+               ((not (= (logand xc (- xc)) xc))
+                None)
+               ;; now xc is a power of 2; e is its exponent
+               (let ((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))))

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

-        # additional restriction for decimal: the modulus must be less
-        # than 10**prec in absolute value
-            return context._raise_error(InvalidOperation,
-                                        'insufficient precision: pow() 3rd '
-                                        'argument must not have more than '
-                                        'precision digits')
+               ;; 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)

-        # define 0**0 == NaN, for consistency with two-argument pow
-        # (even though it hurts!)
-        if not other and not self:
-            return context._raise_error(InvalidOperation,
-                                        'at least one of pow() 1st argument '
-                                        'and 2nd argument must be nonzero ;'
-                                        '0**0 is not defined')
+               ((> e emax)
+                None)

-        # compute sign of result
-        if other._iseven():
-            sign = 0
-        else:
-            sign = self._sign
+               (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))
+                      (xc        (quotient q xz))
+                      (remainder (modulo   q xc))))
+
+                ((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))))

-        # convert modulo to a Python integer, and self and other to
-        # Decimal integers (i.e. force their exponents to be >= 0)
-        modulo = abs(int(modulo))
-        base = _WorkRep(self.to_integral_value())
-        exponent = _WorkRep(other.to_integral_value())
+               ((>= ye (len (str emax)))
+                None)

-        # compute result using integer pow()
-        base = (base.int % modulo * pow(10, base.exp, modulo)) % modulo
-        for i in range(exponent.exp):
-            base = pow(base, 10, modulo)
-        base = pow(base, exponent.int, modulo)
+               (let ()
+                 (set! e  (_decimal_lshift_exact (* e  yc) ye))
+                 (set! xe (_decimal_lshift_exact (* xe yc) ye)))

-        return _dec_from_triple(sign, str(base), 0)
+                ((or (eq? e None= (eq? xe None))) it
+                None)

-    def _power_exact(self, other, p):
-        """Attempt to compute self**other exactly.
+                ((> e emax)
+                None)

-        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.
+               (begin
+                 (set! xc (expt 2 e))
+                 #f)))

-        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).
-
-        x = _WorkRep(self)
-        xc, xe = x.int, x.exp
-        while xc % 10 == 0:
-            xc //= 10
-            xe += 1
-
-        y = _WorkRep(other)
-        yc, ye = y.int, y.exp
-        while yc % 10 == 0:
-            yc //= 10
-            ye += 1
-
-        # case where xc == 1: result is 10**(xe*y), with xe*y
-        # required to be an integer
-        if xc == 1:
-            xe *= yc
-            # result is now 10**(xe * 10**ye);  xe * 10**ye must be integral
-            while xe % 10 == 0:
-                xe //= 10
-                ye += 1
-            if ye < 0:
-                return None
-            exponent = xe * 10**ye
-            if y.sign == 1:
-                exponent = -exponent
-            # if other is a nonnegative integer, use ideal exponent
-            if other._isinteger() and other._sign == 0:
-                ideal_exponent = self._exp*int(other)
-                zeros = min(exponent-ideal_exponent, p-1)
-            else:
-                zeros = 0
-            return _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.
-        if y.sign == 1:
-            last_digit = xc % 10
-            if last_digit in (2,4,6,8):
-                # quick test for power of 2
-                if xc & -xc != xc:
-                    return None
-                # now xc is a power of 2; e is its exponent
-                e = _nbits(xc)-1
-
-                # 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.
-                emax = p*93//65
-                if ye >= len(str(emax)):
-                    return None
-
-                # Find -e*y and -xe*y; both must be integers
-                e = _decimal_lshift_exact(e * yc, ye)
-                xe = _decimal_lshift_exact(xe * yc, ye)
-                if e is None or xe is None:
-                    return None
-
-                if e > emax:
-                    return None
-                xc = 5**e
-
-            elif last_digit == 5:
-                # e >= log_5(xc) if xc is a power of 5; we have
-                # equality all the way up to xc=5**2658
-                e = _nbits(xc)*28//65
-                xc, remainder = divmod(5**e, xc)
-                if remainder:
-                    return None
-                while xc % 5 == 0:
-                    xc //= 5
-                    e -= 1
-
-                # 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).
-                emax = p*10//3
-                if ye >= len(str(emax)):
-                    return None
-
-                e = _decimal_lshift_exact(e * yc, ye)
-                xe = _decimal_lshift_exact(xe * yc, ye)
-                if e is None or xe is None:
-                    return None
-
-                if e > emax:
-                    return None
-                xc = 2**e
-            else:
-                return None
+             (else
+              None)) it it)

-            if xc >= 10**p:
-                return None
-            xe = -e-xe
+           ((>= xc (expt 10 p)) it it)
+
+            (begin
+             (set! xe (+ (- e) (- xe)))
return _dec_from_triple(0, str(xc), xe)

# now y is positive; find m and n such that y = m/n
context._raise_error(exception)

else:
-            ans = ans._fix(context)
+          ans = ans._fix(context)

return ans

@@ -39,20 +39,28 @@ explicitly tell it to not update etc.
(define (pk-obj o)
(pk 'start-pk-obj)
(let ((h (slot-ref o 'h)))
-    (hash-for-each (lambda (k v) (pk k)) h)
+    (hash-for-each (lambda (k v)
+                    (if (member k '(__name__ __qualname__))
+                        (pk k v)
+                        (pk k))) h)
+
(pk 'finished-obj)
-    (aif cl (hash-ref h '__class__)
-         (if (is-a? cl <p>)
-             (if (hash-table? (slot-ref cl 'h))
-                 (hash-for-each (lambda (k v)
-                                  (if (member k '(__name__ __qualname__))
-                                      (pk k v)
-                                      (pk k)))
-                                (slot-ref cl 'h))
-                 (pk 'no-hash-table))
-             (pk 'no-class))
-         (pk 'false-class)))
-  (pk 'end-pk-obj))
+
+    (let lp ((l (ref o '__mro__ '())))
+      (if (pair? l)
+         (let ((cl (car l)))
+           (if (is-a? cl <p>)
+               (if (hash-table? (slot-ref cl 'h))
+                   (hash-for-each (lambda (k v)
+                                    (if (member k '(__name__ __qualname__))
+                                        (pk k v)
+                                        (pk k)))
+                                  (slot-ref cl 'h))
+                   (pk 'no-hash-table))
+               (pk 'no-class))
+           (lp (cdr l)))))
+
+    (pk 'end-pk-obj)))

(define fail (cons 'fail '()))

@@ -200,14 +208,15 @@ explicitly tell it to not update etc.

(define (hashforeach a b) (values))

+
(define (new-class0 meta name parents dict . kw)
+  (pk 'new-class0)
(let* ((goops (pylist-ref dict '__goops__))
(p     (kwclass->class kw meta))
(class (make-p p)))
-    (pk 'new-class0)
(slot-set! class 'procedure
(lambda x
-                (create-object class meta goops x)))
+                (create-object class x)))

(if (hash-table? dict)
(hash-for-each
@@ -217,7 +226,7 @@ explicitly tell it to not update etc.
(lambda (k v) k (set class k v))
dict))

-    (let((mro (ref class '__mro__)))
+    (let ((mro (ref class '__mro__)))
(if (pair? mro)
(let ((p (car mro)))
(aif it (ref p '__zub_classes__)
@@ -227,15 +236,18 @@ explicitly tell it to not update etc.
(aif it (ref p '__init_subclass__)
(apply it class p #f kw)
#f))))
-    (set class '__mro__ (cons class (ref class '__mro__)))
+    (set class '__mro__ (cons class (find-in-class-and-parents
+                                    class '__mro__ '())))
class))

(define (new-class meta name parents dict kw)
+  (pk 'new-class)
(aif it (and meta (ficap meta '__new__ #f))
(apply it meta name parents dict kw)
(apply new-class0 meta name parents dict kw)))

(define (type- meta name parents dict keys)
+  (pk 'type-)
(let ((class (new-class meta name parents dict keys)))
(aif it (and meta (find-in-class meta '__init__ #f))
(it class name parents dict keys)
@@ -244,11 +256,13 @@ explicitly tell it to not update etc.

(define (the-create-object class x)
+  (pk 'the-create-object)
(let* ((meta  (ref class '__class__))
(goops (ref class '__goops__))
-         (obj   (aif it (ficap class '__new__ #f)
-                     (it)
+         (obj   (aif it (pk '__new__ (ficap class '__new__ #f))
+                     (begin (pk-obj class) (apply it class x))
(make-object class meta goops))))
+
(aif it (ref obj '__init__)
(apply it x)
#f)
@@ -261,15 +275,20 @@ explicitly tell it to not update etc.

obj))

-(define (create-object class meta goops x)
-  (with-fluids ((*make-class* #t))
-    (aif it (ficap meta '__call__ #f)
-         (apply it class x)
-         (the-create-object class x))))
+(define (create-object class x)
+  (pk 'create-object)
+  (if (pk 'type? (pytype? class))
+      (apply type-call class x)
+      (let ((meta (find-in-class class '__class__ #f)))
+       (with-fluids ((*make-class* #t))
+         (aif it (ficap meta '__call__ #f)
+              (apply it class x)
+              (the-create-object class x))))))

(define type-call
(lambda (class . l)
-    (if (pytype? class)
+    (pk 'type-call)
+    (if (pk 'type? (pytype? class))
(apply (case-lambda
((meta obj)
(ref obj '__class__ 'None))
@@ -284,6 +303,7 @@ explicitly tell it to not update etc.
(make-hash-table)))

(define (create-class meta name parents gen-methods keys)
+  (pk 'create-class)
(let ((dict (gen-methods (get-dict meta name parents))))
(aif it (ref meta '__class__)
(aif it (find-in-class it '__call__ #f)
@@ -292,6 +312,7 @@ explicitly tell it to not update etc.
(type- meta name parents dict keys))))

(define (make-object class meta goops)
+  (pk 'make-object)
(let ((obj (make-p goops)))
(set obj '__class__ class)
obj))