small steps furter regarding decimal
[software/python-on-guile.git] / modules / language / python / module / decimal.scm
1 (define-module (language python module decimal)
2 #:use-module ((language python module collections) #:select (namedtuple))
3 #:export ())
4
5 (define guile:modulo (@ (guile) moduolo))
6
7 (define __name__ "decimal")
8 (define __xname__ __name__)
9 (define __version__ "1.70")
10 ;; Highest version of the spec this complies with
11 ;; See http://speleotrove.com/decimal/
12
13
14 (define DecimalTuple (namedtuple "DecimalTuple" "sign digits exponent"))
15
16 ;; Rounding
17 (define ROUND_DOWN 'ROUND_DOWN)
18 (define ROUND_HALF_UP 'ROUND_HALF_UP)
19 (define ROUND_HALF_EVEN 'ROUND_HALF_EVEN)
20 (define ROUND_CEILING 'ROUND_CEILING)
21 (define ROUND_FLOOR 'ROUND_FLOOR)
22 (define ROUND_UP 'ROUND_UP)
23 (define ROUND_HALF_DOWN 'ROUND_HALF_DOWN)
24 (define ROUND_05UP 'ROUND_05UP)
25
26 ;; Compatibility with the C version
27 (define MAX_PREC 425000000)
28 (define MAX_EMAX 425000000)
29 (define MIN_EMIN -425000000)
30
31 (if (= sys:maxsize (- (ash 1 63) 1))
32 (begin
33 (set! MAX_PREC 999999999999999999)
34 (set! MAX_EMAX 999999999999999999)
35 (set! MIN_EMIN -999999999999999999)))
36
37 (define MIN_ETINY (- MIN_EMIN (- MAX_PREC 1)))
38
39 ;; Context
40 (define (cx-prec x) (vector-ref x 0))
41 (define (cx-emax x) (vector-ref x 1))
42 (define (cx-raise x) (vector-ref x 2))
43 (define (cx-error x) (vector-ref x 3))
44 (define (cx-capitals x) (vector-ref x 4))
45 (define (cx-rounding x) (vector-ref x 5))
46 ;; Errors
47
48 (define-python-class DecimalException (ArithmeticError)
49 "Base exception class.
50
51 Used exceptions derive from this.
52 If an exception derives from another exception besides this (such as
53 Underflow (Inexact, Rounded, Subnormal) that indicates that it is only
54 called if the others are present. This isn't actually used for
55 anything, though.
56
57 handle -- Called when context._raise_error is called and the
58 trap_enabler is not set. First argument is self, second is the
59 context. More arguments can be given, those being after
60 the explanation in _raise_error (For example,
61 context._raise_error(NewError, '(-x)!', self._sign) would
62 call NewError().handle(context, self._sign).)
63
64 To define a new exception, it should be sufficient to have it derive
65 from DecimalException.
66 "
67
68 (define handle
69 (lambda (self context . args)
70 (values))))
71
72
73 (define-python-class Clamped (DecimalException)
74 """Exponent of a 0 changed to fit bounds.
75
76 This occurs and signals clamped if the exponent of a result has been
77 altered in order to fit the constraints of a specific concrete
78 representation. This may occur when the exponent of a zero result would
79 be outside the bounds of a representation, or when a large normal
80 number would have an encoded exponent that cannot be represented. In
81 this latter case, the exponent is reduced to fit and the corresponding
82 number of zero digits are appended to the coefficient ("fold-down").
83 """)
84
85 (define-python-class InvalidOperation (DecimalException)
86 "An invalid operation was performed.
87
88 Various bad things cause this:
89
90 Something creates a signaling NaN
91 -INF + INF
92 0 * (+-)INF
93 (+-)INF / (+-)INF
94 x % 0
95 (+-)INF % x
96 x._rescale( non-integer )
97 sqrt(-x) , x > 0
98 0 ** 0
99 x ** (non-integer)
100 x ** (+-)INF
101 An operand is invalid
102
103 The result of the operation after these is a quiet positive NaN,
104 except when the cause is a signaling NaN, in which case the result is
105 also a quiet NaN, but with the original sign, and an optional
106 diagnostic information.
107 "
108 (define handle
109 (lambda (self context . args)
110 (if (bool args)
111 (let ((ans (_dec_from_triple
112 (ref (car args) '_sign)
113 (ref (car args) '_int)
114 "n" #t)))
115 ((ref ans '_fix_nan) context))
116 _NaN))))
117
118 (define-python-class ConversionSyntax (InvalidOperation)
119 "Trying to convert badly formed string.
120
121 This occurs and signals invalid-operation if a string is being
122 converted to a number and it does not conform to the numeric string
123 syntax. The result is [0,qNaN].
124 "
125 (define handle
126 (lambda x _NaN)))
127
128 (define-python-class DivisionByZero (DecimalException ZeroDivisionError)
129 "Division by 0.
130
131 This occurs and signals division-by-zero if division of a finite number
132 by zero was attempted (during a divide-integer or divide operation, or a
133 power operation with negative right-hand operand), and the dividend was
134 not zero.
135
136 The result of the operation is [sign,inf], where sign is the exclusive
137 or of the signs of the operands for divide, or is 1 for an odd power of
138 -0, for power.
139 "
140
141 (define handle
142 (lambda (self context sign . args)
143 (pylist-ref _SignedInfinity sign))))
144
145 (define-python-class DivisionImpossible (InvalidOperation)
146 "Cannot perform the division adequately.
147
148 This occurs and signals invalid-operation if the integer result of a
149 divide-integer or remainder operation had too many digits (would be
150 longer than precision). The result is [0,qNaN].
151 "
152
153 (define handle
154 (lambda x _NaN)))
155
156 (define-python-class DivisionUndefined (InvalidOperation ZeroDivisionError)
157 "Undefined result of division.
158
159 This occurs and signals invalid-operation if division by zero was
160 attempted (during a divide-integer, divide, or remainder operation), and
161 the dividend is also zero. The result is [0,qNaN].
162 "
163
164 (define handle
165 (lambda x _NaN)))
166
167 (define-python-class Inexact (DecimalException)
168 "Had to round, losing information.
169
170 This occurs and signals inexact whenever the result of an operation is
171 not exact (that is, it needed to be rounded and any discarded digits
172 were non-zero), or if an overflow or underflow condition occurs. The
173 result in all cases is unchanged.
174
175 The inexact signal may be tested (or trapped) to determine if a given
176 operation (or sequence of operations) was inexact.
177 ")
178
179 (define-python-class InvalidContext (InvalidOperation)
180 "Invalid context. Unknown rounding, for example.
181
182 This occurs and signals invalid-operation if an invalid context was
183 detected during an operation. This can occur if contexts are not checked
184 on creation and either the precision exceeds the capability of the
185 underlying concrete representation or an unknown or unsupported rounding
186 was specified. These aspects of the context need only be checked when
187 the values are required to be used. The result is [0,qNaN].
188 "
189
190 (define handle
191 (lambda x _NaN)))
192
193 (define-python-class Rounded (DecimalException)
194 "Number got rounded (not necessarily changed during rounding).
195
196 This occurs and signals rounded whenever the result of an operation is
197 rounded (that is, some zero or non-zero digits were discarded from the
198 coefficient), or if an overflow or underflow condition occurs. The
199 result in all cases is unchanged.
200
201 The rounded signal may be tested (or trapped) to determine if a given
202 operation (or sequence of operations) caused a loss of precision.
203 ")
204
205 (define-python-class Subnormal (DecimalException)
206 "Exponent < Emin before rounding.
207
208 This occurs and signals subnormal whenever the result of a conversion or
209 operation is subnormal (that is, its adjusted exponent is less than
210 Emin, before any rounding). The result in all cases is unchanged.
211
212 The subnormal signal may be tested (or trapped) to determine if a given
213 or operation (or sequence of operations) yielded a subnormal result.
214 ")
215
216 (define-python-class Overflow (Inexact Rounded)
217 "Numerical overflow.
218
219 This occurs and signals overflow if the adjusted exponent of a result
220 (from a conversion or from an operation that is not an attempt to divide
221 by zero), after rounding, would be greater than the largest value that
222 can be handled by the implementation (the value Emax).
223
224 The result depends on the rounding mode:
225
226 For round-half-up and round-half-even (and for round-half-down and
227 round-up, if implemented), the result of the operation is [sign,inf],
228 where sign is the sign of the intermediate result. For round-down, the
229 result is the largest finite number that can be represented in the
230 current precision, with the sign of the intermediate result. For
231 round-ceiling, the result is the same as for round-down if the sign of
232 the intermediate result is 1, or is [0,inf] otherwise. For round-floor,
233 the result is the same as for round-down if the sign of the intermediate
234 result is 0, or is [1,inf] otherwise. In all cases, Inexact and Rounded
235 will also be raised.
236 "
237
238 (define handle
239 (let ((l (list ROUND_HALF_UP ROUND_HALF_EVEN
240 ROUND_HALF_DOWN ROUND_U)))
241 (lambda (self context sign . args)
242 (let/ec ret
243 (if (memq (ref context 'rounding) l)
244 (ret (pylist-ref _SignedInfinity sign)))
245
246 (if (= sign 0)
247 (if (eq? (ref context 'rounding) ROUND_CEILING)
248 (ret (pylist-ref _SignedInfinity sign))
249 (ret (_dec_from_triple
250 sign
251 (* "9" (cx-prec context))
252 (+ (- (cx-emax context) (cx-prec context)) 1)))))
253
254 (if (= sign 1)
255 (if (eq? (ref context 'rounding) ROUND_FLOOR)
256 (ret (pylist-ref _SignedInfinity sign))
257 (ret (_dec_from_triple
258 sign
259 (* "9" (cx-prec context))
260 (+ (- (cx-emax context) (cx-prec context)) 1))))))))))
261
262
263 (define-python-class Underflow (Inexact Rounded Subnormal)
264 "Numerical underflow with result rounded to 0.
265
266 This occurs and signals underflow if a result is inexact and the
267 adjusted exponent of the result would be smaller (more negative) than
268 the smallest value that can be handled by the implementation (the value
269 Emin). That is, the result is both inexact and subnormal.
270
271 The result after an underflow will be a subnormal number rounded, if
272 necessary, so that its exponent is not less than Etiny. This may result
273 in 0 with the sign of the intermediate result and an exponent of Etiny.
274
275 In all cases, Inexact, Rounded, and Subnormal will also be raised.
276 ")
277
278 (define-python-class FloatOperation (DecimalException TypeError)
279 """Enable stricter semantics for mixing floats and Decimals.
280
281 If the signal is not trapped (default), mixing floats and Decimals is
282 permitted in the Decimal() constructor, context.create_decimal() and
283 all comparison operators. Both conversion and comparisons are exact.
284 Any occurrence of a mixed operation is silently recorded by setting
285 FloatOperation in the context flags. Explicit conversions with
286 Decimal.from_float() or context.create_decimal_from_float() do not
287 set the flag.
288
289 Otherwise (the signal is trapped), only equality comparisons and explicit
290 conversions are silent. All other mixed operations raise FloatOperation.
291 """)
292
293 ;; List of public traps and flags
294 (define _signals
295 (vector Clamped DivisionByZero Inexact Overflow Rounded,
296 Underflow InvalidOperation Subnormal FloatOperation))
297
298 ;; Map conditions (per the spec) to signals
299 (define _condition_map
300 `((,ConversionSyntax . ,InvalidOperation)
301 (,DivisionImpossible . ,InvalidOperation)
302 (,DivisionUndefined . ,InvalidOperation)
303 (,InvalidContext . ,InvalidOperation)))
304
305 ;; Valid rounding modes
306 (define _rounding_modes
307 (list ROUND_DOWN ROUND_HALF_UP ROUND_HALF_EVEN ROUND_CEILING,
308 ROUND_FLOOR ROUND_UP ROUND_HALF_DOWN ROUND_05UP))
309
310 ;; ##### Context Functions ##################################################
311
312 ;; The getcontext() and setcontext() function manage access to a thread-local
313 ;; current context.
314 (define *context* (make-fluid #f))
315 (define (getcontext)
316 (fluid-ref *context*))
317 (define (setcontext context)
318 (fluid-set! *context* context))
319
320 ;; ##### Decimal class #######################################################
321
322 ;; Do not subclass Decimal from numbers.Real and do not register it as such
323 ;; (because Decimals are not interoperable with floats). See the notes in
324 ;; numbers.py for more detail.
325
326 (define _dec_from_triple
327 (lam (sign coefficient exponent (= special #f))
328 "Create a decimal instance directly, without any validation,
329 normalization (e.g. removal of leading zeros) or argument
330 conversion.
331
332 This function is for *internal use only*.
333 "
334 (Decimal sign coeficient exponent special)))
335
336 (def _mk (self (= value "0") (= context None))
337 "Create a decimal point instance.
338
339 >>> Decimal('3.14') # string input
340 Decimal('3.14')
341 >>> Decimal((0, (3, 1, 4), -2)) # tuple (sign, digit_tuple, exponent)
342 Decimal('3.14')
343 >>> Decimal(314) # int
344 Decimal('314')
345 >>> Decimal(Decimal(314)) # another decimal instance
346 Decimal('314')
347 >>> Decimal(' 3.14 \\n') # leading and trailing whitespace okay
348 Decimal('3.14')
349 "
350
351 ;; Note that the coefficient, self._int, is actually stored as
352 ;; a string rather than as a tuple of digits. This speeds up
353 ;; the "digits to integer" and "integer to digits" conversions
354 ;; that are used in almost every arithmetic operation on
355 ;; Decimals. This is an internal detail: the as_tuple function
356 ;; and the Decimal constructor still deal with tuples of
357 ;; digits.
358
359 ;; From a string
360 ;; REs insist on real strings, so we can too.
361 (cond
362 ((isinstance value str)
363 (let ((m (parser (scm-str str))))
364 (if (not m)
365 (let ((context (if (eq? context None)
366 (getcontext)
367 context)))
368 ((cx-raise context)
369 ConversionSyntax
370 (+ "Invalid literal for Decimal: " value))))
371
372 (let ((sign (get-parsed-sign m))
373 (intpart (get-parsed-int m))
374 (fracpart (get-parsed-frac m))
375 (exp (get-parsed-exp m))
376 (diag (get-parsed-diag m))
377 (signal (get-parsed-sig m)))
378
379 (set self 'sign sign)
380
381 (if (not (eq? intpart None))
382 (begin
383 (set self '_int (str (int (+ intpart fracpart))))
384 (set self '_exp (- exp (len fracpart)))
385 (set self '_is_special False))
386 (begin
387 (if (not (eq? diag None))
388 (begin
389 ;; NaN
390 (set self '_int
391 (py-lstrip (str (int (if (bool diag)
392 diag
393 "0")))
394 "0"))
395 (if signal
396 (set self '_exp "N")
397 (set self '_exp "n")))
398 (begin
399 ;; infinity
400 (set self '_int "0")
401 (set self '_exp "F")))
402 (set self '_is_special #t))))))
403
404 ;; From an integer
405 ((isinstance value int)
406 (if (>= value 0)
407 (set self '_sign 0)
408 (set self '_sign 1))
409 (set self '_exp 0)
410 (set self '_int (str (abs value)))
411 (set self '_is_special #f))
412
413 ;; From another decimal
414 ((isinstance value Decimal)
415 (set self '_exp (ref value '_exp ))
416 (set self '_sign (ref value '_sign ))
417 (set self '_int (ref value '_int ))
418 (set self '_is_special (ref value '_is_special)))
419
420 ;; From an internal working value
421 ((isinstance value _WorkRep)
422 (set self '_exp (int (ref value '_exp)))
423 (set self '_sign (ref value '_sign))
424 (set self '_int (str (ref value 'int)))
425 (set self '_is_special #f))
426
427 ;; tuple/list conversion (possibly from as_tuple())
428 ((isinstance value (list list tuple))
429 (if (not (= (len value) 3))
430 (raise (ValueError
431 (+ "Invalid tuple size in creation of Decimal "
432 "from list or tuple. The list or tuple "
433 "should have exactly three elements."))))
434 ;; # process sign. The isinstance test rejects floats
435 (let ((v0 (pylist-ref value 0))
436 (v1 (pylist-ref value 1))
437 (v2 (pylist-ref value 2)))
438 (if (not (and (isinstance v0 int)
439 (or (= v0 0) (= v0 1))))
440 (raise (ValueError
441 (+ "Invalid sign. The first value in the tuple "
442 "should be an integer; either 0 for a "
443 "positive number or 1 for a negative number."))))
444 (set self '_sign v0)
445 (if (eq? v2 'F)
446 (begin
447 (set self '_int "0")
448 (set self '_exp v2)
449 (set self 'is_special #t))
450 (let ((digits (py-list)))
451 ;; process and validate the digits in value[1]
452 (for ((digit : v1)) ()
453 (if (and (isinstance digit int)
454 (<= 0 digit)
455 (<= digit 9))
456 ;; skip leading zeros
457 (if (or (bool digits) (> digit 0))
458 (pylist-append digits digit))
459 (raise (ValueError
460 (+ "The second value in the tuple must "
461 "be composed of integers in the range "
462 "0 through 9.")))))
463
464 (cond
465 ((or (eq? v2 'n) (eq? v2 'N))
466 (begin
467 ;; NaN: digits form the diagnostic
468 (set self '_int (py-join "" (map str digits)))
469 (set self '_exp v2)
470 (set self '_is_special #t)))
471 ((isinstance v2 int)
472 ;; finite number: digits give the coefficient
473 (set self '_int (py-join "" (map str digits)))
474 (set self '_exp v2)
475 (set self '_is_special #f))
476 (else
477 (raise (ValueError
478 (+ "The third value in the tuple must "
479 "be an integer, or one of the "
480 "strings 'F', 'n', 'N'.")))))))))
481
482 ((isinstance value float)
483 (let ((context (if (eq? context None)
484 (getcontext)
485 context)))
486 ((cx-error context)
487 FloatOperation,
488 (+ "strict semantics for mixing floats and Decimals are "
489 "enabled"))
490
491 (__init__ self ((ref Decimal 'from_float) value))))
492
493 (else
494 (raise (TypeError
495 (format #f "Cannot convert %r to Decimal" value))))))
496
497 (define-inlinable (divmod x y)
498 (values (quotient x y) (modulo x y)))
499
500 (define-syntax twix
501 (syntax-rules (let)
502 ((_ a) a)
503 ((_ (let (a ...)) . l)
504 (a ... (twix - l)))
505 ((_ (a it code ...) . l)
506 (aif it a (begin code ...) (twix - l)))))
507
508 (define-syntax-rule (norm-op op)
509 (begin
510 (set! op ((ref self '_convert_other) op))
511 (if (eq? op NotImplemented)
512 other
513 #f)))
514
515 (define-syntax-rule (get-context context code)
516 (let ((context (if (eq? context None)
517 (getcontext)
518 context)))
519 code))
520
521 (define-syntax-rule (un-special self context)
522 (if ((ref self '_is_special))
523 (let ((ans ((ref self '_check_nans) #:context context)))
524 (if (bool ans)
525 (ret ans)
526 #f))
527 #f))
528
529 (define-syntax-rule (bin-special o1 o2 context)
530 (if (or (ref o1 '_is_special)
531 (ref o2 '_is_special))
532 (or (un-special o1 context) (un-special o2 context))))
533
534 (define-syntax-rule (add-special self other context)
535 (or (bin-special self other context)
536 (if ((ref self '_isinfinity))
537 ;; If both INF, same sign =>
538 ;; same as both, opposite => error.
539 (if (and (not (= (ref self '_sign) (ref other '_sign)))
540 ((ref other '_isinfinity)))
541 ((cx-error context) InvalidOperation "-INF + INF")
542 (Decimal self))
543 (if ((ref other '_isinfinity))
544 (ret (Decimal other)) ; Can't both be infinity here
545 #f))))
546
547 (define-syntax-rule (mul-special self other context)
548 (if (or (ref self '_is_special) (ref other '_is_special))
549 (twix
550 ((bin-special self other context) it it)
551
552 ((if ((ref self '_isinfinity))
553 (if (not (bool other))
554 ((cx-error context) InvalidOperation "(+-)INF * 0")
555 (pylist-ref _SignedInfinity resultsign))
556 #f) it it)
557
558 (if ((ref other '_isinfinity))
559 (if (not (bool self))
560 ((cx-error context) InvalidOperation "(+-)INF * 0")
561 (pylist-ref _SignedInfinity resultsign))
562 #f))
563 #f))
564
565 (define-syntax-rule (div-special self other context)
566 (if (or (ref self '_is_special) (ref other '_is_special))
567 (twix
568 ((bin-special self other context) it it)
569
570 ((and ((ref self '_isinfinity)) ((ref other '_isinfinity))) it
571 ((cx-error context) InvalidOperation "(+-)INF/(+-)INF"))
572
573 (((ref self '_isinfinity)) it
574 (pylist-ref _SignedInfinity sign))
575
576 (((ref other '_isinfinity)) it
577 ((cx-error context) Clamped "Division by infinity")
578 (_dec_from_triple sign "0", (cx-etiny context))))))
579
580
581 (define-python-class Decimal (object)
582 "Floating point class for decimal arithmetic."
583
584
585 ;; Generally, the value of the Decimal instance is given by
586 ;; (-1)**_sign * _int * 10**_exp
587 ;; Special values are signified by _is_special == True
588
589 (define __init__
590 (case-lambda
591 ((self sign coefficient exponent special)
592 (set self '_sign sign)
593 (set self '_int coefficient)
594 (set self '_exp exponent)
595 (set self '_is_special special))
596
597 ((self)
598 (_mk self))
599 ((self a)
600 (_mk self a))
601 ((self a b)
602 (_mk self a b))))
603
604 (define from_float
605 (classmethod
606 (lambda (cls f)
607 "Converts a float to a decimal number, exactly.
608
609 Note that Decimal.from_float(0.1) is not the same as Decimal('0.1').
610 Since 0.1 is not exactly representable in binary floating point, the
611 value is stored as the nearest representable value which is
612 0x1.999999999999ap-4. The exact equivalent of the value in decimal
613 is 0.1000000000000000055511151231257827021181583404541015625.
614
615 >>> Decimal.from_float(0.1)
616 Decimal('0.1000000000000000055511151231257827021181583404541015625')
617 >>> Decimal.from_float(float('nan'))
618 Decimal('NaN')
619 >>> Decimal.from_float(float('inf'))
620 Decimal('Infinity')
621 >>> Decimal.from_float(-float('inf'))
622 Decimal('-Infinity')
623 >>> Decimal.from_float(-0.0)
624 Decimal('-0')
625
626 "
627 (cond
628 ((isinstance f int) ; handle integer inputs
629 (cls f))
630 ((not (isinstance f float))
631 (raise (TypeError "argument must be int or float.")))
632 ((or (inf? f) (nan? f))
633 (cls (cond
634 ((nan? f) "")
635 ((eq? f (inf)) "")
636 (eq? f (- (inf))) "")))
637 (else
638 (let* ((sign (if (>= f 0) 0 1))
639 (me (frexp f))
640 (m (car me))
641 (e (cadr me))
642 (res (_dec_from_triple sign, str(m) e)))
643 (if (eq? cls Decimal)
644 res
645 (cls res))))))))
646
647 (define _isnan
648 (lambda (self)
649 "Returns whether the number is not actually one.
650
651 0 if a number
652 1 if NaN
653 2 if sNaN
654 "
655 (if (ref self '_is_special)
656 (let ((exp (ref self '_exp)))
657 (cond
658 ((eq? exp 'n) 1)
659 ((eq? exp 'N) 2)
660 (else 0)))
661 0)))
662
663 (define _isinfinity
664 (lambda (self)
665 "Returns whether the number is infinite
666
667 0 if finite or not a number
668 1 if +INF
669 -1 if -INF
670 "
671 (if (eq? (ref self '_exp) 'F)
672 (if (eq? (ref self '_sign) 1)
673 -1
674 1)
675 0)))
676
677 (define _check_nans
678 (lam (self (= other None) (= context None))
679 "Returns whether the number is not actually one.
680
681 if self, other are sNaN, signal
682 if self, other are NaN return nan
683 return 0
684
685 Done before operations.
686 "
687
688 (let ((self_is_nan ((ref self '_isnan)))
689 (other_is_nan
690 (if (eq? other None)
691 #f
692 ((ref other '_isnan)))))
693
694 (if (or self_is_nan other_is_nan)
695 (let ((context (if (eq? context None)
696 (getcontext)
697 context)))
698 (cond
699 ((eq? self_is_nan 2)
700 ((cx-error context) InvalidOperation "sNaN" self))
701 ((eq? other_is_nan 2)
702 ((cx-error context) InvalidOperation "sNaN" other))
703 (self_is_nan
704 ((ref self '_fix_nan) context))
705 (else
706 ((ref other '_fix_nan) context))))
707 0))))
708
709
710 (define _compare_check_nans
711 (lambda (self other context)
712 "Version of _check_nans used for the signaling comparisons
713 compare_signal, __le__, __lt__, __ge__, __gt__.
714
715 Signal InvalidOperation if either self or other is a (quiet
716 or signaling) NaN. Signaling NaNs take precedence over quiet
717 NaNs.
718
719 Return 0 if neither operand is a NaN.
720
721 "
722 (let ((context (if (eq? context None)
723 (getcontext)
724 context)))
725
726 (if (or (ref self '_is_special)
727 (ref other '_is_special))
728 (cond
729 (((ref self 'is_snan))
730 ((cx-error context)
731 InvalidOperation
732 "comparison involving sNaN" self))
733
734 (((ref other 'is_snan))
735 ((cx-error context)
736 InvalidOperation
737 "comparison involving sNaN" other))
738
739 (((ref self 'is_qnan))
740 ((cx-error context)
741 InvalidOperation
742 "comparison involving NaN" self))
743
744 (((ref other 'is_qnan))
745 ((cx-error context)
746 InvalidOperation
747 "comparison involving NaN" other))
748
749 (else 0))
750 0))))
751
752 (define __bool__
753 (lambda (self)
754 "Return True if self is nonzero; otherwise return False.
755
756 NaNs and infinities are considered nonzero.
757 "
758 (or (ref self '_is_special) (not (equal (ref self '_int) "0")))))
759
760 (define _cmp
761 (lambda (self other)
762 "Compare the two non-NaN decimal instances self and other.
763
764 Returns -1 if self < other, 0 if self == other and 1
765 if self > other. This routine is for internal use only."
766
767 (let ((self_sign (ref self '_sign))
768 (other_sign (ref other '_sign)))
769 (cond
770 ((or (ref self '_is_special) (ref other '_is_special))
771 (let ((self_inf ((ref self '_isinfinity)))
772 (other_inf ((ref other '_isinfinity))))
773 (cond
774 ((eq? self_inf other_inf) 0)
775 ((< self_inf other_inf) -1)
776 (else 1)))
777
778 ;; check for zeros; Decimal('0') == Decimal('-0')
779 ((not (bool self))
780 (if (not (bool other))
781 0
782 (let ((s (ref other '_sign)))
783 (if (= s 0)
784 -1
785 1))))
786 ((not (bool other))
787 (let ((s (ref self '_sign)))
788 (if (= s 0)
789 1
790 -1)))
791
792 ((< other_sign self_sign)
793 -1)
794 ((< self_sign other_sign)
795 1)
796
797 (else
798 (let ((self_adjusted ((ref self 'adjusted)))
799 (other_adjusted ((ref other 'adjusted)))
800 (self_exp (ref self '_exp))
801 (other_exp (ref other '_exp)))
802 (cond
803 ((= self_adjusted other_adjusted)
804 (let ((self_padded (+ (ref self '_int)
805 (* "0" (- self_exp other_exp))))
806 (other_padded (+ (ref other '_int)
807 (* "0" (- other_exp self_exp)))))
808 (cond
809 ((equal? self_padded other_padded)
810 0)
811 ((< self_padded other_padded)
812 (if (= self_sign 0)
813 -1
814 1))
815 (else
816 (if (= self_sign 0)
817 1
818 -1)))))
819 ((> self_adjusted other_adjusted)
820 (if (= self_sign 0)
821 1
822 -1))
823 (else
824 (if (= self_sign 0)
825 -1
826 1))))))))))
827
828 ;; Note: The Decimal standard doesn't cover rich comparisons for
829 ;; Decimals. In particular, the specification is silent on the
830 ;; subject of what should happen for a comparison involving a NaN.
831 ;; We take the following approach:
832 ;;
833 ;; == comparisons involving a quiet NaN always return False
834 ;; != comparisons involving a quiet NaN always return True
835 ;; == or != comparisons involving a signaling NaN signal
836 ;; InvalidOperation, and return False or True as above if the
837 ;; InvalidOperation is not trapped.
838 ;; <, >, <= and >= comparisons involving a (quiet or signaling)
839 ;; NaN signal InvalidOperation, and return False if the
840 ;; InvalidOperation is not trapped.
841 ;;
842 ;; This behavior is designed to conform as closely as possible to
843 ;; that specified by IEEE 754.
844
845 (define __eq__
846 (lam (self other (= context None))
847 (let ((so (_convert_for_comparisonc self other #:equality_op #t))
848 (self (car so))
849 (other (cadr so)))
850
851 (cond
852 ((eq? other NotImplemented)
853 other)
854 ((bool ((ref self '_check_nans) other context))
855 #f)
856 (else (= ((ref self '_cmp) other) 0))))))
857
858 (define _xlt
859 (lambda (<)
860 (lam (self other (= context None))
861 (let ((so (_convert_for_comparisonc self other #:equality_op #t))
862 (self (car so))
863 (other (cadr so)))
864
865 (cond
866 ((eq? other NotImplemented)
867 other)
868 ((bool ((ref self '_compare_check_nans) other context))
869 #f)
870 (else (< ((ref self '_cmp) other) 0)))))))
871
872 (define __lt__ (_xlt < ))
873 (define __le__ (_xlt <=))
874 (define __gt__ (_xlt > ))
875 (define __ge__ (_xlt >=))
876
877 (define compare
878 (lam (self other (= context None))
879 "Compare self to other. Return a decimal value:
880
881 a or b is a NaN ==> Decimal('NaN')
882 a < b ==> Decimal('-1')
883 a == b ==> Decimal('0')
884 a > b ==> Decimal('1')
885 "
886 (let ((other (_convert_other other #:raiseit #t)))
887 ;; Compare(NaN, NaN) = NaN
888 (if (or (ref self '_is_special)
889 (and (bool other)
890 (ref other '_is_special)))
891 (aif it ((ref self '_check_nans) other context)
892 it
893 (Decimal ((ref self '_cmp) other)))))))
894
895 (define __hash__
896 (lambda (self)
897 "x.__hash__() <==> hash(x)"
898
899 ;; In order to make sure that the hash of a Decimal instance
900 ;; agrees with the hash of a numerically equal integer, float
901 ;; or Fraction, we follow the rules for numeric hashes outlined
902 ;; in the documentation. (See library docs, 'Built-in Types').
903 (cond
904 ((ref self '_is_special)
905 (cond
906 (((ref self 'is_snan))
907 (raise (TypeError "Cannot hash a signaling NaN value.")))
908 (((ref self 'is_snan))
909 (hash (nan)))
910 ((= 1 (ref self '_sign))
911 (hash (- (inf))))
912 (else
913 (hash (inf)))))
914
915 (else
916 (let* ((exp (ref self '_exp))
917 (exp_hash
918 (if (>= exp 0)
919 (expt 10 exp _ pyhash-N)
920 (expt _PyHASH_10INV (- exp) pyhash-N)))
921
922 (hash_
923 (modulus (* (int (ref self '_int)) exp_hash)
924 pyhash-N))
925
926 (ans
927 (if (>= self 0) hash_ (- hash_))))
928 (if (= ans -1) -2 ans))))))
929
930 (define as_tuple
931 (lambda (self)
932 "Represents the number as a triple tuple.
933
934 To show the internals exactly as they are.
935 "
936 (DecimalTuple self._sign
937 (tuple (map int (ref self '_int)))
938 (ref self '_exp))))
939
940 (define as_integer_ratio
941 (lambda (self)
942 "Express a finite Decimal instance in the form n / d.
943
944 Returns a pair (n, d) of integers. When called on an infinity
945 or NaN, raises OverflowError or ValueError respectively.
946
947 >>> Decimal('3.14').as_integer_ratio()
948 (157, 50)
949 >>> Decimal('-123e5').as_integer_ratio()
950 (-12300000, 1)
951 >>> Decimal('0.00').as_integer_ratio()
952 (0, 1)
953 "
954 (if (ref self '_is_special)
955 (if ((ref self 'is_nan))
956 (raise (ValueError
957 "cannot convert NaN to integer ratio"))
958 (raise (OverflowError
959 "cannot convert Infinity to integer ratio"))))
960
961 (if (not (bool self))
962 (values 0 1)
963 (let ((s (ref self '_sign))
964 (n (int (ref self '_int)))
965 (e (ref self '_exp))
966 (x
967 (* n (if (> exp 0)
968 (expt 10 exo)
969 (/ 1 (expt 10 (- expt)))))))
970 (values (numerator x)
971 (denomerator x))))))
972
973 (define __repr__
974 (lambda (self)
975 "Represents the number as an instance of Decimal."
976 ;# Invariant: eval(repr(d)) == d
977 (format #f "Decimal('~a')" (str self))))
978
979 (define __str__
980 (lam (self (= eng #f) (= context None))
981 "Return string representation of the number in scientific notation.
982
983 Captures all of the information in the underlying representation.
984 "
985 (let* ((sign (if (= (reg self '_sign) 0) "" "-"))
986 (exp (ref self '_exp))
987 (i (ref self '_int))
988 (leftdigits (+ exp (len i)))
989 (dotplace #f)
990 (intpart #f)
991 (fracpart #f)
992 (exppart #f))
993
994 (cond
995 ((ref self '_is_special)
996 (cond
997 ((eq? (ref self '_exp) 'F)
998 (+ sign "Infinity"))
999 ((eq? (ref self '_exp) 'n)
1000 (+ sign "NaN" (ref self '_int)))
1001 (else ; self._exp == 'N'
1002 (+ sign "sNaN" (ref self '_int)))))
1003 (else
1004 ;; dotplace is number of digits of self._int to the left of the
1005 ;; decimal point in the mantissa of the output string (that is,
1006 ;; after adjusting the exponent)
1007 (cond
1008 ((and (<= exp 0) (> leftdigits -6))
1009 ;; no exponent required
1010 (set! dotplace leftdigits))
1011
1012 ((not eng)
1013 ;; usual scientific notation: 1 digit on left of the point
1014 (set! dotplace 1))
1015
1016 ((equal? i "0")
1017 ;; engineering notation, zero
1018 (set! dotplace (- (modulo (+ leftdigits 1) 3) 1)))
1019 (else
1020 ;; engineering notation, nonzero
1021 (set! dotplace (- (modulo (+ leftdigits 1) 3) 1))))
1022
1023 (cond
1024 ((<= dotplace 0)
1025 (set! intpart "0")
1026 (set! fracpart (+ "." + (* "0" (- dotplace)) + i)))
1027 ((>= dotplace (len i))
1028 (set! intpart (+ i (* "0" (- dotplace (len i)))))
1029 (set! fracpart ""))
1030 (else
1031 (set! intpart (pylist-slice i None dotplace None))
1032 (set! fracpart (+ '.' (pylist-slice i dotplace None None)))))
1033
1034
1035 (cond
1036 ((= leftdigits dotplace)
1037 (set! exp ""))
1038 (else
1039 (let ((context (if (eq? context None)
1040 (getcontext)
1041 context)))
1042 (set! exp
1043 (+ (pylist-ref (lise "e" "E") (cx-capitals context))
1044 (format #f "%@d" (- leftdigits dotplace)))))))
1045 (+ sign intpart fracpart exp))))))
1046
1047 (define to_eng_string
1048 (lam (self (= context None))
1049 "Convert to a string, using engineering notation if an exponent is needed.
1050 Engineering notation has an exponent which is a multiple of 3. This
1051 can leave up to 3 digits to the left of the decimal place and may
1052 require the addition of either one or two trailing zeros.
1053 "
1054 ((ref self '__str__) #:eng #t #:contect context)))
1055
1056 (define __neg__
1057 (lam (self (= contextNone))
1058 "Returns a copy with the sign switched.
1059
1060 Rounds, if it has reason.
1061 "
1062 (twix
1063 ((un-special self context) it it)
1064 (let* ((context (if (eq? context None)
1065 (getcontext)
1066 context))
1067 (ans (if (and (not (bool self))
1068 (not (eq? (cx-rounding context)
1069 ROUND_FLOOR)))
1070 ;; -Decimal('0') is Decimal('0'),
1071 ;; not Decimal('-0'), except
1072 ;; in ROUND_FLOOR rounding mode.
1073 ((ref self 'copy_abs))
1074 ((ref self 'copy_negate)))))
1075
1076 ((ref ans '_fix) context)))))
1077
1078 (define __pos__
1079 (lam (self (= context None))
1080 "Returns a copy, unless it is a sNaN.
1081
1082 Rounds the number (if more than precision digits)
1083 "
1084 (twix
1085 ((un-special self context) it it)
1086
1087 (let* ((context (if (eq? context None)
1088 (getcontext)
1089 context))
1090 (ans (if (and (not (bool self))
1091 (not (eq? (cx-rounding context)
1092 ROUND_FLOOR)))
1093 ;; -Decimal('0') is Decimal('0'),
1094 ;; not Decimal('-0'), except
1095 ;; in ROUND_FLOOR rounding mode.
1096 ((ref self 'copy_abs))
1097 (Decimal self))))
1098
1099 ((ref ans '_fix) context)))))
1100
1101 (define __abs__
1102 (lam (self (= round #t) (= context None))
1103 "Returns the absolute value of self.
1104
1105 If the keyword argument 'round' is false, do not round. The
1106 expression self.__abs__(round=False) is equivalent to
1107 self.copy_abs().
1108 "
1109 (twix
1110 ((not (bool round))
1111 ((ref self 'copy_abs)))
1112
1113 ((un-special self context) it it)
1114
1115 (if (= (ref self '_sign) 1)
1116 ((ref self '__neg__) #:context context)
1117 ((ref self '__pos__) #:context context)))))
1118
1119 (define __add__
1120 (lam (self other (= context None))
1121 "Returns self + other.
1122
1123 -INF + INF (or the reverse) cause InvalidOperation errors.
1124 "
1125 (twix
1126 ((norm-op other) it it)
1127
1128 (let (get-context context))
1129
1130 ((add-special o1 o2 context) it it)
1131
1132 (let (let* ((negativezero 0)
1133 (self_sign (ref self '_sign))
1134 (other_sign (ref other '_sign))
1135 (self_exp (ref self '_sign))
1136 (other_exp (ref other '_sign))
1137 (prec (cx-prec context))
1138 (exp (min self_exp other_exp))
1139 (sign #f)
1140 (ans #f))
1141
1142 (if (and (eq? (cx-rounding context) ROUND_FLOOR)
1143 (not (= self_sign other_sign)))
1144 ;; If the answer is 0, the sign should be negative,
1145 ;; in this case.
1146 (set! negativezero 1))))
1147
1148 ((if (and (not (bool self)) (not (bool other)))
1149 (begin
1150 (set! sign (min self_sign other_sign))
1151 (if (= negativezero 1)
1152 (set! sign 1))
1153 (set! ans (_dec_from_triple sign "0" exp))
1154 (set! ans ((ref ans '_fix) context))
1155 ans)
1156 #f) it it)
1157
1158 ((if (not (bool self))
1159 (begin
1160 (set! exp (max exp (- other_exp prec 1)))
1161 (set! ans ((ref other '_rescale) exp
1162 (cx-rounding rounding)))
1163 (set! ans ((ref ans '_fix) context))
1164 ans)
1165 #f) it it)
1166
1167 ((if (not (bool other))
1168 (begin
1169 (set! exp (max exp (- self_exp prec 1)))
1170 (set! ans ((ref self '_rescale) exp
1171 (cx-rounding rounding)))
1172 (set! ans ((ref ans '_fix) context))
1173 ans)
1174 #f) it it)
1175
1176
1177 (let (let* ((op1 (_WorkRep self))
1178 (op2 (_WorkRep other))
1179 (ab (_normalize op1 op2 prec))
1180 (op1 (car ab))
1181 (op2 (cadr ab))
1182 (result (_WorkRep)))))
1183
1184 ((cond
1185 ((not (= (ref op1 'sign) (ref op2 'sign)))
1186 ;; Equal and opposite
1187 (twix
1188 ((= op1_i op2_i) it
1189 (set! ans (_dec_from_triple negativezero "0" exp))
1190 (set! ans ((ref ans '_fix) context))
1191 ans)
1192
1193 (begin
1194 (if (< op1_i op2_i)
1195 (let ((t op1))
1196 (set! op1 op2)
1197 (set! op2 t)))
1198
1199 (if (= (ref op1 'sign) 1)
1200 (let ((t (ref op1 'sign)))
1201 (set result 'sign 1)
1202 (set op1 'sign (ref op2 'sign))
1203 (set op2 'sign t))
1204 (set result 'sign 0))
1205 #f)))
1206 ((= (ref op1 'sign) 1)
1207 (set result 'sign 1)
1208 #f)
1209
1210 (begin
1211 (set result 'sign 0)
1212 #f)) it it)
1213
1214 (begin
1215 (if (= (ref op2 'sign) 0)
1216 (set result 'int (+ (ref op1 'int) (ref op2 'int)))
1217 (set result 'int (- (ref op1 'int) (ref op2 'int))))
1218
1219 (set result 'exp (ref op1 'exp))
1220 (set! ans (Decimal result))
1221 ((ref ans '_fix) context)))))
1222
1223 (define __radd__ __add__)
1224
1225 (define __sub__
1226 (lam (self other (= context None))
1227 "Return self - other"
1228 (twix
1229 ((norm-op other) it it)
1230 ((bin-special o1 o2 context) it it)
1231 ((ref self '__add__)
1232 ((ref other 'copy_negate)) #:context context))))
1233
1234 (define __rsub__
1235 (lam (self other (= context None))
1236 "Return other - self"
1237 (twix
1238 ((norm-op other) it it)
1239 ((ref 'other '__sub__) self #:context context))))
1240
1241 (define __mul__
1242 (lam (self other (= context None))
1243 "Return self * other.
1244
1245 (+-) INF * 0 (or its reverse) raise InvalidOperation.
1246 "
1247 (twix
1248 ((norm-op other) it it)
1249 (let (get-context context))
1250
1251 (let (let ((resultsign (logxor (ref self '_sign)
1252 (ref other '_sign))))))
1253
1254 ((mul-special o1 o2 context) it it)
1255
1256 (let (let ((resultexp (+ (ref self '_exp) (ref other '_exp))))))
1257
1258 ;; Special case for multiplying by zero
1259 ((or (not (bool self)) (not (bool other)))
1260 (let ((ans (_dec_from_triple resultsign "0" resultexp)))
1261 ((ref and '_fix) context)))
1262
1263 ;; Special case for multiplying by power of 10
1264 ((equal? (ref self '_int) "1")
1265 (let ((ans (_dec_from_triple resultsign (ref other '_int) resultexp)))
1266 ((ref and '_fix) context)))
1267
1268 ((equal? (ref other '_int) "1")
1269 (let ((ans (_dec_from_triple resultsign (ref self '_int) resultexp)))
1270 ((ref and '_fix) context)))
1271
1272 (let* ((op1 (_WorkRep self))
1273 (op2 (_WorkRep other))
1274 (ans (_dec_from_triple resultsign
1275 (str (* (ref op1 ') (ref op2 'int)))
1276 resultexp)))
1277 ((ref and '_fix) context)))))
1278
1279 (define __rmul__ __mul__)
1280
1281 (define __truediv__
1282 (lam (self other (= context None))
1283 "Return self / other."
1284 (twix
1285 ((norm-op other) it it)
1286 (let (get-context context))
1287
1288 (let (let ((sign (logxor (ref self '_sign)
1289 (ref other '_sign))))))
1290
1291 ((div-special o1 o2 context) it it)
1292
1293 ;; Special cases for zeroes
1294 ((if (not (bool other))
1295 (if (not (bool self))
1296 ((cx-error context) DivisionUndefined "0 / 0")
1297 ((cx-error context) DivisionByZero "x / 0" sign))
1298 #f) it it)
1299
1300 (let ((exp #f)
1301 (coeff #f)
1302 (nself (len (ref self '_int)))
1303 (nother (len (ref other '_int))))
1304 (if (not (bool self))
1305 (begin
1306 (set! exp (- (ref self '_exp) (ref other '_exp)))
1307 (set! coeff 0))
1308 ;; OK, so neither = 0, INF or NaN
1309 (let ((shift (+ nother (- nself) prec 1))
1310 (op1 (_WorkRep self))
1311 (op2 (_WorkRep other)))
1312 (set! exp (- (ref self '_exp) (ref other '_exp) shift))
1313 (call-with-values
1314 (lambda ()
1315 (if (>= shift 0)
1316 (divmod (* (ref op1 'int) (expt 10 shift))
1317 (ref op2 'int))
1318 (divmod (ref op1 'int)
1319 (* (ref op2 'int) (expt 10 shift)))))
1320 (lambda (coeff- remainder)
1321 (set! coeff
1322 (if (not (= remainder 0))
1323 ;; result is not exact adjust to ensure
1324 ;; correct rounding
1325 (if (= (modulus coeff- 5) 0)
1326 (+ coeff- 1)
1327 coeff)
1328 (let (ideal_exp (- (ref self '_exp)
1329 (ref other '_exp)))
1330 (let lp ((coeff- coeff-) (exp- exp))
1331 (if (and (< exp- indeal_exp)
1332 (= (modulo coeff 10) 0))
1333 (lp (/ coeff 10) (+ exp- 1))
1334 (begin
1335 (set exp exp-)
1336 coeff))))))))))
1337
1338
1339 (let ((ans (_dec_from_triple sign, (str coeff) exp)))
1340 ((ref ans '_fix) context))))))
1341
1342 (define _divide
1343 (lambda (self other context)
1344 "Return (self // other, self % other), to context.prec precision.
1345
1346 Assumes that neither self nor other is a NaN, that self is not
1347 infinite and that other is nonzero.
1348 "
1349 (apply values
1350 (twix
1351 (let (let ((sign
1352 (logxor (ref self '_sign)
1353 (ref other '_sign)))
1354 (ideal_exp
1355 (if ((ref other '_isinfinity))
1356 (ref self '_exp)
1357 (min (ref self 'exp) (ref other '_exp))))
1358 (expdiff
1359 (- ((ref self 'adjusted)) ((ref other 'adjusted)))))))
1360
1361 ((or (not (bool self))
1362 ((ref other '_isinfinity))
1363 (<= expdiff -1)) it
1364 (list (_dec_from_tripple sign "0" 0)
1365 ((ref self '_rescale) ideal_exp (cx-rounding context))))
1366
1367 ((if (<= expdiff (cx-prec context))
1368 (let ((op1 (_WorkRep self))
1369 (op2 (_WorkRep other)))
1370 (if (>= (ref op1 'exp) (ref op2 'exp))
1371 (set op1 'int (* (ref op1 'int)
1372 (expt 10 (- (ref op1 'exp)
1373 (ref op2 'exp)))))
1374 (set op2 'int (* (ref op2 'int)
1375 (expt 10 (- (ref op2 'exp)
1376 (ref op1 'exp))))))
1377 (call-with-values (lambda () (divmod (ref op1 'int)
1378 (ref op2 'int)))
1379 (lambda (q r)
1380 (if (< q (expt 10 (cx-prec context)))
1381 (list (_dec_from_triple sign (str q) 0)
1382 (_dec_from_triple (ref self '_sign)
1383 (str r)
1384 ideal_exp))
1385 #f))))
1386 #f) it it)
1387 (begin
1388 ;; Here the quotient is too large to be representable
1389 (let ((ans ((cx-raise context) DivisionImpossible
1390 "quotient too large in //, % or divmod")))
1391 (list ans ans)))))))
1392
1393 (define __rtruediv__
1394 (lam (self other (= context None))
1395 ""Swaps self/other and returns __truediv__.""
1396 (twix
1397 ((norm-op other) it it)
1398 ((ref other '__truediv__) self #:context context))))
1399
1400 (define __divmod__
1401 (lam (self other (= context None))
1402 "
1403 Return (self // other, self % other)
1404 "
1405 (apply values
1406 (twix
1407 ((norm-op other) it it)
1408
1409 (let (get-context context))
1410
1411 ((add-special o1 o2 context) it it)
1412
1413 (((ref self '_check_nans) other context) it
1414 (list it it))
1415
1416 (let (let ((sign
1417 (logxor (ref self '_sign)
1418 (ref other '_sign))))))
1419
1420 (((ref self '_isinfinity)) it
1421 (if ((ref other '_isinfinity))
1422 (let ((ans ((cx-error context) InvalidOperation
1423 "divmod(INF, INF)")))
1424 (list ans ans))
1425 (list (list-ref _SignedInfinity sign)
1426 ((cx-raise context) InvalidOperation, "INF % x"))))
1427
1428 ((not (bool other)) it
1429 (if (not (bool self))
1430 (let ((ans ((cx-error context) DivisionUndefined
1431 "divmod(0, 0)")))
1432 (list ans ans))
1433 (list ((cx-error context) DivisionByZero "x // 0" sign)
1434 ((cx-error context) InvalidOperation "x % 0"))))
1435
1436 (call-with-values (lambda () ((ref self '_divide) other context))
1437 (lambda (quotient remainder)
1438 (let ((remainder ((ref remainder '_fix) context)))
1439 (list quotient remainder))))))))
1440
1441 (define __rdivmod__
1442 (lam (self other (= context None))
1443 "Swaps self/other and returns __divmod__."
1444 (twix
1445 ((norm-op other) it it)
1446 ((ref other '__divmod__) self #:context context))))
1447
1448 (define __mod__
1449 (lam (self other (= context None))
1450 "
1451 self % other
1452 "
1453 (twix
1454 ((norm-op other) it it)
1455
1456 (let (get-context context))
1457
1458 ((bin-special o1 o2 context) it it)
1459
1460 (((ref self '_isinfinity)) it
1461 ((cx-error context) InvalidOperation "INF % x"))
1462
1463 ((not (bool other))
1464 (if (bool self)
1465 ((cx-error context) InvalidOperation "x % 0")
1466 ((cx-error context) DivisionUndefined "0 % 0")))
1467
1468 (let* ((remainder ((ref self '_divide) other context)))
1469 ((ref remainder '_fix) context)))))
1470
1471 (define __rmod__
1472 (lam (self other (= context None))
1473 "Swaps self/other and returns __mod__."
1474 (twix
1475 ((norm-op other) it it)
1476 ((ref other '__mod__) self #:context context))))
1477
1478 (define remainder_near
1479 (lambda (self other (= context None))
1480 "
1481 Remainder nearest to 0- abs(remainder-near) <= other/2
1482 "
1483 (twix
1484 ((norm-op other) it it)
1485
1486 (let (get-context context))
1487
1488 ((bin-special self other context) it it)
1489
1490 ;; self == +/-infinity -> InvalidOperation
1491 (((ref self '_isinfinity)) it
1492 ((cx-error context) InvalidOperation "remainder_near(infinity, x)"))
1493
1494 ;; other == 0 -> either InvalidOperation or DivisionUndefined
1495 ((not (bool other)) it
1496 (if (not (bool self))
1497 ((cx-error context) InvalidOperation "remainder_near(x, 0)")
1498 ((cx-error context) DivisionUndefined "remainder_near(0, 0)")))
1499
1500 ;; other = +/-infinity -> remainder = self
1501 (((ref other '_isinfinity())) it
1502 (let ((ans (Decimal self)))
1503 ((ref ans '_fix) context)))
1504
1505 ;; self = 0 -> remainder = self, with ideal exponent
1506 (let (let ((ideal_exponent (min (ref self '_exp) (ref other '_exp))))))
1507
1508 ((not (bool self)) it
1509 (let ((ans (_dec_from_triple (ref self '_sign) "0" ideal_exponent)))
1510 ((ref ans '_fix) context)))
1511
1512 ;; catch most cases of large or small quotient
1513 (let (let ((expdiff
1514 (- ((ref self 'adjusted)) ((red other 'adjusted)))))))
1515
1516 ((>= expdiff (+ (cx-prec context) 1)) it
1517 ;; expdiff >= prec+1 => abs(self/other) > 10**prec
1518 ((cx-error context) DivisionImpossible))
1519
1520 ((<= expdiff -2) it
1521 ;; expdiff <= -2 => abs(self/other) < 0.1
1522 (let ((ans ((ref self '_rescale)
1523 ideal_exponent (cx-rounding context))))
1524 ((ref ans '_fix) context)))
1525
1526 (let ((op1 (_WorkRep self))
1527 (op2 (_WorkRep other)))
1528
1529 ;; adjust both arguments to have the same exponent, then divide
1530 (if (>= (ref op1 'exp) (ref op2 'exp))
1531 (set op1 'int (* (ref op1 'int)
1532 (expt 10 (- (ref op1 'exp) (ref op2 'exp)))))
1533 (set op2 'int (* (ref op2 'int)
1534 (expt 10 (- (ref op2 'exp) (ref op1 'exp))))))
1535
1536 (call-with-values (lambda () (divmod (ref op1 'int) (ref op2 'int)))
1537 (lambda (q r)
1538
1539 ;; remainder is r*10**ideal_exponent; other is +/-op2.int *
1540 ;; 10**ideal_exponent. Apply correction to ensure that
1541 ;; abs(remainder) <= abs(other)/2
1542 (if (> (+ (* 2 r) + (logand q 1)) (ref op2 'int))
1543 (set! r (- r (ref op2 'int)))
1544 (set! q (+ q 1)))
1545
1546 (if (>= q (expt 10 (cx-prec context)))
1547 ((cx-error context) DivisionImpossible)
1548 (let ((sign (ref self '_sign)))
1549 (if (< r 0)
1550 (set! sign (- 1 sign))
1551 (set! r (- r)))
1552 (let ((ans (_dec_from_triple sign (str r) ideal_exponent)))
1553 ((ref ans '_fix) context))))))))))
1554
1555 (define __floordiv__
1556 (lambda (self other (= context None))
1557 "self // other"
1558 (twix
1559 ((norm-op other) it it)
1560
1561 (let (get-context context))
1562
1563 ((bin-special self other context) it it)
1564
1565 (((ref self '_isinfinity)) it
1566 (if ((ref other '_isinfinity))
1567 ((cx-error context) InvalidOperation "INF // INF")
1568 (pylist-ref _SignedInfinity (logxor (ref self '_sign)
1569 (ref other '_sign)))))
1570
1571 ((not (bool other)) it
1572 (if (bool self)
1573 ((cx-error context) DivisionByZero "x // 0"
1574 (logxor (ref self '_sign) (ref other '_sign)))
1575 ((cx-error context) DivisionUndefined "0 // 0")))
1576
1577 ((ref self '_divide) other context))))
1578
1579 (define __rfloordiv__
1580 (lam (self other (= context None))
1581 "Swaps self/other and returns __floordiv__."
1582 (twix
1583 ((norm-op other) it it)
1584 ((ref other '__floordiv__) self #:context context))))
1585
1586 (define __float__
1587 (lambda (self)
1588 "Float representation."
1589 (if ((ref self '_isnan))
1590 (if ((ref self 'is_snan))
1591 (raise (ValueError "Cannot convert signaling NaN to float"))
1592 (if (= (ref self '_sign))
1593 (- (nan))
1594 (nan)))
1595 (if ((ref self '_isspecial))
1596 (if (= (ref self '_sign))
1597 (- (inf))
1598 (inf)))
1599 (float (str self)))))
1600
1601 (define __int__
1602 (lambda (self)
1603 "Converts self to an int, truncating if necessary."
1604 (if ((ref self '_isnan))
1605 (raise (ValueError "Cannot convert NaN to integer"))
1606 (if ((ref self '_isspecial))
1607 (raise (OverflowError "Cannot convert infinity to integer"))
1608 (let ((s (if (= (ref self '_sign) 1) -1 1)))
1609 (if (>= (ref self '_exp) 0)
1610 (* s (int (ref self '_int)) (expt 10 (ref self '_exp)))
1611 (* s (int (or (bool (py-slice (ref self '_int)
1612 None (ref self '_exp) None))
1613 "0")))))))))
1614
1615 (define __trunc__ __int__)
1616
1617 (define real
1618 (property (lambda (self) self)))
1619
1620 (define imag
1621 (property
1622 (lambda (self)
1623 (Decimal 0))))
1624
1625 (define conjugate
1626 (lambda (self) self))
1627
1628 (define __complex__
1629 (lambda (self)
1630 (complex (float self))))
1631
1632 (define _fix_nan
1633 (lambda (self context)
1634 "Decapitate the payload of a NaN to fit the context"
1635 (let ((payload (ref self '_int))
1636
1637 ;; maximum length of payload is precision if clamp=0,
1638 ;; precision-1 if clamp=1.
1639 (max_payload_len
1640 (- (ref context 'prec)
1641 (ref context 'clamp))))
1642
1643 (if (> (len payload) max_payload_len)
1644 (let ((payload (py-lstrip
1645 (pylist-slice payload
1646 (- (len payload) max_payload_len)
1647 None None) "0")))
1648 (_dec_from_triple (ref self '_sign) payload (ref self '_exp)
1649 #t))
1650 (Decimal self)))))
1651
1652 (define _fix
1653 (lambda (self context)
1654 "Round if it is necessary to keep self within prec precision.
1655
1656 Rounds and fixes the exponent. Does not raise on a sNaN.
1657
1658 Arguments:
1659 self - Decimal instance
1660 context - context used.
1661 "
1662
1663 (twix
1664 (((ref self '_is_special)) it
1665 (if ((ref self '_isnan))
1666 ;; decapitate payload if necessary
1667 ((ref self '_fix_nan) context)
1668
1669 ;; self is +/-Infinity; return unaltered
1670 (Decimal self)))
1671
1672 ;; if self is zero then exponent should be between Etiny and
1673 ;; Emax if clamp==0, and between Etiny and Etop if clamp==1.
1674 (let ((Etiny (cx-etiny context))
1675 (Etop (cx-etop context))))
1676
1677 ((not (bool self)) it
1678 (let ((exp_max (if (= (cx-clamp context) 0)
1679 (cx-emax context)
1680 Etop))
1681 (new_exp (min (max (ref self '_exp) Etiny) exp_max)))
1682 (if (not (= new_exp (ref self '_exp)))
1683 (begin
1684 ((cx-error context) Clamped)
1685 (_dec_from_triple (ref self '_sign) "0" new_exp))
1686 (Decimal self))))
1687
1688 ;; exp_min is the smallest allowable exponent of the result,
1689 ;; equal to max(self.adjusted()-context.prec+1, Etiny)
1690 (let ((exp_min (+ (len (ref self '_int))
1691 (ref self '_exp)
1692 (- (cx-prec context)))))))
1693 ((> exp_min Etop) it
1694 ;; overflow: exp_min > Etop iff self.adjusted() > Emax
1695 (let ((ans ((cx-error context) Overflow "above Emax"
1696 (ref self '_sign))))
1697 ((cx-error context) Inexact)
1698 ((cx-error context) Rounded)
1699 ans))
1700
1701 (let* ((self_is_subnormal (< exp_min Etiny))
1702 (exp_min (if self_is_subnormal Eriny exp_min)))))
1703
1704 ;; round if self has too many digits
1705 ((< self._exp exp_min) it
1706 (let ((digits (+ (len (ref self '_int))
1707 (ref self '_exp)
1708 (- exp_min))))
1709 (if (< digits 0)
1710 (set! self (_dec_from_triple (ref self '_sign)
1711 "1" (- exp_min 1)))
1712 (set! digits 0))
1713
1714 (let* ((ans #f)
1715 (rounding_method (pylist-ref
1716 (ref self '_pick_rounding_function)
1717 (cx-rounding context)))
1718 (changed (rounding_method self digits))
1719 (coeff (or (bool (pylist-slice (ref self '_int)
1720 None digits None)) "0")))
1721 (if (> changed 0)
1722 (begin
1723 (set! coeff (str (+ (int coeff) 1)))
1724 (if (> (len coeff) (cx-prec context))
1725 (begin
1726 (set! coeff (pylist-clice coeff None -1 None))
1727 (set! exp_min (+ exp_min 1))))))
1728
1729 ;; check whether the rounding pushed the exponent out of range
1730 (if (> exp_min Etop)
1731 (set! ans
1732 ((cx-error context) Overflow "above Emax"
1733 (ref self '_sign)))
1734 (set! ans (_dec_from_triple (ref self '_sign) coeff exp_min)))
1735
1736 ;; raise the appropriate signals, taking care to respect
1737 ;; the precedence described in the specification
1738 (if (and changed self_is_subnormal)
1739 ((cx-error context) Underflow))
1740 (if self_is_subnormal
1741 ((cx-error context) Subnormal))
1742 (if changed
1743 ((cx-error context) Inexact))
1744
1745 ((cx-error context) Rounded)
1746
1747 (if (not (bool ans))
1748 ;; raise Clamped on underflow to 0
1749 ((cx-error context) Clamped))
1750
1751 ans)))
1752 (begin
1753 (if self_is_subnormal
1754 ((cx-error context) Subnormal))
1755
1756
1757 ;; fold down if clamp == 1 and self has too few digits
1758 (if (and (= (cx-clamp context) 1) (> (ref self '_exp) Etop))
1759 (begin
1760 ((cx-error context) Clamped)
1761 (let ((self_padded (+ (ref self '_int)
1762 (* "0"
1763 (- (ref self '_exp) Etop)))))
1764 (_dec_from_triple (ref self '_sign) self_padded Etop)))
1765
1766 ;; here self was representable to begin with; return unchanged
1767 (Decimal self))))
1768
1769
1770 ;; for each of the rounding functions below:
1771 ;; self is a finite, nonzero Decimal
1772 ;; prec is an integer satisfying 0 <= prec < len(self._int)
1773 ;;
1774 ;; each function returns either -1, 0, or 1, as follows:
1775 ;; 1 indicates that self should be rounded up (away from zero)
1776 ;; 0 indicates that self should be truncated, and that all the
1777 ;; digits to be truncated are zeros (so the value is unchanged)
1778 ;; -1 indicates that there are nonzero digits to be truncated
1779
1780 (define _round_down
1781 (lambda (self prec)
1782 "Also known as round-towards-0, truncate."
1783 (if (_all_zeros (ref self '_int) prec)
1784 0
1785 -1)))
1786
1787 (define _round_up
1788 (lambda (self prec)
1789 "Rounds away from 0."
1790 (- (_round_down self prec))))
1791
1792 (define _round_half_up
1793 (lambda (self prec)
1794 "Rounds 5 up (away from 0)"
1795 (cond
1796 ((in (pylist-ref (ref self '_int) prec) "56789")
1797 1)
1798 ((_all_zeros (ref self '_int) prec)
1799 0)
1800 (else -1))))
1801
1802 (define _round_half_down
1803 (lambda (self prec)
1804 "Round 5 down"
1805 (if (_exact_half (ref self '_int) prec)
1806 -1
1807 (_round_half_up self prec))))
1808
1809 (define _round_half_even
1810 (lambda (self prec)
1811 "Round 5 to even, rest to nearest."
1812 (if (and (_exact_half (ref self '_int) prec)
1813 (or (= prec 0)
1814 (in (pylist-ref (ref self '_int) (- prec 1)) "02468")))
1815 -1)
1816 (_round_half_up self prec)))
1817
1818 (define _round_ceiling
1819 (lambda (self prec)
1820 "Rounds up (not away from 0 if negative.)"
1821 (if (= (ref self '_sign) 1)
1822 (_round_down self prec)
1823 (- (_round_down self prec)))))
1824
1825 (define _round_floor
1826 (lambda (self prec)
1827 "Rounds down (not towards 0 if negative)"
1828 (if (= (ref self '_sign) 1)
1829 (- (_round_down self prec))
1830 (_round_down self prec))))
1831
1832 (define _round_05up
1833 (lambda (self prec)
1834 "Round down unless digit prec-1 is 0 or 5."
1835 (if (and prec (not (in (pylist-ref (ref self '_int) (- prec 1) "05"))))
1836 (_round_down self prec)
1837 (- (_round_down self prec)))))
1838
1839 (define _pick_rounding_function
1840 (dict `((,ROUND_DOWN . ,_round_down )
1841 (,ROUND_UP . ,_round_up )
1842 (,ROUND_HALF_UP . ,_round_half_up)
1843 (,ROUND_HALF_DOWN . ,_round_half_down)
1844 (,ROUND_HALF_EVEN . ,_round_half_even)
1845 (,ROUND_CEILING . ,_round_ceiling)
1846 (,ROUND_FLOOR . ,_round_floor)
1847 (,ROUND_05UP . ,_round_05up))))
1848
1849 (define __round__
1850 (lam (self (= n None))
1851 "Round self to the nearest integer, or to a given precision.
1852
1853 If only one argument is supplied, round a finite Decimal
1854 instance self to the nearest integer. If self is infinite or
1855 a NaN then a Python exception is raised. If self is finite
1856 and lies exactly halfway between two integers then it is
1857 rounded to the integer with even last digit.
1858
1859 >>> round(Decimal('123.456'))
1860 123
1861 >>> round(Decimal('-456.789'))
1862 -457
1863 >>> round(Decimal('-3.0'))
1864 -3
1865 >>> round(Decimal('2.5'))
1866 2
1867 >>> round(Decimal('3.5'))
1868 4
1869 >>> round(Decimal('Inf'))
1870 Traceback (most recent call last):
1871 ...
1872 OverflowError: cannot round an infinity
1873 >>> round(Decimal('NaN'))
1874 Traceback (most recent call last):
1875 ...
1876 ValueError: cannot round a NaN
1877
1878 If a second argument n is supplied, self is rounded to n
1879 decimal places using the rounding mode for the current
1880 context.
1881
1882 For an integer n, round(self, -n) is exactly equivalent to
1883 self.quantize(Decimal('1En')).
1884
1885 >>> round(Decimal('123.456'), 0)
1886 Decimal('123')
1887 >>> round(Decimal('123.456'), 2)
1888 Decimal('123.46')
1889 >>> round(Decimal('123.456'), -2)
1890 Decimal('1E+2')
1891 >>> round(Decimal('-Infinity'), 37)
1892 Decimal('NaN')
1893 >>> round(Decimal('sNaN123'), 0)
1894 Decimal('NaN123')
1895
1896 "
1897 (if (not (eq? n None))
1898 ;; two-argument form: use the equivalent quantize call
1899 (if (not (isinstance n int))
1900 (raise (TypeError
1901 "Second argument to round should be integral"))
1902 (let ((exp (_dec_from_triple 0, "1", (- n))))
1903 ((ref self 'quantize) exp)))
1904
1905 ;; one-argument form
1906 (if (ref self '_is_special)
1907 (if ((ref self 'is_nan))
1908 (raise (ValueError "cannot round a NaN"))
1909 (raise (OverflowError "cannot round an infinity")))
1910 (int ((ref self '_rescale) 0 ROUND_HALF_EVEN))))))
1911
1912 (define __floor__
1913 (lambda (self)
1914 "Return the floor of self, as an integer.
1915
1916 For a finite Decimal instance self, return the greatest
1917 integer n such that n <= self. If self is infinite or a NaN
1918 then a Python exception is raised.
1919
1920 "
1921 (if (ref self '_is_special)
1922 (if ((ref self 'is_nan))
1923 (raise (ValueError "cannot round a NaN"))
1924 (raise (OverflowError "cannot round an infinity")))
1925 (int ((ref self '_rescale) 0 ROUND_FLOOR)))))
1926
1927 (define __ceil__
1928 (lambda (self)
1929 """Return the ceiling of self, as an integer.
1930
1931 For a finite Decimal instance self, return the least integer n
1932 such that n >= self. If self is infinite or a NaN then a
1933 Python exception is raised.
1934
1935 """
1936 (if (ref self '_is_special)
1937 (if ((ref self 'is_nan))
1938 (raise (ValueError "cannot round a NaN"))
1939 (raise (OverflowError "cannot round an infinity")))
1940 (int ((ref self '_rescale) 0 ROUND_CEILING)))))
1941
1942 (define fma
1943 (lam (self other third (= context None))
1944 "Fused multiply-add.
1945
1946 Returns self*other+third with no rounding of the intermediate
1947 product self*other.
1948
1949 self and other are multiplied together, with no rounding of
1950 the result. The third operand is then added to the result,
1951 and a single final rounding is performed.
1952 "
1953 (twix
1954 (let ((other (_convert_other other #:raiseit #t))
1955 (third (_convert_other third #:raiseit #t))
1956 (fin (lambda (product)
1957 ((ref product '__add__) third context)))))
1958 ;; compute product; raise InvalidOperation if either operand is
1959 ;; a signaling NaN or if the product is zero times infinity.
1960 ((if (or (ref self '_is_special) (ref other '_is_special))
1961 (twix
1962 (let (get-context context))
1963 ((equals? (ref self '_exp) "N") it
1964 ((cx-error context) InvalidOperation "sNaN" self))
1965 ((equals? (ref other '_exp) "N") it
1966 ((cx-error context) InvalidOperation "sNaN" other))
1967 ((equals? (ref self '_exp) "n") it
1968 (fin self))
1969 ((equals? (ref other '_exp) "n") it
1970 (fin other))
1971 ((equals? (ref self '_exp) "F") it
1972 (if (not (bool other))
1973 ((cx-error context) InvalidOperation "INF * 0 in fma")
1974 (pylist-ref _SignedInfinity
1975 (logxor (ref self '_sign)
1976 (ref other '_sign)))))
1977 ((equals? (ref other '_exp) "F") it
1978 (if (not (bool self))
1979 ((cx-error context) InvalidOperation "0 * INF in fma")
1980 (pylist-ref _SignedInfinity
1981 (logxor (ref self '_sign)
1982 (ref other '_sign)))))
1983 #f)) it it)
1984
1985 (fin
1986 (_dec_from_triple (logxor (ref self '_sign) (ref other '_sign))
1987 (str (* (int (ref self '_int))
1988 (int (ref other '_int))))
1989 (+ (ref self '_exp) (ref other '_exp)))))))
1990
1991 (define _power_modulo
1992 (lam (self other modulo (= context None))
1993 "Three argument version of __pow__"
1994 (twix
1995 ((norm-op other ) it it)
1996 ((norm-op modulo) it it)
1997 (let (get-context context))
1998
1999 ;; deal with NaNs: if there are any sNaNs then first one wins,
2000 ;; (i.e. behaviour for NaNs is identical to that of fma)
2001 (let ((self_is_nan (ref self '_isnan))
2002 (other_is_nan (ref other '_isnan))
2003 (modulo_is_nan (ref modulo '_isnan))))
2004
2005 ((or (bool self_is_nan) (bool other_is_nan) (bool modulo_is_nan)) it
2006 (cond
2007 ((= self_is_nan 2)
2008 ((cx-error context) InvalidOperation, "sNaN" self))
2009 ((= other_is_nan 2)
2010 ((cx-error context) InvalidOperation, "sNaN" other))
2011 ((modulo_is_nan 2)
2012 ((cx-error context) InvalidOperation, "sNaN" modulo))
2013 ((bool self_is_nan)
2014 (_fix_nan self context))
2015 ((bool other_is_nan)
2016 (_fix_nan other context))
2017 (else
2018 (_fix_nan modulo context))))
2019
2020 ;;check inputs: we apply same restrictions as Python's pow()
2021 ((not (and ((ref self '_isinteger))
2022 ((ref other '_isinteger))
2023 ((ref modulo '_isinteger)))) it
2024 ((cx-error context) InvalidOperation
2025 (+ "pow() 3rd argument not allowed "
2026 "unless all arguments are integers")))
2027
2028 ((< other 0) it
2029 ((cx-error context) InvalidOperation
2030 (+ "pow() 2nd argument cannot be "
2031 "negative when 3rd argument specified")))
2032
2033 ((not (bool modulo)) it
2034 ((cx-error context) InvalidOperation
2035 "pow() 3rd argument cannot be 0"))
2036
2037 ;; additional restriction for decimal: the modulus must be less
2038 ;; than 10**prec in absolute value
2039 ((>= ((ref modulo 'adjusted)) (cx-prec context)) it
2040 ((cx-error context) InvalidOperation
2041 (+ "insufficient precision: pow() 3rd "
2042 "argument must not have more than "
2043 "precision digits")))
2044
2045 ;; define 0**0 == NaN, for consistency with two-argument pow
2046 ;; (even though it hurts!)
2047 ((and (not (bool other)) (not (bool self)))
2048 ((cx-error context) InvalidOperation
2049 (+ "at least one of pow() 1st argument "
2050 "and 2nd argument must be nonzero ;"
2051 "0**0 is not defined")))
2052
2053 ;; compute sign of result
2054 (let ((sign (if ((ref other '_iseven))
2055 0
2056 (ref self '_sign)))
2057 (base (_WorkRep ((ref self 'to_integral_value))))
2058 (exponent (_WorkRep ((ref other 'to_integral_value)))))
2059
2060
2061 ;; convert modulo to a Python integer, and self and other to
2062 ;; Decimal integers (i.e. force their exponents to be >= 0)
2063 (set! modulo (abs (int modulo)))
2064
2065 ;; compute result using integer pow()
2066 (set! base (guile:modulo
2067 (* (guile:modulo (ref base 'int) modulo)
2068 (modulo-expt 10 (ref base 'exp) modulo))
2069 modulo))
2070
2071 (let lp ((i (ref exponent 'exp)))
2072 (if (> i 0)
2073 (begin
2074 (set! base (modulo-expt base 10 modulo))
2075 (lp (- i 1)))))
2076
2077 (set! base (modulo-expt base (ref exponent 'int) modulo))
2078
2079 (_dec_from_triple sign (str base) 0)))))
2080
2081 (define _power_exact
2082 (lambda (self other p)
2083 "Attempt to compute self**other exactly.
2084
2085 Given Decimals self and other and an integer p, attempt to
2086 compute an exact result for the power self**other, with p
2087 digits of precision. Return None if self**other is not
2088 exactly representable in p digits.
2089
2090 Assumes that elimination of special cases has already been
2091 performed: self and other must both be nonspecial; self must
2092 be positive and not numerically equal to 1; other must be
2093 nonzero. For efficiency, other._exp should not be too large,
2094 so that 10**abs(other._exp) is a feasible calculation."
2095
2096 ;; In the comments below, we write x for the value of self and y for the
2097 ;; value of other. Write x = xc*10**xe and abs(y) = yc*10**ye, with xc
2098 ;; and yc positive integers not divisible by 10.
2099
2100 ;; The main purpose of this method is to identify the *failure*
2101 ;; of x**y to be exactly representable with as little effort as
2102 ;; possible. So we look for cheap and easy tests that
2103 ;; eliminate the possibility of x**y being exact. Only if all
2104 ;; these tests are passed do we go on to actually compute x**y.
2105
2106 ;; Here's the main idea. Express y as a rational number m/n, with m and
2107 ;; n relatively prime and n>0. Then for x**y to be exactly
2108 ;; representable (at *any* precision), xc must be the nth power of a
2109 ;; positive integer and xe must be divisible by n. If y is negative
2110 ;; then additionally xc must be a power of either 2 or 5, hence a power
2111 ;; of 2**n or 5**n.
2112 ;;
2113 ;; There's a limit to how small |y| can be: if y=m/n as above
2114 ;; then:
2115 ;;
2116 ;; (1) if xc != 1 then for the result to be representable we
2117 ;; need xc**(1/n) >= 2, and hence also xc**|y| >= 2. So
2118 ;; if |y| <= 1/nbits(xc) then xc < 2**nbits(xc) <=
2119 ;; 2**(1/|y|), hence xc**|y| < 2 and the result is not
2120 ;; representable.
2121 ;;
2122 ;; (2) if xe != 0, |xe|*(1/n) >= 1, so |xe|*|y| >= 1. Hence if
2123 ;; |y| < 1/|xe| then the result is not representable.
2124 ;;
2125 ;; Note that since x is not equal to 1, at least one of (1) and
2126 ;; (2) must apply. Now |y| < 1/nbits(xc) iff |yc|*nbits(xc) <
2127 ;; 10**-ye iff len(str(|yc|*nbits(xc)) <= -ye.
2128 ;;
2129 ;; There's also a limit to how large y can be, at least if it's
2130 ;; positive: the normalized result will have coefficient xc**y,
2131 ;; so if it's representable then xc**y < 10**p, and y <
2132 ;; p/log10(xc). Hence if y*log10(xc) >= p then the result is
2133 ;; not exactly representable.
2134
2135 ;; if len(str(abs(yc*xe)) <= -ye then abs(yc*xe) < 10**-ye,
2136 ;; so |y| < 1/xe and the result is not representable.
2137 ;; Similarly, len(str(abs(yc)*xc_bits)) <= -ye implies |y|
2138 ;; < 1/nbits(xc).
2139
2140 (twix
2141 (let ()
2142 (define-syntax-rule (clean xc xe n +)
2143 (let lp ()
2144 (if (= (modulo xc n) 0)
2145 (begin
2146 (set! xc (/ xc n))
2147 (set! xe (+ xe 1))
2148 (lp)))))
2149
2150 (let* ((x (_WorkRep self))
2151 (xc (ref x 'int))
2152 (xe (ref x 'exp)))
2153 (clean xc xe 10 +))
2154
2155 (let* ((y (_WorkRep other))
2156 (yc (ref y 'int))
2157 (ye (ref y 'exp)))
2158 (clean yc ye 10 +))
2159
2160 ;; case where xc == 1: result is 10**(xe*y), with xe*y
2161 ;; required to be an integer
2162 ((= xc 1) it
2163 (set! xe (* xe yc))
2164
2165 ;; result is now 10**(xe * 10**ye); xe * 10**ye must be integral
2166 (clean xe ye 10 +)
2167
2168 (if (< ye 0)
2169 None
2170 (let ((exponent (* xe (expt 10 ye)))
2171 (zeros #f))
2172 (if (= (ref y 'sign) 1)
2173 (set! exponent (- exponent)))
2174
2175 ;;if other is a nonnegative integer, use ideal exponent
2176 (if (and ((ref other '_isinteger)) (= (ref other '_sign) 0))
2177 (begin
2178 (let ((ideal_exponent (* (ref self '_exp) (int other))))
2179 (set! zeros (min (- exponent ideal_exponent) (- p 1)))))
2180 (set! zeros 0))
2181
2182 (_dec_from_triple 0 (+ "1" (* "0" zeros)) exponent-zeros))))
2183
2184 ;; case where y is negative: xc must be either a power
2185 ;; of 2 or a power of 5.
2186 ((= (ref y 'sign) 1) it
2187 (let ((last_digit (modulo xc 10)))
2188 (twix
2189 ((cond
2190 ((= (modulo last_digit 2) 0)
2191 ;; quick test for power of 2
2192 (twix
2193 ((not (= (logand xc (- xc)) xc))
2194 None)
2195 ;; now xc is a power of 2; e is its exponent
2196 (let ((e (- (_nbits xc) 1))))
2197
2198 ;; We now have:
2199 ;;
2200 ;; x = 2**e * 10**xe, e > 0, and y < 0.
2201 ;;
2202 ;; The exact result is:
2203 ;;
2204 ;; x**y = 5**(-e*y) * 10**(e*y + xe*y)
2205 ;;
2206 ;; provided that both e*y and xe*y are integers.
2207 ;; Note that if
2208 ;; 5**(-e*y) >= 10**p, then the result can't be expressed
2209 ;; exactly with p digits of precision.
2210 ;;
2211 ;; Using the above, we can guard against large values of ye.
2212 ;; 93/65 is an upper bound for log(10)/log(5), so if
2213 ;;
2214 ;; ye >= len(str(93*p//65))
2215 ;;
2216 ;; then
2217 ;;
2218 ;; -e*y >= -y >= 10**ye > 93*p/65 > p*log(10)/log(5),
2219 ;;
2220 ;; so 5**(-e*y) >= 10**p, and the coefficient of the result
2221 ;; can't be expressed in p digits.
2222
2223 ;; emax >= largest e such that 5**e < 10**p.
2224 (let ((emax (quotient (* p 93) 65))))
2225
2226 ((>= ye (len (str emax))) it
2227 None)
2228
2229 ;; Find -e*y and -xe*y; both must be integers
2230 (let ()
2231 (set! e (_decimal_lshift_exact (* e yc) ye))
2232 (set! xe (_decimal_lshift_exact (* xe yc) ye)))
2233
2234 ((or (eq? e None) (eq? xe None)) it
2235 None)
2236
2237 ((> e emax)
2238 None)
2239
2240 (begin
2241 (set! xc (expt 5 e))
2242 #f)))
2243
2244 ((= last_digit 5)
2245 (twix
2246 ;; e >= log_5(xc) if xc is a power of 5; we have
2247 ;; equality all the way up to xc=5**2658
2248 (let* ((e (quotient (* (_nbits xc) 28) 65))
2249 (q (expt 5 e))
2250 (xc (quotient q xz))
2251 (remainder (modulo q xc))))
2252
2253 ((not (= remainder 0)) it
2254 None)
2255
2256 (let () (clean xc e 5 -))
2257
2258 ;; Guard against large values of ye, using the same logic as in
2259 ;; the 'xc is a power of 2' branch. 10/3 is an upper bound for
2260 ;; log(10)/log(2).
2261 (let ((emax (quotient (* p 10) 3))))
2262
2263 ((>= ye (len (str emax)))
2264 None)
2265
2266 (let ()
2267 (set! e (_decimal_lshift_exact (* e yc) ye))
2268 (set! xe (_decimal_lshift_exact (* xe yc) ye)))
2269
2270 ((or (eq? e None= (eq? xe None))) it
2271 None)
2272
2273 ((> e emax)
2274 None)
2275
2276 (begin
2277 (set! xc (expt 2 e))
2278 #f)))
2279
2280 (else
2281 None)) it it)
2282
2283 ((>= xc (expt 10 p)) it it)
2284
2285 (begin
2286 (set! xe (+ (- e) (- xe)))
2287 return _dec_from_triple(0, str(xc), xe)
2288
2289 # now y is positive; find m and n such that y = m/n
2290 if ye >= 0:
2291 m, n = yc*10**ye, 1
2292 else:
2293 if xe != 0 and len(str(abs(yc*xe))) <= -ye:
2294 return None
2295 xc_bits = _nbits(xc)
2296 if xc != 1 and len(str(abs(yc)*xc_bits)) <= -ye:
2297 return None
2298 m, n = yc, 10**(-ye)
2299 while m % 2 == n % 2 == 0:
2300 m //= 2
2301 n //= 2
2302 while m % 5 == n % 5 == 0:
2303 m //= 5
2304 n //= 5
2305
2306 # compute nth root of xc*10**xe
2307 if n > 1:
2308 # if 1 < xc < 2**n then xc isn't an nth power
2309 if xc != 1 and xc_bits <= n:
2310 return None
2311
2312 xe, rem = divmod(xe, n)
2313 if rem != 0:
2314 return None
2315
2316 # compute nth root of xc using Newton's method
2317 a = 1 << -(-_nbits(xc)//n) # initial estimate
2318 while True:
2319 q, r = divmod(xc, a**(n-1))
2320 if a <= q:
2321 break
2322 else:
2323 a = (a*(n-1) + q)//n
2324 if not (a == q and r == 0):
2325 return None
2326 xc = a
2327
2328 # now xc*10**xe is the nth root of the original xc*10**xe
2329 # compute mth power of xc*10**xe
2330
2331 # if m > p*100//_log10_lb(xc) then m > p/log10(xc), hence xc**m >
2332 # 10**p and the result is not representable.
2333 if xc > 1 and m > p*100//_log10_lb(xc):
2334 return None
2335 xc = xc**m
2336 xe *= m
2337 if xc > 10**p:
2338 return None
2339
2340 # by this point the result *is* exactly representable
2341 # adjust the exponent to get as close as possible to the ideal
2342 # exponent, if necessary
2343 str_xc = str(xc)
2344 if other._isinteger() and other._sign == 0:
2345 ideal_exponent = self._exp*int(other)
2346 zeros = min(xe-ideal_exponent, p-len(str_xc))
2347 else:
2348 zeros = 0
2349 return _dec_from_triple(0, str_xc+'0'*zeros, xe-zeros)
2350
2351 def __pow__(self, other, modulo=None, context=None):
2352 """Return self ** other [ % modulo].
2353
2354 With two arguments, compute self**other.
2355
2356 With three arguments, compute (self**other) % modulo. For the
2357 three argument form, the following restrictions on the
2358 arguments hold:
2359
2360 - all three arguments must be integral
2361 - other must be nonnegative
2362 - either self or other (or both) must be nonzero
2363 - modulo must be nonzero and must have at most p digits,
2364 where p is the context precision.
2365
2366 If any of these restrictions is violated the InvalidOperation
2367 flag is raised.
2368
2369 The result of pow(self, other, modulo) is identical to the
2370 result that would be obtained by computing (self**other) %
2371 modulo with unbounded precision, but is computed more
2372 efficiently. It is always exact.
2373 """
2374
2375 if modulo is not None:
2376 return self._power_modulo(other, modulo, context)
2377
2378 other = _convert_other(other)
2379 if other is NotImplemented:
2380 return other
2381
2382 if context is None:
2383 context = getcontext()
2384
2385 # either argument is a NaN => result is NaN
2386 ans = self._check_nans(other, context)
2387 if ans:
2388 return ans
2389
2390 # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
2391 if not other:
2392 if not self:
2393 return context._raise_error(InvalidOperation, '0 ** 0')
2394 else:
2395 return _One
2396
2397 # result has sign 1 iff self._sign is 1 and other is an odd integer
2398 result_sign = 0
2399 if self._sign == 1:
2400 if other._isinteger():
2401 if not other._iseven():
2402 result_sign = 1
2403 else:
2404 # -ve**noninteger = NaN
2405 # (-0)**noninteger = 0**noninteger
2406 if self:
2407 return context._raise_error(InvalidOperation,
2408 'x ** y with x negative and y not an integer')
2409 # negate self, without doing any unwanted rounding
2410 self = self.copy_negate()
2411
2412 # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
2413 if not self:
2414 if other._sign == 0:
2415 return _dec_from_triple(result_sign, '0', 0)
2416 else:
2417 return _SignedInfinity[result_sign]
2418
2419 # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
2420 if self._isinfinity():
2421 if other._sign == 0:
2422 return _SignedInfinity[result_sign]
2423 else:
2424 return _dec_from_triple(result_sign, '0', 0)
2425
2426 # 1**other = 1, but the choice of exponent and the flags
2427 # depend on the exponent of self, and on whether other is a
2428 # positive integer, a negative integer, or neither
2429 if self == _One:
2430 if other._isinteger():
2431 # exp = max(self._exp*max(int(other), 0),
2432 # 1-context.prec) but evaluating int(other) directly
2433 # is dangerous until we know other is small (other
2434 # could be 1e999999999)
2435 if other._sign == 1:
2436 multiplier = 0
2437 elif other > context.prec:
2438 multiplier = context.prec
2439 else:
2440 multiplier = int(other)
2441
2442 exp = self._exp * multiplier
2443 if exp < 1-context.prec:
2444 exp = 1-context.prec
2445 context._raise_error(Rounded)
2446 else:
2447 context._raise_error(Inexact)
2448 context._raise_error(Rounded)
2449 exp = 1-context.prec
2450
2451 return _dec_from_triple(result_sign, '1'+'0'*-exp, exp)
2452
2453 # compute adjusted exponent of self
2454 self_adj = self.adjusted()
2455
2456 # self ** infinity is infinity if self > 1, 0 if self < 1
2457 # self ** -infinity is infinity if self < 1, 0 if self > 1
2458 if other._isinfinity():
2459 if (other._sign == 0) == (self_adj < 0):
2460 return _dec_from_triple(result_sign, '0', 0)
2461 else:
2462 return _SignedInfinity[result_sign]
2463
2464 # from here on, the result always goes through the call
2465 # to _fix at the end of this function.
2466 ans = None
2467 exact = False
2468
2469 # crude test to catch cases of extreme overflow/underflow. If
2470 # log10(self)*other >= 10**bound and bound >= len(str(Emax))
2471 # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
2472 # self**other >= 10**(Emax+1), so overflow occurs. The test
2473 # for underflow is similar.
2474 bound = self._log10_exp_bound() + other.adjusted()
2475 if (self_adj >= 0) == (other._sign == 0):
2476 # self > 1 and other +ve, or self < 1 and other -ve
2477 # possibility of overflow
2478 if bound >= len(str(context.Emax)):
2479 ans = _dec_from_triple(result_sign, '1', context.Emax+1)
2480 else:
2481 # self > 1 and other -ve, or self < 1 and other +ve
2482 # possibility of underflow to 0
2483 Etiny = context.Etiny()
2484 if bound >= len(str(-Etiny)):
2485 ans = _dec_from_triple(result_sign, '1', Etiny-1)
2486
2487 # try for an exact result with precision +1
2488 if ans is None:
2489 ans = self._power_exact(other, context.prec + 1)
2490 if ans is not None:
2491 if result_sign == 1:
2492 ans = _dec_from_triple(1, ans._int, ans._exp)
2493 exact = True
2494
2495 # usual case: inexact result, x**y computed directly as exp(y*log(x))
2496 if ans is None:
2497 p = context.prec
2498 x = _WorkRep(self)
2499 xc, xe = x.int, x.exp
2500 y = _WorkRep(other)
2501 yc, ye = y.int, y.exp
2502 if y.sign == 1:
2503 yc = -yc
2504
2505 # compute correctly rounded result: start with precision +3,
2506 # then increase precision until result is unambiguously roundable
2507 extra = 3
2508 while True:
2509 coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
2510 if coeff % (5*10**(len(str(coeff))-p-1)):
2511 break
2512 extra += 3
2513
2514 ans = _dec_from_triple(result_sign, str(coeff), exp)
2515
2516 # unlike exp, ln and log10, the power function respects the
2517 # rounding mode; no need to switch to ROUND_HALF_EVEN here
2518
2519 # There's a difficulty here when 'other' is not an integer and
2520 # the result is exact. In this case, the specification
2521 # requires that the Inexact flag be raised (in spite of
2522 # exactness), but since the result is exact _fix won't do this
2523 # for us. (Correspondingly, the Underflow signal should also
2524 # be raised for subnormal results.) We can't directly raise
2525 # these signals either before or after calling _fix, since
2526 # that would violate the precedence for signals. So we wrap
2527 # the ._fix call in a temporary context, and reraise
2528 # afterwards.
2529 if exact and not other._isinteger():
2530 # pad with zeros up to length context.prec+1 if necessary; this
2531 # ensures that the Rounded signal will be raised.
2532 if len(ans._int) <= context.prec:
2533 expdiff = context.prec + 1 - len(ans._int)
2534 ans = _dec_from_triple(ans._sign, ans._int+'0'*expdiff,
2535 ans._exp-expdiff)
2536
2537 # create a copy of the current context, with cleared flags/traps
2538 newcontext = context.copy()
2539 newcontext.clear_flags()
2540 for exception in _signals:
2541 newcontext.traps[exception] = 0
2542
2543 # round in the new context
2544 ans = ans._fix(newcontext)
2545
2546 # raise Inexact, and if necessary, Underflow
2547 newcontext._raise_error(Inexact)
2548 if newcontext.flags[Subnormal]:
2549 newcontext._raise_error(Underflow)
2550
2551 # propagate signals to the original context; _fix could
2552 # have raised any of Overflow, Underflow, Subnormal,
2553 # Inexact, Rounded, Clamped. Overflow needs the correct
2554 # arguments. Note that the order of the exceptions is
2555 # important here.
2556 if newcontext.flags[Overflow]:
2557 context._raise_error(Overflow, 'above Emax', ans._sign)
2558 for exception in Underflow, Subnormal, Inexact, Rounded, Clamped:
2559 if newcontext.flags[exception]:
2560 context._raise_error(exception)
2561
2562 else:
2563 ans = ans._fix(context)
2564
2565 return ans
2566
2567 def __rpow__(self, other, context=None):
2568 """Swaps self/other and returns __pow__."""
2569 other = _convert_other(other)
2570 if other is NotImplemented:
2571 return other
2572 return other.__pow__(self, context=context)
2573
2574 def normalize(self, context=None):
2575 """Normalize- strip trailing 0s, change anything equal to 0 to 0e0"""
2576
2577 if context is None:
2578 context = getcontext()
2579
2580 if self._is_special:
2581 ans = self._check_nans(context=context)
2582 if ans:
2583 return ans
2584
2585 dup = self._fix(context)
2586 if dup._isinfinity():
2587 return dup
2588
2589 if not dup:
2590 return _dec_from_triple(dup._sign, '0', 0)
2591 exp_max = [context.Emax, context.Etop()][context.clamp]
2592 end = len(dup._int)
2593 exp = dup._exp
2594 while dup._int[end-1] == '0' and exp < exp_max:
2595 exp += 1
2596 end -= 1
2597 return _dec_from_triple(dup._sign, dup._int[:end], exp)
2598
2599 def quantize(self, exp, rounding=None, context=None):
2600 """Quantize self so its exponent is the same as that of exp.
2601
2602 Similar to self._rescale(exp._exp) but with error checking.
2603 """
2604 exp = _convert_other(exp, raiseit=True)
2605
2606 if context is None:
2607 context = getcontext()
2608 if rounding is None:
2609 rounding = context.rounding
2610
2611 if self._is_special or exp._is_special:
2612 ans = self._check_nans(exp, context)
2613 if ans:
2614 return ans
2615
2616 if exp._isinfinity() or self._isinfinity():
2617 if exp._isinfinity() and self._isinfinity():
2618 return Decimal(self) # if both are inf, it is OK
2619 return context._raise_error(InvalidOperation,
2620 'quantize with one INF')
2621
2622 # exp._exp should be between Etiny and Emax
2623 if not (context.Etiny() <= exp._exp <= context.Emax):
2624 return context._raise_error(InvalidOperation,
2625 'target exponent out of bounds in quantize')
2626
2627 if not self:
2628 ans = _dec_from_triple(self._sign, '0', exp._exp)
2629 return ans._fix(context)
2630
2631 self_adjusted = self.adjusted()
2632 if self_adjusted > context.Emax:
2633 return context._raise_error(InvalidOperation,
2634 'exponent of quantize result too large for current context')
2635 if self_adjusted - exp._exp + 1 > context.prec:
2636 return context._raise_error(InvalidOperation,
2637 'quantize result has too many digits for current context')
2638
2639 ans = self._rescale(exp._exp, rounding)
2640 if ans.adjusted() > context.Emax:
2641 return context._raise_error(InvalidOperation,
2642 'exponent of quantize result too large for current context')
2643 if len(ans._int) > context.prec:
2644 return context._raise_error(InvalidOperation,
2645 'quantize result has too many digits for current context')
2646
2647 # raise appropriate flags
2648 if ans and ans.adjusted() < context.Emin:
2649 context._raise_error(Subnormal)
2650 if ans._exp > self._exp:
2651 if ans != self:
2652 context._raise_error(Inexact)
2653 context._raise_error(Rounded)
2654
2655 # call to fix takes care of any necessary folddown, and
2656 # signals Clamped if necessary
2657 ans = ans._fix(context)
2658 return ans
2659
2660 def same_quantum(self, other, context=None):
2661 """Return True if self and other have the same exponent; otherwise
2662 return False.
2663
2664 If either operand is a special value, the following rules are used:
2665 * return True if both operands are infinities
2666 * return True if both operands are NaNs
2667 * otherwise, return False.
2668 """
2669 other = _convert_other(other, raiseit=True)
2670 if self._is_special or other._is_special:
2671 return (self.is_nan() and other.is_nan() or
2672 self.is_infinite() and other.is_infinite())
2673 return self._exp == other._exp
2674
2675 def _rescale(self, exp, rounding):
2676 """Rescale self so that the exponent is exp, either by padding with zeros
2677 or by truncating digits, using the given rounding mode.
2678
2679 Specials are returned without change. This operation is
2680 quiet: it raises no flags, and uses no information from the
2681 context.
2682
2683 exp = exp to scale to (an integer)
2684 rounding = rounding mode
2685 """
2686 if self._is_special:
2687 return Decimal(self)
2688 if not self:
2689 return _dec_from_triple(self._sign, '0', exp)
2690
2691 if self._exp >= exp:
2692 # pad answer with zeros if necessary
2693 return _dec_from_triple(self._sign,
2694 self._int + '0'*(self._exp - exp), exp)
2695
2696 # too many digits; round and lose data. If self.adjusted() <
2697 # exp-1, replace self by 10**(exp-1) before rounding
2698 digits = len(self._int) + self._exp - exp
2699 if digits < 0:
2700 self = _dec_from_triple(self._sign, '1', exp-1)
2701 digits = 0
2702 this_function = self._pick_rounding_function[rounding]
2703 changed = this_function(self, digits)
2704 coeff = self._int[:digits] or '0'
2705 if changed == 1:
2706 coeff = str(int(coeff)+1)
2707 return _dec_from_triple(self._sign, coeff, exp)
2708
2709 def _round(self, places, rounding):
2710 """Round a nonzero, nonspecial Decimal to a fixed number of
2711 significant figures, using the given rounding mode.
2712
2713 Infinities, NaNs and zeros are returned unaltered.
2714
2715 This operation is quiet: it raises no flags, and uses no
2716 information from the context.
2717
2718 """
2719 if places <= 0:
2720 raise ValueError("argument should be at least 1 in _round")
2721 if self._is_special or not self:
2722 return Decimal(self)
2723 ans = self._rescale(self.adjusted()+1-places, rounding)
2724 # it can happen that the rescale alters the adjusted exponent;
2725 # for example when rounding 99.97 to 3 significant figures.
2726 # When this happens we end up with an extra 0 at the end of
2727 # the number; a second rescale fixes this.
2728 if ans.adjusted() != self.adjusted():
2729 ans = ans._rescale(ans.adjusted()+1-places, rounding)
2730 return ans
2731
2732 def to_integral_exact(self, rounding=None, context=None):
2733 """Rounds to a nearby integer.
2734
2735 If no rounding mode is specified, take the rounding mode from
2736 the context. This method raises the Rounded and Inexact flags
2737 when appropriate.
2738
2739 See also: to_integral_value, which does exactly the same as
2740 this method except that it doesn't raise Inexact or Rounded.
2741 """
2742 if self._is_special:
2743 ans = self._check_nans(context=context)
2744 if ans:
2745 return ans
2746 return Decimal(self)
2747 if self._exp >= 0:
2748 return Decimal(self)
2749 if not self:
2750 return _dec_from_triple(self._sign, '0', 0)
2751 if context is None:
2752 context = getcontext()
2753 if rounding is None:
2754 rounding = context.rounding
2755 ans = self._rescale(0, rounding)
2756 if ans != self:
2757 context._raise_error(Inexact)
2758 context._raise_error(Rounded)
2759 return ans
2760
2761 def to_integral_value(self, rounding=None, context=None):
2762 """Rounds to the nearest integer, without raising inexact, rounded."""
2763 if context is None:
2764 context = getcontext()
2765 if rounding is None:
2766 rounding = context.rounding
2767 if self._is_special:
2768 ans = self._check_nans(context=context)
2769 if ans:
2770 return ans
2771 return Decimal(self)
2772 if self._exp >= 0:
2773 return Decimal(self)
2774 else:
2775 return self._rescale(0, rounding)
2776
2777 # the method name changed, but we provide also the old one, for compatibility
2778 to_integral = to_integral_value
2779
2780 def sqrt(self, context=None):
2781 """Return the square root of self."""
2782 if context is None:
2783 context = getcontext()
2784
2785 if self._is_special:
2786 ans = self._check_nans(context=context)
2787 if ans:
2788 return ans
2789
2790 if self._isinfinity() and self._sign == 0:
2791 return Decimal(self)
2792
2793 if not self:
2794 # exponent = self._exp // 2. sqrt(-0) = -0
2795 ans = _dec_from_triple(self._sign, '0', self._exp // 2)
2796 return ans._fix(context)
2797
2798 if self._sign == 1:
2799 return context._raise_error(InvalidOperation, 'sqrt(-x), x > 0')
2800
2801 # At this point self represents a positive number. Let p be
2802 # the desired precision and express self in the form c*100**e
2803 # with c a positive real number and e an integer, c and e
2804 # being chosen so that 100**(p-1) <= c < 100**p. Then the
2805 # (exact) square root of self is sqrt(c)*10**e, and 10**(p-1)
2806 # <= sqrt(c) < 10**p, so the closest representable Decimal at
2807 # precision p is n*10**e where n = round_half_even(sqrt(c)),
2808 # the closest integer to sqrt(c) with the even integer chosen
2809 # in the case of a tie.
2810 #
2811 # To ensure correct rounding in all cases, we use the
2812 # following trick: we compute the square root to an extra
2813 # place (precision p+1 instead of precision p), rounding down.
2814 # Then, if the result is inexact and its last digit is 0 or 5,
2815 # we increase the last digit to 1 or 6 respectively; if it's
2816 # e