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