From 9f09b5ef3c86e26d0b5462e267633f25d0240df5 Mon Sep 17 00:00:00 2001 From: Stefan Israelsson Tampe Date: Mon, 16 Apr 2018 15:55:30 +0200 Subject: small steps furter regarding decimal --- modules/language/python/module/decimal.scm | 830 +++++++++++++++-------------- 1 file changed, 439 insertions(+), 391 deletions(-) (limited to 'modules/language/python') diff --git a/modules/language/python/module/decimal.scm b/modules/language/python/module/decimal.scm index 1411be8..2cf1ce5 100644 --- a/modules/language/python/module/decimal.scm +++ b/modules/language/python/module/decimal.scm @@ -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") @@ -1631,7 +1633,7 @@ (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 @@ -1765,85 +1767,88 @@ (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 @@ -1888,38 +1893,39 @@ >>> 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 @@ -1927,15 +1933,15 @@ 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): - """Fused multiply-add. + (define fma + (lam (self other third (= context None)) + "Fused multiply-add. Returns self*other+third with no rounding of the intermediate product self*other. @@ -1943,299 +1949,341 @@ 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. - return product.__add__(third, context) + 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 - if modulo.adjusted() >= context.prec: - 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 @@ -2512,7 +2560,7 @@ context._raise_error(exception) else: - ans = ans._fix(context) + ans = ans._fix(context) return ans -- cgit v1.2.3