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.scm1881
1 files changed, 1371 insertions, 510 deletions
diff --git a/modules/language/python/module/decimal.scm b/modules/language/python/module/decimal.scm
index 2cf1ce5..a4c31c7 100644
--- a/modules/language/python/module/decimal.scm
+++ b/modules/language/python/module/decimal.scm
@@ -2145,7 +2145,7 @@
(begin
(set! xc (/ xc n))
(set! xe (+ xe 1))
- (lp)))))
+ (lp))))))
(let* ((x (_WorkRep self))
(xc (ref x 'int))
@@ -2284,72 +2284,102 @@
(begin
(set! xe (+ (- e) (- xe)))
- return _dec_from_triple(0, str(xc), xe)
+ (_dec_from_triple 0 (str xc) xe)))))
- # now y is positive; find m and n such that y = m/n
- if ye >= 0:
- m, n = yc*10**ye, 1
- else:
- if xe != 0 and len(str(abs(yc*xe))) <= -ye:
- return None
- xc_bits = _nbits(xc)
- if xc != 1 and len(str(abs(yc)*xc_bits)) <= -ye:
- return None
- m, n = yc, 10**(-ye)
- while m % 2 == n % 2 == 0:
- m //= 2
- n //= 2
- while m % 5 == n % 5 == 0:
- m //= 5
- n //= 5
-
- # compute nth root of xc*10**xe
- if n > 1:
- # if 1 < xc < 2**n then xc isn't an nth power
- if xc != 1 and xc_bits <= n:
- return None
-
- xe, rem = divmod(xe, n)
- if rem != 0:
- return None
-
- # compute nth root of xc using Newton's method
- a = 1 << -(-_nbits(xc)//n) # initial estimate
- while True:
- q, r = divmod(xc, a**(n-1))
- if a <= q:
- break
- else:
- a = (a*(n-1) + q)//n
- if not (a == q and r == 0):
- return None
- xc = a
-
- # now xc*10**xe is the nth root of the original xc*10**xe
- # compute mth power of xc*10**xe
-
- # if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
- # 10**p and the result is not representable.
- if xc > 1 and m > p*100//_log10_lb(xc):
- return None
- xc = xc**m
- xe *= m
- if xc > 10**p:
- return None
-
- # by this point the result *is* exactly representable
- # adjust the exponent to get as close as possible to the ideal
- # exponent, if necessary
- str_xc = str(xc)
- if other._isinteger() and other._sign == 0:
- ideal_exponent = self._exp*int(other)
- zeros = min(xe-ideal_exponent, p-len(str_xc))
- else:
- zeros = 0
- return _dec_from_triple(0, str_xc+'0'*zeros, xe-zeros)
+ ;; now y is positive; find m and n such that y = m/n
+ (let ((m #f) (n #f) (xc_bits (_nbits xc))))
+ ((if (>= ye 0) it
+ (begin
+ (set! m (* yc (expt 10 ye)))
+ (set! n 1)
+ #f)
+ (twix
+ ((and (not (= xe 0)) (<= (len (str (abs (* yc xe)))) (- ye))) it
+ None)
+
+ ((and (not (= xc 1))
+ (<= (len (str (* (abs yc) xc_bits))) (- ye))) it
+ None)
+
+ (begin
+ (set! m yc)
+ (set! n (expt 10 (- ye)))
+
+ (let lp()
+ (if (and (= (modulo m 2) 0) (= (modulo n 2) 0))
+ (begin
+ (set! m (quotient m 2))
+ (set! n (quotient n 2)))))
+ (let lp()
+ (if (and (= (modulo m 5) 0) (= (modulo n 5) 0))
+ (begin
+ (set! m (quotient m 5))
+ (set! n (quotient n 5)))))
+ #f))) it it)
+
+ ;; compute nth root of xc*10**xe
+ ((if (> n 1)
+ (begin
+ (twix
+ ;; if 1 < xc < 2**n then xc isn't an nth power
+ ((and (not (= xc 1)) (<= xc_bits n)) it
+ None)
+
+ ((not (= (modulo xe n) 0)) it
+ None)
+
+ (begin
+ (let ((a (ash 1 (- (quotient (- xc_bits) n)))))
+ (set! xe (quotient xe n))
+
+ ;; compute nth root of xc using Newton's method
+ (let lp ()
+ (let* ((x (expt a (- n 1)))
+ (q (quotient xc x))
+ (r (modulo xc x)))
+ (if (<= a q)
+ (if (not (and (= a q) (= r 0)))
+ None
+ (begin
+ (set! xc a)
+ #f))
+ (begin
+ (set! a (quotient (+ (* a (- n 1)) q) n))
+ (lp)))))))))
+ #f) it it)
+
+ ;; now xc*10**xe is the nth root of the original xc*10**xe
+ ;; compute mth power of xc*10**xe
- def __pow__(self, other, modulo=None, context=None):
- """Return self ** other [ % modulo].
+ ;; if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
+ ;; 10**p and the result is not representable.
+ ((and (> xc 1) (> m (quotient (* p 100) (_log10_lb xc)))) it
+ None)
+
+ (let ()
+ (set! xc (expt xc m))
+ (set! xe (xe * m)))
+
+ ((> xc (expt 10 p)) it
+ None)
+
+ (begin
+ ;; by this point the result *is* exactly representable
+ ;; adjust the exponent to get as close as possible to the ideal
+ ;; exponent, if necessary
+ (let* ((str_xc (str xc))
+ (zeros (if (and ((ref other '_isinteger))
+ (= (ref other '_sign) 0))
+ (let ((ideal_exponent
+ (* (ref self '_exp) (int other))))
+ (min (- xe ideal_exponent)
+ (- p (len str_xc))))
+ 0)))
+ (_dec_from_triple 0 (+ str_xc (* '0' zeros)) (- xe zeros)))))))
+
+ (define __pow__
+ (lam (self other (= modulo None) (= context None))
+ "Return self ** other [ % modulo].
With two arguments, compute self**other.
@@ -2370,310 +2400,357 @@
result that would be obtained by computing (self**other) %
modulo with unbounded precision, but is computed more
efficiently. It is always exact.
- """
-
- if modulo is not None:
- return self._power_modulo(other, modulo, context)
-
- other = _convert_other(other)
- if other is NotImplemented:
- return other
-
- if context is None:
- context = getcontext()
-
- # either argument is a NaN => result is NaN
- ans = self._check_nans(other, context)
- if ans:
- return ans
-
- # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
- if not other:
- if not self:
- return context._raise_error(InvalidOperation, '0 ** 0')
- else:
- return _One
-
- # result has sign 1 iff self._sign is 1 and other is an odd integer
- result_sign = 0
- if self._sign == 1:
- if other._isinteger():
- if not other._iseven():
- result_sign = 1
- else:
- # -ve**noninteger = NaN
- # (-0)**noninteger = 0**noninteger
- if self:
- return context._raise_error(InvalidOperation,
- 'x ** y with x negative and y not an integer')
- # negate self, without doing any unwanted rounding
- self = self.copy_negate()
-
- # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
- if not self:
- if other._sign == 0:
- return _dec_from_triple(result_sign, '0', 0)
- else:
- return _SignedInfinity[result_sign]
-
- # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
- if self._isinfinity():
- if other._sign == 0:
- return _SignedInfinity[result_sign]
- else:
- return _dec_from_triple(result_sign, '0', 0)
-
- # 1**other = 1, but the choice of exponent and the flags
- # depend on the exponent of self, and on whether other is a
- # positive integer, a negative integer, or neither
- if self == _One:
- if other._isinteger():
- # exp = max(self._exp*max(int(other), 0),
- # 1-context.prec) but evaluating int(other) directly
- # is dangerous until we know other is small (other
- # could be 1e999999999)
- if other._sign == 1:
- multiplier = 0
- elif other > context.prec:
- multiplier = context.prec
- else:
- multiplier = int(other)
-
- exp = self._exp * multiplier
- if exp < 1-context.prec:
- exp = 1-context.prec
- context._raise_error(Rounded)
- else:
- context._raise_error(Inexact)
- context._raise_error(Rounded)
- exp = 1-context.prec
-
- return _dec_from_triple(result_sign, '1'+'0'*-exp, exp)
-
- # compute adjusted exponent of self
- self_adj = self.adjusted()
+ "
- # self ** infinity is infinity if self > 1, 0 if self < 1
- # self ** -infinity is infinity if self < 1, 0 if self > 1
- if other._isinfinity():
- if (other._sign == 0) == (self_adj < 0):
- return _dec_from_triple(result_sign, '0', 0)
- else:
- return _SignedInfinity[result_sign]
-
- # from here on, the result always goes through the call
- # to _fix at the end of this function.
- ans = None
- exact = False
-
- # crude test to catch cases of extreme overflow/underflow. If
- # log10(self)*other >= 10**bound and bound >= len(str(Emax))
- # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
- # self**other >= 10**(Emax+1), so overflow occurs. The test
- # for underflow is similar.
- bound = self._log10_exp_bound() + other.adjusted()
- if (self_adj >= 0) == (other._sign == 0):
- # self > 1 and other +ve, or self < 1 and other -ve
- # possibility of overflow
- if bound >= len(str(context.Emax)):
- ans = _dec_from_triple(result_sign, '1', context.Emax+1)
- else:
- # self > 1 and other -ve, or self < 1 and other +ve
- # possibility of underflow to 0
- Etiny = context.Etiny()
- if bound >= len(str(-Etiny)):
- ans = _dec_from_triple(result_sign, '1', Etiny-1)
-
- # try for an exact result with precision +1
- if ans is None:
- ans = self._power_exact(other, context.prec + 1)
- if ans is not None:
- if result_sign == 1:
- ans = _dec_from_triple(1, ans._int, ans._exp)
- exact = True
-
- # usual case: inexact result, x**y computed directly as exp(y*log(x))
- if ans is None:
- p = context.prec
- x = _WorkRep(self)
- xc, xe = x.int, x.exp
- y = _WorkRep(other)
- yc, ye = y.int, y.exp
- if y.sign == 1:
- yc = -yc
-
- # compute correctly rounded result: start with precision +3,
- # then increase precision until result is unambiguously roundable
- extra = 3
- while True:
- coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
- if coeff % (5*10**(len(str(coeff))-p-1)):
- break
- extra += 3
+ (twix
+ ((not (eq= modulo None)) it
+ ((ref self '_power_modulo) other modulo context))
- ans = _dec_from_triple(result_sign, str(coeff), exp)
-
- # unlike exp, ln and log10, the power function respects the
- # rounding mode; no need to switch to ROUND_HALF_EVEN here
-
- # There's a difficulty here when 'other' is not an integer and
- # the result is exact. In this case, the specification
- # requires that the Inexact flag be raised (in spite of
- # exactness), but since the result is exact _fix won't do this
- # for us. (Correspondingly, the Underflow signal should also
- # be raised for subnormal results.) We can't directly raise
- # these signals either before or after calling _fix, since
- # that would violate the precedence for signals. So we wrap
- # the ._fix call in a temporary context, and reraise
- # afterwards.
- if exact and not other._isinteger():
- # pad with zeros up to length context.prec+1 if necessary; this
- # ensures that the Rounded signal will be raised.
- if len(ans._int) <= context.prec:
- expdiff = context.prec + 1 - len(ans._int)
- ans = _dec_from_triple(ans._sign, ans._int+'0'*expdiff,
- ans._exp-expdiff)
-
- # create a copy of the current context, with cleared flags/traps
- newcontext = context.copy()
- newcontext.clear_flags()
- for exception in _signals:
- newcontext.traps[exception] = 0
-
- # round in the new context
- ans = ans._fix(newcontext)
-
- # raise Inexact, and if necessary, Underflow
- newcontext._raise_error(Inexact)
- if newcontext.flags[Subnormal]:
- newcontext._raise_error(Underflow)
-
- # propagate signals to the original context; _fix could
- # have raised any of Overflow, Underflow, Subnormal,
- # Inexact, Rounded, Clamped. Overflow needs the correct
- # arguments. Note that the order of the exceptions is
- # important here.
- if newcontext.flags[Overflow]:
- context._raise_error(Overflow, 'above Emax', ans._sign)
- for exception in Underflow, Subnormal, Inexact, Rounded, Clamped:
- if newcontext.flags[exception]:
- context._raise_error(exception)
+ ((norm-op other ) it it)
+ (let (get-context context))
- else:
- ans = ans._fix(context)
+ ;; either argument is a NaN => result is NaN
+ (bin-special self other context)
- return ans
+ ;; 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
+ ((not (bool other))
+ (if (not (bool self))
+ ((cx-error context) InvalidOperation "0 ** 0")
+ _One))
+
+ ;; result has sign 1 iff self._sign is 1 and other is an odd integer
+ (let ((result_sign 0)))
+
+ ((if (= (ref self '_sign) 1)
+ (twix
+ ((if ((ref other '_isinteger))
+ (if (not ((ref other '_iseven)))
+ (begin
+ (set! result_sign 1)
+ #f)
+ #f)
+ ;; -ve**noninteger = NaN
+ ;; (-0)**noninteger = 0**noninteger
+ (if (bool self)
+ ((cx-error context) InvalidOperation
+ "x ** y with x negative and y not an integer")
+ #f)) it it)
+ (begin
+ ;; negate self, without doing any unwanted rounding
+ (set! self ((ref self 'copy_negate)))
+ #f))
+ #f) it it)
- def __rpow__(self, other, context=None):
- """Swaps self/other and returns __pow__."""
- other = _convert_other(other)
- if other is NotImplemented:
- return other
- return other.__pow__(self, context=context)
+ ;; 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
+ ((not (bool self)) it
+ (if (= (ref other '_sign) 0)
+ (_dec_from_triple result_sign "0" 0)
+ (pylist-ref _SignedInfinity result_sign)))
+
+ ;; Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
+ (((self '_isinfinity)) it
+ (if (= (ref other '_sign) 0)
+ (pylist-ref _SignedInfinity result_sign)
+ (_dec_from_triple result_sign "0" 0)))
+
+ ;; 1**other = 1, but the choice of exponent and the flags
+ ;; depend on the exponent of self, and on whether other is a
+ ;; positive integer, a negative integer, or neither
+ (let ((prec (cx-prec context))))
+
+ ((equal? self _One) it
+ (let ((exp #f))
+ (if ((ref other '_isinteger))
+ ;; exp = max(self._exp*max(int(other), 0),
+ ;; 1-context.prec) but evaluating int(other) directly
+ ;; is dangerous until we know other is small (other
+ ;; could be 1e999999999)
+ (let ((multiplier
+ (cond
+ ((= (ref other '_sign) 1)
+ 0)
+ ((> other prec)
+ prec)
+ (else
+ (int other)))))
+
+ (set! exp (* (ref self '_exp) multiplier))
+ (if (< exp (- 1 prec))
+ (begin
+ (set! exp (- 1 prec))
+ ((cx-error context) Rounded))))
+ (begin
+ ((cx-error context) Inexact)
+ ((cx-error context) Rounded)
+ (set! exp (- 1 prec))))
+
+ (_dec_from_triple result_sign (+ "1" (* "0" (- exp)) exp))))
+
+ ;; compute adjusted exponent of self
+ (let ((self_adj ((ref self 'adjusted)))))
+
+ ;; self ** infinity is infinity if self > 1, 0 if self < 1
+ ;; self ** -infinity is infinity if self < 1, 0 if self > 1
+ (((ref other '_isinfinity))
+ (if (eq? (= (ref other '_sign) 0)
+ (< self_adj 0))
+ (_dec_from_triple result_sign "0" 0)
+ (pylist-ref _SignedInfinity result_sign)))
+
+ ;; from here on, the result always goes through the call
+ ;; to _fix at the end of this function.
+ (let ((ans None)
+ (exact False)
+ (bound (+ ((ref self '_log10_exp_bound))
+ ((ref other 'adjusted)))))
+
+ ;; crude test to catch cases of extreme overflow/underflow. If
+ ;; log10(self)*other >= 10**bound and bound >= len(str(Emax))
+ ;; then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
+ ;; self**other >= 10**(Emax+1), so overflow occurs. The test
+ ;; for underflow is similar.
+
+ (if (eq? (>= self_adj 0) (= (ref other '_sign) 0))
+ ;; self > 1 and other +ve, or self < 1 and other -ve
+ ;; possibility of overflow
+ (if (>= bound (len (str (cx-emax context))))
+ (set! ans
+ (_dec_from_triple result_sign "1"
+ (+ (cx-emax context) 1))))
+
+ ;; self > 1 and other -ve, or self < 1 and other +ve
+ ;; possibility of underflow to 0
+ (let ((Etiny (cx-etiny context)))
+ (if (>= bound (len (str (- Etiny))))
+ (set! ans (_dec_from_triple result_sign "1" (- Etiny 1))))))
+
+ ;; try for an exact result with precision +1
+ (when (eq? ans None)
+ (set! ans ((ref self '_power_exact) other (+ prec 1)))
+ (when (not (eq? ans None))
+ (if (= result_sign 1)
+ (set! ans (_dec_from_triple 1 (ref ans '_int)
+ (ref ans '_exp))))
+ (set! exact #t)))
+
+ ;; usual case: inexact result, x**y computed directly as exp(y*log(x))
+ (when (eq? ans None)
+ (let* ((p prec)
+
+ (x (_WorkRep self))
+ (xc (ref x 'int))
+ (xe (ref x 'exp))
+
+ (y (_WorkRep other))
+ (yc (ref y 'int))
+ (ye (ref y 'exp)))
+
+ (if (= (ref y 'sign) 1)
+ (set! yc (- yc)))
+
+ ;; compute correctly rounded result: start with precision +3,
+ ;; then increase precision until result is unambiguously roundable
+ (call-with-values
+ (lambda ()
+ (let lp ((extra 3))
+ (call-with-values
+ (lambda () (_dpower xc xe yc ye (+ p extra)))
+ (lambda (coeff exp)
+ (if (modulo coeff
+ (* 5 (expt 10 (- (len (str coeff)) p 1))))
+ (values coeff exp)
+ (lp (+ extra 3)))))))
+ (lambda (coeff exp)
+ (set! ans (_dec_from_triple result_sign (strcoeff) exp))))))
+
+ ;; unlike exp, ln and log10, the power function respects the
+ ;; rounding mode; no need to switch to ROUND_HALF_EVEN here
+
+ ;; There's a difficulty here when 'other' is not an integer and
+ ;; the result is exact. In this case, the specification
+ ;; requires that the Inexact flag be raised (in spite of
+ ;; exactness), but since the result is exact _fix won't do this
+ ;; for us. (Correspondingly, the Underflow signal should also
+ ;; be raised for subnormal results.) We can't directly raise
+ ;; these signals either before or after calling _fix, since
+ ;; that would violate the precedence for signals. So we wrap
+ ;; the ._fix call in a temporary context, and reraise
+ ;; afterwards.
+ (if (and exact (not ((ref other '_isinteger))))
+ (begin
+ ;; pad with zeros up to length context.prec+1 if necessary; this
+ ;; ensures that the Rounded signal will be raised.
+ (if (<= (len (ref ans '_int)) prec)
+ (let ((expdiff (+ prec 1 (- (len (ref ans '_int))))))
+ (set! ans (_dec_from_triple (ref ans '_sign)
+ (+ (ref ans '_int)
+ (* "0" expdiff))
+ (- (ref ans '_exp) expdiff)))))
+
+ ;; create a copy of the current context, with cleared flags/traps
+ (let ((newcontext (cx-copy context)))
+ (cx-clear_flags newcontext))
+
+ (for ((exception : _signals)) ()
+ (pylist-set! (cx-traps newcontext) exception 0)
+ (values))
+
+ ;; round in the new context
+ (set! ans ((ref ans '_fix) newcontext))
+
+ ;; raise Inexact, and if necessary, Underflow
+ ((cx-error newcontext) Inexact)
+ (if (bool (pylist-ref (cx-flags newcontext) Subnormal))
+ ((cx-error newcontext) Underflow))
+
+ ;; propagate signals to the original context; _fix could
+ ;; have raised any of Overflow, Underflow, Subnormal,
+ ;; Inexact, Rounded, Clamped. Overflow needs the correct
+ ;; arguments. Note that the order of the exceptions is
+ ;; important here.
+ (if (bool (pylist-ref (cx-flags newcontext) Overflow))
+ ((cx-error newcontext)
+ Overflow "above Emax" (ref ans '_sign)))
+
+ (for ((exception : (list Underflow Subnormal
+ Inexact Rounded Clamped))) ()
+ (if (bool (pylist-ref (cx-flags newcontext) exception))
+ ((cx-error newcontext) exception)
+ (values))))
+
+ (set! ans ((ref ans '_fix) context)))
- def normalize(self, context=None):
- """Normalize- strip trailing 0s, change anything equal to 0 to 0e0"""
+ ans))))
- if context is None:
- context = getcontext()
+ (define __rpow__
+ (lam (self other (= context None))
+ "Swaps self/other and returns __pow__."
+ (twix
+ ((norm-op other) it it)
+ ((ref 'other '__pow__) self #:context context))))
- if self._is_special:
- ans = self._check_nans(context=context)
- if ans:
- return ans
-
- dup = self._fix(context)
- if dup._isinfinity():
- return dup
-
- if not dup:
- return _dec_from_triple(dup._sign, '0', 0)
- exp_max = [context.Emax, context.Etop()][context.clamp]
- end = len(dup._int)
- exp = dup._exp
- while dup._int[end-1] == '0' and exp < exp_max:
- exp += 1
- end -= 1
- return _dec_from_triple(dup._sign, dup._int[:end], exp)
-
- def quantize(self, exp, rounding=None, context=None):
- """Quantize self so its exponent is the same as that of exp.
+ (define normalize
+ (lam (self (= context None))
+ "Normalize- strip trailing 0s, change anything equal to 0 to 0e0"
+
+ (twix
+ (let (get-context context))
+ (un-special self context)
+ (let ((dup ((ref self _fix) context))))
+ (((dup '_isinfinity)) it dup)
+ ((not (bool dup))
+ (_dec_from_triple (reg dup '_sign) "0" 0))
+
+ (let* ((_int (ref dup '_int))
+ (exp_max (let ((i (cx-clamp context)))
+ (if (= i 0)
+ (cx-emax context)
+ (cx-etop context))))
+ (let lp ((end (len _int)) (exp (ref dup '_exp)))
+ (if (and (equal? (pylist-ref _int (- end-1))
+ "0")
+ (< exp exp_max))
+ (lp (- end 1) (+ exp 1))
+ (_dec_from_triple
+ (ref dup '_sign)
+ (pylist-slice _int None end None)
+ exp))))))))
+
+ (define quantize
+ (lam (self exp (= rounding None) (= context None))
+ "Quantize self so its exponent is the same as that of exp.
Similar to self._rescale(exp._exp) but with error checking.
- """
- exp = _convert_other(exp, raiseit=True)
-
- if context is None:
- context = getcontext()
- if rounding is None:
- rounding = context.rounding
-
- if self._is_special or exp._is_special:
- ans = self._check_nans(exp, context)
- if ans:
- return ans
-
- if exp._isinfinity() or self._isinfinity():
- if exp._isinfinity() and self._isinfinity():
- return Decimal(self) # if both are inf, it is OK
- return context._raise_error(InvalidOperation,
- 'quantize with one INF')
-
- # exp._exp should be between Etiny and Emax
- if not (context.Etiny() <= exp._exp <= context.Emax):
- return context._raise_error(InvalidOperation,
- 'target exponent out of bounds in quantize')
-
- if not self:
- ans = _dec_from_triple(self._sign, '0', exp._exp)
- return ans._fix(context)
-
- self_adjusted = self.adjusted()
- if self_adjusted > context.Emax:
- return context._raise_error(InvalidOperation,
- 'exponent of quantize result too large for current context')
- if self_adjusted - exp._exp + 1 > context.prec:
- return context._raise_error(InvalidOperation,
- 'quantize result has too many digits for current context')
+ "
+ (twix
+ (let* ((exp (_convert_other exp #:raiseit #t))
+ (context (if (eq? context None) (getcontext) context))
+ (rounding (if (eq? rounding None)
+ (cx-rounding context)
+ rounding))))
+
+ ((if (or ((self '_is_special)) ((exp '_is_special)))
+ (let ((ans ((ref self '_check_nans) exp context)))
+ (cond
+ ((bool ans)
+ ans)
+ ((or ((ref exp '_isinfinity)) ((ref self '_isinfinity)))
+ (if (and ((ref exp '_isinfinity))
+ ((ref self '_isinfinity)))
+ (Decimal self)) ; if both are inf, it is OK
+ ((cx-error context) InvalidOperation "quantize with one INF"))
+ (else
+ #f)))
+ #f) it it)
+
+ ;; exp._exp should be between Etiny and Emax
+ (let ((_eexp (ref exp '_exp))
+ (Emax (cx-emax context))))
+
+ ((not (and (<= (cx-etiny context) eexp) (<= eexp Emax))) it
+ ((cx-error context) InvalidOperation,
+ "target exponent out of bounds in quantize"))
- ans = self._rescale(exp._exp, rounding)
- if ans.adjusted() > context.Emax:
- return context._raise_error(InvalidOperation,
- 'exponent of quantize result too large for current context')
- if len(ans._int) > context.prec:
- return context._raise_error(InvalidOperation,
- 'quantize result has too many digits for current context')
+ ((not (bool self)) it
+ (let ((ans (_dec_from_triple (ref self '_sign) "0" _eexp)))
+ ((ref ans '_fix) context)))
- # raise appropriate flags
- if ans and ans.adjusted() < context.Emin:
- context._raise_error(Subnormal)
- if ans._exp > self._exp:
- if ans != self:
- context._raise_error(Inexact)
- context._raise_error(Rounded)
+ (let ((self_adjusted ((ref self 'adjusted)))
+ (prec (cx-prec context))))
- # call to fix takes care of any necessary folddown, and
- # signals Clamped if necessary
- ans = ans._fix(context)
- return ans
+ ((> self_adjusted (cx-emax context)) it
+ ((cx-error context) InvalidOperation,
+ "exponent of quantize result too large for current context"))
+
+ ((> (+ self_adjusted (- _eexp) 1) prec) it
+ ((cx-error context) InvalidOperation,
+ "quantize result has too many digits for current context"))
- def same_quantum(self, other, context=None):
- """Return True if self and other have the same exponent; otherwise
+ (let ((ans ((ref self '_rescale) _eexp rounding))))
+
+ (if (> ((ref ans 'adjusted)) Emax) it
+ ((cx-error context) InvalidOperation,
+ "exponent of quantize result too large for current context"))
+
+ ((> (len (ref ans '_int)) prec)
+ ((cx-error context) InvalidOperation,
+ "quantize result has too many digits for current context"))
+
+
+ (begin
+ ;; raise appropriate flags
+ (if (and (bool ans) (< ((ref ans 'adjusted)) (cx-emin context)))
+ ((cx-error context) Subnormal))
+
+ (when (> (reg ans '_exp) (ref self '_exp))
+ (if (not (equal? ans self))
+ ((cx-error context) Inexact))
+ ((cx-error context) Rounded))
+
+ ;; call to fix takes care of any necessary folddown, and
+ ;; signals Clamped if necessary
+ ((ref ans '_fix) context)))))
+
+ (define same_quantum
+ (lam (self other (= context None))
+ "Return True if self and other have the same exponent; otherwise
return False.
If either operand is a special value, the following rules are used:
* return True if both operands are infinities
* return True if both operands are NaNs
* otherwise, return False.
- """
- other = _convert_other(other, raiseit=True)
- if self._is_special or other._is_special:
- return (self.is_nan() and other.is_nan() or
- self.is_infinite() and other.is_infinite())
- return self._exp == other._exp
+ "
+ (let ((other (_convert_other other #raiseit #t)))
+ (if (or (ref self '_is_special) (ref other '_is_special))
+ (or (and ((ref self 'is_nan)) ((ref other 'is_nan)))
+ (and ((ref self 'is_infinite)) ((ref other 'is_infinite))))))
+
+ (= (ref self '_exp) (ref other '_exp))))
- def _rescale(self, exp, rounding):
- """Rescale self so that the exponent is exp, either by padding with zeros
+ (define _rescale
+ (lam (self exp rounding)
+ "Rescale self so that the exponent is exp, either by padding with zeros
or by truncating digits, using the given rounding mode.
Specials are returned without change. This operation is
@@ -2682,32 +2759,45 @@
exp = exp to scale to (an integer)
rounding = rounding mode
- """
- if self._is_special:
- return Decimal(self)
- if not self:
- return _dec_from_triple(self._sign, '0', exp)
-
- if self._exp >= exp:
- # pad answer with zeros if necessary
- return _dec_from_triple(self._sign,
- self._int + '0'*(self._exp - exp), exp)
-
- # too many digits; round and lose data. If self.adjusted() <
- # exp-1, replace self by 10**(exp-1) before rounding
- digits = len(self._int) + self._exp - exp
- if digits < 0:
- self = _dec_from_triple(self._sign, '1', exp-1)
- digits = 0
- this_function = self._pick_rounding_function[rounding]
- changed = this_function(self, digits)
- coeff = self._int[:digits] or '0'
- if changed == 1:
- coeff = str(int(coeff)+1)
- return _dec_from_triple(self._sign, coeff, exp)
-
- def _round(self, places, rounding):
- """Round a nonzero, nonspecial Decimal to a fixed number of
+ "
+
+ (cond
+ ((ref self '_is_special) it
+ (Decimal self))
+
+ ((not (bool self)) it
+ (_dec_from_triple (ref self '_sign) "0" exp))
+
+ (let ((_exp (ref self '_exp))
+ (_sign (ref self '_sign))
+ (_int (ref self '_int))))
+
+ ((>= _exp exp)
+ ;; pad answer with zeros if necessary
+ (_dec_from_triple _sign (+ _int (* "0" (- _exp exp))) exp))
+
+
+ ;; too many digits; round and lose data. If self.adjusted() <
+ ;; exp-1, replace self by 10**(exp-1) before rounding
+ (let ((digits (+ (len _int) _exp (- exp))))
+ (if (< digits 0)
+ (set! self (_dec_from_triple _sign "1" (- exp 1)))
+ (set! digits 0))
+
+ (let* ((this_function (pylist-ref (ref self '_pick_rounding_function)
+ rounding))
+ (changed (this_function self digits))
+ (coeff (or (bool
+ (pylist-slice _int None digits None))
+ "0")))
+ (if (= changed 1)
+ (set! coeff (str (+ (int coeff) 1))))
+
+ (_dec_from_triple _sign coeff exp))))))
+
+ (define _round
+ (lam (self places rounding)
+ "Round a nonzero, nonspecial Decimal to a fixed number of
significant figures, using the given rounding mode.
Infinities, NaNs and zeros are returned unaltered.
@@ -2715,22 +2805,28 @@
This operation is quiet: it raises no flags, and uses no
information from the context.
- """
- if places <= 0:
- raise ValueError("argument should be at least 1 in _round")
- if self._is_special or not self:
- return Decimal(self)
- ans = self._rescale(self.adjusted()+1-places, rounding)
- # it can happen that the rescale alters the adjusted exponent;
- # for example when rounding 99.97 to 3 significant figures.
- # When this happens we end up with an extra 0 at the end of
- # the number; a second rescale fixes this.
- if ans.adjusted() != self.adjusted():
- ans = ans._rescale(ans.adjusted()+1-places, rounding)
- return ans
-
- def to_integral_exact(self, rounding=None, context=None):
- """Rounds to a nearby integer.
+ "
+ (cond
+ ((<= places 0)
+ (raise (ValueError "argument should be at least 1 in _round")))
+ ((or (ref self '_is_special) (not (bool self)))
+ (Decimal self))
+ (else
+ (let ((ans ((ref self '_rescale) (+ ((self 'adjusted)) 1 (- places))
+ rounding)))
+ ;; it can happen that the rescale alters the adjusted exponent;
+ ;; for example when rounding 99.97 to 3 significant figures.
+ ;; When this happens we end up with an extra 0 at the end of
+ ;; the number; a second rescale fixes this.
+ (if (not (= ((ref ans 'adjusted)) ((ref self 'adjusted))))
+ (set! ans ((ref ans '_rescale) (+ ((ans 'adjusted)) 1
+ (- places))
+ rounding)))
+ ans)))))
+
+ (define to_integral_exact
+ (lam (self (= rounding None) (= context None))
+ "Rounds to a nearby integer.
If no rounding mode is specified, take the rounding mode from
the context. This method raises the Rounded and Inexact flags
@@ -2738,133 +2834,155 @@
See also: to_integral_value, which does exactly the same as
this method except that it doesn't raise Inexact or Rounded.
- """
- if self._is_special:
- ans = self._check_nans(context=context)
- if ans:
- return ans
- return Decimal(self)
- if self._exp >= 0:
- return Decimal(self)
- if not self:
- return _dec_from_triple(self._sign, '0', 0)
- if context is None:
- context = getcontext()
- if rounding is None:
- rounding = context.rounding
- ans = self._rescale(0, rounding)
- if ans != self:
- context._raise_error(Inexact)
- context._raise_error(Rounded)
- return ans
+ "
+ (cond
+ ((ref self '_is_special)
+ (let ((ans ((ref self '_check_nans) #:context context)))
+ (if (bool ans)
+ ans
+ (Decimal self))))
+
+ ((>= (ref self '_exp) 0)
+ (Decimal self))
+ ((not (boool self))
+ (_dec_from_triple (ref self '_sign) "0" 0))
+ (else
+ (let* ((context (if (eq? context None) (getcontext) context))
+ (rounding (if (eq? rounding None)
+ (cx-rounding context)
+ rounding))
+ (ans ((ref self '_rescale) 0 rounding)))
+
+ (if (not (equal? ans self))
+ ((cx-error context) Inexact))
+
+ ((cx-error context) Rounded)
- def to_integral_value(self, rounding=None, context=None):
- """Rounds to the nearest integer, without raising inexact, rounded."""
- if context is None:
- context = getcontext()
- if rounding is None:
- rounding = context.rounding
- if self._is_special:
- ans = self._check_nans(context=context)
- if ans:
- return ans
- return Decimal(self)
- if self._exp >= 0:
- return Decimal(self)
- else:
- return self._rescale(0, rounding)
+ ans)))))
- # the method name changed, but we provide also the old one, for compatibility
- to_integral = to_integral_value
+ (define to_integral_value
+ (lam (self (= rounding None) (= context None))
+ "Rounds to the nearest integer, without raising inexact, rounded."
+ (let* ((context (if (eq? context None) (getcontext) context))
+ (rounding (if (eq? rounding None)
+ (cx-rounding context)
+ rounding)))
+
+ (cond
+ ((ref self '_is_special)
+ (let ((ans ((ref self '_check_nans) #:context context)))
+ (if (bool ans)
+ ans
+ (Decimal self))))
+
+ ((>= (ref self '_exp) 0)
+ (Decimal self))
+
+ (else
+ ((ref self '_rescale) 0 rounding))))))
- def sqrt(self, context=None):
- """Return the square root of self."""
- if context is None:
- context = getcontext()
+ ;; the method name changed, but we provide also the old one,
+ ;; for compatibility
+ (define to_integral to_integral_value)
- if self._is_special:
- ans = self._check_nans(context=context)
- if ans:
- return ans
+ (define sqrt
+ (lam (self (= context None))
+ "Return the square root of self."
+ (twix
+ (let (get-context context))
- if self._isinfinity() and self._sign == 0:
- return Decimal(self)
+ ((if (ref self '_is_special)
+ (let ((ans ((ref self '_check_nans) #:context context)))
+ (if (bool ans)
+ ans
+ (if (and ((self '_isinfinity)) (= (ref self '_sign) 0))
+ (Decimal self)
+ #f)))
- if not self:
- # exponent = self._exp // 2. sqrt(-0) = -0
- ans = _dec_from_triple(self._sign, '0', self._exp // 2)
- return ans._fix(context)
+ #f) it it)
- if self._sign == 1:
- return context._raise_error(InvalidOperation, 'sqrt(-x), x > 0')
-
- # At this point self represents a positive number. Let p be
- # the desired precision and express self in the form c*100**e
- # with c a positive real number and e an integer, c and e
- # being chosen so that 100**(p-1) <= c < 100**p. Then the
- # (exact) square root of self is sqrt(c)*10**e, and 10**(p-1)
- # <= sqrt(c) < 10**p, so the closest representable Decimal at
- # precision p is n*10**e where n = round_half_even(sqrt(c)),
- # the closest integer to sqrt(c) with the even integer chosen
- # in the case of a tie.
- #
- # To ensure correct rounding in all cases, we use the
- # following trick: we compute the square root to an extra
- # place (precision p+1 instead of precision p), rounding down.
- # Then, if the result is inexact and its last digit is 0 or 5,
- # we increase the last digit to 1 or 6 respectively; if it's
- # exact we leave the last digit alone. Now the final round to
- # p places (or fewer in the case of underflow) will round
- # correctly and raise the appropriate flags.
-
- # use an extra digit of precision
- prec = context.prec+1
-
- # write argument in the form c*100**e where e = self._exp//2
- # is the 'ideal' exponent, to be used if the square root is
- # exactly representable. l is the number of 'digits' of c in
- # base 100, so that 100**(l-1) <= c < 100**l.
- op = _WorkRep(self)
- e = op.exp >> 1
- if op.exp & 1:
- c = op.int * 10
- l = (len(self._int) >> 1) + 1
- else:
- c = op.int
- l = len(self._int)+1 >> 1
-
- # rescale so that c has exactly prec base 100 'digits'
- shift = prec-l
- if shift >= 0:
- c *= 100**shift
- exact = True
- else:
- c, remainder = divmod(c, 100**-shift)
- exact = not remainder
- e -= shift
+ ((not (bool self)) it
+ ;; exponent = self._exp // 2. sqrt(-0) = -0
+ (let ((ans (_dec_from_triple (ref self '_sign) "0"
+ (quotient (ref self '_exp) 2))))
+ ((ref ans '_fix) context)))
+
+ ((= (ref self '_sign) 1)
+ ((cx-error context) InvalidOperation "sqrt(-x), x > 0"))
+
+ ;; At this point self represents a positive number. Let p be
+ ;; the desired precision and express self in the form c*100**e
+ ;; with c a positive real number and e an integer, c and e
+ ;; being chosen so that 100**(p-1) <= c < 100**p. Then the
+ ;; (exact) square root of self is sqrt(c)*10**e, and 10**(p-1)
+ ;; <= sqrt(c) < 10**p, so the closest representable Decimal at
+ ;; precision p is n*10**e where n = round_half_even(sqrt(c)),
+ ;; the closest integer to sqrt(c) with the even integer chosen
+ ;; in the case of a tie.
+ ;;
+ ;; To ensure correct rounding in all cases, we use the
+ ;; following trick: we compute the square root to an extra
+ ;; place (precision p+1 instead of precision p), rounding down.
+ ;; Then, if the result is inexact and its last digit is 0 or 5,
+ ;; we increase the last digit to 1 or 6 respectively; if it's
+ ;; exact we leave the last digit alone. Now the final round to
+ ;; p places (or fewer in the case of underflow) will round
+ ;; correctly and raise the appropriate flags.
+
+ ;; use an extra digit of precision
+ (let* ((prec (+ (cx-prec context) 1))
+ (op (_WorkRep self))
+ (e (ash (ref op 'exp) -1))
+ (c #f)
+ (l #f)
+ (shift #f)
+ (exact #f))
+
+ ;; write argument in the form c*100**e where e = self._exp//2
+ ;; is the 'ideal' exponent, to be used if the square root is
+ ;; exactly representable. l is the number of 'digits' of c in
+ ;; base 100, so that 100**(l-1) <= c < 100**l.
+ (if (= (logand (ref op 'exp) 1) 1)
+ (begin
+ (set! c (* (ref op 'int) 10))
+ (set! l (+ (ash (len (ref self '_int)) -1) 1)))
+ (begin
+ (set! c (ref op 'int))
+ (set! l (ash (+ (len (ref self '_int)) 1) -1))))
- # find n = floor(sqrt(c)) using Newton's method
- n = 10**prec
- while True:
- q = c//n
- if n <= q:
- break
- else:
- n = n + q >> 1
- exact = exact and n*n == c
-
- if exact:
- # result is exact; rescale to use ideal exponent e
- if shift >= 0:
- # assert n % 10**shift == 0
- n //= 10**shift
- else:
- n *= 10**-shift
- e += shift
- else:
- # result is not exact; fix last digit as described above
- if n % 5 == 0:
- n += 1
+ ;; rescale so that c has exactly prec base 100 'digits'
+ (set! shift (- prec l))
+ (if (>= shift 0)
+ (begin
+ (set! c (* c (expt 100 shift)))
+ (set! exact #t))
+ (let ((x (expt 100 (- shift))))
+ (set! c (quotient c x))
+ (let ((remainder (modulo c x)))
+ (set! exact (= remainder 0)))))
+
+ (set! e (- e shift))
+
+ ;; find n = floor(sqrt(c)) using Newton's method
+ (let ((n (let lp ((n (expt 10 prec)))
+ (let ((q (quotient c n)))
+ (if (<= n q)
+ n
+ (lp (ash (+ n q) -1)))))))
+
+ (set! exact (and exact (= (* n n) c)))
+
+ (if exact
+ ;; result is exact; rescale to use ideal exponent e
+ (begin
+ (if (>= shift 0)
+ ;; assert n % 10**shift == 0
+ (set! n (quotient n (expt 10 shift)))
+ (set! n (* n (expt 10 (- shift)))))
+ (set! e (+ e shift)))
+ ;; result is not exact; fix last digit as described above
+ (if (= (modulo n 5) 0)
+ (set! n (+ n 1))))
ans = _dec_from_triple(0, str(n), e)
@@ -2918,7 +3036,7 @@
return ans._fix(context)
- def min(self, other, context=None):
+ def min(self, other, context=None):
"""Returns the smaller value.
Like min(self, other) except if one is not a number, returns
@@ -5731,3 +5849,746 @@ def _normalize(op1, op2, prec = 0):
return op1, op2
+##### Integer arithmetic functions used by ln, log10, exp and __pow__ #####
+
+_nbits = int.bit_length
+
+def _decimal_lshift_exact(n, e):
+ """ Given integers n and e, return n * 10**e if it's an integer, else None.
+
+ The computation is designed to avoid computing large powers of 10
+ unnecessarily.
+
+ >>> _decimal_lshift_exact(3, 4)
+ 30000
+ >>> _decimal_lshift_exact(300, -999999999) # returns None
+
+ """
+ if n == 0:
+ return 0
+ elif e >= 0:
+ return n * 10**e
+ else:
+ # val_n = largest power of 10 dividing n.
+ str_n = str(abs(n))
+ val_n = len(str_n) - len(str_n.rstrip('0'))
+ return None if val_n < -e else n // 10**-e
+
+def _sqrt_nearest(n, a):
+ """Closest integer to the square root of the positive integer n. a is
+ an initial approximation to the square root. Any positive integer
+ will do for a, but the closer a is to the square root of n the
+ faster convergence will be.
+
+ """
+ if n <= 0 or a <= 0:
+ raise ValueError("Both arguments to _sqrt_nearest should be positive.")
+
+ b=0
+ while a != b:
+ b, a = a, a--n//a>>1
+ return a
+
+def _rshift_nearest(x, shift):
+ """Given an integer x and a nonnegative integer shift, return closest
+ integer to x / 2**shift; use round-to-even in case of a tie.
+
+ """
+ b, q = 1 << shift, x >> shift
+ return q + (2*(x & (b-1)) + (q&1) > b)
+
+def _div_nearest(a, b):
+ """Closest integer to a/b, a and b positive integers; rounds to even
+ in the case of a tie.
+
+ """
+ q, r = divmod(a, b)
+ return q + (2*r + (q&1) > b)
+
+def _ilog(x, M, L = 8):
+ """Integer approximation to M*log(x/M), with absolute error boundable
+ in terms only of x/M.
+
+ Given positive integers x and M, return an integer approximation to
+ M * log(x/M). For L = 8 and 0.1 <= x/M <= 10 the difference
+ between the approximation and the exact result is at most 22. For
+ L = 8 and 1.0 <= x/M <= 10.0 the difference is at most 15. In
+ both cases these are upper bounds on the error; it will usually be
+ much smaller."""
+
+ # The basic algorithm is the following: let log1p be the function
+ # log1p(x) = log(1+x). Then log(x/M) = log1p((x-M)/M). We use
+ # the reduction
+ #
+ # log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
+ #
+ # repeatedly until the argument to log1p is small (< 2**-L in
+ # absolute value). For small y we can use the Taylor series
+ # expansion
+ #
+ # log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
+ #
+ # truncating at T such that y**T is small enough. The whole
+ # computation is carried out in a form of fixed-point arithmetic,
+ # with a real number z being represented by an integer
+ # approximation to z*M. To avoid loss of precision, the y below
+ # is actually an integer approximation to 2**R*y*M, where R is the
+ # number of reductions performed so far.
+
+ y = x-M
+ # argument reduction; R = number of reductions performed
+ R = 0
+ while (R <= L and abs(y) << L-R >= M or
+ R > L and abs(y) >> R-L >= M):
+ y = _div_nearest((M*y) << 1,
+ M + _sqrt_nearest(M*(M+_rshift_nearest(y, R)), M))
+ R += 1
+
+ # Taylor series with T terms
+ T = -int(-10*len(str(M))//(3*L))
+ yshift = _rshift_nearest(y, R)
+ w = _div_nearest(M, T)
+ for k in range(T-1, 0, -1):
+ w = _div_nearest(M, k) - _div_nearest(yshift*w, M)
+
+ return _div_nearest(w*y, M)
+
+def _dlog10(c, e, p):
+ """Given integers c, e and p with c > 0, p >= 0, compute an integer
+ approximation to 10**p * log10(c*10**e), with an absolute error of
+ at most 1. Assumes that c*10**e is not exactly 1."""
+
+ # increase precision by 2; compensate for this by dividing
+ # final result by 100
+ p += 2
+
+ # write c*10**e as d*10**f with either:
+ # f >= 0 and 1 <= d <= 10, or
+ # f <= 0 and 0.1 <= d <= 1.
+ # Thus for c*10**e close to 1, f = 0
+ l = len(str(c))
+ f = e+l - (e+l >= 1)
+
+ if p > 0:
+ M = 10**p
+ k = e+p-f
+ if k >= 0:
+ c *= 10**k
+ else:
+ c = _div_nearest(c, 10**-k)
+
+ log_d = _ilog(c, M) # error < 5 + 22 = 27
+ log_10 = _log10_digits(p) # error < 1
+ log_d = _div_nearest(log_d*M, log_10)
+ log_tenpower = f*M # exact
+ else:
+ log_d = 0 # error < 2.31
+ log_tenpower = _div_nearest(f, 10**-p) # error < 0.5
+
+ return _div_nearest(log_tenpower+log_d, 100)
+
+def _dlog(c, e, p):
+ """Given integers c, e and p with c > 0, compute an integer
+ approximation to 10**p * log(c*10**e), with an absolute error of
+ at most 1. Assumes that c*10**e is not exactly 1."""
+
+ # Increase precision by 2. The precision increase is compensated
+ # for at the end with a division by 100.
+ p += 2
+
+ # rewrite c*10**e as d*10**f with either f >= 0 and 1 <= d <= 10,
+ # or f <= 0 and 0.1 <= d <= 1. Then we can compute 10**p * log(c*10**e)
+ # as 10**p * log(d) + 10**p*f * log(10).
+ l = len(str(c))
+ f = e+l - (e+l >= 1)
+
+ # compute approximation to 10**p*log(d), with error < 27
+ if p > 0:
+ k = e+p-f
+ if k >= 0:
+ c *= 10**k
+ else:
+ c = _div_nearest(c, 10**-k) # error of <= 0.5 in c
+
+ # _ilog magnifies existing error in c by a factor of at most 10
+ log_d = _ilog(c, 10**p) # error < 5 + 22 = 27
+ else:
+ # p <= 0: just approximate the whole thing by 0; error < 2.31
+ log_d = 0
+
+ # compute approximation to f*10**p*log(10), with error < 11.
+ if f:
+ extra = len(str(abs(f)))-1
+ if p + extra >= 0:
+ # error in f * _log10_digits(p+extra) < |f| * 1 = |f|
+ # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
+ f_log_ten = _div_nearest(f*_log10_digits(p+extra), 10**extra)
+ else:
+ f_log_ten = 0
+ else:
+ f_log_ten = 0
+
+ # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1
+ return _div_nearest(f_log_ten + log_d, 100)
+
+class _Log10Memoize(object):
+ """Class to compute, store, and allow retrieval of, digits of the
+ constant log(10) = 2.302585.... This constant is needed by
+ Decimal.ln, Decimal.log10, Decimal.exp and Decimal.__pow__."""
+ def __init__(self):
+ self.digits = "23025850929940456840179914546843642076011014886"
+
+ def getdigits(self, p):
+ """Given an integer p >= 0, return floor(10**p)*log(10).
+
+ For example, self.getdigits(3) returns 2302.
+ """
+ # digits are stored as a string, for quick conversion to
+ # integer in the case that we've already computed enough
+ # digits; the stored digits should always be correct
+ # (truncated, not rounded to nearest).
+ if p < 0:
+ raise ValueError("p should be nonnegative")
+
+ if p >= len(self.digits):
+ # compute p+3, p+6, p+9, ... digits; continue until at
+ # least one of the extra digits is nonzero
+ extra = 3
+ while True:
+ # compute p+extra digits, correct to within 1ulp
+ M = 10**(p+extra+2)
+ digits = str(_div_nearest(_ilog(10*M, M), 100))
+ if digits[-extra:] != '0'*extra:
+ break
+ extra += 3
+ # keep all reliable digits so far; remove trailing zeros
+ # and next nonzero digit
+ self.digits = digits.rstrip('0')[:-1]
+ return int(self.digits[:p+1])
+
+_log10_digits = _Log10Memoize().getdigits
+
+def _iexp(x, M, L=8):
+ """Given integers x and M, M > 0, such that x/M is small in absolute
+ value, compute an integer approximation to M*exp(x/M). For 0 <=
+ x/M <= 2.4, the absolute error in the result is bounded by 60 (and
+ is usually much smaller)."""
+
+ # Algorithm: to compute exp(z) for a real number z, first divide z
+ # by a suitable power R of 2 so that |z/2**R| < 2**-L. Then
+ # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
+ # series
+ #
+ # expm1(x) = x + x**2/2! + x**3/3! + ...
+ #
+ # Now use the identity
+ #
+ # expm1(2x) = expm1(x)*(expm1(x)+2)
+ #
+ # R times to compute the sequence expm1(z/2**R),
+ # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).
+
+ # Find R such that x/2**R/M <= 2**-L
+ R = _nbits((x<<L)//M)
+
+ # Taylor series. (2**L)**T > M
+ T = -int(-10*len(str(M))//(3*L))
+ y = _div_nearest(x, T)
+ Mshift = M<<R
+ for i in range(T-1, 0, -1):
+ y = _div_nearest(x*(Mshift + y), Mshift * i)
+
+ # Expansion
+ for k in range(R-1, -1, -1):
+ Mshift = M<<(k+2)
+ y = _div_nearest(y*(y+Mshift), Mshift)
+
+ return M+y
+
+def _dexp(c, e, p):
+ """Compute an approximation to exp(c*10**e), with p decimal places of
+ precision.
+
+ Returns integers d, f such that:
+
+ 10**(p-1) <= d <= 10**p, and
+ (d-1)*10**f < exp(c*10**e) < (d+1)*10**f
+
+ In other words, d*10**f is an approximation to exp(c*10**e) with p
+ digits of precision, and with an error in d of at most 1. This is
+ almost, but not quite, the same as the error being < 1ulp: when d
+ = 10**(p-1) the error could be up to 10 ulp."""
+
+ # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision
+ p += 2
+
+ # compute log(10) with extra precision = adjusted exponent of c*10**e
+ extra = max(0, e + len(str(c)) - 1)
+ q = p + extra
+
+ # compute quotient c*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q),
+ # rounding down
+ shift = e+q
+ if shift >= 0:
+ cshift = c*10**shift
+ else:
+ cshift = c//10**-shift
+ quot, rem = divmod(cshift, _log10_digits(q))
+
+ # reduce remainder back to original precision
+ rem = _div_nearest(rem, 10**extra)
+
+ # error in result of _iexp < 120; error after division < 0.62
+ return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
+
+def _dpower(xc, xe, yc, ye, p):
+ """Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and
+ y = yc*10**ye, compute x**y. Returns a pair of integers (c, e) such that:
+
+ 10**(p-1) <= c <= 10**p, and
+ (c-1)*10**e < x**y < (c+1)*10**e
+
+ in other words, c*10**e is an approximation to x**y with p digits
+ of precision, and with an error in c of at most 1. (This is
+ almost, but not quite, the same as the error being < 1ulp: when c
+ == 10**(p-1) we can only guarantee error < 10ulp.)
+
+ We assume that: x is positive and not equal to 1, and y is nonzero.
+ """
+
+ # Find b such that 10**(b-1) <= |y| <= 10**b
+ b = len(str(abs(yc))) + ye
+
+ # log(x) = lxc*10**(-p-b-1), to p+b+1 places after the decimal point
+ lxc = _dlog(xc, xe, p+b+1)
+
+ # compute product y*log(x) = yc*lxc*10**(-p-b-1+ye) = pc*10**(-p-1)
+ shift = ye-b
+ if shift >= 0:
+ pc = lxc*yc*10**shift
+ else:
+ pc = _div_nearest(lxc*yc, 10**-shift)
+
+ if pc == 0:
+ # we prefer a result that isn't exactly 1; this makes it
+ # easier to compute a correctly rounded result in __pow__
+ if ((len(str(xc)) + xe >= 1) == (yc > 0)): # if x**y > 1:
+ coeff, exp = 10**(p-1)+1, 1-p
+ else:
+ coeff, exp = 10**p-1, -p
+ else:
+ coeff, exp = _dexp(pc, -(p+1), p+1)
+ coeff = _div_nearest(coeff, 10)
+ exp += 1
+
+ return coeff, exp
+
+def _log10_lb(c, correction = {
+ '1': 100, '2': 70, '3': 53, '4': 40, '5': 31,
+ '6': 23, '7': 16, '8': 10, '9': 5}):
+ """Compute a lower bound for 100*log10(c) for a positive integer c."""
+ if c <= 0:
+ raise ValueError("The argument to _log10_lb should be nonnegative.")
+ str_c = str(c)
+ return 100*len(str_c) - correction[str_c[0]]
+
+##### Helper Functions ####################################################
+
+def _convert_other(other, raiseit=False, allow_float=False):
+ """Convert other to Decimal.
+
+ Verifies that it's ok to use in an implicit construction.
+ If allow_float is true, allow conversion from float; this
+ is used in the comparison methods (__eq__ and friends).
+
+ """
+ if isinstance(other, Decimal):
+ return other
+ if isinstance(other, int):
+ return Decimal(other)
+ if allow_float and isinstance(other, float):
+ return Decimal.from_float(other)
+
+ if raiseit:
+ raise TypeError("Unable to convert %s to Decimal" % other)
+ return NotImplemented
+
+def _convert_for_comparison(self, other, equality_op=False):
+ """Given a Decimal instance self and a Python object other, return
+ a pair (s, o) of Decimal instances such that "s op o" is
+ equivalent to "self op other" for any of the 6 comparison
+ operators "op".
+
+ """
+ if isinstance(other, Decimal):
+ return self, other
+
+ # Comparison with a Rational instance (also includes integers):
+ # self op n/d <=> self*d op n (for n and d integers, d positive).
+ # A NaN or infinity can be left unchanged without affecting the
+ # comparison result.
+ if isinstance(other, _numbers.Rational):
+ if not self._is_special:
+ self = _dec_from_triple(self._sign,
+ str(int(self._int) * other.denominator),
+ self._exp)
+ return self, Decimal(other.numerator)
+
+ # Comparisons with float and complex types. == and != comparisons
+ # with complex numbers should succeed, returning either True or False
+ # as appropriate. Other comparisons return NotImplemented.
+ if equality_op and isinstance(other, _numbers.Complex) and other.imag == 0:
+ other = other.real
+ if isinstance(other, float):
+ context = getcontext()
+ if equality_op:
+ context.flags[FloatOperation] = 1
+ else:
+ context._raise_error(FloatOperation,
+ "strict semantics for mixing floats and Decimals are enabled")
+ return self, Decimal.from_float(other)
+ return NotImplemented, NotImplemented
+
+
+##### Setup Specific Contexts ############################################
+
+# The default context prototype used by Context()
+# Is mutable, so that new contexts can have different default values
+
+DefaultContext = Context(
+ prec=28, rounding=ROUND_HALF_EVEN,
+ traps=[DivisionByZero, Overflow, InvalidOperation],
+ flags=[],
+ Emax=999999,
+ Emin=-999999,
+ capitals=1,
+ clamp=0
+)
+
+# Pre-made alternate contexts offered by the specification
+# Don't change these; the user should be able to select these
+# contexts and be able to reproduce results from other implementations
+# of the spec.
+
+BasicContext = Context(
+ prec=9, rounding=ROUND_HALF_UP,
+ traps=[DivisionByZero, Overflow, InvalidOperation, Clamped, Underflow],
+ flags=[],
+)
+
+ExtendedContext = Context(
+ prec=9, rounding=ROUND_HALF_EVEN,
+ traps=[],
+ flags=[],
+)
+
+
+##### crud for parsing strings #############################################
+#
+# Regular expression used for parsing numeric strings. Additional
+# comments:
+#
+# 1. Uncomment the two '\s*' lines to allow leading and/or trailing
+# whitespace. But note that the specification disallows whitespace in
+# a numeric string.
+#
+# 2. For finite numbers (not infinities and NaNs) the body of the
+# number between the optional sign and the optional exponent must have
+# at least one decimal digit, possibly after the decimal point. The
+# lookahead expression '(?=\d|\.\d)' checks this.
+
+import re
+_parser = re.compile(r""" # A numeric string consists of:
+# \s*
+ (?P<sign>[-+])? # an optional sign, followed by either...
+ (
+ (?=\d|\.\d) # ...a number (with at least one digit)
+ (?P<int>\d*) # having a (possibly empty) integer part
+ (\.(?P<frac>\d*))? # followed by an optional fractional part
+ (E(?P<exp>[-+]?\d+))? # followed by an optional exponent, or...
+ |
+ Inf(inity)? # ...an infinity, or...
+ |
+ (?P<signal>s)? # ...an (optionally signaling)
+ NaN # NaN
+ (?P<diag>\d*) # with (possibly empty) diagnostic info.
+ )
+# \s*
+ \Z
+""", re.VERBOSE | re.IGNORECASE).match
+
+_all_zeros = re.compile('0*$').match
+_exact_half = re.compile('50*$').match
+
+##### PEP3101 support functions ##############################################
+# The functions in this section have little to do with the Decimal
+# class, and could potentially be reused or adapted for other pure
+# Python numeric classes that want to implement __format__
+#
+# A format specifier for Decimal looks like:
+#
+# [[fill]align][sign][#][0][minimumwidth][,][.precision][type]
+
+_parse_format_specifier_regex = re.compile(r"""\A
+(?:
+ (?P<fill>.)?
+ (?P<align>[<>=^])
+)?
+(?P<sign>[-+ ])?
+(?P<alt>\#)?
+(?P<zeropad>0)?
+(?P<minimumwidth>(?!0)\d+)?
+(?P<thousands_sep>,)?
+(?:\.(?P<precision>0|(?!0)\d+))?
+(?P<type>[eEfFgGn%])?
+\Z
+""", re.VERBOSE|re.DOTALL)
+
+del re
+
+# The locale module is only needed for the 'n' format specifier. The
+# rest of the PEP 3101 code functions quite happily without it, so we
+# don't care too much if locale isn't present.
+try:
+ import locale as _locale
+except ImportError:
+ pass
+
+def _parse_format_specifier(format_spec, _localeconv=None):
+ """Parse and validate a format specifier.
+
+ Turns a standard numeric format specifier into a dict, with the
+ following entries:
+
+ fill: fill character to pad field to minimum width
+ align: alignment type, either '<', '>', '=' or '^'
+ sign: either '+', '-' or ' '
+ minimumwidth: nonnegative integer giving minimum width
+ zeropad: boolean, indicating whether to pad with zeros
+ thousands_sep: string to use as thousands separator, or ''
+ grouping: grouping for thousands separators, in format
+ used by localeconv
+ decimal_point: string to use for decimal point
+ precision: nonnegative integer giving precision, or None
+ type: one of the characters 'eEfFgG%', or None
+
+ """
+ m = _parse_format_specifier_regex.match(format_spec)
+ if m is None:
+ raise ValueError("Invalid format specifier: " + format_spec)
+
+ # get the dictionary
+ format_dict = m.groupdict()
+
+ # zeropad; defaults for fill and alignment. If zero padding
+ # is requested, the fill and align fields should be absent.
+ fill = format_dict['fill']
+ align = format_dict['align']
+ format_dict['zeropad'] = (format_dict['zeropad'] is not None)
+ if format_dict['zeropad']:
+ if fill is not None:
+ raise ValueError("Fill character conflicts with '0'"
+ " in format specifier: " + format_spec)
+ if align is not None:
+ raise ValueError("Alignment conflicts with '0' in "
+ "format specifier: " + format_spec)
+ format_dict['fill'] = fill or ' '
+ # PEP 3101 originally specified that the default alignment should
+ # be left; it was later agreed that right-aligned makes more sense
+ # for numeric types. See http://bugs.python.org/issue6857.
+ format_dict['align'] = align or '>'
+
+ # default sign handling: '-' for negative, '' for positive
+ if format_dict['sign'] is None:
+ format_dict['sign'] = '-'
+
+ # minimumwidth defaults to 0; precision remains None if not given
+ format_dict['minimumwidth'] = int(format_dict['minimumwidth'] or '0')
+ if format_dict['precision'] is not None:
+ format_dict['precision'] = int(format_dict['precision'])
+
+ # if format type is 'g' or 'G' then a precision of 0 makes little
+ # sense; convert it to 1. Same if format type is unspecified.
+ if format_dict['precision'] == 0:
+ if format_dict['type'] is None or format_dict['type'] in 'gGn':
+ format_dict['precision'] = 1
+
+ # determine thousands separator, grouping, and decimal separator, and
+ # add appropriate entries to format_dict
+ if format_dict['type'] == 'n':
+ # apart from separators, 'n' behaves just like 'g'
+ format_dict['type'] = 'g'
+ if _localeconv is None:
+ _localeconv = _locale.localeconv()
+ if format_dict['thousands_sep'] is not None:
+ raise ValueError("Explicit thousands separator conflicts with "
+ "'n' type in format specifier: " + format_spec)
+ format_dict['thousands_sep'] = _localeconv['thousands_sep']
+ format_dict['grouping'] = _localeconv['grouping']
+ format_dict['decimal_point'] = _localeconv['decimal_point']
+ else:
+ if format_dict['thousands_sep'] is None:
+ format_dict['thousands_sep'] = ''
+ format_dict['grouping'] = [3, 0]
+ format_dict['decimal_point'] = '.'
+
+ return format_dict
+
+def _format_align(sign, body, spec):
+ """Given an unpadded, non-aligned numeric string 'body' and sign
+ string 'sign', add padding and alignment conforming to the given
+ format specifier dictionary 'spec' (as produced by
+ parse_format_specifier).
+
+ """
+ # how much extra space do we have to play with?
+ minimumwidth = spec['minimumwidth']
+ fill = spec['fill']
+ padding = fill*(minimumwidth - len(sign) - len(body))
+
+ align = spec['align']
+ if align == '<':
+ result = sign + body + padding
+ elif align == '>':
+ result = padding + sign + body
+ elif align == '=':
+ result = sign + padding + body
+ elif align == '^':
+ half = len(padding)//2
+ result = padding[:half] + sign + body + padding[half:]
+ else:
+ raise ValueError('Unrecognised alignment field')
+
+ return result
+
+def _group_lengths(grouping):
+ """Convert a localeconv-style grouping into a (possibly infinite)
+ iterable of integers representing group lengths.
+
+ """
+ # The result from localeconv()['grouping'], and the input to this
+ # function, should be a list of integers in one of the
+ # following three forms:
+ #
+ # (1) an empty list, or
+ # (2) nonempty list of positive integers + [0]
+ # (3) list of positive integers + [locale.CHAR_MAX], or
+
+ from itertools import chain, repeat
+ if not grouping:
+ return []
+ elif grouping[-1] == 0 and len(grouping) >= 2:
+ return chain(grouping[:-1], repeat(grouping[-2]))
+ elif grouping[-1] == _locale.CHAR_MAX:
+ return grouping[:-1]
+ else:
+ raise ValueError('unrecognised format for grouping')
+
+def _insert_thousands_sep(digits, spec, min_width=1):
+ """Insert thousands separators into a digit string.
+
+ spec is a dictionary whose keys should include 'thousands_sep' and
+ 'grouping'; typically it's the result of parsing the format
+ specifier using _parse_format_specifier.
+
+ The min_width keyword argument gives the minimum length of the
+ result, which will be padded on the left with zeros if necessary.
+
+ If necessary, the zero padding adds an extra '0' on the left to
+ avoid a leading thousands separator. For example, inserting
+ commas every three digits in '123456', with min_width=8, gives
+ '0,123,456', even though that has length 9.
+
+ """
+
+ sep = spec['thousands_sep']
+ grouping = spec['grouping']
+
+ groups = []
+ for l in _group_lengths(grouping):
+ if l <= 0:
+ raise ValueError("group length should be positive")
+ # max(..., 1) forces at least 1 digit to the left of a separator
+ l = min(max(len(digits), min_width, 1), l)
+ groups.append('0'*(l - len(digits)) + digits[-l:])
+ digits = digits[:-l]
+ min_width -= l
+ if not digits and min_width <= 0:
+ break
+ min_width -= len(sep)
+ else:
+ l = max(len(digits), min_width, 1)
+ groups.append('0'*(l - len(digits)) + digits[-l:])
+ return sep.join(reversed(groups))
+
+def _format_sign(is_negative, spec):
+ """Determine sign character."""
+
+ if is_negative:
+ return '-'
+ elif spec['sign'] in ' +':
+ return spec['sign']
+ else:
+ return ''
+
+def _format_number(is_negative, intpart, fracpart, exp, spec):
+ """Format a number, given the following data:
+
+ is_negative: true if the number is negative, else false
+ intpart: string of digits that must appear before the decimal point
+ fracpart: string of digits that must come after the point
+ exp: exponent, as an integer
+ spec: dictionary resulting from parsing the format specifier
+
+ This function uses the information in spec to:
+ insert separators (decimal separator and thousands separators)
+ format the sign
+ format the exponent
+ add trailing '%' for the '%' type
+ zero-pad if necessary
+ fill and align if necessary
+ """
+
+ sign = _format_sign(is_negative, spec)
+
+ if fracpart or spec['alt']:
+ fracpart = spec['decimal_point'] + fracpart
+
+ if exp != 0 or spec['type'] in 'eE':
+ echar = {'E': 'E', 'e': 'e', 'G': 'E', 'g': 'e'}[spec['type']]
+ fracpart += "{0}{1:+}".format(echar, exp)
+ if spec['type'] == '%':
+ fracpart += '%'
+
+ if spec['zeropad']:
+ min_width = spec['minimumwidth'] - len(fracpart) - len(sign)
+ else:
+ min_width = 0
+ intpart = _insert_thousands_sep(intpart, spec, min_width)
+
+ return _format_align(sign, intpart+fracpart, spec)
+
+
+##### Useful Constants (internal use only) ################################
+
+# Reusable defaults
+_Infinity = Decimal('Inf')
+_NegativeInfinity = Decimal('-Inf')
+_NaN = Decimal('NaN')
+_Zero = Decimal(0)
+_One = Decimal(1)
+_NegativeOne = Decimal(-1)
+
+# _SignedInfinity[sign] is infinity w/ that sign
+_SignedInfinity = (_Infinity, _NegativeInfinity)
+
+# Constants related to the hash implementation; hash(x) is based
+# on the reduction of x modulo _PyHASH_MODULUS
+_PyHASH_MODULUS = sys.hash_info.modulus
+# hash values to use for positive and negative infinities, and nans
+_PyHASH_INF = sys.hash_info.inf
+_PyHASH_NAN = sys.hash_info.nan
+
+# _PyHASH_10INV is the inverse of 10 modulo the prime _PyHASH_MODULUS
+_PyHASH_10INV = pow(10, _PyHASH_MODULUS - 2, _PyHASH_MODULUS)
+del sys
generated by cgit v1.2.3 (git 2.46.0) at 2024-10-19 02:23:29 +0000