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