From be1c8b5d06e50eed7df45c29ba98bc0c913b4207 Mon Sep 17 00:00:00 2001 From: Stefan Israelsson Tampe Date: Fri, 13 Apr 2018 15:06:25 +0200 Subject: small small steps ... --- modules/language/python/module/decimal.scm | 708 ++++++++++++++++------------- 1 file changed, 383 insertions(+), 325 deletions(-) (limited to 'modules/language') diff --git a/modules/language/python/module/decimal.scm b/modules/language/python/module/decimal.scm index 55a6675..1411be8 100644 --- a/modules/language/python/module/decimal.scm +++ b/modules/language/python/module/decimal.scm @@ -1337,375 +1337,433 @@ (let ((ans (_dec_from_triple sign, (str coeff) exp))) ((ref ans '_fix) context)))))) - def _divide(self, other, context): - """Return (self // other, self % other), to context.prec precision. + (define _divide + (lambda (self other context) + "Return (self // other, self % other), to context.prec precision. Assumes that neither self nor other is a NaN, that self is not infinite and that other is nonzero. - """ - sign = self._sign ^ other._sign - if other._isinfinity(): - ideal_exp = self._exp - else: - ideal_exp = min(self._exp, other._exp) - - expdiff = self.adjusted() - other.adjusted() - if not self or other._isinfinity() or expdiff <= -2: - return (_dec_from_triple(sign, '0', 0), - self._rescale(ideal_exp, context.rounding)) - if expdiff <= context.prec: - op1 = _WorkRep(self) - op2 = _WorkRep(other) - if op1.exp >= op2.exp: - op1.int *= 10**(op1.exp - op2.exp) - else: - op2.int *= 10**(op2.exp - op1.exp) - q, r = divmod(op1.int, op2.int) - if q < 10**context.prec: - return (_dec_from_triple(sign, str(q), 0), - _dec_from_triple(self._sign, str(r), ideal_exp)) - - # Here the quotient is too large to be representable - ans = context._raise_error(DivisionImpossible, - 'quotient too large in //, % or divmod') - return ans, ans - - def __rtruediv__(self, other, context=None): - """Swaps self/other and returns __truediv__.""" - other = _convert_other(other) - if other is NotImplemented: - return other - return other.__truediv__(self, context=context) + " + (apply values + (twix + (let (let ((sign + (logxor (ref self '_sign) + (ref other '_sign))) + (ideal_exp + (if ((ref other '_isinfinity)) + (ref self '_exp) + (min (ref self 'exp) (ref other '_exp)))) + (expdiff + (- ((ref self 'adjusted)) ((ref other 'adjusted))))))) + + ((or (not (bool self)) + ((ref other '_isinfinity)) + (<= expdiff -1)) it + (list (_dec_from_tripple sign "0" 0) + ((ref self '_rescale) ideal_exp (cx-rounding context)))) + + ((if (<= expdiff (cx-prec context)) + (let ((op1 (_WorkRep self)) + (op2 (_WorkRep other))) + (if (>= (ref op1 'exp) (ref op2 'exp)) + (set op1 'int (* (ref op1 'int) + (expt 10 (- (ref op1 'exp) + (ref op2 'exp))))) + (set op2 'int (* (ref op2 'int) + (expt 10 (- (ref op2 'exp) + (ref op1 'exp)))))) + (call-with-values (lambda () (divmod (ref op1 'int) + (ref op2 'int))) + (lambda (q r) + (if (< q (expt 10 (cx-prec context))) + (list (_dec_from_triple sign (str q) 0) + (_dec_from_triple (ref self '_sign) + (str r) + ideal_exp)) + #f)))) + #f) it it) + (begin + ;; Here the quotient is too large to be representable + (let ((ans ((cx-raise context) DivisionImpossible + "quotient too large in //, % or divmod"))) + (list ans ans))))))) - def __divmod__(self, other, context=None): - """ + (define __rtruediv__ + (lam (self other (= context None)) + ""Swaps self/other and returns __truediv__."" + (twix + ((norm-op other) it it) + ((ref other '__truediv__) self #:context context)))) + + (define __divmod__ + (lam (self other (= context None)) + " Return (self // other, self % other) - """ - other = _convert_other(other) - if other is NotImplemented: - return other - - if context is None: - context = getcontext() + " + (apply values + (twix + ((norm-op other) it it) - ans = self._check_nans(other, context) - if ans: - return (ans, ans) + (let (get-context context)) + + ((add-special o1 o2 context) it it) + + (((ref self '_check_nans) other context) it + (list it it)) - sign = self._sign ^ other._sign - if self._isinfinity(): - if other._isinfinity(): - ans = context._raise_error(InvalidOperation, 'divmod(INF, INF)') - return ans, ans - else: - return (_SignedInfinity[sign], - context._raise_error(InvalidOperation, 'INF % x')) + (let (let ((sign + (logxor (ref self '_sign) + (ref other '_sign)))))) - if not other: - if not self: - ans = context._raise_error(DivisionUndefined, 'divmod(0, 0)') - return ans, ans - else: - return (context._raise_error(DivisionByZero, 'x // 0', sign), - context._raise_error(InvalidOperation, 'x % 0')) + (((ref self '_isinfinity)) it + (if ((ref other '_isinfinity)) + (let ((ans ((cx-error context) InvalidOperation + "divmod(INF, INF)"))) + (list ans ans)) + (list (list-ref _SignedInfinity sign) + ((cx-raise context) InvalidOperation, "INF % x")))) - quotient, remainder = self._divide(other, context) - remainder = remainder._fix(context) - return quotient, remainder + ((not (bool other)) it + (if (not (bool self)) + (let ((ans ((cx-error context) DivisionUndefined + "divmod(0, 0)"))) + (list ans ans)) + (list ((cx-error context) DivisionByZero "x // 0" sign) + ((cx-error context) InvalidOperation "x % 0")))) + + (call-with-values (lambda () ((ref self '_divide) other context)) + (lambda (quotient remainder) + (let ((remainder ((ref remainder '_fix) context))) + (list quotient remainder)))))))) - def __rdivmod__(self, other, context=None): - """Swaps self/other and returns __divmod__.""" - other = _convert_other(other) - if other is NotImplemented: - return other - return other.__divmod__(self, context=context) + (define __rdivmod__ + (lam (self other (= context None)) + "Swaps self/other and returns __divmod__." + (twix + ((norm-op other) it it) + ((ref other '__divmod__) self #:context context)))) - def __mod__(self, other, context=None): - """ + (define __mod__ + (lam (self other (= context None)) + " self % other - """ - other = _convert_other(other) - if other is NotImplemented: - return other - - if context is None: - context = getcontext() + " + (twix + ((norm-op other) it it) - ans = self._check_nans(other, context) - if ans: - return ans + (let (get-context context)) - if self._isinfinity(): - return context._raise_error(InvalidOperation, 'INF % x') - elif not other: - if self: - return context._raise_error(InvalidOperation, 'x % 0') - else: - return context._raise_error(DivisionUndefined, '0 % 0') + ((bin-special o1 o2 context) it it) - remainder = self._divide(other, context)[1] - remainder = remainder._fix(context) - return remainder + (((ref self '_isinfinity)) it + ((cx-error context) InvalidOperation "INF % x")) + + ((not (bool other)) + (if (bool self) + ((cx-error context) InvalidOperation "x % 0") + ((cx-error context) DivisionUndefined "0 % 0"))) - def __rmod__(self, other, context=None): - """Swaps self/other and returns __mod__.""" - other = _convert_other(other) - if other is NotImplemented: - return other - return other.__mod__(self, context=context) + (let* ((remainder ((ref self '_divide) other context))) + ((ref remainder '_fix) context))))) + + (define __rmod__ + (lam (self other (= context None)) + "Swaps self/other and returns __mod__." + (twix + ((norm-op other) it it) + ((ref other '__mod__) self #:context context)))) - def remainder_near(self, other, context=None): - """ + (define remainder_near + (lambda (self other (= context None)) + " Remainder nearest to 0- abs(remainder-near) <= other/2 - """ - if context is None: - context = getcontext() - - other = _convert_other(other, raiseit=True) - - ans = self._check_nans(other, context) - if ans: - return ans - - # self == +/-infinity -> InvalidOperation - if self._isinfinity(): - return context._raise_error(InvalidOperation, - 'remainder_near(infinity, x)') - - # other == 0 -> either InvalidOperation or DivisionUndefined - if not other: - if self: - return context._raise_error(InvalidOperation, - 'remainder_near(x, 0)') - else: - return context._raise_error(DivisionUndefined, - 'remainder_near(0, 0)') - - # other = +/-infinity -> remainder = self - if other._isinfinity(): - ans = Decimal(self) - return ans._fix(context) - - # self = 0 -> remainder = self, with ideal exponent - ideal_exponent = min(self._exp, other._exp) - if not self: - ans = _dec_from_triple(self._sign, '0', ideal_exponent) - return ans._fix(context) - - # catch most cases of large or small quotient - expdiff = self.adjusted() - other.adjusted() - if expdiff >= context.prec + 1: - # expdiff >= prec+1 => abs(self/other) > 10**prec - return context._raise_error(DivisionImpossible) - if expdiff <= -2: - # expdiff <= -2 => abs(self/other) < 0.1 - ans = self._rescale(ideal_exponent, context.rounding) - return ans._fix(context) - - # adjust both arguments to have the same exponent, then divide - op1 = _WorkRep(self) - op2 = _WorkRep(other) - if op1.exp >= op2.exp: - op1.int *= 10**(op1.exp - op2.exp) - else: - op2.int *= 10**(op2.exp - op1.exp) - q, r = divmod(op1.int, op2.int) - # remainder is r*10**ideal_exponent; other is +/-op2.int * - # 10**ideal_exponent. Apply correction to ensure that - # abs(remainder) <= abs(other)/2 - if 2*r + (q&1) > op2.int: - r -= op2.int - q += 1 - - if q >= 10**context.prec: - return context._raise_error(DivisionImpossible) - - # result has same sign as self unless r is negative - sign = self._sign - if r < 0: - sign = 1-sign - r = -r - - ans = _dec_from_triple(sign, str(r), ideal_exponent) - return ans._fix(context) + " + (twix + ((norm-op other) it it) - def __floordiv__(self, other, context=None): - """self // other""" - other = _convert_other(other) - if other is NotImplemented: - return other + (let (get-context context)) - if context is None: - context = getcontext() + ((bin-special self other context) it it) - ans = self._check_nans(other, context) - if ans: - return ans + ;; self == +/-infinity -> InvalidOperation + (((ref self '_isinfinity)) it + ((cx-error context) InvalidOperation "remainder_near(infinity, x)")) - if self._isinfinity(): - if other._isinfinity(): - return context._raise_error(InvalidOperation, 'INF // INF') - else: - return _SignedInfinity[self._sign ^ other._sign] + ;; other == 0 -> either InvalidOperation or DivisionUndefined + ((not (bool other)) it + (if (not (bool self)) + ((cx-error context) InvalidOperation "remainder_near(x, 0)") + ((cx-error context) DivisionUndefined "remainder_near(0, 0)"))) - if not other: - if self: - return context._raise_error(DivisionByZero, 'x // 0', - self._sign ^ other._sign) - else: - return context._raise_error(DivisionUndefined, '0 // 0') + ;; other = +/-infinity -> remainder = self + (((ref other '_isinfinity())) it + (let ((ans (Decimal self))) + ((ref ans '_fix) context))) - return self._divide(other, context)[0] + ;; self = 0 -> remainder = self, with ideal exponent + (let (let ((ideal_exponent (min (ref self '_exp) (ref other '_exp)))))) - def __rfloordiv__(self, other, context=None): - """Swaps self/other and returns __floordiv__.""" - other = _convert_other(other) - if other is NotImplemented: - return other - return other.__floordiv__(self, context=context) - - def __float__(self): - """Float representation.""" - if self._isnan(): - if self.is_snan(): - raise ValueError("Cannot convert signaling NaN to float") - s = "-nan" if self._sign else "nan" - else: - s = str(self) - return float(s) + ((not (bool self)) it + (let ((ans (_dec_from_triple (ref self '_sign) "0" ideal_exponent))) + ((ref ans '_fix) context))) - def __int__(self): - """Converts self to an int, truncating if necessary.""" - if self._is_special: - if self._isnan(): - raise ValueError("Cannot convert NaN to integer") - elif self._isinfinity(): - raise OverflowError("Cannot convert infinity to integer") - s = (-1)**self._sign - if self._exp >= 0: - return s*int(self._int)*10**self._exp - else: - return s*int(self._int[:self._exp] or '0') + ;; catch most cases of large or small quotient + (let (let ((expdiff + (- ((ref self 'adjusted)) ((red other 'adjusted))))))) + + ((>= expdiff (+ (cx-prec context) 1)) it + ;; expdiff >= prec+1 => abs(self/other) > 10**prec + ((cx-error context) DivisionImpossible)) + + ((<= expdiff -2) it + ;; expdiff <= -2 => abs(self/other) < 0.1 + (let ((ans ((ref self '_rescale) + ideal_exponent (cx-rounding context)))) + ((ref ans '_fix) context))) + + (let ((op1 (_WorkRep self)) + (op2 (_WorkRep other))) + + ;; adjust both arguments to have the same exponent, then divide + (if (>= (ref op1 'exp) (ref op2 'exp)) + (set op1 'int (* (ref op1 'int) + (expt 10 (- (ref op1 'exp) (ref op2 'exp))))) + (set op2 'int (* (ref op2 'int) + (expt 10 (- (ref op2 'exp) (ref op1 'exp)))))) + + (call-with-values (lambda () (divmod (ref op1 'int) (ref op2 'int))) + (lambda (q r) + + ;; remainder is r*10**ideal_exponent; other is +/-op2.int * + ;; 10**ideal_exponent. Apply correction to ensure that + ;; abs(remainder) <= abs(other)/2 + (if (> (+ (* 2 r) + (logand q 1)) (ref op2 'int)) + (set! r (- r (ref op2 'int))) + (set! q (+ q 1))) + + (if (>= q (expt 10 (cx-prec context))) + ((cx-error context) DivisionImpossible) + (let ((sign (ref self '_sign))) + (if (< r 0) + (set! sign (- 1 sign)) + (set! r (- r))) + (let ((ans (_dec_from_triple sign (str r) ideal_exponent))) + ((ref ans '_fix) context)))))))))) + + (define __floordiv__ + (lambda (self other (= context None)) + "self // other" + (twix + ((norm-op other) it it) - __trunc__ = __int__ + (let (get-context context)) - def real(self): - return self - real = property(real) + ((bin-special self other context) it it) - def imag(self): - return Decimal(0) - imag = property(imag) + (((ref self '_isinfinity)) it + (if ((ref other '_isinfinity)) + ((cx-error context) InvalidOperation "INF // INF") + (pylist-ref _SignedInfinity (logxor (ref self '_sign) + (ref other '_sign))))) - def conjugate(self): - return self + ((not (bool other)) it + (if (bool self) + ((cx-error context) DivisionByZero "x // 0" + (logxor (ref self '_sign) (ref other '_sign))) + ((cx-error context) DivisionUndefined "0 // 0"))) - def __complex__(self): - return complex(float(self)) + ((ref self '_divide) other context)))) - def _fix_nan(self, context): - """Decapitate the payload of a NaN to fit the context""" - payload = self._int + (define __rfloordiv__ + (lam (self other (= context None)) + "Swaps self/other and returns __floordiv__." + (twix + ((norm-op other) it it) + ((ref other '__floordiv__) self #:context context)))) - # maximum length of payload is precision if clamp=0, - # precision-1 if clamp=1. - max_payload_len = context.prec - context.clamp - if len(payload) > max_payload_len: - payload = payload[len(payload)-max_payload_len:].lstrip('0') - return _dec_from_triple(self._sign, payload, self._exp, True) - return Decimal(self) + (define __float__ + (lambda (self) + "Float representation." + (if ((ref self '_isnan)) + (if ((ref self 'is_snan)) + (raise (ValueError "Cannot convert signaling NaN to float")) + (if (= (ref self '_sign)) + (- (nan)) + (nan))) + (if ((ref self '_isspecial)) + (if (= (ref self '_sign)) + (- (inf)) + (inf))) + (float (str self))))) + + (define __int__ + (lambda (self) + "Converts self to an int, truncating if necessary." + (if ((ref self '_isnan)) + (raise (ValueError "Cannot convert NaN to integer")) + (if ((ref self '_isspecial)) + (raise (OverflowError "Cannot convert infinity to integer")) + (let ((s (if (= (ref self '_sign) 1) -1 1))) + (if (>= (ref self '_exp) 0) + (* s (int (ref self '_int)) (expt 10 (ref self '_exp))) + (* s (int (or (bool (py-slice (ref self '_int) + None (ref self '_exp) None)) + "0"))))))))) + + (define __trunc__ __int__) + + (define real + (property (lambda (self) self))) + + (define imag + (property + (lambda (self) + (Decimal 0)))) - def _fix(self, context): - """Round if it is necessary to keep self within prec precision. + (define conjugate + (lambda (self) self)) + + (define __complex__ + (lambda (self) + (complex (float self)))) + + (define _fix_nan + (lambda (self context) + "Decapitate the payload of a NaN to fit the context" + (let ((payload (ref self '_int)) + + ;; maximum length of payload is precision if clamp=0, + ;; precision-1 if clamp=1. + (max_payload_len + (- (ref context 'prec) + (ref context 'clamp)))) + + (if (> (len payload) max_payload_len) + (let ((payload (py-lstrip + (pylist-slice payload + (- (len payload) max_payload_len) + None None) "0"))) + (_dec_from_triple (ref self '_sign) payload (ref self '_exp) + #t)) + (Decimal self))))) + + (define _fix + (lambda (self context) + "Round if it is necessary to keep self within prec precision. Rounds and fixes the exponent. Does not raise on a sNaN. Arguments: self - Decimal instance context - context used. - """ - - if self._is_special: - if self._isnan(): - # decapitate payload if necessary - return self._fix_nan(context) - else: - # self is +/-Infinity; return unaltered - return Decimal(self) - - # if self is zero then exponent should be between Etiny and - # Emax if clamp==0, and between Etiny and Etop if clamp==1. - Etiny = context.Etiny() - Etop = context.Etop() - if not self: - exp_max = [context.Emax, Etop][context.clamp] - new_exp = min(max(self._exp, Etiny), exp_max) - if new_exp != self._exp: - context._raise_error(Clamped) - return _dec_from_triple(self._sign, '0', new_exp) - else: - return Decimal(self) - - # exp_min is the smallest allowable exponent of the result, - # equal to max(self.adjusted()-context.prec+1, Etiny) - exp_min = len(self._int) + self._exp - context.prec - if exp_min > Etop: - # overflow: exp_min > Etop iff self.adjusted() > Emax - ans = context._raise_error(Overflow, 'above Emax', self._sign) - context._raise_error(Inexact) - context._raise_error(Rounded) - return ans + " - self_is_subnormal = exp_min < Etiny - if self_is_subnormal: - exp_min = Etiny - - # round if self has too many digits - if self._exp < exp_min: - digits = len(self._int) + self._exp - exp_min - if digits < 0: - self = _dec_from_triple(self._sign, '1', exp_min-1) - digits = 0 - rounding_method = self._pick_rounding_function[context.rounding] - changed = rounding_method(self, digits) - coeff = self._int[:digits] or '0' - if changed > 0: - coeff = str(int(coeff)+1) - if len(coeff) > context.prec: - coeff = coeff[:-1] - exp_min += 1 - - # check whether the rounding pushed the exponent out of range - if exp_min > Etop: - ans = context._raise_error(Overflow, 'above Emax', self._sign) - else: - ans = _dec_from_triple(self._sign, coeff, exp_min) - - # raise the appropriate signals, taking care to respect - # the precedence described in the specification - if changed and self_is_subnormal: - context._raise_error(Underflow) - if self_is_subnormal: - context._raise_error(Subnormal) - if changed: - context._raise_error(Inexact) - context._raise_error(Rounded) - if not ans: - # raise Clamped on underflow to 0 - context._raise_error(Clamped) - return ans + (twix + (((ref self '_is_special)) it + (if ((ref self '_isnan)) + ;; decapitate payload if necessary + ((ref self '_fix_nan) context) + + ;; self is +/-Infinity; return unaltered + (Decimal self))) + + ;; if self is zero then exponent should be between Etiny and + ;; Emax if clamp==0, and between Etiny and Etop if clamp==1. + (let ((Etiny (cx-etiny context)) + (Etop (cx-etop context)))) + + ((not (bool self)) it + (let ((exp_max (if (= (cx-clamp context) 0) + (cx-emax context) + Etop)) + (new_exp (min (max (ref self '_exp) Etiny) exp_max))) + (if (not (= new_exp (ref self '_exp))) + (begin + ((cx-error context) Clamped) + (_dec_from_triple (ref self '_sign) "0" new_exp)) + (Decimal self)))) + + ;; exp_min is the smallest allowable exponent of the result, + ;; equal to max(self.adjusted()-context.prec+1, Etiny) + (let ((exp_min (+ (len (ref self '_int)) + (ref self '_exp) + (- (cx-prec context))))))) + ((> exp_min Etop) it + ;; overflow: exp_min > Etop iff self.adjusted() > Emax + (let ((ans ((cx-error context) Overflow "above Emax" + (ref self '_sign)))) + ((cx-error context) Inexact) + ((cx-error context) Rounded) + ans)) + + (let* ((self_is_subnormal (< exp_min Etiny)) + (exp_min (if self_is_subnormal Eriny exp_min))))) + + ;; round if self has too many digits + ((< self._exp exp_min) it + (let ((digits (+ (len (ref self '_int)) + (ref self '_exp) + (- exp_min)))) + (if (< digits 0) + (set! self (_dec_from_triple (ref self '_sign) + "1" (- exp_min 1))) + (set! digits 0)) + + (let* ((ans #f) + (rounding_method (pylist-ref + (ref self '_pick_rounding_function) + (cx-rounding context))) + (changed (rounding_method self digits)) + (coeff (or (bool (pylist-slice (ref self '_int) + None digits None)) "0"))) + (if (> changed 0) + (begin + (set! coeff (str (+ (int coeff) 1))) + (if (> (len coeff) (cx-prec context)) + (begin + (set! coeff (pylist-clice coeff None -1 None)) + (set! exp_min (+ exp_min 1)))))) + + ;; check whether the rounding pushed the exponent out of range + (if (> exp_min Etop) + (set! ans + ((cx-error context) Overflow "above Emax" + (ref self '_sign))) + (set! ans (_dec_from_triple (ref self '_sign) coeff exp_min))) + + ;; raise the appropriate signals, taking care to respect + ;; the precedence described in the specification + (if (and changed self_is_subnormal) + ((cx-error context) Underflow)) + (if self_is_subnormal + ((cx-error context) Subnormal)) + (if changed + ((cx-error context) Inexact)) + + ((cx-error context) Rounded) + + (if (not (bool ans)) + ;; raise Clamped on underflow to 0 + ((cx-error context) Clamped)) + + ans))) + (begin + (if self_is_subnormal + ((cx-error context) Subnormal)) - if self_is_subnormal: - context._raise_error(Subnormal) - # fold down if clamp == 1 and self has too few digits - if context.clamp == 1 and self._exp > Etop: - context._raise_error(Clamped) - self_padded = self._int + '0'*(self._exp - Etop) - return _dec_from_triple(self._sign, self_padded, Etop) + ;; fold down if clamp == 1 and self has too few digits + (if (and (= (cx-clamp context) 1) (> (ref self '_exp) Etop)) + (begin + ((cx-error context) Clamped) + (let ((self_padded (+ (ref self '_int) + (* "0" + (- (ref self '_exp) Etop))))) + (_dec_from_triple (ref self '_sign) self_padded Etop))) + + ;; here self was representable to begin with; return unchanged + (Decimal self)))) - # here self was representable to begin with; return unchanged - return Decimal(self) # for each of the rounding functions below: # self is a finite, nonzero Decimal -- cgit v1.2.3