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