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