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