1 (define-module (language python module decimal
)
2 #:use-module
((language python module collections
) #:select
(namedtuple))
3 #:use-module
((language python module itertools
) #:select
(chain repeat
))
4 #:use-module
((language python module sys
) #:select
(maxsize hash_info
))
5 #:use-module
(language python module
)
6 #:use-module
((language python module python
) #:select
7 (isinstance str float int tuple classmethod pow property
8 complex range reversed
(format . str-format
)))
9 #:use-module
(language python list
)
10 #:use-module
(language python number
)
11 #:use-module
(language python string
)
12 #:use-module
(language python for
)
13 #:use-module
(language python try
)
14 #:use-module
(language python hash
)
15 #:use-module
(language python dict
)
16 #:use-module
(language python def
)
17 #:use-module
(language python exceptions
)
18 #:use-module
(language python bool
)
19 #:use-module
(language python module
)
20 #:use-module
(oop pf-objects
)
21 #:use-module
(language python module re
)
22 #:use-module
(ice-9 control
)
23 #:use-module
((ice-9 match
) #:select
((match . ice
:match
)))
25 ( ;; Two major classes
28 ;; Named tuple representation
32 DefaultContext BasicContext ExtendedContext
35 DecimalException Clamped InvalidOperation DivisionByZero
36 Inexact Rounded Subnormal Overflow Underflow
39 ;; Exceptional conditions that trigger InvalidOperation
40 DivisionImpossible InvalidContext ConversionSyntax DivisionUndefined
42 ;; Constants for use in setting up contexts
43 ROUND_DOWN ROUND_HALF_UP ROUND_HALF_EVEN ROUND_CEILING
44 ROUND_FLOOR ROUND_UP ROUND_HALF_DOWN ROUND_05UP
46 ;; Functions for manipulating contexts
47 setcontext getcontext localcontext
49 ;; Limits for the C version for compatibility
50 MAX_PREC MAX_EMAX MIN_EMIN MIN_ETINY
))
52 (define-syntax-rule (aif it p . l
) (let ((it p
)) (if it . l
)))
55 This is the copyright information of the file ported over to scheme
56 # Copyright
(c) 2004 Python Software Foundation.
57 # All rights reserved.
59 # Written by Eric Price
<eprice at tjhsst.edu
>
60 # and Facundo Batista
<facundo at taniquetil.com.ar
>
61 # and Raymond Hettinger
<python at rcn.com
>
62 # and Aahz
<aahz at pobox.com
>
65 # This module should be kept in sync with the latest updates of the
66 # IBM specification as it evolves. Those updates will be treated
67 # as bug fixes
(deviation from the spec is a compatibility
, usability
68 # bug
) and will be backported. At this point the spec is stabilizing
69 # and the updates are becoming fewer
, smaller
, and less significant.
72 (define guile
:modulo
(@ (guile) modulo
))
74 (define __name__
"decimal")
75 (define __xname__ __name__
)
76 (define __version__
"1.70")
77 ;; Highest version of the spec this complies with
78 ;; See http://speleotrove.com/decimal/
80 (define DecimalTuple
(namedtuple "DecimalTuple" "sign,digits,exponent"))
83 (define ROUND_DOWN
'ROUND_DOWN
)
84 (define ROUND_HALF_UP
'ROUND_HALF_UP
)
85 (define ROUND_HALF_EVEN
'ROUND_HALF_EVEN
)
86 (define ROUND_CEILING
'ROUND_CEILING
)
87 (define ROUND_FLOOR
'ROUND_FLOOR
)
88 (define ROUND_UP
'ROUND_UP
)
89 (define ROUND_HALF_DOWN
'ROUND_HALF_DOWN
)
90 (define ROUND_05UP
'ROUND_05UP
)
92 ;; Compatibility with the C version
93 (define MAX_PREC
425000000)
94 (define MAX_EMAX
425000000)
95 (define MIN_EMIN -
425000000)
97 (if (= maxsize
(- (ash 1 63) 1))
99 (set! MAX_PREC
999999999999999999)
100 (set! MAX_EMAX
999999999999999999)
101 (set! MIN_EMIN -
999999999999999999)))
103 (define MIN_ETINY
(- MIN_EMIN
(- MAX_PREC
1)))
106 (define-inlinable (cx-prec x
) (rawref x
'prec
))
107 (define-inlinable (cx-emax x
) (rawref x
'Emax
))
108 (define-inlinable (cx-emin x
) (rawref x
'Emin
))
109 (define-inlinable (cx-etiny x
) ((ref x
'Etiny
)))
110 (define-inlinable (cx-etop x
) ((ref x
'Etop
)))
111 (define-inlinable (cx-copy x
) ((ref x
'copy
)))
112 (define-inlinable (cx-clear_flags x
) ((ref x
'clear_flags
)))
113 (define-inlinable (cx-raise x
) (ref x
'_raise_error
))
114 (define-inlinable (cx-error x
) (ref x
'_raise_error
))
115 (define-inlinable (cx-capitals x
) (rawref x
'capitals
))
116 (define-inlinable (cx-cap x
) (rawref x
'capitals
))
117 (define-inlinable (cx-rounding x
) (rawref x
'rounding
))
118 (define-inlinable (cx-clamp x
) (rawref x
'clamp
))
119 (define-inlinable (cx-traps x
) (rawref x
'traps
))
120 (define-inlinable (cx-flags x
) (rawref x
'flags
))
124 (define-python-class DecimalException
(ArithmeticError)
125 "Base exception class.
127 Used exceptions derive from this.
128 If an exception derives from another exception besides this (such as
129 Underflow (Inexact, Rounded, Subnormal) that indicates that it is only
130 called if the others are present. This isn't actually used for
133 handle -- Called when context._raise_error is called and the
134 trap_enabler is not set. First argument is self, second is the
135 context. More arguments can be given, those being after
136 the explanation in _raise_error (For example,
137 context._raise_error(NewError, '(-x)!', self._sign) would
138 call NewError().handle(context, self._sign).)
140 To define a new exception, it should be sufficient to have it derive
141 from DecimalException.
145 (lambda (self context . args
)
149 (define-python-class Clamped
(DecimalException)
150 "Exponent of a 0 changed to fit bounds.
152 This occurs and signals clamped if the exponent of a result has been
153 altered in order to fit the constraints of a specific concrete
154 representation. This may occur when the exponent of a zero result would
155 be outside the bounds of a representation, or when a large normal
156 number would have an encoded exponent that cannot be represented. In
157 this latter case, the exponent is reduced to fit and the corresponding
158 number of zero digits are appended to the coefficient ('fold-down').
161 (define-python-class InvalidOperation
(DecimalException)
162 "An invalid operation was performed.
164 Various bad things cause this:
166 Something creates a signaling NaN
172 x._rescale( non-integer )
177 An operand is invalid
179 The result of the operation after these is a quiet positive NaN,
180 except when the cause is a signaling NaN, in which case the result is
181 also a quiet NaN, but with the original sign, and an optional
182 diagnostic information.
185 (lambda (self context . args
)
187 (let ((ans (_dec_from_triple
188 (ref (car args
) '_sign
)
189 (ref (car args
) '_int
)
191 ((ref ans
'_fix_nan
) context
))
194 (define-python-class ConversionSyntax
(InvalidOperation)
195 "Trying to convert badly formed string.
197 This occurs and signals invalid-operation if a string is being
198 converted to a number and it does not conform to the numeric string
199 syntax. The result is [0,qNaN].
204 (define-python-class DivisionByZero
(DecimalException ZeroDivisionError
)
207 This occurs and signals division-by-zero if division of a finite number
208 by zero was attempted (during a divide-integer or divide operation, or a
209 power operation with negative right-hand operand), and the dividend was
212 The result of the operation is [sign,inf], where sign is the exclusive
213 or of the signs of the operands for divide, or is 1 for an odd power of
218 (lambda (self context sign . args
)
219 (pylist-ref _SignedInfinity sign
))))
221 (define-python-class DivisionImpossible
(InvalidOperation)
222 "Cannot perform the division adequately.
224 This occurs and signals invalid-operation if the integer result of a
225 divide-integer or remainder operation had too many digits (would be
226 longer than precision). The result is [0,qNaN].
232 (define-python-class DivisionUndefined
(InvalidOperation ZeroDivisionError
)
233 "Undefined result of division.
235 This occurs and signals invalid-operation if division by zero was
236 attempted (during a divide-integer, divide, or remainder operation), and
237 the dividend is also zero. The result is [0,qNaN].
243 (define-python-class Inexact
(DecimalException)
244 "Had to round, losing information.
246 This occurs and signals inexact whenever the result of an operation is
247 not exact (that is, it needed to be rounded and any discarded digits
248 were non-zero), or if an overflow or underflow condition occurs. The
249 result in all cases is unchanged.
251 The inexact signal may be tested (or trapped) to determine if a given
252 operation (or sequence of operations) was inexact.
255 (define-python-class InvalidContext
(InvalidOperation)
256 "Invalid context. Unknown rounding, for example.
258 This occurs and signals invalid-operation if an invalid context was
259 detected during an operation. This can occur if contexts are not checked
260 on creation and either the precision exceeds the capability of the
261 underlying concrete representation or an unknown or unsupported rounding
262 was specified. These aspects of the context need only be checked when
263 the values are required to be used. The result is [0,qNaN].
269 (define-python-class Rounded
(DecimalException)
270 "Number got rounded (not necessarily changed during rounding).
272 This occurs and signals rounded whenever the result of an operation is
273 rounded (that is, some zero or non-zero digits were discarded from the
274 coefficient), or if an overflow or underflow condition occurs. The
275 result in all cases is unchanged.
277 The rounded signal may be tested (or trapped) to determine if a given
278 operation (or sequence of operations) caused a loss of precision.
281 (define-python-class Subnormal
(DecimalException)
282 "Exponent < Emin before rounding.
284 This occurs and signals subnormal whenever the result of a conversion or
285 operation is subnormal (that is, its adjusted exponent is less than
286 Emin, before any rounding). The result in all cases is unchanged.
288 The subnormal signal may be tested (or trapped) to determine if a given
289 or operation (or sequence of operations) yielded a subnormal result.
292 (define-python-class Overflow
(Inexact Rounded
)
295 This occurs and signals overflow if the adjusted exponent of a result
296 (from a conversion or from an operation that is not an attempt to divide
297 by zero), after rounding, would be greater than the largest value that
298 can be handled by the implementation (the value Emax).
300 The result depends on the rounding mode:
302 For round-half-up and round-half-even (and for round-half-down and
303 round-up, if implemented), the result of the operation is [sign,inf],
304 where sign is the sign of the intermediate result. For round-down, the
305 result is the largest finite number that can be represented in the
306 current precision, with the sign of the intermediate result. For
307 round-ceiling, the result is the same as for round-down if the sign of
308 the intermediate result is 1, or is [0,inf] otherwise. For round-floor,
309 the result is the same as for round-down if the sign of the intermediate
310 result is 0, or is [1,inf] otherwise. In all cases, Inexact and Rounded
315 (let ((l (list ROUND_HALF_UP ROUND_HALF_EVEN
316 ROUND_HALF_DOWN ROUND_UP
)))
317 (lambda (self context sign . args
)
319 (if (memq (ref context
'rounding
) l
)
320 (ret (pylist-ref _SignedInfinity sign
)))
323 (if (eq?
(ref context
'rounding
) ROUND_CEILING
)
324 (ret (pylist-ref _SignedInfinity sign
))
325 (ret (_dec_from_triple
327 (* "9" (cx-prec context
))
328 (+ (- (cx-emax context
) (cx-prec context
)) 1)))))
331 (if (eq?
(ref context
'rounding
) ROUND_FLOOR
)
332 (ret (pylist-ref _SignedInfinity sign
))
333 (ret (_dec_from_triple
335 (* "9" (cx-prec context
))
336 (+ (- (cx-emax context
) (cx-prec context
)) 1))))))))))
339 (define-python-class Underflow
(Inexact Rounded Subnormal
)
340 "Numerical underflow with result rounded to 0.
342 This occurs and signals underflow if a result is inexact and the
343 adjusted exponent of the result would be smaller (more negative) than
344 the smallest value that can be handled by the implementation (the value
345 Emin). That is, the result is both inexact and subnormal.
347 The result after an underflow will be a subnormal number rounded, if
348 necessary, so that its exponent is not less than Etiny. This may result
349 in 0 with the sign of the intermediate result and an exponent of Etiny.
351 In all cases, Inexact, Rounded, and Subnormal will also be raised.
354 (define-python-class FloatOperation
(DecimalException TypeError
)
355 "Enable stricter semantics for mixing floats and Decimals.
357 If the signal is not trapped (default), mixing floats and Decimals is
358 permitted in the Decimal() constructor, context.create_decimal() and
359 all comparison operators. Both conversion and comparisons are exact.
360 Any occurrence of a mixed operation is silently recorded by setting
361 FloatOperation in the context flags. Explicit conversions with
362 Decimal.from_float() or context.create_decimal_from_float() do not
365 Otherwise (the signal is trapped), only equality comparisons and explicit
366 conversions are silent. All other mixed operations raise FloatOperation.
369 ;; List of public traps and flags
371 (list Clamped DivisionByZero Inexact Overflow Rounded
372 Underflow InvalidOperation Subnormal FloatOperation
))
374 ;; Map conditions (per the spec) to signals
375 (define _condition_map
376 `((,ConversionSyntax .
,InvalidOperation
)
377 (,DivisionImpossible .
,InvalidOperation
)
378 (,DivisionUndefined .
,InvalidOperation
)
379 (,InvalidContext .
,InvalidOperation
)))
381 ;; Valid rounding modes
382 (define _rounding_modes
383 (list ROUND_DOWN ROUND_HALF_UP ROUND_HALF_EVEN ROUND_CEILING
384 ROUND_FLOOR ROUND_UP ROUND_HALF_DOWN ROUND_05UP
))
386 ;; ##### Context Functions ##################################################
388 ;; The getcontext() and setcontext() function manage access to a thread-local
390 (define *context
* (make-fluid #f
))
392 (fluid-ref *context
*))
393 (define (setcontext context
)
394 (fluid-set! *context
* context
))
396 ;; ##### Decimal class #######################################################
398 ;; Do not subclass Decimal from numbers.Real and do not register it as such
399 ;; (because Decimals are not interoperable with floats). See the notes in
400 ;; numbers.py for more detail.
402 (define _dec_from_triple
403 (lam (sign coefficient exponent
(= special
#f
))
404 "Create a decimal instance directly, without any validation,
405 normalization (e.g. removal of leading zeros) or argument
408 This function is for *internal use only*.
410 (Decimal sign coefficient exponent special
)))
412 (define (get-parsed-sign m
)
413 (if (equal?
((ref m
'group
) "sign") "-")
417 (define (get-parsed-int m
)
418 ((ref m
'group
) "int"))
420 (define (get-parsed-frac m
)
421 ((ref m
'group
) "frac"))
423 (define (get-parsed-exp m
)
424 ((ref m
'group
) "exp"))
426 (define (get-parsed-diag m
)
427 ((ref m
'group
) "diag"))
429 (define (get-parsed-sig m
)
430 ((ref m
'group
) "signal"))
432 (def (_mk self __init__
(= value
"0") (= context None
))
433 "Create a decimal point instance.
435 >>> Decimal('3.14') # string input
437 >>> Decimal((0, (3, 1, 4), -2)) # tuple (sign, digit_tuple, exponent)
439 >>> Decimal(314) # int
441 >>> Decimal(Decimal(314)) # another decimal instance
443 >>> Decimal(' 3.14 \\n') # leading and trailing whitespace okay
447 ;; Note that the coefficient, self._int, is actually stored as
448 ;; a string rather than as a tuple of digits. This speeds up
449 ;; the "digits to integer" and "integer to digits" conversions
450 ;; that are used in almost every arithmetic operation on
451 ;; Decimals. This is an internal detail: the as_tuple function
452 ;; and the Decimal constructor still deal with tuples of
456 ;; REs insist on real strings, so we can too.
458 ((isinstance value str
)
459 (let ((m (_parser (scm-str str
))))
461 (let ((context (if (eq? context None
)
466 (+ "Invalid literal for Decimal: " value
))))
468 (let ((sign (get-parsed-sign m
))
469 (intpart (get-parsed-int m
))
470 (fracpart (get-parsed-frac m
))
471 (exp (get-parsed-exp m
))
472 (diag (get-parsed-diag m
))
473 (signal (get-parsed-sig m
)))
475 (set self
'sign sign
)
477 (if (not (eq? intpart None
))
479 (set self
'_int
(str (int (+ intpart fracpart
))))
480 (set self
'_exp
(- exp
(len fracpart
)))
481 (set self
'_is_special
#f
))
483 (if (not (eq? diag None
))
487 (py-lstrip (str (int (if (bool diag
)
493 (set self
'_exp
"n")))
497 (set self
'_exp
"F")))
498 (set self
'_is_special
#t
))))))
501 ((isinstance value int
)
506 (set self
'_int
(str (abs value
)))
507 (set self
'_is_special
#f
))
509 ;; From another decimal
510 ((isinstance value Decimal
)
511 (set self
'_exp
(ref value
'_exp
))
512 (set self
'_sign
(ref value
'_sign
))
513 (set self
'_int
(ref value
'_int
))
514 (set self
'_is_special
(ref value
'_is_special
)))
516 ;; From an internal working value
517 ((isinstance value _WorkRep
)
518 (set self
'_exp
(int (ref value
'_exp
)))
519 (set self
'_sign
(ref value
'_sign
))
520 (set self
'_int
(str (ref value
'int
)))
521 (set self
'_is_special
#f
))
523 ;; tuple/list conversion (possibly from as_tuple())
524 ((isinstance value
(list list tuple
))
525 (if (not (= (len value
) 3))
527 (+ "Invalid tuple size in creation of Decimal "
528 "from list or tuple. The list or tuple "
529 "should have exactly three elements."))))
530 ;; # process sign. The isinstance test rejects floats
531 (let ((v0 (pylist-ref value
0))
532 (v1 (pylist-ref value
1))
533 (v2 (pylist-ref value
2)))
534 (if (not (and (isinstance v0 int
)
535 (or (= v0
0) (= v0
1))))
537 (+ "Invalid sign. The first value in the tuple "
538 "should be an integer; either 0 for a "
539 "positive number or 1 for a negative number."))))
545 (set self
'is_special
#t
))
546 (let ((digits (py-list)))
547 ;; process and validate the digits in value[1]
548 (for ((digit : v1
)) ()
549 (if (and (isinstance digit int
)
552 ;; skip leading zeros
553 (if (or (bool digits
) (> digit
0))
554 (pylist-append! digits digit
))
556 (+ "The second value in the tuple must "
557 "be composed of integers in the range "
561 ((or (eq? v2
'n
) (eq? v2
'N
))
563 ;; NaN: digits form the diagnostic
564 (set self
'_int
(py-join "" (map str digits
)))
566 (set self
'_is_special
#t
)))
568 ;; finite number: digits give the coefficient
569 (set self
'_int
(py-join "" (map str digits
)))
571 (set self
'_is_special
#f
))
574 (+ "The third value in the tuple must "
575 "be an integer, or one of the "
576 "strings 'F', 'n', 'N'.")))))))))
578 ((isinstance value float
)
579 (let ((context (if (eq? context None
)
584 (+ "strict semantics for mixing floats and Decimals are "
587 (__init__ self
((ref Decimal
'from_float
) value
))))
591 (format #f
"Cannot convert ~a to Decimal" value
))))))
593 (define-inlinable (divmod x y
)
594 (values (quotient x y
) (modulo x y
)))
597 (syntax-rules (when let let
* if
)
599 ((_ (let () a ...
) . l
)
600 (begin a ...
(twix . l
)))
602 ((_ (let ((a aa
) ...
) b ...
) . l
)
603 (let ((a aa
) ...
) b ...
(twix . l
)))
604 ((_ (let (a ...
)) . l
)
606 ((_ (let* (a ...
) b ...
) . l
)
607 (let* (a ...
) b ...
(twix . l
)))
609 (begin (if . u
) (twix . l
)))
611 (begin (when . u
) (twix . l
)))
612 ((_ (a it code ...
) . l
)
613 (aif it a
(begin code ...
) (twix . l
)))))
615 (define-syntax-rule (norm-op self op
)
617 (set! op
((ref self
'_convert_other
) op
))
618 (if (eq? op NotImplemented
)
622 (define-syntax-rule (get-context context code
)
623 (let ((context (if (eq? context None
)
628 (define-syntax-rule (un-special self context
)
629 (if ((ref self
'_is_special
))
630 (let ((ans ((ref self
'_check_nans
) #:context context
)))
636 (define-syntax-rule (bin-special o1 o2 context
)
637 (if (or (ref o1
'_is_special
)
638 (ref o2
'_is_special
))
639 (or (un-special o1 context
) (un-special o2 context
))))
641 (define-syntax-rule (add-special self other context
)
642 (or (bin-special self other context
)
643 (if ((ref self
'_isinfinity
))
644 ;; If both INF, same sign =>
645 ;; same as both, opposite => error.
646 (if (and (not (= (ref self
'_sign
) (ref other
'_sign
)))
647 ((ref other
'_isinfinity
)))
648 ((cx-error context
) InvalidOperation
"-INF + INF")
650 (if ((ref other
'_isinfinity
))
651 (Decimal other
) ; Can't both be infinity here
654 (define-syntax-rule (mul-special self other context resultsign
)
655 (if (or (ref self
'_is_special
) (ref other
'_is_special
))
657 ((bin-special self other context
) it it
)
659 ((if ((ref self
'_isinfinity
))
660 (if (not (bool other
))
661 ((cx-error context
) InvalidOperation
"(+-)INF * 0")
662 (pylist-ref _SignedInfinity resultsign
))
665 (if ((ref other
'_isinfinity
))
666 (if (not (bool self
))
667 ((cx-error context
) InvalidOperation
"(+-)INF * 0")
668 (pylist-ref _SignedInfinity resultsign
))
672 (define-syntax-rule (div-special self other context sign
)
673 (if (or (ref self
'_is_special
) (ref other
'_is_special
))
675 ((bin-special self other context
) it it
)
677 ((and ((ref self
'_isinfinity
)) ((ref other
'_isinfinity
))) it
678 ((cx-error context
) InvalidOperation
"(+-)INF/(+-)INF"))
680 (((ref self
'_isinfinity
)) it
681 (pylist-ref _SignedInfinity sign
))
683 (((ref other
'_isinfinity
)) it
684 ((cx-error context
) Clamped
"Division by infinity")
685 (_dec_from_triple sign
"0" (cx-etiny context
))))))
688 (define-python-class Decimal
(object)
689 "Floating point class for decimal arithmetic."
692 ;; Generally, the value of the Decimal instance is given by
693 ;; (-1)**_sign * _int * 10**_exp
694 ;; Special values are signified by _is_special == True
698 ((self sign coefficient exponent special
)
699 (set self
'_sign sign
)
700 (set self
'_int coefficient
)
701 (set self
'_exp exponent
)
702 (set self
'_is_special special
))
707 (_mk self __init__ a
))
709 (_mk self __init__ a b
))))
714 "Converts a float to a decimal number, exactly.
716 Note that Decimal.from_float(0.1) is not the same as Decimal('0.1').
717 Since 0.1 is not exactly representable in binary floating point, the
718 value is stored as the nearest representable value which is
719 0x1.999999999999ap-4. The exact equivalent of the value in decimal
720 is 0.1000000000000000055511151231257827021181583404541015625.
722 >>> Decimal.from_float(0.1)
723 Decimal('0.1000000000000000055511151231257827021181583404541015625')
724 >>> Decimal.from_float(float('nan'))
726 >>> Decimal.from_float(float('inf'))
728 >>> Decimal.from_float(-float('inf'))
730 >>> Decimal.from_float(-0.0)
735 (if (< x
0) (set! x
(- x
)))
737 (let lp
((l (string->list
(format #f
"~e" x
))) (r1 '()))
740 (let lp
((l l
) (r2 '()))
743 (let* ((n (length r1
))
744 (m (list->string
(append (reverse r1
) (reverse r2
))))
745 (e (+ (- n
1) (string->number
(list->string l
)))))
748 (lp l
(cons x r2
))))))
751 (lp l
(cons x r1
))))))
755 ((isinstance f int
) ; handle integer inputs
757 ((not (isinstance f float
))
758 (raise (TypeError "argument must be int or float.")))
759 ((or (inf? f
) (nan? f
))
763 ((eq? f
(- (inf))) ""))))
765 (let* ((sign (if (>= f
0) 0 1))
769 (res (_dec_from_triple sign m e
)))
770 (if (eq? cls Decimal
)
776 "Returns whether the number is not actually one.
782 (if (ref self
'_is_special
)
783 (let ((exp (ref self
'_exp
)))
792 "Returns whether the number is infinite
794 0 if finite or not a number
798 (if (eq?
(ref self
'_exp
) 'F
)
799 (if (eq?
(ref self
'_sign
) 1)
805 (lam (self (= other None
) (= context None
))
806 "Returns whether the number is not actually one.
808 if self, other are sNaN, signal
809 if self, other are NaN return nan
812 Done before operations.
815 (let ((self_is_nan ((ref self
'_isnan
)))
819 ((ref other
'_isnan
)))))
821 (if (or self_is_nan other_is_nan
)
822 (let ((context (if (eq? context None
)
827 ((cx-error context
) InvalidOperation
"sNaN" self
))
828 ((eq? other_is_nan
2)
829 ((cx-error context
) InvalidOperation
"sNaN" other
))
831 ((ref self
'_fix_nan
) context
))
833 ((ref other
'_fix_nan
) context
))))
837 (define _compare_check_nans
838 (lambda (self other context
)
839 "Version of _check_nans used for the signaling comparisons
840 compare_signal, __le__, __lt__, __ge__, __gt__.
842 Signal InvalidOperation if either self or other is a (quiet
843 or signaling) NaN. Signaling NaNs take precedence over quiet
846 Return 0 if neither operand is a NaN.
849 (let ((context (if (eq? context None
)
853 (if (or (ref self
'_is_special
)
854 (ref other
'_is_special
))
856 (((ref self
'is_snan
))
859 "comparison involving sNaN" self
))
861 (((ref other
'is_snan
))
864 "comparison involving sNaN" other
))
866 (((ref self
'is_qnan
))
869 "comparison involving NaN" self
))
871 (((ref other
'is_qnan
))
874 "comparison involving NaN" other
))
881 "Return True if self is nonzero; otherwise return False.
883 NaNs and infinities are considered nonzero.
885 (or (ref self
'_is_special
) (not (equal?
(ref self
'_int
) "0")))))
889 "Compare the two non-NaN decimal instances self and other.
891 Returns -1 if self < other, 0 if self == other and 1
892 if self > other. This routine is for internal use only."
894 (let ((self_sign (ref self
'_sign
))
895 (other_sign (ref other
'_sign
)))
897 ((or (ref self
'_is_special
) (ref other
'_is_special
))
898 (let ((self_inf ((ref self
'_isinfinity
)))
899 (other_inf ((ref other
'_isinfinity
))))
901 ((eq? self_inf other_inf
) 0)
902 ((< self_inf other_inf
) -
1)
905 ;; check for zeros; Decimal('0') == Decimal('-0')
907 (if (not (bool other
))
909 (let ((s (ref other
'_sign
)))
914 (let ((s (ref self
'_sign
)))
919 ((< other_sign self_sign
)
921 ((< self_sign other_sign
)
925 (let ((self_adjusted ((ref self
'adjusted
)))
926 (other_adjusted ((ref other
'adjusted
)))
927 (self_exp (ref self
'_exp
))
928 (other_exp (ref other
'_exp
)))
930 ((= self_adjusted other_adjusted
)
931 (let ((self_padded (+ (ref self
'_int
)
932 (* "0" (- self_exp other_exp
))))
933 (other_padded (+ (ref other
'_int
)
934 (* "0" (- other_exp self_exp
)))))
936 ((equal? self_padded other_padded
)
938 ((< self_padded other_padded
)
946 ((> self_adjusted other_adjusted
)
955 ;; Note: The Decimal standard doesn't cover rich comparisons for
956 ;; Decimals. In particular, the specification is silent on the
957 ;; subject of what should happen for a comparison involving a NaN.
958 ;; We take the following approach:
960 ;; == comparisons involving a quiet NaN always return False
961 ;; != comparisons involving a quiet NaN always return True
962 ;; == or != comparisons involving a signaling NaN signal
963 ;; InvalidOperation, and return False or True as above if the
964 ;; InvalidOperation is not trapped.
965 ;; <, >, <= and >= comparisons involving a (quiet or signaling)
966 ;; NaN signal InvalidOperation, and return False if the
967 ;; InvalidOperation is not trapped.
969 ;; This behavior is designed to conform as closely as possible to
970 ;; that specified by IEEE 754.
973 (lam (self other
(= context None
))
974 (let* ((so (_convert_for_comparison self other
#:equality_op
#t
))
979 ((eq? other NotImplemented
)
981 ((bool ((ref self
'_check_nans
) other context
))
983 (else (= ((ref self
'_cmp
) other
) 0))))))
987 (lam (self other
(= context None
))
988 (let* ((so (_convert_for_comparison self other
#:equality_op
#t
))
993 ((eq? other NotImplemented
)
995 ((bool ((ref self
'_compare_check_nans
) other context
))
997 (else (< ((ref self
'_cmp
) other
) 0)))))))
999 (define __lt__
(lambda x
(apply (_xlt < ) x
)))
1000 (define __le__
(lambda x
(apply (_xlt <= ) x
)))
1001 (define __gt__
(lambda x
(apply (_xlt > ) x
)))
1002 (define __ge__
(lambda x
(apply (_xlt >= ) x
)))
1005 (lam (self other
(= context None
))
1006 "Compare self to other. Return a decimal value:
1008 a or b is a NaN ==> Decimal('NaN')
1009 a < b ==> Decimal('-1')
1010 a == b ==> Decimal('0')
1011 a > b ==> Decimal('1')
1013 (let ((other (_convert_other other
#:raiseit
#t
)))
1014 ;; Compare(NaN, NaN) = NaN
1015 (if (or (ref self
'_is_special
)
1017 (ref other
'_is_special
)))
1018 (aif it
((ref self
'_check_nans
) other context
)
1020 (Decimal ((ref self
'_cmp
) other
)))))))
1024 "x.__hash__() <==> hash(x)"
1026 ;; In order to make sure that the hash of a Decimal instance
1027 ;; agrees with the hash of a numerically equal integer, float
1028 ;; or Fraction, we follow the rules for numeric hashes outlined
1029 ;; in the documentation. (See library docs, 'Built-in Types').
1031 ((ref self
'_is_special
)
1033 (((ref self
'is_snan
))
1034 (raise (TypeError "Cannot hash a signaling NaN value.")))
1035 (((ref self
'is_snan
))
1036 (hash (nan) pyhash-N
))
1037 ((= 1 (ref self
'_sign
))
1038 (hash (- (inf)) pyhash-N
))
1040 (hash (inf) pyhash-N
))))
1043 (let* ((exp (ref self
'_exp
))
1046 (pow 10 exp pyhash-N
)
1047 (pow _PyHASH_10INV
(- exp
) pyhash-N
)))
1050 (modulo (* (int (ref self
'_int
)) exp_hash
)
1054 (if (>= self
0) hash_
(- hash_
))))
1055 (if (= ans -
1) -
2 ans
))))))
1059 "Represents the number as a triple tuple.
1061 To show the internals exactly as they are.
1063 (DecimalTuple (ref self
'_sign
)
1064 (tuple (map int
(ref self
'_int
)))
1067 (define as_integer_ratio
1069 "Express a finite Decimal instance in the form n / d.
1071 Returns a pair (n, d) of integers. When called on an infinity
1072 or NaN, raises OverflowError or ValueError respectively.
1074 >>> Decimal('3.14').as_integer_ratio()
1076 >>> Decimal('-123e5').as_integer_ratio()
1078 >>> Decimal('0.00').as_integer_ratio()
1081 (if (ref self
'_is_special
)
1082 (if ((ref self
'is_nan
))
1084 "cannot convert NaN to integer ratio"))
1085 (raise (OverflowError
1086 "cannot convert Infinity to integer ratio"))))
1088 (if (not (bool self
))
1090 (let* ((s (ref self
'_sign
))
1091 (n (int (ref self
'_int
)))
1092 (e (ref self
'_exp
))
1096 (/ 1 (expt 10 (- expt
)))))))
1097 (values (numerator x
)
1098 (denominator x
))))))
1102 "Represents the number as an instance of Decimal."
1103 ;# Invariant: eval(repr(d)) == d
1104 (format #f
"Decimal('~a')" (str self
))))
1107 (lam (self (= eng
#f
) (= context None
))
1108 "Return string representation of the number in scientific notation.
1110 Captures all of the information in the underlying representation.
1112 (let* ((sign (if (= (ref self
'_sign
) 0) "" "-"))
1113 (exp (ref self
'_exp
))
1114 (i (ref self
'_int
))
1115 (leftdigits (+ exp
(len i
)))
1122 ((ref self
'_is_special
)
1124 ((eq?
(ref self
'_exp
) 'F
)
1125 (+ sign
"Infinity"))
1126 ((eq?
(ref self
'_exp
) 'n
)
1127 (+ sign
"NaN" (ref self
'_int
)))
1128 (else ; self._exp == 'N'
1129 (+ sign
"sNaN" (ref self
'_int
)))))
1131 ;; dotplace is number of digits of self._int to the left of the
1132 ;; decimal point in the mantissa of the output string (that is,
1133 ;; after adjusting the exponent)
1135 ((and (<= exp
0) (> leftdigits -
6))
1136 ;; no exponent required
1137 (set! dotplace leftdigits
))
1140 ;; usual scientific notation: 1 digit on left of the point
1144 ;; engineering notation, zero
1145 (set! dotplace
(- (modulo (+ leftdigits
1) 3) 1)))
1147 ;; engineering notation, nonzero
1148 (set! dotplace
(- (modulo (+ leftdigits
1) 3) 1))))
1153 (set! fracpart
(+ "." + (* "0" (- dotplace
)) + i
)))
1154 ((>= dotplace
(len i
))
1155 (set! intpart
(+ i
(* "0" (- dotplace
(len i
)))))
1158 (set! intpart
(pylist-slice i None dotplace None
))
1159 (set! fracpart
(+ '.
' (pylist-slice i dotplace None None
)))))
1163 ((= leftdigits dotplace
)
1166 (let ((context (if (eq? context None
)
1170 (+ (pylist-ref '("e" "E") (cx-capitals context
))
1171 (format #f
"~@d" (- leftdigits dotplace
)))))))
1172 (+ sign intpart fracpart exp
))))))
1174 (define to_eng_string
1175 (lam (self (= context None
))
1176 "Convert to a string, using engineering notation if an exponent is needed.
1177 Engineering notation has an exponent which is a multiple of 3. This
1178 can leave up to 3 digits to the left of the decimal place and may
1179 require the addition of either one or two trailing zeros.
1181 ((ref self
'__str__
) #:eng
#t
#:contect context
)))
1184 (lam (self (= context None
))
1185 "Returns a copy with the sign switched.
1187 Rounds, if it has reason.
1190 ((un-special self context
) it it
)
1191 (let* ((context (if (eq? context None
)
1194 (ans (if (and (not (bool self
))
1195 (not (eq?
(cx-rounding context
)
1197 ;; -Decimal('0') is Decimal('0'),
1198 ;; not Decimal('-0'), except
1199 ;; in ROUND_FLOOR rounding mode.
1200 ((ref self
'copy_abs
))
1201 ((ref self
'copy_negate
)))))
1203 ((ref ans
'_fix
) context
)))))
1206 (lam (self (= context None
))
1207 "Returns a copy, unless it is a sNaN.
1209 Rounds the number (if more than precision digits)
1212 ((un-special self context
) it it
)
1214 (let* ((context (if (eq? context None
)
1217 (ans (if (and (not (bool self
))
1218 (not (eq?
(cx-rounding context
)
1220 ;; -Decimal('0') is Decimal('0'),
1221 ;; not Decimal('-0'), except
1222 ;; in ROUND_FLOOR rounding mode.
1223 ((ref self
'copy_abs
))
1226 ((ref ans
'_fix
) context
)))))
1229 (lam (self (= round
#t
) (= context None
))
1230 "Returns the absolute value of self.
1232 If the keyword argument 'round' is false, do not round. The
1233 expression self.__abs__(round=False) is equivalent to
1237 ((not (bool round
)) it
1238 ((ref self
'copy_abs
)))
1240 ((un-special self context
) it it
)
1242 (if (= (ref self
'_sign
) 1)
1243 ((ref self
'__neg__
) #:context context
)
1244 ((ref self
'__pos__
) #:context context
)))))
1248 (lam (self other
(= context None
))
1249 "Returns self + other.
1251 -INF + INF (or the reverse) cause InvalidOperation errors.
1254 ((norm-op self other
) it it
)
1256 (let (get-context context
))
1258 ((add-special self other context
) it it
)
1260 (let* ((negativezero 0)
1261 (self_sign (ref self
'_sign
))
1262 (other_sign (ref other
'_sign
))
1263 (self_exp (ref self
'_sign
))
1264 (other_exp (ref other
'_sign
))
1265 (prec (cx-prec context
))
1266 (exp (min self_exp other_exp
))
1270 (if (and (eq?
(cx-rounding context
) ROUND_FLOOR
)
1271 (not (= self_sign other_sign
)))
1272 ;; If the answer is 0, the sign should be negative,
1274 (set! negativezero
1)))
1276 ((if (and (not (bool self
)) (not (bool other
)))
1278 (set! sign
(min self_sign other_sign
))
1279 (if (= negativezero
1)
1281 (set! ans
(_dec_from_triple sign
"0" exp
))
1282 (set! ans
((ref ans
'_fix
) context
))
1286 ((if (not (bool self
))
1288 (set! exp
(max exp
(- other_exp prec
1)))
1289 (set! ans
((ref other
'_rescale
) exp
1290 (cx-rounding context
)))
1291 (set! ans
((ref ans
'_fix
) context
))
1295 ((if (not (bool other
))
1297 (set! exp
(max exp
(- self_exp prec
1)))
1298 (set! ans
((ref self
'_rescale
) exp
1299 (cx-rounding context
)))
1300 (set! ans
((ref ans
'_fix
) context
))
1305 (let* ((op1 (_WorkRep self
))
1306 (op2 (_WorkRep other
))
1307 (ab (_normalize op1 op2 prec
))
1310 (result (_WorkRep))))
1313 ((not (= (ref op1
'sign
) (ref op2
'sign
)))
1314 ;; Equal and opposite
1317 (set! ans
(_dec_from_triple negativezero
"0" exp
))
1318 (set! ans
((ref ans
'_fix
) context
))
1327 (if (= (ref op1
'sign
) 1)
1328 (let ((t (ref op1
'sign
)))
1329 (set result
'sign
1)
1330 (set op1
'sign
(ref op2
'sign
))
1332 (set result
'sign
0))
1334 ((= (ref op1
'sign
) 1)
1335 (set result
'sign
1)
1339 (set result
'sign
0)
1343 (if (= (ref op2
'sign
) 0)
1344 (set result
'int
(+ (ref op1
'int
) (ref op2
'int
)))
1345 (set result
'int
(- (ref op1
'int
) (ref op2
'int
))))
1347 (set result
'exp
(ref op1
'exp
))
1348 (set! ans
(Decimal result
))
1349 ((ref ans
'_fix
) context
)))))
1351 (define __radd__
(lambda x
(apply __add__ x
)))
1354 (lam (self other
(= context None
))
1355 "Return self - other"
1357 ((norm-op self other
) it it
)
1358 ((bin-special self other context
) it it
)
1359 ((ref self
'__add__
)
1360 ((ref other
'copy_negate
)) #:context context
))))
1363 (lam (self other
(= context None
))
1364 "Return other - self"
1366 ((norm-op self other
) it it
)
1367 ((ref 'other
'__sub__
) self
#:context context
))))
1370 (lam (self other
(= context None
))
1371 "Return self * other.
1373 (+-) INF * 0 (or its reverse) raise InvalidOperation.
1376 ((norm-op self other
) it it
)
1377 (let (get-context context
))
1379 (let ((resultsign (logxor (ref self
'_sign
)
1380 (ref other
'_sign
)))))
1382 ((mul-special self other context resultsign
) it it
)
1384 (let ((resultexp (+ (ref self
'_exp
) (ref other
'_exp
)))))
1386 ;; Special case for multiplying by zero
1387 ((or (not (bool self
)) (not (bool other
))) it
1388 (let ((ans (_dec_from_triple resultsign
"0" resultexp
)))
1389 ((ref ans
'_fix
) context
)))
1391 ;; Special case for multiplying by power of 10
1392 ((equal?
(ref self
'_int
) "1") it
1393 (let ((ans (_dec_from_triple resultsign
(ref other
'_int
) resultexp
)))
1394 ((ref ans
'_fix
) context
)))
1396 ((equal?
(ref other
'_int
) "1") it
1397 (let ((ans (_dec_from_triple resultsign
(ref self
'_int
) resultexp
)))
1398 ((ref ans
'_fix
) context
)))
1400 (let* ((op1 (_WorkRep self
))
1401 (op2 (_WorkRep other
))
1402 (ans (_dec_from_triple resultsign
1403 (str (* (ref op1
'int
) (ref op2
'int
)))
1405 ((ref ans
'_fix
) context
)))))
1407 (define __rmul__ __mul__
)
1410 (lam (self other
(= context None
))
1411 "Return self / other."
1413 ((norm-op self other
) it it
)
1414 (let (get-context context
))
1416 (let ((sign (logxor (ref self
'_sign
)
1417 (ref other
'_sign
)))))
1419 ((div-special self other context sign
) it it
)
1421 ;; Special cases for zeroes
1422 ((if (not (bool other
))
1423 (if (not (bool self
))
1424 ((cx-error context
) DivisionUndefined
"0 / 0")
1425 ((cx-error context
) DivisionByZero
"x / 0" sign
))
1430 (prec (cx-prec context
))
1431 (nself (len (ref self
'_int
)))
1432 (nother (len (ref other
'_int
))))
1433 (if (not (bool self
))
1435 (set! exp
(- (ref self
'_exp
) (ref other
'_exp
)))
1437 ;; OK, so neither = 0, INF or NaN
1438 (let ((shift (+ nother
(- nself
) prec
1))
1439 (op1 (_WorkRep self
))
1440 (op2 (_WorkRep other
)))
1441 (set! exp
(- (ref self
'_exp
) (ref other
'_exp
) shift
))
1445 (divmod (* (ref op1
'int
) (expt 10 shift
))
1447 (divmod (ref op1
'int
)
1448 (* (ref op2
'int
) (expt 10 shift
)))))
1449 (lambda (coeff- remainder
)
1451 (if (not (= remainder
0))
1452 ;; result is not exact adjust to ensure
1454 (if (= (modulo coeff-
5) 0)
1457 (let ((ideal_exp (- (ref self
'_exp
)
1458 (ref other
'_exp
))))
1459 (let lp
((coeff- coeff-
) (exp- exp
))
1460 (if (and (< exp- ideal_exp
)
1461 (= (modulo coeff
10) 0))
1462 (lp (/ coeff
10) (+ exp-
1))
1468 (let ((ans (_dec_from_triple sign
(str coeff
) exp
)))
1469 ((ref ans
'_fix
) context
))))))
1472 (lambda (self other context
)
1473 "Return (self // other, self % other), to context.prec precision.
1475 Assumes that neither self nor other is a NaN, that self is not
1476 infinite and that other is nonzero.
1481 (logxor (ref self
'_sign
)
1482 (ref other
'_sign
)))
1484 (if ((ref other
'_isinfinity
))
1486 (min (ref self
'exp
) (ref other
'_exp
))))
1488 (- ((ref self
'adjusted
)) ((ref other
'adjusted
)))))))
1490 ((or (not (bool self
))
1491 ((ref other
'_isinfinity
))
1493 (list (_dec_from_triple sign
"0" 0)
1494 ((ref self
'_rescale
) ideal_exp
(cx-rounding context
))))
1496 ((if (<= expdiff
(cx-prec context
))
1497 (let ((op1 (_WorkRep self
))
1498 (op2 (_WorkRep other
)))
1499 (if (>= (ref op1
'exp
) (ref op2
'exp
))
1500 (set op1
'int
(* (ref op1
'int
)
1501 (expt 10 (- (ref op1
'exp
)
1503 (set op2
'int
(* (ref op2
'int
)
1504 (expt 10 (- (ref op2
'exp
)
1506 (call-with-values (lambda () (divmod (ref op1
'int
)
1509 (if (< q
(expt 10 (cx-prec context
)))
1510 (list (_dec_from_triple sign
(str q
) 0)
1511 (_dec_from_triple (ref self
'_sign
)
1517 ;; Here the quotient is too large to be representable
1518 (let ((ans ((cx-raise context
) DivisionImpossible
1519 "quotient too large in //, % or divmod")))
1520 (list ans ans
)))))))
1522 (define __rtruediv__
1523 (lam (self other
(= context None
))
1524 "Swaps self/other and returns __truediv__."
1526 ((norm-op self other
) it it
)
1527 ((ref other
'__truediv__
) self
#:context context
))))
1530 (lam (self other
(= context None
))
1532 Return (self // other, self % other)
1536 ((norm-op self other
) it it
)
1538 (let (get-context context
))
1540 ((add-special self other context
) it it
)
1542 (((ref self
'_check_nans
) other context
) it
1546 (logxor (ref self
'_sign
)
1547 (ref other
'_sign
))))))
1549 (((ref self
'_isinfinity
)) it
1550 (if ((ref other
'_isinfinity
))
1551 (let ((ans ((cx-error context
) InvalidOperation
1552 "divmod(INF, INF)")))
1554 (list (list-ref _SignedInfinity sign
)
1555 ((cx-raise context
) InvalidOperation
"INF % x"))))
1557 ((not (bool other
)) it
1558 (if (not (bool self
))
1559 (let ((ans ((cx-error context
) DivisionUndefined
1562 (list ((cx-error context
) DivisionByZero
"x // 0" sign
)
1563 ((cx-error context
) InvalidOperation
"x % 0"))))
1565 (call-with-values (lambda () ((ref self
'_divide
) other context
))
1566 (lambda (quotient remainder
)
1567 (let ((remainder ((ref remainder
'_fix
) context
)))
1568 (list quotient remainder
))))))))
1571 (lam (self other
(= context None
))
1572 "Swaps self/other and returns __divmod__."
1574 ((norm-op self other
) it it
)
1575 ((ref other
'__divmod__
) self
#:context context
))))
1578 (lam (self other
(= context None
))
1583 ((norm-op self other
) it it
)
1585 (let (get-context context
))
1587 ((bin-special self other context
) it it
)
1589 (((ref self
'_isinfinity
)) it
1590 ((cx-error context
) InvalidOperation
"INF % x"))
1592 ((not (bool other
)) it
1594 ((cx-error context
) InvalidOperation
"x % 0")
1595 ((cx-error context
) DivisionUndefined
"0 % 0")))
1597 (let* ((remainder ((ref self
'_divide
) other context
)))
1598 ((ref remainder
'_fix
) context
)))))
1601 (lam (self other
(= context None
))
1602 "Swaps self/other and returns __mod__."
1604 ((norm-op self other
) it it
)
1605 ((ref other
'__mod__
) self
#:context context
))))
1607 (define remainder_near
1608 (lam (self other
(= context None
))
1610 Remainder nearest to 0- abs(remainder-near) <= other/2
1613 ((norm-op self other
) it it
)
1615 (let (get-context context
))
1617 ((bin-special self other context
) it it
)
1619 ;; self == +/-infinity -> InvalidOperation
1620 (((ref self
'_isinfinity
)) it
1621 ((cx-error context
) InvalidOperation
"remainder_near(infinity, x)"))
1623 ;; other == 0 -> either InvalidOperation or DivisionUndefined
1624 ((not (bool other
)) it
1625 (if (not (bool self
))
1626 ((cx-error context
) InvalidOperation
"remainder_near(x, 0)")
1627 ((cx-error context
) DivisionUndefined
"remainder_near(0, 0)")))
1629 ;; other = +/-infinity -> remainder = self
1630 (((ref other
'_isinfinity
)) it
1631 (let ((ans (Decimal self
)))
1632 ((ref ans
'_fix
) context
)))
1634 ;; self = 0 -> remainder = self, with ideal exponent
1635 (let (let ((ideal_exponent (min (ref self
'_exp
) (ref other
'_exp
))))))
1637 ((not (bool self
)) it
1638 (let ((ans (_dec_from_triple (ref self
'_sign
) "0" ideal_exponent
)))
1639 ((ref ans
'_fix
) context
)))
1641 ;; catch most cases of large or small quotient
1643 (- ((ref self
'adjusted
)) ((ref other
'adjusted
)))))))
1645 ((>= expdiff
(+ (cx-prec context
) 1)) it
1646 ;; expdiff >= prec+1 => abs(self/other) > 10**prec
1647 ((cx-error context
) DivisionImpossible
))
1650 ;; expdiff <= -2 => abs(self/other) < 0.1
1651 (let ((ans ((ref self
'_rescale
)
1652 ideal_exponent
(cx-rounding context
))))
1653 ((ref ans
'_fix
) context
)))
1655 (let ((op1 (_WorkRep self
))
1656 (op2 (_WorkRep other
)))
1658 ;; adjust both arguments to have the same exponent, then divide
1659 (if (>= (ref op1
'exp
) (ref op2
'exp
))
1660 (set op1
'int
(* (ref op1
'int
)
1661 (expt 10 (- (ref op1
'exp
) (ref op2
'exp
)))))
1662 (set op2
'int
(* (ref op2
'int
)
1663 (expt 10 (- (ref op2
'exp
) (ref op1
'exp
))))))
1665 (call-with-values (lambda () (divmod (ref op1
'int
) (ref op2
'int
)))
1668 ;; remainder is r*10**ideal_exponent; other is +/-op2.int *
1669 ;; 10**ideal_exponent. Apply correction to ensure that
1670 ;; abs(remainder) <= abs(other)/2
1671 (if (> (+ (* 2 r
) + (logand q
1)) (ref op2
'int
))
1672 (set! r
(- r
(ref op2
'int
)))
1675 (if (>= q
(expt 10 (cx-prec context
)))
1676 ((cx-error context
) DivisionImpossible
)
1677 (let ((sign (ref self
'_sign
)))
1679 (set! sign
(- 1 sign
))
1681 (let ((ans (_dec_from_triple sign
(str r
) ideal_exponent
)))
1682 ((ref ans
'_fix
) context
))))))))))
1684 (define __floordiv__
1685 (lam (self other
(= context None
))
1688 ((norm-op self other
) it it
)
1690 (let (get-context context
))
1692 ((bin-special self other context
) it it
)
1694 (((ref self
'_isinfinity
)) it
1695 (if ((ref other
'_isinfinity
))
1696 ((cx-error context
) InvalidOperation
"INF // INF")
1697 (pylist-ref _SignedInfinity
(logxor (ref self
'_sign
)
1698 (ref other
'_sign
)))))
1700 ((not (bool other
)) it
1702 ((cx-error context
) DivisionByZero
"x // 0"
1703 (logxor (ref self
'_sign
) (ref other
'_sign
)))
1704 ((cx-error context
) DivisionUndefined
"0 // 0")))
1706 ((ref self
'_divide
) other context
))))
1708 (define __rfloordiv__
1709 (lam (self other
(= context None
))
1710 "Swaps self/other and returns __floordiv__."
1712 ((norm-op self other
) it it
)
1713 ((ref other
'__floordiv__
) self
#:context context
))))
1717 "Float representation."
1718 (if ((ref self
'_isnan
))
1719 (if ((ref self
'is_snan
))
1720 (raise (ValueError "Cannot convert signaling NaN to float"))
1721 (if (= (ref self
'_sign
))
1724 (if ((ref self
'_isspecial
))
1725 (if (= (ref self
'_sign
))
1728 (float (str self
))))))
1732 "Converts self to an int, truncating if necessary."
1733 (if ((ref self
'_isnan
))
1734 (raise (ValueError "Cannot convert NaN to integer"))
1735 (if ((ref self
'_isspecial
))
1736 (raise (OverflowError "Cannot convert infinity to integer"))
1737 (let ((s (if (= (ref self
'_sign
) 1) -
1 1)))
1738 (if (>= (ref self
'_exp
) 0)
1739 (* s
(int (ref self
'_int
)) (expt 10 (ref self
'_exp
)))
1740 (* s
(int (or (bool (pylist-slice (ref self
'_int
)
1741 None
(ref self
'_exp
) None
))
1744 (define __trunc__ __int__
)
1747 (property (lambda (self) self
)))
1755 (lambda (self) self
))
1759 (complex (float self
))))
1762 (lambda (self context
)
1763 "Decapitate the payload of a NaN to fit the context"
1764 (let ((payload (ref self
'_int
))
1766 ;; maximum length of payload is precision if clamp=0,
1767 ;; precision-1 if clamp=1.
1769 (- (ref context
'prec
)
1770 (ref context
'clamp
))))
1772 (if (> (len payload
) max_payload_len
)
1773 (let ((payload (py-lstrip
1774 (pylist-slice payload
1775 (- (len payload
) max_payload_len
)
1777 (_dec_from_triple (ref self
'_sign
) payload
(ref self
'_exp
)
1782 (lambda (self context
)
1783 "Round if it is necessary to keep self within prec precision.
1785 Rounds and fixes the exponent. Does not raise on a sNaN.
1788 self - Decimal instance
1789 context - context used.
1793 (((ref self
'_is_special
)) it
1794 (if ((ref self
'_isnan
))
1795 ;; decapitate payload if necessary
1796 ((ref self
'_fix_nan
) context
)
1798 ;; self is +/-Infinity; return unaltered
1801 ;; if self is zero then exponent should be between Etiny and
1802 ;; Emax if clamp==0, and between Etiny and Etop if clamp==1.
1803 (let ((Etiny (cx-etiny context
))
1804 (Etop (cx-etop context
))))
1806 ((not (bool self
)) it
1807 (let* ((exp_max (if (= (cx-clamp context
) 0)
1810 (new_exp (min (max (ref self
'_exp
) Etiny
) exp_max
)))
1811 (if (not (= new_exp
(ref self
'_exp
)))
1813 ((cx-error context
) Clamped
)
1814 (_dec_from_triple (ref self
'_sign
) "0" new_exp
))
1817 ;; exp_min is the smallest allowable exponent of the result,
1818 ;; equal to max(self.adjusted()-context.prec+1, Etiny)
1819 (let ((exp_min (+ (len (ref self
'_int
))
1821 (- (cx-prec context
))))))
1822 ((> exp_min Etop
) it
1823 ;; overflow: exp_min > Etop iff self.adjusted() > Emax
1824 (let ((ans ((cx-error context
) Overflow
"above Emax"
1825 (ref self
'_sign
))))
1826 ((cx-error context
) Inexact
)
1827 ((cx-error context
) Rounded
)
1830 (let* ((self_is_subnormal (< exp_min Etiny
))
1831 (exp_min (if self_is_subnormal Etiny exp_min
))))
1833 ;; round if self has too many digits
1834 ((< (ref self
'_exp
) exp_min
) it
1835 (let ((digits (+ (len (ref self
'_int
))
1839 (set! self
(_dec_from_triple (ref self
'_sign
)
1844 (rounding_method (pylist-ref
1845 (ref self
'_pick_rounding_function
)
1846 (cx-rounding context
)))
1847 (changed (rounding_method self digits
))
1848 (coeff (or (bool (pylist-slice (ref self
'_int
)
1849 None digits None
)) "0")))
1852 (set! coeff
(str (+ (int coeff
) 1)))
1853 (if (> (len coeff
) (cx-prec context
))
1855 (set! coeff
(pylist-slice coeff None -
1 None
))
1856 (set! exp_min
(+ exp_min
1))))))
1858 ;; check whether the rounding pushed the exponent out of range
1859 (if (> exp_min Etop
)
1861 ((cx-error context
) Overflow
"above Emax"
1863 (set! ans
(_dec_from_triple (ref self
'_sign
) coeff exp_min
)))
1865 ;; raise the appropriate signals, taking care to respect
1866 ;; the precedence described in the specification
1867 (if (and changed self_is_subnormal
)
1868 ((cx-error context
) Underflow
))
1869 (if self_is_subnormal
1870 ((cx-error context
) Subnormal
))
1872 ((cx-error context
) Inexact
))
1874 ((cx-error context
) Rounded
)
1876 (if (not (bool ans
))
1877 ;; raise Clamped on underflow to 0
1878 ((cx-error context
) Clamped
))
1882 (if self_is_subnormal
1883 ((cx-error context
) Subnormal
))
1886 ;; fold down if clamp == 1 and self has too few digits
1887 (if (and (= (cx-clamp context
) 1) (> (ref self
'_exp
) Etop
))
1889 ((cx-error context
) Clamped
)
1890 (let ((self_padded (+ (ref self
'_int
)
1892 (- (ref self
'_exp
) Etop
)))))
1893 (_dec_from_triple (ref self
'_sign
) self_padded Etop
)))
1895 ;; here self was representable to begin with; return unchanged
1900 ;; for each of the rounding functions below:
1901 ;; self is a finite, nonzero Decimal
1902 ;; prec is an integer satisfying 0 <= prec < len(self._int)
1904 ;; each function returns either -1, 0, or 1, as follows:
1905 ;; 1 indicates that self should be rounded up (away from zero)
1906 ;; 0 indicates that self should be truncated, and that all the
1907 ;; digits to be truncated are zeros (so the value is unchanged)
1908 ;; -1 indicates that there are nonzero digits to be truncated
1912 "Also known as round-towards-0, truncate."
1913 (if (_all_zeros (ref self
'_int
) prec
)
1919 "Rounds away from 0."
1920 (- (_round_down self prec
))))
1922 (define _round_half_up
1924 "Rounds 5 up (away from 0)"
1926 ((in (pylist-ref (ref self
'_int
) prec
) "56789")
1928 ((_all_zeros (ref self
'_int
) prec
)
1932 (define _round_half_down
1935 (if (_exact_half (ref self
'_int
) prec
)
1937 (_round_half_up self prec
))))
1939 (define _round_half_even
1941 "Round 5 to even, rest to nearest."
1942 (if (and (_exact_half (ref self
'_int
) prec
)
1944 (in (pylist-ref (ref self
'_int
) (- prec
1)) "02468")))
1946 (_round_half_up self prec
)))
1948 (define _round_ceiling
1950 "Rounds up (not away from 0 if negative.)"
1951 (if (= (ref self
'_sign
) 1)
1952 (_round_down self prec
)
1953 (- (_round_down self prec
)))))
1955 (define _round_floor
1957 "Rounds down (not towards 0 if negative)"
1958 (if (= (ref self
'_sign
) 1)
1959 (- (_round_down self prec
))
1960 (_round_down self prec
))))
1964 "Round down unless digit prec-1 is 0 or 5."
1965 (if (and prec
(not (in (pylist-ref (ref self
'_int
) (- prec
1) "05"))))
1966 (_round_down self prec
)
1967 (- (_round_down self prec
)))))
1969 (define _pick_rounding_function
1970 (dict `((,ROUND_DOWN .
,_round_down
)
1971 (,ROUND_UP .
,_round_up
)
1972 (,ROUND_HALF_UP .
,_round_half_up
)
1973 (,ROUND_HALF_DOWN .
,_round_half_down
)
1974 (,ROUND_HALF_EVEN .
,_round_half_even
)
1975 (,ROUND_CEILING .
,_round_ceiling
)
1976 (,ROUND_FLOOR .
,_round_floor
)
1977 (,ROUND_05UP .
,_round_05up
))))
1980 (lam (self (= n None
))
1981 "Round self to the nearest integer, or to a given precision.
1983 If only one argument is supplied, round a finite Decimal
1984 instance self to the nearest integer. If self is infinite or
1985 a NaN then a Python exception is raised. If self is finite
1986 and lies exactly halfway between two integers then it is
1987 rounded to the integer with even last digit.
1989 >>> round(Decimal('123.456'))
1991 >>> round(Decimal('-456.789'))
1993 >>> round(Decimal('-3.0'))
1995 >>> round(Decimal('2.5'))
1997 >>> round(Decimal('3.5'))
1999 >>> round(Decimal('Inf'))
2000 Traceback (most recent call last):
2002 OverflowError: cannot round an infinity
2003 >>> round(Decimal('NaN'))
2004 Traceback (most recent call last):
2006 ValueError: cannot round a NaN
2008 If a second argument n is supplied, self is rounded to n
2009 decimal places using the rounding mode for the current
2012 For an integer n, round(self, -n) is exactly equivalent to
2013 self.quantize(Decimal('1En')).
2015 >>> round(Decimal('123.456'), 0)
2017 >>> round(Decimal('123.456'), 2)
2019 >>> round(Decimal('123.456'), -2)
2021 >>> round(Decimal('-Infinity'), 37)
2023 >>> round(Decimal('sNaN123'), 0)
2027 (if (not (eq? n None
))
2028 ;; two-argument form: use the equivalent quantize call
2029 (if (not (isinstance n int
))
2031 "Second argument to round should be integral"))
2032 (let ((exp (_dec_from_triple 0 "1" (- n
))))
2033 ((ref self
'quantize
) exp
)))
2035 ;; one-argument form§x
2036 (if (ref self
'_is_special
)
2037 (if ((ref self
'is_nan
))
2038 (raise (ValueError "cannot round a NaN"))
2039 (raise (OverflowError "cannot round an infinity")))
2040 (int ((ref self
'_rescale
) 0 ROUND_HALF_EVEN
))))))
2044 "Return the floor of self, as an integer.
2046 For a finite Decimal instance self, return the greatest
2047 integer n such that n <= self. If self is infinite or a NaN
2048 then a Python exception is raised.
2051 (if (ref self
'_is_special
)
2052 (if ((ref self
'is_nan
))
2053 (raise (ValueError "cannot round a NaN"))
2054 (raise (OverflowError "cannot round an infinity")))
2055 (int ((ref self
'_rescale
) 0 ROUND_FLOOR
)))))
2059 """Return the ceiling of self, as an integer.
2061 For a finite Decimal instance self, return the least integer n
2062 such that n >= self. If self is infinite or a NaN then a
2063 Python exception is raised.
2066 (if (ref self
'_is_special
)
2067 (if ((ref self
'is_nan
))
2068 (raise (ValueError "cannot round a NaN"))
2069 (raise (OverflowError "cannot round an infinity")))
2070 (int ((ref self
'_rescale
) 0 ROUND_CEILING
)))))
2073 (lam (self other third
(= context None
))
2074 "Fused multiply-add.
2076 Returns self*other+third with no rounding of the intermediate
2079 self and other are multiplied together, with no rounding of
2080 the result. The third operand is then added to the result,
2081 and a single final rounding is performed.
2084 (let ((other (_convert_other other
#:raiseit
#t
))
2085 (third (_convert_other third
#:raiseit
#t
))
2086 (fin (lambda (product)
2087 ((ref product
'__add__
) third context
)))))
2088 ;; compute product; raise InvalidOperation if either operand is
2089 ;; a signaling NaN or if the product is zero times infinity.
2090 ((if (or (ref self
'_is_special
) (ref other
'_is_special
))
2092 (let (get-context context
))
2093 ((equal?
(ref self
'_exp
) "N") it
2094 ((cx-error context
) InvalidOperation
"sNaN" self
))
2095 ((equal?
(ref other
'_exp
) "N") it
2096 ((cx-error context
) InvalidOperation
"sNaN" other
))
2097 ((equal?
(ref self
'_exp
) "n") it
2099 ((equal?
(ref other
'_exp
) "n") it
2101 ((equal?
(ref self
'_exp
) "F") it
2102 (if (not (bool other
))
2103 ((cx-error context
) InvalidOperation
"INF * 0 in fma")
2104 (pylist-ref _SignedInfinity
2105 (logxor (ref self
'_sign
)
2106 (ref other
'_sign
)))))
2107 ((equal?
(ref other
'_exp
) "F") it
2108 (if (not (bool self
))
2109 ((cx-error context
) InvalidOperation
"0 * INF in fma")
2110 (pylist-ref _SignedInfinity
2111 (logxor (ref self
'_sign
)
2112 (ref other
'_sign
)))))
2116 (_dec_from_triple (logxor (ref self
'_sign
) (ref other
'_sign
))
2117 (str (* (int (ref self
'_int
))
2118 (int (ref other
'_int
))))
2119 (+ (ref self
'_exp
) (ref other
'_exp
)))))))
2121 (define _power_modulo
2122 (lam (self other modulo
(= context None
))
2123 "Three argument version of __pow__"
2125 ((norm-op self other
) it it
)
2126 ((norm-op self modulo
) it it
)
2127 (let (get-context context
))
2129 ;; deal with NaNs: if there are any sNaNs then first one wins,
2130 ;; (i.e. behaviour for NaNs is identical to that of fma)
2131 (let ((self_is_nan (ref self
'_isnan
))
2132 (other_is_nan (ref other
'_isnan
))
2133 (modulo_is_nan (ref modulo
'_isnan
))))
2135 ((or (bool self_is_nan
) (bool other_is_nan
) (bool modulo_is_nan
)) it
2138 ((cx-error context
) InvalidOperation
"sNaN" self
))
2140 ((cx-error context
) InvalidOperation
"sNaN" other
))
2142 ((cx-error context
) InvalidOperation
"sNaN" modulo
))
2144 (_fix_nan self context
))
2145 ((bool other_is_nan
)
2146 (_fix_nan other context
))
2148 (_fix_nan modulo context
))))
2150 ;;check inputs: we apply same restrictions as Python's pow()
2151 ((not (and ((ref self
'_isinteger
))
2152 ((ref other
'_isinteger
))
2153 ((ref modulo
'_isinteger
)))) it
2154 ((cx-error context
) InvalidOperation
2155 (+ "pow() 3rd argument not allowed "
2156 "unless all arguments are integers")))
2159 ((cx-error context
) InvalidOperation
2160 (+ "pow() 2nd argument cannot be "
2161 "negative when 3rd argument specified")))
2163 ((not (bool modulo
)) it
2164 ((cx-error context
) InvalidOperation
2165 "pow() 3rd argument cannot be 0"))
2167 ;; additional restriction for decimal: the modulus must be less
2168 ;; than 10**prec in absolute value
2169 ((>= ((ref modulo
'adjusted
)) (cx-prec context
)) it
2170 ((cx-error context
) InvalidOperation
2171 (+ "insufficient precision: pow() 3rd "
2172 "argument must not have more than "
2173 "precision digits")))
2175 ;; define 0**0 == NaN, for consistency with two-argument pow
2176 ;; (even though it hurts!)
2177 ((and (not (bool other
)) (not (bool self
))) it
2178 ((cx-error context
) InvalidOperation
2179 (+ "at least one of pow() 1st argument "
2180 "and 2nd argument must be nonzero ;"
2181 "0**0 is not defined")))
2183 ;; compute sign of result
2184 (let ((sign (if ((ref other
'_iseven
))
2187 (base (_WorkRep ((ref self
'to_integral_value
))))
2188 (exponent (_WorkRep ((ref other
'to_integral_value
)))))
2191 ;; convert modulo to a Python integer, and self and other to
2192 ;; Decimal integers (i.e. force their exponents to be >= 0)
2193 (set! modulo
(abs (int modulo
)))
2195 ;; compute result using integer pow()
2196 (set! base
(guile:modulo
2197 (* (guile:modulo
(ref base
'int
) modulo
)
2198 (modulo-expt 10 (ref base
'exp
) modulo
))
2201 (let lp
((i (ref exponent
'exp
)))
2204 (set! base
(modulo-expt base
10 modulo
))
2207 (set! base
(modulo-expt base
(ref exponent
'int
) modulo
))
2209 (_dec_from_triple sign
(str base
) 0)))))
2211 (define _power_exact
2212 (lambda (self other p
)
2213 "Attempt to compute self**other exactly.
2215 Given Decimals self and other and an integer p, attempt to
2216 compute an exact result for the power self**other, with p
2217 digits of precision. Return None if self**other is not
2218 exactly representable in p digits.
2220 Assumes that elimination of special cases has already been
2221 performed: self and other must both be nonspecial; self must
2222 be positive and not numerically equal to 1; other must be
2223 nonzero. For efficiency, other._exp should not be too large,
2224 so that 10**abs(other._exp) is a feasible calculation."
2226 ;; In the comments below, we write x for the value of self and y for the
2227 ;; value of other. Write x = xc*10**xe and abs(y) = yc*10**ye, with xc
2228 ;; and yc positive integers not divisible by 10.
2230 ;; The main purpose of this method is to identify the *failure*
2231 ;; of x**y to be exactly representable with as little effort as
2232 ;; possible. So we look for cheap and easy tests that
2233 ;; eliminate the possibility of x**y being exact. Only if all
2234 ;; these tests are passed do we go on to actually compute x**y.
2236 ;; Here's the main idea. Express y as a rational number m/n, with m and
2237 ;; n relatively prime and n>0. Then for x**y to be exactly
2238 ;; representable (at *any* precision), xc must be the nth power of a
2239 ;; positive integer and xe must be divisible by n. If y is negative
2240 ;; then additionally xc must be a power of either 2 or 5, hence a power
2243 ;; There's a limit to how small |y| can be: if y=m/n as above
2246 ;; (1) if xc != 1 then for the result to be representable we
2247 ;; need xc**(1/n) >= 2, and hence also xc**|y| >= 2. So
2248 ;; if |y| <= 1/nbits(xc) then xc < 2**nbits(xc) <=
2249 ;; 2**(1/|y|), hence xc**|y| < 2 and the result is not
2252 ;; (2) if xe != 0, |xe|*(1/n) >= 1, so |xe|*|y| >= 1. Hence if
2253 ;; |y| < 1/|xe| then the result is not representable.
2255 ;; Note that since x is not equal to 1, at least one of (1) and
2256 ;; (2) must apply. Now |y| < 1/nbits(xc) iff |yc|*nbits(xc) <
2257 ;; 10**-ye iff len(str(|yc|*nbits(xc)) <= -ye.
2259 ;; There's also a limit to how large y can be, at least if it's
2260 ;; positive: the normalized result will have coefficient xc**y,
2261 ;; so if it's representable then xc**y < 10**p, and y <
2262 ;; p/log10(xc). Hence if y*log10(xc) >= p then the result is
2263 ;; not exactly representable.
2265 ;; if len(str(abs(yc*xe)) <= -ye then abs(yc*xe) < 10**-ye,
2266 ;; so |y| < 1/xe and the result is not representable.
2267 ;; Similarly, len(str(abs(yc)*xc_bits)) <= -ye implies |y|
2272 (define-syntax-rule (clean xc xe n
+)
2274 (if (= (modulo xc n
) 0)
2280 (let* ((x (_WorkRep self
))
2285 (let* ((y (_WorkRep other
))
2290 ;; case where xc == 1: result is 10**(xe*y), with xe*y
2291 ;; required to be an integer
2295 ;; result is now 10**(xe * 10**ye); xe * 10**ye must be integral
2300 (let ((exponent (* xe
(expt 10 ye
)))
2302 (if (= (ref y
'sign
) 1)
2303 (set! exponent
(- exponent
)))
2305 ;;if other is a nonnegative integer, use ideal exponent
2306 (if (and ((ref other
'_isinteger
)) (= (ref other
'_sign
) 0))
2308 (let ((ideal_exponent (* (ref self
'_exp
) (int other
))))
2309 (set! zeros
(min (- exponent ideal_exponent
) (- p
1)))))
2312 (_dec_from_triple 0 (+ "1" (* "0" zeros
)) (- exponent zeros
)))))
2314 ;; case where y is negative: xc must be either a power
2315 ;; of 2 or a power of 5.
2316 ((= (ref y
'sign
) 1) it
2317 (let ((last_digit (modulo xc
10)))
2321 ((= (modulo last_digit
2) 0)
2322 ;; quick test for power of 2
2324 ((not (= (logand xc
(- xc
)) xc
)) it
2326 ;; now xc is a power of 2; e is its exponent
2327 (let () (set! e
(- (_nbits xc
) 1)))
2331 ;; x = 2**e * 10**xe, e > 0, and y < 0.
2333 ;; The exact result is:
2335 ;; x**y = 5**(-e*y) * 10**(e*y + xe*y)
2337 ;; provided that both e*y and xe*y are integers.
2339 ;; 5**(-e*y) >= 10**p, then the result can't be expressed
2340 ;; exactly with p digits of precision.
2342 ;; Using the above, we can guard against large values of ye.
2343 ;; 93/65 is an upper bound for log(10)/log(5), so if
2345 ;; ye >= len(str(93*p//65))
2349 ;; -e*y >= -y >= 10**ye > 93*p/65 > p*log(10)/log(5),
2351 ;; so 5**(-e*y) >= 10**p, and the coefficient of the result
2352 ;; can't be expressed in p digits.
2354 ;; emax >= largest e such that 5**e < 10**p.
2355 (let ((emax (quotient (* p
93) 65))))
2357 ((>= ye
(len (str emax
))) it
2360 ;; Find -e*y and -xe*y; both must be integers
2362 (set! e
(_decimal_lshift_exact (* e yc
) ye
))
2363 (set! xe
(_decimal_lshift_exact (* xe yc
) ye
)))
2365 ((or (eq? e None
) (eq? xe None
)) it
2372 (set! xc
(expt 5 e
))
2377 ;; e >= log_5(xc) if xc is a power of 5; we have
2378 ;; equality all the way up to xc=5**2658
2379 (let* ((e (quotient (* (_nbits xc
) 28) 65))
2382 (xc (quotient q xz
))
2383 (remainder (modulo q xz
))))
2385 ((not (= remainder
0)) it
2388 (let () (clean xc e
5 -
))
2390 ;; Guard against large values of ye, using the same logic as in
2391 ;; the 'xc is a power of 2' branch. 10/3 is an upper bound for
2393 (let ((emax (quotient (* p
10) 3))))
2395 ((>= ye
(len (str emax
))) it
2399 (set! e
(_decimal_lshift_exact (* e yc
) ye
))
2400 (set! xe
(_decimal_lshift_exact (* xe yc
) ye
)))
2402 ((or (eq? e None
) (eq? xe None
)) it
2409 (set! xc
(expt 2 e
))
2415 ((>= xc
(expt 10 p
)) it it
)
2418 (set! xe
(+ (- e
) (- xe
)))
2419 (_dec_from_triple 0 (str xc
) xe
)))))
2421 ;; now y is positive; find m and n such that y = m/n
2422 (let ((m #f
) (n #f
) (xc_bits (_nbits xc
))))
2425 (set! m
(* yc
(expt 10 ye
)))
2429 ((and (not (= xe
0)) (<= (len (str (abs (* yc xe
)))) (- ye
))) it
2432 ((and (not (= xc
1))
2433 (<= (len (str (* (abs yc
) xc_bits
))) (- ye
))) it
2438 (set! n
(expt 10 (- ye
)))
2441 (if (and (= (modulo m
2) 0) (= (modulo n
2) 0))
2443 (set! m
(quotient m
2))
2444 (set! n
(quotient n
2)))))
2446 (if (and (= (modulo m
5) 0) (= (modulo n
5) 0))
2448 (set! m
(quotient m
5))
2449 (set! n
(quotient n
5)))))
2452 ;; compute nth root of xc*10**xe
2456 ;; if 1 < xc < 2**n then xc isn't an nth power
2457 ((and (not (= xc
1)) (<= xc_bits n
)) it
2460 ((not (= (modulo xe n
) 0)) it
2464 (let ((a (ash 1 (- (quotient (- xc_bits
) n
)))))
2465 (set! xe
(quotient xe n
))
2467 ;; compute nth root of xc using Newton's method
2469 (let* ((x (expt a
(- n
1)))
2473 (if (not (and (= a q
) (= r
0)))
2479 (set! a
(quotient (+ (* a
(- n
1)) q
) n
))
2483 ;; now xc*10**xe is the nth root of the original xc*10**xe
2484 ;; compute mth power of xc*10**xe
2486 ;; if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
2487 ;; 10**p and the result is not representable.
2488 ((and (> xc
1) (> m
(quotient (* p
100) (_log10_lb xc
)))) it
2492 (set! xc
(expt xc m
))
2495 ((> xc
(expt 10 p
)) it
2499 ;; by this point the result *is* exactly representable
2500 ;; adjust the exponent to get as close as possible to the ideal
2501 ;; exponent, if necessary
2502 (let* ((str_xc (str xc
))
2503 (zeros (if (and ((ref other
'_isinteger
))
2504 (= (ref other
'_sign
) 0))
2505 (let ((ideal_exponent
2506 (* (ref self
'_exp
) (int other
))))
2507 (min (- xe ideal_exponent
)
2508 (- p
(len str_xc
))))
2510 (_dec_from_triple 0 (+ str_xc
(* '0' zeros
)) (- xe zeros
)))))))
2513 (lam (self other
(= modulo None
) (= context None
))
2514 "Return self ** other [ % modulo].
2516 With two arguments, compute self**other.
2518 With three arguments, compute (self**other) % modulo. For the
2519 three argument form, the following restrictions on the
2522 - all three arguments must be integral
2523 - other must be nonnegative
2524 - either self or other (or both) must be nonzero
2525 - modulo must be nonzero and must have at most p digits,
2526 where p is the context precision.
2528 If any of these restrictions is violated the InvalidOperation
2531 The result of pow(self, other, modulo) is identical to the
2532 result that would be obtained by computing (self**other) %
2533 modulo with unbounded precision but is computed more
2534 efficiently. It is always exact.
2538 ((not (eq? modulo None
)) it
2539 ((ref self
'_power_modulo
) other modulo context
))
2541 ((norm-op self other
) it it
)
2542 (let (get-context context
))
2544 ;; either argument is a NaN => result is NaN
2545 ((bin-special self other context
) it it
)
2547 ;; 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
2548 ((not (bool other
)) it
2549 (if (not (bool self
))
2550 ((cx-error context
) InvalidOperation
"0 ** 0")
2553 ;; result has sign 1 iff self._sign is 1 and other is an odd integer
2554 (let ((result_sign 0)))
2556 ((if (= (ref self
'_sign
) 1)
2558 ((if ((ref other
'_isinteger
))
2559 (if (not ((ref other
'_iseven
)))
2561 (set! result_sign
1)
2564 ;; -ve**noninteger = NaN
2565 ;; (-0)**noninteger = 0**noninteger
2567 ((cx-error context
) InvalidOperation
2568 "x ** y with x negative and y not an integer")
2571 ;; negate self, without doing any unwanted rounding
2572 (set! self
((ref self
'copy_negate
)))
2576 ;; 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
2577 ((not (bool self
)) it
2578 (if (= (ref other
'_sign
) 0)
2579 (_dec_from_triple result_sign
"0" 0)
2580 (pylist-ref _SignedInfinity result_sign
)))
2582 ;; Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
2583 (((self '_isinfinity
)) it
2584 (if (= (ref other
'_sign
) 0)
2585 (pylist-ref _SignedInfinity result_sign
)
2586 (_dec_from_triple result_sign
"0" 0)))
2588 ;; 1**other = 1, but the choice of exponent and the flags
2589 ;; depend on the exponent of self, and on whether other is a
2590 ;; positive integer, a negative integer, or neither
2591 (let ((prec (cx-prec context
))))
2593 ((equal? self _One
) it
2595 (if ((ref other
'_isinteger
))
2596 ;; exp = max(self._exp*max(int(other), 0),
2597 ;; 1-context.prec) but evaluating int(other) directly
2598 ;; is dangerous until we know other is small (other
2599 ;; could be 1e999999999)
2602 ((= (ref other
'_sign
) 1)
2609 (set! exp
(* (ref self
'_exp
) multiplier
))
2610 (if (< exp
(- 1 prec
))
2612 (set! exp
(- 1 prec
))
2613 ((cx-error context
) Rounded
))))
2615 ((cx-error context
) Inexact
)
2616 ((cx-error context
) Rounded
)
2617 (set! exp
(- 1 prec
))))
2619 (_dec_from_triple result_sign
(+ "1" (* "0" (- exp
)) exp
))))
2621 ;; compute adjusted exponent of self
2622 (let ((self_adj ((ref self
'adjusted
)))))
2624 ;; self ** infinity is infinity if self > 1, 0 if self < 1
2625 ;; self ** -infinity is infinity if self < 1, 0 if self > 1
2626 (((ref other
'_isinfinity
)) it
2627 (if (eq?
(= (ref other
'_sign
) 0)
2629 (_dec_from_triple result_sign
"0" 0)
2630 (pylist-ref _SignedInfinity result_sign
)))
2632 ;; from here on, the result always goes through the call
2633 ;; to _fix at the end of this function.
2636 (bound (+ ((ref self
'_log10_exp_bound
))
2637 ((ref other
'adjusted
)))))
2639 ;; crude test to catch cases of extreme overflow/underflow. If
2640 ;; log10(self)*other >= 10**bound and bound >= len(str(Emax))
2641 ;; then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
2642 ;; self**other >= 10**(Emax+1), so overflow occurs. The test
2643 ;; for underflow is similar.
2645 (if (eq?
(>= self_adj
0) (= (ref other
'_sign
) 0))
2646 ;; self > 1 and other +ve, or self < 1 and other -ve
2647 ;; possibility of overflow
2648 (if (>= bound
(len (str (cx-emax context
))))
2650 (_dec_from_triple result_sign
"1"
2651 (+ (cx-emax context
) 1))))
2653 ;; self > 1 and other -ve, or self < 1 and other +ve
2654 ;; possibility of underflow to 0
2655 (let ((Etiny (cx-etiny context
)))
2656 (if (>= bound
(len (str (- Etiny
))))
2657 (set! ans
(_dec_from_triple result_sign
"1" (- Etiny
1))))))
2659 ;; try for an exact result with precision +1
2660 (when (eq? ans None
)
2661 (set! ans
((ref self
'_power_exact
) other
(+ prec
1)))
2662 (when (not (eq? ans None
))
2663 (if (= result_sign
1)
2664 (set! ans
(_dec_from_triple 1 (ref ans
'_int
)
2668 ;; usual case: inexact result, x**y computed directly as exp(y*log(x))
2669 (when (eq? ans None
)
2676 (y (_WorkRep other
))
2680 (if (= (ref y
'sign
) 1)
2683 ;; compute correctly rounded result: start with precision +3,
2684 ;; then increase precision until result is unambiguously roundable
2689 (lambda () (_dpower xc xe yc ye
(+ p extra
)))
2692 (* 5 (expt 10 (- (len (str coeff
)) p
1))))
2694 (lp (+ extra
3)))))))
2696 (set! ans
(_dec_from_triple result_sign
(str coeff
) exp
))))))
2698 ;; unlike exp, ln and log10, the power function respects the
2699 ;; rounding mode; no need to switch to ROUND_HALF_EVEN here
2701 ;; There's a difficulty here when 'other' is not an integer and
2702 ;; the result is exact. In this case, the specification
2703 ;; requires that the Inexact flag be raised (in spite of
2704 ;; exactness), but since the result is exact _fix won't do this
2705 ;; for us. (Correspondingly, the Underflow signal should also
2706 ;; be raised for subnormal results.) We can't directly raise
2707 ;; these signals either before or after calling _fix, since
2708 ;; that would violate the precedence for signals. So we wrap
2709 ;; the ._fix call in a temporary context, and reraise
2711 (if (and exact
(not ((ref other
'_isinteger
))))
2713 ;; pad with zeros up to length context.prec+1 if necessary; this
2714 ;; ensures that the Rounded signal will be raised.
2715 (if (<= (len (ref ans
'_int
)) prec
)
2716 (let ((expdiff (+ prec
1 (- (len (ref ans
'_int
))))))
2717 (set! ans
(_dec_from_triple (ref ans
'_sign
)
2720 (- (ref ans
'_exp
) expdiff
)))))
2722 ;; create a copy of the current context, with cleared flags/traps
2723 (let ((newcontext (cx-copy context
)))
2724 (cx-clear_flags newcontext
)
2726 (for ((exception : _signals
)) ()
2727 (pylist-set! (cx-traps newcontext
) exception
0)
2730 ;; round in the new context
2731 (set! ans
((ref ans
'_fix
) newcontext
))
2733 ;; raise Inexact, and if necessary, Underflow
2734 ((cx-error newcontext
) Inexact
)
2735 (if (bool (pylist-ref (cx-flags newcontext
) Subnormal
))
2736 ((cx-error newcontext
) Underflow
))
2738 ;; propagate signals to the original context; _fix could
2739 ;; have raised any of Overflow, Underflow, Subnormal,
2740 ;; Inexact, Rounded, Clamped. Overflow needs the correct
2741 ;; arguments. Note that the order of the exceptions is
2743 (if (bool (pylist-ref (cx-flags newcontext
) Overflow
))
2744 ((cx-error newcontext
)
2745 Overflow
"above Emax" (ref ans
'_sign
)))
2747 (for ((exception : (list Underflow Subnormal
2748 Inexact Rounded Clamped
))) ()
2749 (if (bool (pylist-ref (cx-flags newcontext
) exception
))
2750 ((cx-error newcontext
) exception
)
2753 (set! ans
((ref ans
'_fix
) context
)))
2758 (lam (self other
(= context None
))
2759 "Swaps self/other and returns __pow__."
2761 ((norm-op self other
) it it
)
2762 ((ref 'other
'__pow__
) self
#:context context
))))
2765 (lam (self (= context None
))
2766 "Normalize- strip trailing 0s, change anything equal to 0 to 0e0"
2769 (let (get-context context
))
2770 ((un-special self context
) it it
)
2771 (let ((dup ((ref self _fix
) context
))))
2772 (((dup '_isinfinity
)) it dup
)
2773 ((not (bool dup
)) it
2774 (_dec_from_triple (ref dup
'_sign
) "0" 0))
2776 (let* ((_int (ref dup
'_int
))
2777 (exp_max (let ((i (cx-clamp context
)))
2780 (cx-etop context
)))))