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