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