summaryrefslogtreecommitdiff
path: root/modules/language/python/module
diff options
context:
space:
mode:
Diffstat (limited to 'modules/language/python/module')
-rw-r--r--modules/language/python/module/decimal.scm830
1 files changed, 439 insertions, 391 deletions
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
generated by cgit v1.2.3 (git 2.46.0) at 2024-10-19 04:40:22 +0000