1 /* Decimal Number module for the decNumber C Library
2 Copyright (C) 2005 Free Software Foundation, Inc.
3 Contributed by IBM Corporation. Author Mike Cowlishaw.
5 This file is part of GCC.
7 GCC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU General Public License as published by the Free
9 Software Foundation; either version 2, or (at your option) any later
12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 You should have received a copy of the GNU General Public License
18 along with GCC; see the file COPYING. If not, write to the Free
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
22 /* ------------------------------------------------------------------ */
23 /* This module comprises the routines for Standard Decimal Arithmetic */
24 /* as defined in the specification which may be found on the */
25 /* http://www2.hursley.ibm.com/decimal web pages. It implements both */
26 /* the full ('extended') arithmetic and the simpler ('subset') */
31 /* 1. This code is ANSI C89 except: */
33 /* a) Line comments (double forward slash) are used. (Most C */
34 /* compilers accept these. If yours does not, a simple script */
35 /* can be used to convert them to ANSI C comments.) */
37 /* b) Types from C99 stdint.h are used. If you do not have this */
38 /* header file, see the User's Guide section of the decNumber */
39 /* documentation; this lists the necessary definitions. */
41 /* c) If DECDPUN>4, non-ANSI 64-bit 'long long' types are used. */
42 /* To avoid these, set DECDPUN <= 4 (see documentation). */
44 /* 2. The decNumber format which this library uses is optimized for */
45 /* efficient processing of relatively short numbers; in particular */
46 /* it allows the use of fixed sized structures and minimizes copy */
47 /* and move operations. It does, however, support arbitrary */
48 /* precision (up to 999,999,999 digits) and arbitrary exponent */
49 /* range (Emax in the range 0 through 999,999,999 and Emin in the */
50 /* range -999,999,999 through 0). */
52 /* 3. Operands to operator functions are never modified unless they */
53 /* are also specified to be the result number (which is always */
54 /* permitted). Other than that case, operands may not overlap. */
56 /* 4. Error handling: the type of the error is ORed into the status */
57 /* flags in the current context (decContext structure). The */
58 /* SIGFPE signal is then raised if the corresponding trap-enabler */
59 /* flag in the decContext is set (is 1). */
61 /* It is the responsibility of the caller to clear the status */
62 /* flags as required. */
64 /* The result of any routine which returns a number will always */
65 /* be a valid number (which may be a special value, such as an */
66 /* Infinity or NaN). */
68 /* 5. The decNumber format is not an exchangeable concrete */
69 /* representation as it comprises fields which may be machine- */
70 /* dependent (big-endian or little-endian, for example). */
71 /* Canonical conversions to and from strings are provided; other */
72 /* conversions are available in separate modules. */
74 /* 6. Normally, input operands are assumed to be valid. Set DECCHECK */
75 /* to 1 for extended operand checking (including NULL operands). */
76 /* Results are undefined if a badly-formed structure (or a NULL */
77 /* NULL pointer to a structure) is provided, though with DECCHECK */
78 /* enabled the operator routines are protected against exceptions. */
79 /* (Except if the result pointer is NULL, which is unrecoverable.) */
81 /* However, the routines will never cause exceptions if they are */
82 /* given well-formed operands, even if the value of the operands */
83 /* is inappropriate for the operation and DECCHECK is not set. */
85 /* 7. Subset arithmetic is available only if DECSUBSET is set to 1. */
86 /* ------------------------------------------------------------------ */
87 /* Implementation notes for maintenance of this module: */
89 /* 1. Storage leak protection: Routines which use malloc are not */
90 /* permitted to use return for fastpath or error exits (i.e., */
91 /* they follow strict structured programming conventions). */
92 /* Instead they have a do{}while(0); construct surrounding the */
93 /* code which is protected -- break may be used from this. */
94 /* Other routines are allowed to use the return statement inline. */
96 /* Storage leak accounting can be enabled using DECALLOC. */
98 /* 2. All loops use the for(;;) construct. Any do construct is for */
99 /* protection as just described. */
101 /* 3. Setting status in the context must always be the very last */
102 /* action in a routine, as non-0 status may raise a trap and hence */
103 /* the call to set status may not return (if the handler uses long */
104 /* jump). Therefore all cleanup must be done first. In general, */
105 /* to achieve this we accumulate status and only finally apply it */
106 /* by calling decContextSetStatus (via decStatus). */
108 /* Routines which allocate storage cannot, therefore, use the */
109 /* 'top level' routines which could cause a non-returning */
110 /* transfer of control. The decXxxxOp routines are safe (do not */
111 /* call decStatus even if traps are set in the context) and should */
112 /* be used instead (they are also a little faster). */
114 /* 4. Exponent checking is minimized by allowing the exponent to */
115 /* grow outside its limits during calculations, provided that */
116 /* the decFinalize function is called later. Multiplication and */
117 /* division, and intermediate calculations in exponentiation, */
118 /* require more careful checks because of the risk of 31-bit */
119 /* overflow (the most negative valid exponent is -1999999997, for */
120 /* a 999999999-digit number with adjusted exponent of -999999999). */
122 /* 5. Rounding is deferred until finalization of results, with any */
123 /* 'off to the right' data being represented as a single digit */
124 /* residue (in the range -1 through 9). This avoids any double- */
125 /* rounding when more than one shortening takes place (for */
126 /* example, when a result is subnormal). */
128 /* 6. The digits count is allowed to rise to a multiple of DECDPUN */
129 /* during many operations, so whole Units are handled and exact */
130 /* accounting of digits is not needed. The correct digits value */
131 /* is found by decGetDigits, which accounts for leading zeros. */
132 /* This must be called before any rounding if the number of digits */
133 /* is not known exactly. */
135 /* 7. We use the multiply-by-reciprocal 'trick' for partitioning */
136 /* numbers up to four digits, using appropriate constants. This */
137 /* is not useful for longer numbers because overflow of 32 bits */
138 /* would lead to 4 multiplies, which is almost as expensive as */
139 /* a divide (unless we assumed floating-point multiply available). */
141 /* 8. Unusual abbreviations possibly used in the commentary: */
142 /* lhs -- left hand side (operand, of an operation) */
143 /* lsd -- least significant digit (of coefficient) */
144 /* lsu -- least significant Unit (of coefficient) */
145 /* msd -- most significant digit (of coefficient) */
146 /* msu -- most significant Unit (of coefficient) */
147 /* rhs -- right hand side (operand, of an operation) */
148 /* +ve -- positive */
149 /* -ve -- negative */
150 /* ------------------------------------------------------------------ */
152 #include <stdlib.h> /* for malloc, free, etc. */
153 #include <stdio.h> /* for printf [if needed] */
154 #include <string.h> /* for strcpy */
155 #include <ctype.h> /* for lower */
157 #include "decNumber.h" /* base number library */
158 #include "decNumberLocal.h" /* decNumber local types, etc. */
161 /* Public constant array: powers of ten (powers[n]==10**n) */
162 const uInt powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
163 10000000, 100000000, 1000000000
166 /* Local constants */
167 #define DIVIDE 0x80 /* Divide operators */
168 #define REMAINDER 0x40 /* .. */
169 #define DIVIDEINT 0x20 /* .. */
170 #define REMNEAR 0x10 /* .. */
171 #define COMPARE 0x01 /* Compare operators */
172 #define COMPMAX 0x02 /* .. */
173 #define COMPMIN 0x03 /* .. */
174 #define COMPNAN 0x04 /* .. [NaN processing] */
176 #define DEC_sNaN 0x40000000 /* local status: sNaN signal */
177 #define BADINT (Int)0x80000000 /* most-negative Int; error indicator */
179 static Unit one[] = { 1 }; /* Unit array of 1, used for incrementing */
181 /* Granularity-dependent code */
183 #define eInt Int /* extended integer */
184 #define ueInt uInt /* unsigned extended integer */
185 /* Constant multipliers for divide-by-power-of five using reciprocal */
186 /* multiply, after removing powers of 2 by shifting, and final shift */
187 /* of 17 [we only need up to **4] */
188 static const uInt multies[] = { 131073, 26215, 5243, 1049, 210 };
190 /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
191 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
193 /* For DECDPUN>4 we currently use non-ANSI 64-bit types. These could */
194 /* be replaced by subroutine calls later. */
198 typedef signed long long Long;
199 typedef unsigned long long uLong;
200 #define eInt Long /* extended integer */
201 #define ueInt uLong /* unsigned extended integer */
205 static decNumber *decAddOp (decNumber *, decNumber *, decNumber *,
206 decContext *, uByte, uInt *);
207 static void decApplyRound (decNumber *, decContext *, Int, uInt *);
208 static Int decCompare (decNumber * lhs, decNumber * rhs);
209 static decNumber *decCompareOp (decNumber *, decNumber *, decNumber *,
210 decContext *, Flag, uInt *);
211 static void decCopyFit (decNumber *, decNumber *, decContext *,
213 static decNumber *decDivideOp (decNumber *, decNumber *, decNumber *,
214 decContext *, Flag, uInt *);
215 static void decFinalize (decNumber *, decContext *, Int *, uInt *);
216 static Int decGetDigits (Unit *, Int);
218 static Int decGetInt (decNumber *, decContext *);
220 static Int decGetInt (decNumber *);
222 static decNumber *decMultiplyOp (decNumber *, decNumber *, decNumber *,
223 decContext *, uInt *);
224 static decNumber *decNaNs (decNumber *, decNumber *, decNumber *, uInt *);
225 static decNumber *decQuantizeOp (decNumber *, decNumber *, decNumber *,
226 decContext *, Flag, uInt *);
227 static void decSetCoeff (decNumber *, decContext *, Unit *,
229 static void decSetOverflow (decNumber *, decContext *, uInt *);
230 static void decSetSubnormal (decNumber *, decContext *, Int *, uInt *);
231 static Int decShiftToLeast (Unit *, Int, Int);
232 static Int decShiftToMost (Unit *, Int, Int);
233 static void decStatus (decNumber *, uInt, decContext *);
234 static Flag decStrEq (const char *, const char *);
235 static void decToString (decNumber *, char[], Flag);
236 static decNumber *decTrim (decNumber *, Flag, Int *);
237 static Int decUnitAddSub (Unit *, Int, Unit *, Int, Int, Unit *, Int);
238 static Int decUnitCompare (Unit *, Int, Unit *, Int, Int);
241 /* decFinish == decFinalize when no subset arithmetic needed */
242 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
244 static void decFinish (decNumber *, decContext *, Int *, uInt *);
245 static decNumber *decRoundOperand (decNumber *, decContext *, uInt *);
248 /* Diagnostic macros, etc. */
250 /* Handle malloc/free accounting. If enabled, our accountable routines */
251 /* are used; otherwise the code just goes straight to the system malloc */
252 /* and free routines. */
253 #define malloc(a) decMalloc(a)
254 #define free(a) decFree(a)
255 #define DECFENCE 0x5a /* corruption detector */
256 /* 'Our' malloc and free: */
257 static void *decMalloc (size_t);
258 static void decFree (void *);
259 uInt decAllocBytes = 0; /* count of bytes allocated */
260 /* Note that DECALLOC code only checks for storage buffer overflow. */
261 /* To check for memory leaks, the decAllocBytes variable should be */
262 /* checked to be 0 at appropriate times (e.g., after the test */
263 /* harness completes a set of tests). This checking may be unreliable */
264 /* if the testing is done in a multi-thread environment. */
268 /* Optional operand checking routines. Enabling these means that */
269 /* decNumber and decContext operands to operator routines are checked */
270 /* for correctness. This roughly doubles the execution time of the */
271 /* fastest routines (and adds 600+ bytes), so should not normally be */
272 /* used in 'production'. */
273 #define DECUNUSED (void *)(0xffffffff)
274 static Flag decCheckOperands (decNumber *, decNumber *, decNumber *,
276 static Flag decCheckNumber (decNumber *, decContext *);
279 #if DECTRACE || DECCHECK
280 /* Optional trace/debugging routines. */
281 void decNumberShow (decNumber *); /* displays the components of a number */
282 static void decDumpAr (char, Unit *, Int);
285 /* ================================================================== */
287 /* ================================================================== */
289 /* ------------------------------------------------------------------ */
290 /* to-scientific-string -- conversion to numeric string */
291 /* to-engineering-string -- conversion to numeric string */
293 /* decNumberToString(dn, string); */
294 /* decNumberToEngString(dn, string); */
296 /* dn is the decNumber to convert */
297 /* string is the string where the result will be laid out */
299 /* string must be at least dn->digits+14 characters long */
301 /* No error is possible, and no status can be set. */
302 /* ------------------------------------------------------------------ */
304 decNumberToString (decNumber * dn, char *string)
306 decToString (dn, string, 0);
311 decNumberToEngString (decNumber * dn, char *string)
313 decToString (dn, string, 1);
317 /* ------------------------------------------------------------------ */
318 /* to-number -- conversion from numeric string */
320 /* decNumberFromString -- convert string to decNumber */
321 /* dn -- the number structure to fill */
322 /* chars[] -- the string to convert ('\0' terminated) */
323 /* set -- the context used for processing any error, */
324 /* determining the maximum precision available */
325 /* (set.digits), determining the maximum and minimum */
326 /* exponent (set.emax and set.emin), determining if */
327 /* extended values are allowed, and checking the */
328 /* rounding mode if overflow occurs or rounding is */
331 /* The length of the coefficient and the size of the exponent are */
332 /* checked by this routine, so the correct error (Underflow or */
333 /* Overflow) can be reported or rounding applied, as necessary. */
335 /* If bad syntax is detected, the result will be a quiet NaN. */
336 /* ------------------------------------------------------------------ */
338 decNumberFromString (decNumber * dn, char chars[], decContext * set)
340 Int exponent = 0; /* working exponent [assume 0] */
341 uByte bits = 0; /* working flags [assume +ve] */
342 Unit *res; /* where result will be built */
343 Unit resbuff[D2U (DECBUFFER + 1)]; /* local buffer in case need temporary */
344 Unit *allocres = NULL; /* -> allocated result, iff allocated */
345 Int need; /* units needed for result */
346 Int d = 0; /* count of digits found in decimal part */
347 char *dotchar = NULL; /* where dot was found */
348 char *cfirst; /* -> first character of decimal part */
349 char *last = NULL; /* -> last digit of decimal part */
350 char *firstexp; /* -> first significant exponent digit */
356 Int residue = 0; /* rounding residue */
357 uInt status = 0; /* error code */
360 if (decCheckOperands (DECUNUSED, DECUNUSED, DECUNUSED, set))
361 return decNumberZero (dn);
365 { /* status & malloc protection */
366 c = chars; /* -> input character */
368 { /* handle leading '-' */
373 c++; /* step over leading '+' */
374 /* We're at the start of the number [we think] */
375 cfirst = c; /* save */
378 if (*c >= '0' && *c <= '9')
379 { /* test for Arabic digit */
381 d++; /* count of real digits */
382 continue; /* still in decimal part */
385 break; /* done with decimal part */
386 /* dot: record, check, and ignore */
389 last = NULL; /* indicate bad */
391 } /* .. and go report */
392 dotchar = c; /* offset into decimal part */
396 { /* no decimal digits, or >1 . */
398 /* If subset then infinities and NaNs are not allowed */
401 status = DEC_Conversion_syntax;
402 break; /* all done */
407 /* Infinities and NaNs are possible, here */
408 decNumberZero (dn); /* be optimistic */
409 if (decStrEq (c, "Infinity") || decStrEq (c, "Inf"))
411 dn->bits = bits | DECINF;
412 break; /* all done */
415 { /* a NaN expected */
416 /* 2003.09.10 NaNs are now permitted to have a sign */
417 status = DEC_Conversion_syntax; /* assume the worst */
418 dn->bits = bits | DECNAN; /* assume simple NaN */
419 if (*c == 's' || *c == 'S')
420 { /* looks like an` sNaN */
422 dn->bits = bits | DECSNAN;
424 if (*c != 'n' && *c != 'N')
425 break; /* check caseless "NaN" */
427 if (*c != 'a' && *c != 'A')
430 if (*c != 'n' && *c != 'N')
433 /* now nothing, or nnnn, expected */
434 /* -> start of integer and skip leading 0s [including plain 0] */
435 for (cfirst = c; *cfirst == '0';)
438 { /* "NaN" or "sNaN", maybe with all 0s */
439 status = 0; /* it's good */
442 /* something other than 0s; setup last and d as usual [no dots] */
443 for (c = cfirst;; c++, d++)
445 if (*c < '0' || *c > '9')
446 break; /* test for Arabic digit */
450 break; /* not all digits */
452 break; /* too many digits */
453 /* good; drop through and convert the integer */
455 bits = dn->bits; /* for copy-back */
463 { /* more there; exponent expected... */
464 Flag nege = 0; /* 1=negative exponent */
465 if (*c != 'e' && *c != 'E')
467 status = DEC_Conversion_syntax;
471 /* Found 'e' or 'E' -- now process explicit exponent */
472 /* 1998.07.11: sign no longer required */
473 c++; /* to (expected) sign */
483 status = DEC_Conversion_syntax;
487 for (; *c == '0' && *(c + 1) != '\0';)
488 c++; /* strip insignificant zeros */
489 firstexp = c; /* save exponent digit place */
492 if (*c < '0' || *c > '9')
493 break; /* not a digit */
494 exponent = X10 (exponent) + (Int) * c - (Int) '0';
496 /* if we didn't end on '\0' must not be a digit */
499 status = DEC_Conversion_syntax;
503 /* (this next test must be after the syntax check) */
504 /* if it was too long the exponent may have wrapped, so check */
505 /* carefully and set it to a certain overflow if wrap possible */
506 if (c >= firstexp + 9 + 1)
508 if (c > firstexp + 9 + 1 || *firstexp > '1')
509 exponent = DECNUMMAXE * 2;
510 /* [up to 1999999999 is OK, for example 1E-1000000998] */
513 exponent = -exponent; /* was negative */
515 /* Here when all inspected; syntax is good */
517 /* Handle decimal point... */
518 if (dotchar != NULL && dotchar < last) /* embedded . found, so */
519 exponent = exponent - (last - dotchar); /* .. adjust exponent */
520 /* [we can now ignore the .] */
522 /* strip leading zeros/dot (leave final if all 0's) */
523 for (c = cfirst; c < last; c++)
526 d--; /* 0 stripped */
529 cfirst++; /* step past leader */
533 /* We can now make a rapid exit for zeros if !extended */
534 if (*cfirst == '0' && !set->extended)
536 decNumberZero (dn); /* clean result */
537 break; /* [could be return] */
541 /* OK, the digits string is good. Copy to the decNumber, or to
542 a temporary decNumber if rounding is needed */
543 if (d <= set->digits)
544 res = dn->lsu; /* fits into given decNumber */
546 { /* rounding needed */
547 need = D2U (d); /* units needed */
548 res = resbuff; /* assume use local buffer */
549 if (need * sizeof (Unit) > sizeof (resbuff))
550 { /* too big for local */
551 allocres = (Unit *) malloc (need * sizeof (Unit));
552 if (allocres == NULL)
554 status |= DEC_Insufficient_storage;
560 /* res now -> number lsu, buffer, or allocated storage for Unit array */
562 /* Place the coefficient into the selected Unit array */
564 i = d % DECDPUN; /* digits in top unit */
567 up = res + D2U (d) - 1; /* -> msu */
569 for (c = cfirst;; c++)
570 { /* along the digits */
572 { /* ignore . [don't decrement i] */
577 *up = (Unit) (X10 (*up) + (Int) * c - (Int) '0');
580 continue; /* more for this unit */
582 break; /* just filled the last unit */
589 up = res; /* -> lsu */
590 for (c = last; c >= cfirst; c--)
591 { /* over each character, from least */
593 continue; /* ignore . [don't step b] */
594 *up = (Unit) ((Int) * c - (Int) '0');
600 dn->exponent = exponent;
603 /* if not in number (too long) shorten into the number */
605 decSetCoeff (dn, set, res, d, &residue, &status);
607 /* Finally check for overflow or subnormal and round as needed */
608 decFinalize (dn, set, &residue, &status);
609 /* decNumberShow(dn); */
611 while (0); /* [for break] */
613 if (allocres != NULL)
614 free (allocres); /* drop any storage we used */
616 decStatus (dn, status, set);
620 /* ================================================================== */
622 /* ================================================================== */
624 /* ------------------------------------------------------------------ */
625 /* decNumberAbs -- absolute value operator */
627 /* This computes C = abs(A) */
629 /* res is C, the result. C may be A */
631 /* set is the context */
633 /* C must have space for set->digits digits. */
634 /* ------------------------------------------------------------------ */
635 /* This has the same effect as decNumberPlus unless A is negative, */
636 /* in which case it has the same effect as decNumberMinus. */
637 /* ------------------------------------------------------------------ */
639 decNumberAbs (decNumber * res, decNumber * rhs, decContext * set)
641 decNumber dzero; /* for 0 */
642 uInt status = 0; /* accumulator */
645 if (decCheckOperands (res, DECUNUSED, rhs, set))
649 decNumberZero (&dzero); /* set 0 */
650 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
651 decAddOp (res, &dzero, rhs, set, (uByte) (rhs->bits & DECNEG), &status);
653 decStatus (res, status, set);
657 /* ------------------------------------------------------------------ */
658 /* decNumberAdd -- add two Numbers */
660 /* This computes C = A + B */
662 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
665 /* set is the context */
667 /* C must have space for set->digits digits. */
668 /* ------------------------------------------------------------------ */
669 /* This just calls the routine shared with Subtract */
671 decNumberAdd (decNumber * res, decNumber * lhs, decNumber * rhs,
674 uInt status = 0; /* accumulator */
675 decAddOp (res, lhs, rhs, set, 0, &status);
677 decStatus (res, status, set);
681 /* ------------------------------------------------------------------ */
682 /* decNumberCompare -- compare two Numbers */
684 /* This computes C = A ? B */
686 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
689 /* set is the context */
691 /* C must have space for one digit. */
692 /* ------------------------------------------------------------------ */
694 decNumberCompare (decNumber * res, decNumber * lhs, decNumber * rhs,
697 uInt status = 0; /* accumulator */
698 decCompareOp (res, lhs, rhs, set, COMPARE, &status);
700 decStatus (res, status, set);
704 /* ------------------------------------------------------------------ */
705 /* decNumberDivide -- divide one number by another */
707 /* This computes C = A / B */
709 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
712 /* set is the context */
714 /* C must have space for set->digits digits. */
715 /* ------------------------------------------------------------------ */
717 decNumberDivide (decNumber * res, decNumber * lhs,
718 decNumber * rhs, decContext * set)
720 uInt status = 0; /* accumulator */
721 decDivideOp (res, lhs, rhs, set, DIVIDE, &status);
723 decStatus (res, status, set);
727 /* ------------------------------------------------------------------ */
728 /* decNumberDivideInteger -- divide and return integer quotient */
730 /* This computes C = A # B, where # is the integer divide operator */
732 /* res is C, the result. C may be A and/or B (e.g., X=X#X) */
735 /* set is the context */
737 /* C must have space for set->digits digits. */
738 /* ------------------------------------------------------------------ */
740 decNumberDivideInteger (decNumber * res, decNumber * lhs,
741 decNumber * rhs, decContext * set)
743 uInt status = 0; /* accumulator */
744 decDivideOp (res, lhs, rhs, set, DIVIDEINT, &status);
746 decStatus (res, status, set);
750 /* ------------------------------------------------------------------ */
751 /* decNumberMax -- compare two Numbers and return the maximum */
753 /* This computes C = A ? B, returning the maximum or A if equal */
755 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
758 /* set is the context */
760 /* C must have space for set->digits digits. */
761 /* ------------------------------------------------------------------ */
763 decNumberMax (decNumber * res, decNumber * lhs, decNumber * rhs,
766 uInt status = 0; /* accumulator */
767 decCompareOp (res, lhs, rhs, set, COMPMAX, &status);
769 decStatus (res, status, set);
773 /* ------------------------------------------------------------------ */
774 /* decNumberMin -- compare two Numbers and return the minimum */
776 /* This computes C = A ? B, returning the minimum or A if equal */
778 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
781 /* set is the context */
783 /* C must have space for set->digits digits. */
784 /* ------------------------------------------------------------------ */
786 decNumberMin (decNumber * res, decNumber * lhs, decNumber * rhs,
789 uInt status = 0; /* accumulator */
790 decCompareOp (res, lhs, rhs, set, COMPMIN, &status);
792 decStatus (res, status, set);
796 /* ------------------------------------------------------------------ */
797 /* decNumberMinus -- prefix minus operator */
799 /* This computes C = 0 - A */
801 /* res is C, the result. C may be A */
803 /* set is the context */
805 /* C must have space for set->digits digits. */
806 /* ------------------------------------------------------------------ */
807 /* We simply use AddOp for the subtract, which will do the necessary. */
808 /* ------------------------------------------------------------------ */
810 decNumberMinus (decNumber * res, decNumber * rhs, decContext * set)
813 uInt status = 0; /* accumulator */
816 if (decCheckOperands (res, DECUNUSED, rhs, set))
820 decNumberZero (&dzero); /* make 0 */
821 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
822 decAddOp (res, &dzero, rhs, set, DECNEG, &status);
824 decStatus (res, status, set);
828 /* ------------------------------------------------------------------ */
829 /* decNumberPlus -- prefix plus operator */
831 /* This computes C = 0 + A */
833 /* res is C, the result. C may be A */
835 /* set is the context */
837 /* C must have space for set->digits digits. */
838 /* ------------------------------------------------------------------ */
839 /* We simply use AddOp; Add will take fast path after preparing A. */
840 /* Performance is a concern here, as this routine is often used to */
841 /* check operands and apply rounding and overflow/underflow testing. */
842 /* ------------------------------------------------------------------ */
844 decNumberPlus (decNumber * res, decNumber * rhs, decContext * set)
847 uInt status = 0; /* accumulator */
850 if (decCheckOperands (res, DECUNUSED, rhs, set))
854 decNumberZero (&dzero); /* make 0 */
855 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
856 decAddOp (res, &dzero, rhs, set, 0, &status);
858 decStatus (res, status, set);
862 /* ------------------------------------------------------------------ */
863 /* decNumberMultiply -- multiply two Numbers */
865 /* This computes C = A x B */
867 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
870 /* set is the context */
872 /* C must have space for set->digits digits. */
873 /* ------------------------------------------------------------------ */
875 decNumberMultiply (decNumber * res, decNumber * lhs,
876 decNumber * rhs, decContext * set)
878 uInt status = 0; /* accumulator */
879 decMultiplyOp (res, lhs, rhs, set, &status);
881 decStatus (res, status, set);
885 /* ------------------------------------------------------------------ */
886 /* decNumberNormalize -- remove trailing zeros */
888 /* This computes C = 0 + A, and normalizes the result */
890 /* res is C, the result. C may be A */
892 /* set is the context */
894 /* C must have space for set->digits digits. */
895 /* ------------------------------------------------------------------ */
897 decNumberNormalize (decNumber * res, decNumber * rhs, decContext * set)
899 decNumber *allocrhs = NULL; /* non-NULL if rounded rhs allocated */
900 uInt status = 0; /* as usual */
901 Int residue = 0; /* as usual */
902 Int dropped; /* work */
905 if (decCheckOperands (res, DECUNUSED, rhs, set))
910 { /* protect allocated storage */
914 /* reduce operand and set lostDigits status, as needed */
915 if (rhs->digits > set->digits)
917 allocrhs = decRoundOperand (rhs, set, &status);
918 if (allocrhs == NULL)
924 /* [following code does not require input rounding] */
926 /* specials copy through, except NaNs need care */
927 if (decNumberIsNaN (rhs))
929 decNaNs (res, rhs, NULL, &status);
933 /* reduce result to the requested length and copy to result */
934 decCopyFit (res, rhs, set, &residue, &status); /* copy & round */
935 decFinish (res, set, &residue, &status); /* cleanup/set flags */
936 decTrim (res, 1, &dropped); /* normalize in place */
938 while (0); /* end protected */
940 if (allocrhs != NULL)
941 free (allocrhs); /* .. */
943 decStatus (res, status, set); /* then report status */
947 /* ------------------------------------------------------------------ */
948 /* decNumberPower -- raise a number to an integer power */
950 /* This computes C = A ** B */
952 /* res is C, the result. C may be A and/or B (e.g., X=X**X) */
955 /* set is the context */
957 /* C must have space for set->digits digits. */
959 /* Specification restriction: abs(n) must be <=999999999 */
960 /* ------------------------------------------------------------------ */
962 decNumberPower (decNumber * res, decNumber * lhs,
963 decNumber * rhs, decContext * set)
965 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
966 decNumber *allocrhs = NULL; /* .., rhs */
967 decNumber *allocdac = NULL; /* -> allocated acc buffer, iff used */
968 decNumber *inrhs = rhs; /* save original rhs */
969 Int reqdigits = set->digits; /* requested DIGITS */
970 Int n; /* RHS in binary */
973 Int dropped; /* .. */
975 uInt needbytes; /* buffer size needed */
976 Flag seenbit; /* seen a bit while powering */
977 Int residue = 0; /* rounding residue */
978 uInt status = 0; /* accumulator */
979 uByte bits = 0; /* result sign if errors */
980 decContext workset; /* working context */
981 decNumber dnOne; /* work value 1... */
982 /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
983 uByte dacbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
984 /* same again for possible 1/lhs calculation */
985 uByte lhsbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
986 decNumber *dac = (decNumber *) dacbuff; /* -> result accumulator */
989 if (decCheckOperands (res, lhs, rhs, set))
994 { /* protect allocated storage */
998 /* reduce operands and set lostDigits status, as needed */
999 if (lhs->digits > reqdigits)
1001 alloclhs = decRoundOperand (lhs, set, &status);
1002 if (alloclhs == NULL)
1006 /* rounding won't affect the result, but we might signal lostDigits */
1007 /* as well as the error for non-integer [x**y would need this too] */
1008 if (rhs->digits > reqdigits)
1010 allocrhs = decRoundOperand (rhs, set, &status);
1011 if (allocrhs == NULL)
1017 /* [following code does not require input rounding] */
1019 /* handle rhs Infinity */
1020 if (decNumberIsInfinite (rhs))
1022 status |= DEC_Invalid_operation; /* bad */
1026 if ((lhs->bits | rhs->bits) & (DECNAN | DECSNAN))
1028 decNaNs (res, lhs, rhs, &status);
1032 /* Original rhs must be an integer that fits and is in range */
1034 n = decGetInt (inrhs, set);
1036 n = decGetInt (inrhs);
1038 if (n == BADINT || n > 999999999 || n < -999999999)
1040 status |= DEC_Invalid_operation;
1045 n = -n; /* use the absolute value */
1047 if (decNumberIsNegative (lhs) /* -x .. */
1048 && (n & 0x00000001))
1049 bits = DECNEG; /* .. to an odd power */
1051 /* handle LHS infinity */
1052 if (decNumberIsInfinite (lhs))
1053 { /* [NaNs already handled] */
1054 uByte rbits = rhs->bits; /* save */
1055 decNumberZero (res);
1057 *res->lsu = 1; /* [-]Inf**0 => 1 */
1060 if (!(rbits & DECNEG))
1061 bits |= DECINF; /* was not a **-n */
1062 /* [otherwise will be 0 or -0] */
1068 /* clone the context */
1069 workset = *set; /* copy all fields */
1070 /* calculate the working DIGITS */
1071 workset.digits = reqdigits + (inrhs->digits + inrhs->exponent) + 1;
1072 /* it's an error if this is more than we can handle */
1073 if (workset.digits > DECNUMMAXP)
1075 status |= DEC_Invalid_operation;
1079 /* workset.digits is the count of digits for the accumulator we need */
1080 /* if accumulator is too long for local storage, then allocate */
1082 sizeof (decNumber) + (D2U (workset.digits) - 1) * sizeof (Unit);
1083 /* [needbytes also used below if 1/lhs needed] */
1084 if (needbytes > sizeof (dacbuff))
1086 allocdac = (decNumber *) malloc (needbytes);
1087 if (allocdac == NULL)
1088 { /* hopeless -- abandon */
1089 status |= DEC_Insufficient_storage;
1092 dac = allocdac; /* use the allocated space */
1094 decNumberZero (dac); /* acc=1 */
1095 *dac->lsu = 1; /* .. */
1098 { /* x**0 is usually 1 */
1099 /* 0**0 is bad unless subset, when it becomes 1 */
1105 status |= DEC_Invalid_operation;
1107 decNumberCopy (res, dac); /* copy the 1 */
1111 /* if a negative power we'll need the constant 1, and if not subset */
1112 /* we'll invert the lhs now rather than inverting the result later */
1113 if (decNumberIsNegative (rhs))
1114 { /* was a **-n [hence digits>0] */
1115 decNumberCopy (&dnOne, dac); /* dnOne=1; [needed now or later] */
1118 { /* need to calculate 1/lhs */
1120 /* divide lhs into 1, putting result in dac [dac=1/dac] */
1121 decDivideOp (dac, &dnOne, lhs, &workset, DIVIDE, &status);
1122 if (alloclhs != NULL)
1124 free (alloclhs); /* done with intermediate */
1125 alloclhs = NULL; /* indicate freed */
1127 /* now locate or allocate space for the inverted lhs */
1128 if (needbytes > sizeof (lhsbuff))
1130 alloclhs = (decNumber *) malloc (needbytes);
1131 if (alloclhs == NULL)
1132 { /* hopeless -- abandon */
1133 status |= DEC_Insufficient_storage;
1136 lhs = alloclhs; /* use the allocated space */
1139 lhs = (decNumber *) lhsbuff; /* use stack storage */
1140 /* [lhs now points to buffer or allocated storage] */
1141 decNumberCopy (lhs, dac); /* copy the 1/lhs */
1142 decNumberCopy (dac, &dnOne); /* restore acc=1 */
1148 /* Raise-to-the-power loop... */
1149 seenbit = 0; /* set once we've seen a 1-bit */
1151 { /* for each bit [top bit ignored] */
1152 /* abandon if we have had overflow or terminal underflow */
1153 if (status & (DEC_Overflow | DEC_Underflow))
1154 { /* interesting? */
1155 if (status & DEC_Overflow || ISZERO (dac))
1158 /* [the following two lines revealed an optimizer bug in a C++ */
1159 /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
1160 n = n << 1; /* move next bit to testable position */
1162 { /* top bit is set */
1163 seenbit = 1; /* OK, we're off */
1164 decMultiplyOp (dac, dac, lhs, &workset, &status); /* dac=dac*x */
1167 break; /* that was the last bit */
1169 continue; /* we don't have to square 1 */
1170 decMultiplyOp (dac, dac, dac, &workset, &status); /* dac=dac*dac [square] */
1171 } /*i *//* 32 bits */
1173 /* complete internal overflow or underflow processing */
1174 if (status & (DEC_Overflow | DEC_Subnormal))
1177 /* If subset, and power was negative, reverse the kind of -erflow */
1178 /* [1/x not yet done] */
1179 if (!set->extended && decNumberIsNegative (rhs))
1181 if (status & DEC_Overflow)
1182 status ^= DEC_Overflow | DEC_Underflow | DEC_Subnormal;
1184 { /* trickier -- Underflow may or may not be set */
1185 status &= ~(DEC_Underflow | DEC_Subnormal); /* [one or both] */
1186 status |= DEC_Overflow;
1190 dac->bits = (dac->bits & ~DECNEG) | bits; /* force correct sign */
1191 /* round subnormals [to set.digits rather than workset.digits] */
1192 /* or set overflow result similarly as required */
1193 decFinalize (dac, set, &residue, &status);
1194 decNumberCopy (res, dac); /* copy to result (is now OK length) */
1199 if (!set->extended && /* subset math */
1200 decNumberIsNegative (rhs))
1201 { /* was a **-n [hence digits>0] */
1202 /* so divide result into 1 [dac=1/dac] */
1203 decDivideOp (dac, &dnOne, dac, &workset, DIVIDE, &status);
1207 /* reduce result to the requested length and copy to result */
1208 decCopyFit (res, dac, set, &residue, &status);
1209 decFinish (res, set, &residue, &status); /* final cleanup */
1212 decTrim (res, 0, &dropped); /* trailing zeros */
1215 while (0); /* end protected */
1217 if (allocdac != NULL)
1218 free (allocdac); /* drop any storage we used */
1219 if (allocrhs != NULL)
1220 free (allocrhs); /* .. */
1221 if (alloclhs != NULL)
1222 free (alloclhs); /* .. */
1224 decStatus (res, status, set);
1228 /* ------------------------------------------------------------------ */
1229 /* decNumberQuantize -- force exponent to requested value */
1231 /* This computes C = op(A, B), where op adjusts the coefficient */
1232 /* of C (by rounding or shifting) such that the exponent (-scale) */
1233 /* of C has exponent of B. The numerical value of C will equal A, */
1234 /* except for the effects of any rounding that occurred. */
1236 /* res is C, the result. C may be A or B */
1237 /* lhs is A, the number to adjust */
1238 /* rhs is B, the number with exponent to match */
1239 /* set is the context */
1241 /* C must have space for set->digits digits. */
1243 /* Unless there is an error or the result is infinite, the exponent */
1244 /* after the operation is guaranteed to be equal to that of B. */
1245 /* ------------------------------------------------------------------ */
1247 decNumberQuantize (decNumber * res, decNumber * lhs,
1248 decNumber * rhs, decContext * set)
1250 uInt status = 0; /* accumulator */
1251 decQuantizeOp (res, lhs, rhs, set, 1, &status);
1253 decStatus (res, status, set);
1257 /* ------------------------------------------------------------------ */
1258 /* decNumberRescale -- force exponent to requested value */
1260 /* This computes C = op(A, B), where op adjusts the coefficient */
1261 /* of C (by rounding or shifting) such that the exponent (-scale) */
1262 /* of C has the value B. The numerical value of C will equal A, */
1263 /* except for the effects of any rounding that occurred. */
1265 /* res is C, the result. C may be A or B */
1266 /* lhs is A, the number to adjust */
1267 /* rhs is B, the requested exponent */
1268 /* set is the context */
1270 /* C must have space for set->digits digits. */
1272 /* Unless there is an error or the result is infinite, the exponent */
1273 /* after the operation is guaranteed to be equal to B. */
1274 /* ------------------------------------------------------------------ */
1276 decNumberRescale (decNumber * res, decNumber * lhs,
1277 decNumber * rhs, decContext * set)
1279 uInt status = 0; /* accumulator */
1280 decQuantizeOp (res, lhs, rhs, set, 0, &status);
1282 decStatus (res, status, set);
1286 /* ------------------------------------------------------------------ */
1287 /* decNumberRemainder -- divide and return remainder */
1289 /* This computes C = A % B */
1291 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1294 /* set is the context */
1296 /* C must have space for set->digits digits. */
1297 /* ------------------------------------------------------------------ */
1299 decNumberRemainder (decNumber * res, decNumber * lhs,
1300 decNumber * rhs, decContext * set)
1302 uInt status = 0; /* accumulator */
1303 decDivideOp (res, lhs, rhs, set, REMAINDER, &status);
1305 decStatus (res, status, set);
1309 /* ------------------------------------------------------------------ */
1310 /* decNumberRemainderNear -- divide and return remainder from nearest */
1312 /* This computes C = A % B, where % is the IEEE remainder operator */
1314 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1317 /* set is the context */
1319 /* C must have space for set->digits digits. */
1320 /* ------------------------------------------------------------------ */
1322 decNumberRemainderNear (decNumber * res, decNumber * lhs,
1323 decNumber * rhs, decContext * set)
1325 uInt status = 0; /* accumulator */
1326 decDivideOp (res, lhs, rhs, set, REMNEAR, &status);
1328 decStatus (res, status, set);
1332 /* ------------------------------------------------------------------ */
1333 /* decNumberSameQuantum -- test for equal exponents */
1335 /* res is the result number, which will contain either 0 or 1 */
1336 /* lhs is a number to test */
1337 /* rhs is the second (usually a pattern) */
1339 /* No errors are possible and no context is needed. */
1340 /* ------------------------------------------------------------------ */
1342 decNumberSameQuantum (decNumber * res, decNumber * lhs, decNumber * rhs)
1344 uByte merged; /* merged flags */
1345 Unit ret = 0; /* return value */
1348 if (decCheckOperands (res, lhs, rhs, DECUNUSED))
1352 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
1355 if (decNumberIsNaN (lhs) && decNumberIsNaN (rhs))
1357 else if (decNumberIsInfinite (lhs) && decNumberIsInfinite (rhs))
1359 /* [anything else with a special gives 0] */
1361 else if (lhs->exponent == rhs->exponent)
1364 decNumberZero (res); /* OK to overwrite an operand */
1369 /* ------------------------------------------------------------------ */
1370 /* decNumberSquareRoot -- square root operator */
1372 /* This computes C = squareroot(A) */
1374 /* res is C, the result. C may be A */
1376 /* set is the context; note that rounding mode has no effect */
1378 /* C must have space for set->digits digits. */
1379 /* ------------------------------------------------------------------ */
1380 /* This uses the following varying-precision algorithm in: */
1382 /* Properly Rounded Variable Precision Square Root, T. E. Hull and */
1383 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
1384 /* pp229-237, ACM, September 1985. */
1386 /* % [Reformatted original Numerical Turing source code follows.] */
1387 /* function sqrt(x : real) : real */
1388 /* % sqrt(x) returns the properly rounded approximation to the square */
1389 /* % root of x, in the precision of the calling environment, or it */
1390 /* % fails if x < 0. */
1391 /* % t e hull and a abrham, august, 1984 */
1392 /* if x <= 0 then */
1399 /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
1400 /* var e := getexp(x) % exponent part of x */
1401 /* var approx : real */
1402 /* if e mod 2 = 0 then */
1403 /* approx := .259 + .819 * f % approx to root of f */
1405 /* f := f/l0 % adjustments */
1406 /* e := e + 1 % for odd */
1407 /* approx := .0819 + 2.59 * f % exponent */
1411 /* const maxp := currentprecision + 2 */
1413 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
1415 /* approx := .5 * (approx + f/approx) */
1416 /* exit when p = maxp */
1419 /* % approx is now within 1 ulp of the properly rounded square root */
1420 /* % of f; to ensure proper rounding, compare squares of (approx - */
1421 /* % l/2 ulp) and (approx + l/2 ulp) with f. */
1422 /* p := currentprecision */
1424 /* precision p + 2 */
1425 /* const approxsubhalf := approx - setexp(.5, -p) */
1426 /* if mulru(approxsubhalf, approxsubhalf) > f then */
1427 /* approx := approx - setexp(.l, -p + 1) */
1429 /* const approxaddhalf := approx + setexp(.5, -p) */
1430 /* if mulrd(approxaddhalf, approxaddhalf) < f then */
1431 /* approx := approx + setexp(.l, -p + 1) */
1435 /* result setexp(approx, e div 2) % fix exponent */
1437 /* ------------------------------------------------------------------ */
1439 decNumberSquareRoot (decNumber * res, decNumber * rhs, decContext * set)
1441 decContext workset, approxset; /* work contexts */
1442 decNumber dzero; /* used for constant zero */
1443 Int maxp = set->digits + 2; /* largest working precision */
1444 Int residue = 0; /* rounding residue */
1445 uInt status = 0, ignore = 0; /* status accumulators */
1446 Int exp; /* working exponent */
1447 Int ideal; /* ideal (preferred) exponent */
1448 uInt needbytes; /* work */
1449 Int dropped; /* .. */
1451 decNumber *allocrhs = NULL; /* non-NULL if rounded rhs allocated */
1452 /* buffer for f [needs +1 in case DECBUFFER 0] */
1453 uByte buff[sizeof (decNumber) + (D2U (DECBUFFER + 1) - 1) * sizeof (Unit)];
1454 /* buffer for a [needs +2 to match maxp] */
1455 uByte bufa[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1456 /* buffer for temporary, b [must be same size as a] */
1457 uByte bufb[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1458 decNumber *allocbuff = NULL; /* -> allocated buff, iff allocated */
1459 decNumber *allocbufa = NULL; /* -> allocated bufa, iff allocated */
1460 decNumber *allocbufb = NULL; /* -> allocated bufb, iff allocated */
1461 decNumber *f = (decNumber *) buff; /* reduced fraction */
1462 decNumber *a = (decNumber *) bufa; /* approximation to result */
1463 decNumber *b = (decNumber *) bufb; /* intermediate result */
1464 /* buffer for temporary variable, up to 3 digits */
1465 uByte buft[sizeof (decNumber) + (D2U (3) - 1) * sizeof (Unit)];
1466 decNumber *t = (decNumber *) buft; /* up-to-3-digit constant or work */
1469 if (decCheckOperands (res, DECUNUSED, rhs, set))
1474 { /* protect allocated storage */
1478 /* reduce operand and set lostDigits status, as needed */
1479 if (rhs->digits > set->digits)
1481 allocrhs = decRoundOperand (rhs, set, &status);
1482 if (allocrhs == NULL)
1484 /* [Note: 'f' allocation below could reuse this buffer if */
1485 /* used, but as this is rare we keep them separate for clarity.] */
1490 /* [following code does not require input rounding] */
1492 /* handle infinities and NaNs */
1493 if (rhs->bits & DECSPECIAL)
1495 if (decNumberIsInfinite (rhs))
1497 if (decNumberIsNegative (rhs))
1498 status |= DEC_Invalid_operation;
1500 decNumberCopy (res, rhs); /* +Infinity */
1503 decNaNs (res, rhs, NULL, &status); /* a NaN */
1507 /* calculate the ideal (preferred) exponent [floor(exp/2)] */
1508 /* [We would like to write: ideal=rhs->exponent>>1, but this */
1509 /* generates a compiler warning. Generated code is the same.] */
1510 ideal = (rhs->exponent & ~1) / 2; /* target */
1515 decNumberCopy (res, rhs); /* could be 0 or -0 */
1516 res->exponent = ideal; /* use the ideal [safe] */
1520 /* any other -x is an oops */
1521 if (decNumberIsNegative (rhs))
1523 status |= DEC_Invalid_operation;
1527 /* we need space for three working variables */
1528 /* f -- the same precision as the RHS, reduced to 0.01->0.99... */
1529 /* a -- Hull's approx -- precision, when assigned, is */
1530 /* currentprecision (we allow +2 for use as temporary) */
1531 /* b -- intermediate temporary result */
1532 /* if any is too long for local storage, then allocate */
1534 sizeof (decNumber) + (D2U (rhs->digits) - 1) * sizeof (Unit);
1535 if (needbytes > sizeof (buff))
1537 allocbuff = (decNumber *) malloc (needbytes);
1538 if (allocbuff == NULL)
1539 { /* hopeless -- abandon */
1540 status |= DEC_Insufficient_storage;
1543 f = allocbuff; /* use the allocated space */
1545 /* a and b both need to be able to hold a maxp-length number */
1546 needbytes = sizeof (decNumber) + (D2U (maxp) - 1) * sizeof (Unit);
1547 if (needbytes > sizeof (bufa))
1548 { /* [same applies to b] */
1549 allocbufa = (decNumber *) malloc (needbytes);
1550 allocbufb = (decNumber *) malloc (needbytes);
1551 if (allocbufa == NULL || allocbufb == NULL)
1553 status |= DEC_Insufficient_storage;
1556 a = allocbufa; /* use the allocated space */
1557 b = allocbufb; /* .. */
1560 /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
1561 decNumberCopy (f, rhs);
1562 exp = f->exponent + f->digits; /* adjusted to Hull rules */
1563 f->exponent = -(f->digits); /* to range */
1565 /* set up working contexts (the second is used for Numerical */
1566 /* Turing assignment) */
1567 decContextDefault (&workset, DEC_INIT_DECIMAL64);
1568 decContextDefault (&approxset, DEC_INIT_DECIMAL64);
1569 approxset.digits = set->digits; /* approx's length */
1571 /* [Until further notice, no error is possible and status bits */
1572 /* (Rounded, etc.) should be ignored, not accumulated.] */
1574 /* Calculate initial approximation, and allow for odd exponent */
1575 workset.digits = set->digits; /* p for initial calculation */
1581 { /* even exponent */
1582 /* Set t=0.259, a=0.819 */
1603 { /* odd exponent */
1604 /* Set t=0.0819, a=2.59 */
1605 f->exponent--; /* f=f/10 */
1626 decMultiplyOp (a, a, f, &workset, &ignore); /* a=a*f */
1627 decAddOp (a, a, t, &workset, 0, &ignore); /* ..+t */
1628 /* [a is now the initial approximation for sqrt(f), calculated with */
1629 /* currentprecision, which is also a's precision.] */
1631 /* the main calculation loop */
1632 decNumberZero (&dzero); /* make 0 */
1633 decNumberZero (t); /* set t = 0.5 */
1634 t->lsu[0] = 5; /* .. */
1635 t->exponent = -1; /* .. */
1636 workset.digits = 3; /* initial p */
1639 /* set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp] */
1640 workset.digits = workset.digits * 2 - 2;
1641 if (workset.digits > maxp)
1642 workset.digits = maxp;
1643 /* a = 0.5 * (a + f/a) */
1644 /* [calculated at p then rounded to currentprecision] */
1645 decDivideOp (b, f, a, &workset, DIVIDE, &ignore); /* b=f/a */
1646 decAddOp (b, b, a, &workset, 0, &ignore); /* b=b+a */
1647 decMultiplyOp (a, b, t, &workset, &ignore); /* a=b*0.5 */
1648 /* assign to approx [round to length] */
1649 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1650 if (workset.digits == maxp)
1651 break; /* just did final */
1654 /* a is now at currentprecision and within 1 ulp of the properly */
1655 /* rounded square root of f; to ensure proper rounding, compare */
1656 /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
1657 /* Here workset.digits=maxp and t=0.5 */
1658 workset.digits--; /* maxp-1 is OK now */
1659 t->exponent = -set->digits - 1; /* make 0.5 ulp */
1660 decNumberCopy (b, a);
1661 decAddOp (b, b, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */
1662 workset.round = DEC_ROUND_UP;
1663 decMultiplyOp (b, b, b, &workset, &ignore); /* b = mulru(b, b) */
1664 decCompareOp (b, f, b, &workset, COMPARE, &ignore); /* b ? f, reversed */
1665 if (decNumberIsNegative (b))
1666 { /* f < b [i.e., b > f] */
1667 /* this is the more common adjustment, though both are rare */
1668 t->exponent++; /* make 1.0 ulp */
1669 t->lsu[0] = 1; /* .. */
1670 decAddOp (a, a, t, &workset, DECNEG, &ignore); /* a = a - 1 ulp */
1671 /* assign to approx [round to length] */
1672 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1676 decNumberCopy (b, a);
1677 decAddOp (b, b, t, &workset, 0, &ignore); /* b = a + 0.5 ulp */
1678 workset.round = DEC_ROUND_DOWN;
1679 decMultiplyOp (b, b, b, &workset, &ignore); /* b = mulrd(b, b) */
1680 decCompareOp (b, b, f, &workset, COMPARE, &ignore); /* b ? f */
1681 if (decNumberIsNegative (b))
1683 t->exponent++; /* make 1.0 ulp */
1684 t->lsu[0] = 1; /* .. */
1685 decAddOp (a, a, t, &workset, 0, &ignore); /* a = a + 1 ulp */
1686 /* assign to approx [round to length] */
1687 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1690 /* [no errors are possible in the above, and rounding/inexact during */
1691 /* estimation are irrelevant, so status was not accumulated] */
1693 /* Here, 0.1 <= a < 1 [Hull] */
1694 a->exponent += exp / 2; /* set correct exponent */
1696 /* Process Subnormals */
1697 decFinalize (a, set, &residue, &status);
1699 /* count dropable zeros [after any subnormal rounding] */
1700 decNumberCopy (b, a);
1701 decTrim (b, 1, &dropped); /* [drops trailing zeros] */
1703 /* Finally set Inexact and Rounded. The answer can only be exact if */
1704 /* it is short enough so that squaring it could fit in set->digits, */
1705 /* so this is the only (relatively rare) time we have to check */
1707 if (b->digits * 2 - 1 > set->digits)
1709 status |= DEC_Inexact | DEC_Rounded;
1712 { /* could be exact/unrounded */
1713 uInt mstatus = 0; /* local status */
1714 decMultiplyOp (b, b, b, &workset, &mstatus); /* try the multiply */
1716 { /* result won't fit */
1717 status |= DEC_Inexact | DEC_Rounded;
1721 decCompareOp (t, b, rhs, &workset, COMPARE, &mstatus); /* b ? rhs */
1724 status |= DEC_Inexact | DEC_Rounded;
1728 /* here, dropped is the count of trailing zeros in 'a' */
1729 /* use closest exponent to ideal... */
1730 Int todrop = ideal - a->exponent; /* most we can drop */
1733 { /* ideally would add 0s */
1734 status |= DEC_Rounded;
1738 if (dropped < todrop)
1739 todrop = dropped; /* clamp to those available */
1741 { /* OK, some to drop */
1742 decShiftToLeast (a->lsu, D2U (a->digits), todrop);
1743 a->exponent += todrop; /* maintain numerical value */
1744 a->digits -= todrop; /* new length */
1750 decNumberCopy (res, a); /* assume this is the result */
1752 while (0); /* end protected */
1754 if (allocbuff != NULL)
1755 free (allocbuff); /* drop any storage we used */
1756 if (allocbufa != NULL)
1757 free (allocbufa); /* .. */
1758 if (allocbufb != NULL)
1759 free (allocbufb); /* .. */
1760 if (allocrhs != NULL)
1761 free (allocrhs); /* .. */
1763 decStatus (res, status, set); /* then report status */
1767 /* ------------------------------------------------------------------ */
1768 /* decNumberSubtract -- subtract two Numbers */
1770 /* This computes C = A - B */
1772 /* res is C, the result. C may be A and/or B (e.g., X=X-X) */
1775 /* set is the context */
1777 /* C must have space for set->digits digits. */
1778 /* ------------------------------------------------------------------ */
1780 decNumberSubtract (decNumber * res, decNumber * lhs,
1781 decNumber * rhs, decContext * set)
1783 uInt status = 0; /* accumulator */
1785 decAddOp (res, lhs, rhs, set, DECNEG, &status);
1787 decStatus (res, status, set);
1791 /* ------------------------------------------------------------------ */
1792 /* decNumberToIntegralValue -- round-to-integral-value */
1794 /* res is the result */
1795 /* rhs is input number */
1796 /* set is the context */
1798 /* res must have space for any value of rhs. */
1800 /* This implements the IEEE special operator and therefore treats */
1801 /* special values as valid, and also never sets Inexact. For finite */
1802 /* numbers it returns rescale(rhs, 0) if rhs->exponent is <0. */
1803 /* Otherwise the result is rhs (so no error is possible). */
1805 /* The context is used for rounding mode and status after sNaN, but */
1806 /* the digits setting is ignored. */
1807 /* ------------------------------------------------------------------ */
1809 decNumberToIntegralValue (decNumber * res, decNumber * rhs, decContext * set)
1812 decContext workset; /* working context */
1815 if (decCheckOperands (res, DECUNUSED, rhs, set))
1819 /* handle infinities and NaNs */
1820 if (rhs->bits & DECSPECIAL)
1823 if (decNumberIsInfinite (rhs))
1824 decNumberCopy (res, rhs); /* an Infinity */
1826 decNaNs (res, rhs, NULL, &status); /* a NaN */
1828 decStatus (res, status, set);
1832 /* we have a finite number; no error possible */
1833 if (rhs->exponent >= 0)
1834 return decNumberCopy (res, rhs);
1835 /* that was easy, but if negative exponent we have work to do... */
1836 workset = *set; /* clone rounding, etc. */
1837 workset.digits = rhs->digits; /* no length rounding */
1838 workset.traps = 0; /* no traps */
1839 decNumberZero (&dn); /* make a number with exponent 0 */
1840 return decNumberQuantize (res, rhs, &dn, &workset);
1843 /* ================================================================== */
1844 /* Utility routines */
1845 /* ================================================================== */
1847 /* ------------------------------------------------------------------ */
1848 /* decNumberCopy -- copy a number */
1850 /* dest is the target decNumber */
1851 /* src is the source decNumber */
1854 /* (dest==src is allowed and is a no-op) */
1855 /* All fields are updated as required. This is a utility operation, */
1856 /* so special values are unchanged and no error is possible. */
1857 /* ------------------------------------------------------------------ */
1859 decNumberCopy (decNumber * dest, decNumber * src)
1864 return decNumberZero (dest);
1868 return dest; /* no copy required */
1870 /* We use explicit assignments here as structure assignment can copy */
1871 /* more than just the lsu (for small DECDPUN). This would not affect */
1872 /* the value of the results, but would disturb test harness spill */
1874 dest->bits = src->bits;
1875 dest->exponent = src->exponent;
1876 dest->digits = src->digits;
1877 dest->lsu[0] = src->lsu[0];
1878 if (src->digits > DECDPUN)
1879 { /* more Units to come */
1880 Unit *s, *d, *smsup; /* work */
1881 /* memcpy for the remaining Units would be safe as they cannot */
1882 /* overlap. However, this explicit loop is faster in short cases. */
1883 d = dest->lsu + 1; /* -> first destination */
1884 smsup = src->lsu + D2U (src->digits); /* -> source msu+1 */
1885 for (s = src->lsu + 1; s < smsup; s++, d++)
1891 /* ------------------------------------------------------------------ */
1892 /* decNumberTrim -- remove insignificant zeros */
1894 /* dn is the number to trim */
1897 /* All fields are updated as required. This is a utility operation, */
1898 /* so special values are unchanged and no error is possible. */
1899 /* ------------------------------------------------------------------ */
1901 decNumberTrim (decNumber * dn)
1903 Int dropped; /* work */
1904 return decTrim (dn, 0, &dropped);
1907 /* ------------------------------------------------------------------ */
1908 /* decNumberVersion -- return the name and version of this module */
1910 /* No error is possible. */
1911 /* ------------------------------------------------------------------ */
1913 decNumberVersion (void)
1918 /* ------------------------------------------------------------------ */
1919 /* decNumberZero -- set a number to 0 */
1921 /* dn is the number to set, with space for one digit */
1924 /* No error is possible. */
1925 /* ------------------------------------------------------------------ */
1926 /* Memset is not used as it is much slower in some environments. */
1928 decNumberZero (decNumber * dn)
1932 if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
1943 /* ================================================================== */
1944 /* Local routines */
1945 /* ================================================================== */
1947 /* ------------------------------------------------------------------ */
1948 /* decToString -- lay out a number into a string */
1950 /* dn is the number to lay out */
1951 /* string is where to lay out the number */
1952 /* eng is 1 if Engineering, 0 if Scientific */
1954 /* str must be at least dn->digits+14 characters long */
1955 /* No error is possible. */
1957 /* Note that this routine can generate a -0 or 0.000. These are */
1958 /* never generated in subset to-number or arithmetic, but can occur */
1959 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234). */
1960 /* ------------------------------------------------------------------ */
1961 /* If DECCHECK is enabled the string "?" is returned if a number is */
1964 /* TODIGIT -- macro to remove the leading digit from the unsigned */
1965 /* integer u at column cut (counting from the right, LSD=0) and place */
1966 /* it as an ASCII character into the character pointed to by c. Note */
1967 /* that cut must be <= 9, and the maximum value for u is 2,000,000,000 */
1968 /* (as is needed for negative exponents of subnormals). The unsigned */
1969 /* integer pow is used as a temporary variable. */
1970 #define TODIGIT(u, cut, c) { \
1972 pow=powers[cut]*2; \
1975 if ((u)>=pow) {(u)-=pow; *(c)+=8;} \
1977 if ((u)>=pow) {(u)-=pow; *(c)+=4;} \
1980 if ((u)>=pow) {(u)-=pow; *(c)+=2;} \
1982 if ((u)>=pow) {(u)-=pow; *(c)+=1;} \
1986 decToString (decNumber * dn, char *string, Flag eng)
1988 Int exp = dn->exponent; /* local copy */
1989 Int e; /* E-part value */
1990 Int pre; /* digits before the '.' */
1991 Int cut; /* for counting digits in a Unit */
1992 char *c = string; /* work [output pointer] */
1993 Unit *up = dn->lsu + D2U (dn->digits) - 1; /* -> msu [input pointer] */
1994 uInt u, pow; /* work */
1997 if (decCheckOperands (DECUNUSED, dn, DECUNUSED, DECUNUSED))
1999 strcpy (string, "?");
2004 if (decNumberIsNegative (dn))
2005 { /* Negatives get a minus (except */
2006 *c = '-'; /* NaNs, which remove the '-' below) */
2009 if (dn->bits & DECSPECIAL)
2010 { /* Is a special value */
2011 if (decNumberIsInfinite (dn))
2013 strcpy (c, "Infinity");
2017 if (dn->bits & DECSNAN)
2018 { /* signalling NaN */
2023 c += 3; /* step past */
2024 /* if not a clean non-zero coefficient, that's all we have in a */
2026 if (exp != 0 || (*dn->lsu == 0 && dn->digits == 1))
2028 /* [drop through to add integer] */
2031 /* calculate how many digits in msu, and hence first cut */
2032 cut = dn->digits % DECDPUN;
2034 cut = DECDPUN; /* msu is full */
2035 cut--; /* power of ten for digit */
2038 { /* simple integer [common fastpath, */
2039 /* used for NaNs, too] */
2040 for (; up >= dn->lsu; up--)
2041 { /* each Unit from msu */
2042 u = *up; /* contains DECDPUN digits to lay out */
2043 for (; cut >= 0; c++, cut--)
2044 TODIGIT (u, cut, c);
2045 cut = DECDPUN - 1; /* next Unit has all digits */
2047 *c = '\0'; /* terminate the string */
2051 /* non-0 exponent -- assume plain form */
2052 pre = dn->digits + exp; /* digits before '.' */
2054 if ((exp > 0) || (pre < -5))
2055 { /* need exponential form */
2056 e = exp + dn->digits - 1; /* calculate E value */
2057 pre = 1; /* assume one digit before '.' */
2058 if (eng && (e != 0))
2059 { /* may need to adjust */
2060 Int adj; /* adjustment */
2061 /* The C remainder operator is undefined for negative numbers, so */
2062 /* we must use positive remainder calculation here */
2074 /* if we are dealing with zero we will use exponent which is a */
2075 /* multiple of three, as expected, but there will only be the */
2076 /* one zero before the E, still. Otherwise note the padding. */
2082 { /* 0.00Esnn needed */
2090 /* lay out the digits of the coefficient, adding 0s and . as needed */
2093 { /* xxx.xxx or xx00 (engineering) form */
2094 for (; pre > 0; pre--, c++, cut--)
2097 { /* need new Unit */
2099 break; /* out of input digits (pre>digits) */
2104 TODIGIT (u, cut, c);
2106 if (up > dn->lsu || (up == dn->lsu && cut >= 0))
2107 { /* more to come, after '.' */
2113 { /* need new Unit */
2115 break; /* out of input digits */
2120 TODIGIT (u, cut, c);
2124 for (; pre > 0; pre--, c++)
2125 *c = '0'; /* 0 padding (for engineering) needed */
2128 { /* 0.xxx or 0.000xxx form */
2133 for (; pre < 0; pre++, c++)
2134 *c = '0'; /* add any 0's after '.' */
2138 { /* need new Unit */
2140 break; /* out of input digits */
2145 TODIGIT (u, cut, c);
2149 /* Finally add the E-part, if needed. It will never be 0, has a
2150 base maximum and minimum of +999999999 through -999999999, but
2151 could range down to -1999999998 for subnormal numbers */
2154 Flag had = 0; /* 1=had non-zero */
2158 c++; /* assume positive */
2162 *(c - 1) = '-'; /* oops, need - */
2163 u = -e; /* uInt, please */
2165 /* layout the exponent (_itoa is not ANSI C) */
2166 for (cut = 9; cut >= 0; cut--)
2168 TODIGIT (u, cut, c);
2169 if (*c == '0' && !had)
2170 continue; /* skip leading zeros */
2171 had = 1; /* had non-0 */
2172 c++; /* step for next */
2175 *c = '\0'; /* terminate the string (all paths) */
2179 /* ------------------------------------------------------------------ */
2180 /* decAddOp -- add/subtract operation */
2182 /* This computes C = A + B */
2184 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
2187 /* set is the context */
2188 /* negate is DECNEG if rhs should be negated, or 0 otherwise */
2189 /* status accumulates status for the caller */
2191 /* C must have space for set->digits digits. */
2192 /* ------------------------------------------------------------------ */
2193 /* If possible, we calculate the coefficient directly into C. */
2195 /* -- we need a digits+1 calculation because numbers are unaligned */
2196 /* and span more than set->digits digits */
2197 /* -- a carry to digits+1 digits looks possible */
2198 /* -- C is the same as A or B, and the result would destructively */
2199 /* overlap the A or B coefficient */
2200 /* then we must calculate into a temporary buffer. In this latter */
2201 /* case we use the local (stack) buffer if possible, and only if too */
2202 /* long for that do we resort to malloc. */
2204 /* Misalignment is handled as follows: */
2205 /* Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp. */
2206 /* BPad: Apply the padding by a combination of shifting (whole */
2207 /* units) and multiplication (part units). */
2209 /* Addition, especially x=x+1, is speed-critical, so we take pains */
2210 /* to make returning as fast as possible, by flagging any allocation. */
2211 /* ------------------------------------------------------------------ */
2213 decAddOp (decNumber * res, decNumber * lhs,
2214 decNumber * rhs, decContext * set, uByte negate, uInt * status)
2216 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
2217 decNumber *allocrhs = NULL; /* .., rhs */
2218 Int rhsshift; /* working shift (in Units) */
2219 Int maxdigits; /* longest logical length */
2220 Int mult; /* multiplier */
2221 Int residue; /* rounding accumulator */
2222 uByte bits; /* result bits */
2223 Flag diffsign; /* non-0 if arguments have different sign */
2224 Unit *acc; /* accumulator for result */
2225 Unit accbuff[D2U (DECBUFFER + 1)]; /* local buffer [+1 is for possible */
2226 /* final carry digit or DECBUFFER=0] */
2227 Unit *allocacc = NULL; /* -> allocated acc buffer, iff allocated */
2228 Flag alloced = 0; /* set non-0 if any allocations */
2229 Int reqdigits = set->digits; /* local copy; requested DIGITS */
2230 uByte merged; /* merged flags */
2231 Int padding; /* work */
2234 if (decCheckOperands (res, lhs, rhs, set))
2239 { /* protect allocated storage */
2243 /* reduce operands and set lostDigits status, as needed */
2244 if (lhs->digits > reqdigits)
2246 alloclhs = decRoundOperand (lhs, set, status);
2247 if (alloclhs == NULL)
2252 if (rhs->digits > reqdigits)
2254 allocrhs = decRoundOperand (rhs, set, status);
2255 if (allocrhs == NULL)
2262 /* [following code does not require input rounding] */
2264 /* note whether signs differ */
2265 diffsign = (Flag) ((lhs->bits ^ rhs->bits ^ negate) & DECNEG);
2267 /* handle infinities and NaNs */
2268 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2270 { /* a special bit set */
2271 if (merged & (DECSNAN | DECNAN)) /* a NaN */
2272 decNaNs (res, lhs, rhs, status);
2274 { /* one or two infinities */
2275 if (decNumberIsInfinite (lhs))
2276 { /* LHS is infinity */
2277 /* two infinities with different signs is invalid */
2278 if (decNumberIsInfinite (rhs) && diffsign)
2280 *status |= DEC_Invalid_operation;
2283 bits = lhs->bits & DECNEG; /* get sign from LHS */
2286 bits = (rhs->bits ^ negate) & DECNEG; /* RHS must be Infinity */
2288 decNumberZero (res);
2289 res->bits = bits; /* set +/- infinity */
2294 /* Quick exit for add 0s; return the non-0, modified as need be */
2297 Int adjust; /* work */
2298 Int lexp = lhs->exponent; /* save in case LHS==RES */
2299 bits = lhs->bits; /* .. */
2300 residue = 0; /* clear accumulator */
2301 decCopyFit (res, rhs, set, &residue, status); /* copy (as needed) */
2302 res->bits ^= negate; /* flip if rhs was negated */
2305 { /* exponents on zeros count */
2307 /* exponent will be the lower of the two */
2308 adjust = lexp - res->exponent; /* adjustment needed [if -ve] */
2310 { /* both 0: special IEEE 854 rules */
2312 res->exponent = lexp; /* set exponent */
2313 /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
2316 if (set->round != DEC_ROUND_FLOOR)
2319 res->bits = DECNEG; /* preserve 0 sign */
2325 { /* 0-padding needed */
2326 if ((res->digits - adjust) > set->digits)
2328 adjust = res->digits - set->digits; /* to fit exactly */
2329 *status |= DEC_Rounded; /* [but exact] */
2332 decShiftToMost (res->lsu, res->digits, -adjust);
2333 res->exponent += adjust; /* set the exponent. */
2339 decFinish (res, set, &residue, status); /* clean and finalize */
2344 { /* [lhs is non-zero] */
2345 Int adjust; /* work */
2346 Int rexp = rhs->exponent; /* save in case RHS==RES */
2347 bits = rhs->bits; /* be clean */
2348 residue = 0; /* clear accumulator */
2349 decCopyFit (res, lhs, set, &residue, status); /* copy (as needed) */
2352 { /* exponents on zeros count */
2354 /* exponent will be the lower of the two */
2355 /* [0-0 case handled above] */
2356 adjust = rexp - res->exponent; /* adjustment needed [if -ve] */
2358 { /* 0-padding needed */
2359 if ((res->digits - adjust) > set->digits)
2361 adjust = res->digits - set->digits; /* to fit exactly */
2362 *status |= DEC_Rounded; /* [but exact] */
2365 decShiftToMost (res->lsu, res->digits, -adjust);
2366 res->exponent += adjust; /* set the exponent. */
2371 decFinish (res, set, &residue, status); /* clean and finalize */
2374 /* [both fastpath and mainpath code below assume these cases */
2375 /* (notably 0-0) have already been handled] */
2377 /* calculate the padding needed to align the operands */
2378 padding = rhs->exponent - lhs->exponent;
2380 /* Fastpath cases where the numbers are aligned and normal, the RHS */
2381 /* is all in one unit, no operand rounding is needed, and no carry, */
2382 /* lengthening, or borrow is needed */
2383 if (rhs->digits <= DECDPUN && padding == 0 && rhs->exponent >= set->emin /* [some normals drop through] */
2384 && rhs->digits <= reqdigits && lhs->digits <= reqdigits)
2386 Int partial = *lhs->lsu;
2389 Int maxv = DECDPUNMAX; /* highest no-overflow */
2390 if (lhs->digits < DECDPUN)
2391 maxv = powers[lhs->digits] - 1;
2392 partial += *rhs->lsu;
2393 if (partial <= maxv)
2396 decNumberCopy (res, lhs); /* not in place */
2397 *res->lsu = (Unit) partial; /* [copy could have overwritten RHS] */
2400 /* else drop out for careful add */
2403 { /* signs differ */
2404 partial -= *rhs->lsu;
2406 { /* no borrow needed, and non-0 result */
2408 decNumberCopy (res, lhs); /* not in place */
2409 *res->lsu = (Unit) partial;
2410 /* this could have reduced digits [but result>0] */
2411 res->digits = decGetDigits (res->lsu, D2U (res->digits));
2414 /* else drop out for careful subtract */
2418 /* Now align (pad) the lhs or rhs so we can add or subtract them, as
2419 necessary. If one number is much larger than the other (that is,
2420 if in plain form there is a least one digit between the lowest
2421 digit or one and the highest of the other) we need to pad with up
2422 to DIGITS-1 trailing zeros, and then apply rounding (as exotic
2423 rounding modes may be affected by the residue).
2425 rhsshift = 0; /* rhs shift to left (padding) in Units */
2426 bits = lhs->bits; /* assume sign is that of LHS */
2427 mult = 1; /* likely multiplier */
2429 /* if padding==0 the operands are aligned; no padding needed */
2432 /* some padding needed */
2433 /* We always pad the RHS, as we can then effect any required */
2434 /* padding by a combination of shifts and a multiply */
2437 { /* LHS needs the padding */
2439 padding = -padding; /* will be +ve */
2440 bits = (uByte) (rhs->bits ^ negate); /* assumed sign is now that of RHS */
2447 /* If, after pad, rhs would be longer than lhs by digits+1 or */
2448 /* more then lhs cannot affect the answer, except as a residue, */
2449 /* so we only need to pad up to a length of DIGITS+1. */
2450 if (rhs->digits + padding > lhs->digits + reqdigits + 1)
2452 /* The RHS is sufficient */
2453 /* for residue we use the relative sign indication... */
2454 Int shift = reqdigits - rhs->digits; /* left shift needed */
2455 residue = 1; /* residue for rounding */
2457 residue = -residue; /* signs differ */
2458 /* copy, shortening if necessary */
2459 decCopyFit (res, rhs, set, &residue, status);
2460 /* if it was already shorter, then need to pad with zeros */
2463 res->digits = decShiftToMost (res->lsu, res->digits, shift);
2464 res->exponent -= shift; /* adjust the exponent. */
2466 /* flip the result sign if unswapped and rhs was negated */
2468 res->bits ^= negate;
2469 decFinish (res, set, &residue, status); /* done */
2473 /* LHS digits may affect result */
2474 rhsshift = D2U (padding + 1) - 1; /* this much by Unit shift .. */
2475 mult = powers[padding - (rhsshift * DECDPUN)]; /* .. this by multiplication */
2476 } /* padding needed */
2479 mult = -mult; /* signs differ */
2481 /* determine the longer operand */
2482 maxdigits = rhs->digits + padding; /* virtual length of RHS */
2483 if (lhs->digits > maxdigits)
2484 maxdigits = lhs->digits;
2486 /* Decide on the result buffer to use; if possible place directly */
2488 acc = res->lsu; /* assume build direct */
2489 /* If destructive overlap, or the number is too long, or a carry or */
2490 /* borrow to DIGITS+1 might be possible we must use a buffer. */
2491 /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
2492 if ((maxdigits >= reqdigits) /* is, or could be, too large */
2493 || (res == rhs && rhsshift > 0))
2494 { /* destructive overlap */
2495 /* buffer needed; choose it */
2496 /* we'll need units for maxdigits digits, +1 Unit for carry or borrow */
2497 Int need = D2U (maxdigits) + 1;
2498 acc = accbuff; /* assume use local buffer */
2499 if (need * sizeof (Unit) > sizeof (accbuff))
2501 allocacc = (Unit *) malloc (need * sizeof (Unit));
2502 if (allocacc == NULL)
2503 { /* hopeless -- abandon */
2504 *status |= DEC_Insufficient_storage;
2512 res->bits = (uByte) (bits & DECNEG); /* it's now safe to overwrite.. */
2513 res->exponent = lhs->exponent; /* .. operands (even if aliased) */
2516 decDumpAr ('A', lhs->lsu, D2U (lhs->digits));
2517 decDumpAr ('B', rhs->lsu, D2U (rhs->digits));
2518 printf (" :h: %d %d\n", rhsshift, mult);
2521 /* add [A+B*m] or subtract [A+B*(-m)] */
2522 res->digits = decUnitAddSub (lhs->lsu, D2U (lhs->digits), rhs->lsu, D2U (rhs->digits), rhsshift, acc, mult) * DECDPUN; /* [units -> digits] */
2523 if (res->digits < 0)
2525 res->digits = -res->digits;
2526 res->bits ^= DECNEG; /* flip the sign */
2529 decDumpAr ('+', acc, D2U (res->digits));
2532 /* If we used a buffer we need to copy back, possibly shortening */
2533 /* (If we didn't use buffer it must have fit, so can't need rounding */
2534 /* and residue must be 0.) */
2535 residue = 0; /* clear accumulator */
2536 if (acc != res->lsu)
2540 { /* round from first significant digit */
2542 /* remove leading zeros that we added due to rounding up to */
2543 /* integral Units -- before the test for rounding. */
2544 if (res->digits > reqdigits)
2545 res->digits = decGetDigits (acc, D2U (res->digits));
2546 decSetCoeff (res, set, acc, res->digits, &residue, status);
2550 { /* subset arithmetic rounds from original significant digit */
2551 /* We may have an underestimate. This only occurs when both */
2552 /* numbers fit in DECDPUN digits and we are padding with a */
2553 /* negative multiple (-10, -100...) and the top digit(s) become */
2554 /* 0. (This only matters if we are using X3.274 rules where the */
2555 /* leading zero could be included in the rounding.) */
2556 if (res->digits < maxdigits)
2558 *(acc + D2U (res->digits)) = 0; /* ensure leading 0 is there */
2559 res->digits = maxdigits;
2563 /* remove leading zeros that we added due to rounding up to */
2564 /* integral Units (but only those in excess of the original */
2565 /* maxdigits length, unless extended) before test for rounding. */
2566 if (res->digits > reqdigits)
2568 res->digits = decGetDigits (acc, D2U (res->digits));
2569 if (res->digits < maxdigits)
2570 res->digits = maxdigits;
2573 decSetCoeff (res, set, acc, res->digits, &residue, status);
2574 /* Now apply rounding if needed before removing leading zeros. */
2575 /* This is safe because subnormals are not a possibility */
2578 decApplyRound (res, set, residue, status);
2579 residue = 0; /* we did what we had to do */
2585 /* strip leading zeros [these were left on in case of subset subtract] */
2586 res->digits = decGetDigits (res->lsu, D2U (res->digits));
2588 /* apply checks and rounding */
2589 decFinish (res, set, &residue, status);
2591 /* "When the sum of two operands with opposite signs is exactly */
2592 /* zero, the sign of that sum shall be '+' in all rounding modes */
2593 /* except round toward -Infinity, in which mode that sign shall be */
2594 /* '-'." [Subset zeros also never have '-', set by decFinish.] */
2595 if (ISZERO (res) && diffsign
2599 && (*status & DEC_Inexact) == 0)
2601 if (set->round == DEC_ROUND_FLOOR)
2602 res->bits |= DECNEG; /* sign - */
2604 res->bits &= ~DECNEG; /* sign + */
2607 while (0); /* end protected */
2611 if (allocacc != NULL)
2612 free (allocacc); /* drop any storage we used */
2613 if (allocrhs != NULL)
2614 free (allocrhs); /* .. */
2615 if (alloclhs != NULL)
2616 free (alloclhs); /* .. */
2621 /* ------------------------------------------------------------------ */
2622 /* decDivideOp -- division operation */
2624 /* This routine performs the calculations for all four division */
2625 /* operators (divide, divideInteger, remainder, remainderNear). */
2629 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
2632 /* set is the context */
2633 /* op is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively. */
2634 /* status is the usual accumulator */
2636 /* C must have space for set->digits digits. */
2638 /* ------------------------------------------------------------------ */
2639 /* The underlying algorithm of this routine is the same as in the */
2640 /* 1981 S/370 implementation, that is, non-restoring long division */
2641 /* with bi-unit (rather than bi-digit) estimation for each unit */
2642 /* multiplier. In this pseudocode overview, complications for the */
2643 /* Remainder operators and division residues for exact rounding are */
2644 /* omitted for clarity. */
2646 /* Prepare operands and handle special values */
2647 /* Test for x/0 and then 0/x */
2648 /* Exp =Exp1 - Exp2 */
2649 /* Exp =Exp +len(var1) -len(var2) */
2650 /* Sign=Sign1 * Sign2 */
2651 /* Pad accumulator (Var1) to double-length with 0's (pad1) */
2652 /* Pad Var2 to same length as Var1 */
2653 /* msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round */
2655 /* Do until (have=digits+1 OR residue=0) */
2656 /* if exp<0 then if integer divide/residue then leave */
2659 /* compare numbers */
2660 /* if <0 then leave inner_loop */
2661 /* if =0 then (* quick exit without subtract *) do */
2662 /* this_unit=this_unit+1; output this_unit */
2663 /* leave outer_loop; end */
2664 /* Compare lengths of numbers (mantissae): */
2665 /* If same then tops2=msu2pair -- {units 1&2 of var2} */
2666 /* else tops2=msu2plus -- {0, unit 1 of var2} */
2667 /* tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
2668 /* mult=tops1/tops2 -- Good and safe guess at divisor */
2669 /* if mult=0 then mult=1 */
2670 /* this_unit=this_unit+mult */
2672 /* end inner_loop */
2673 /* if have\=0 | this_unit\=0 then do */
2674 /* output this_unit */
2675 /* have=have+1; end */
2678 /* end outer_loop */
2679 /* exp=exp+1 -- set the proper exponent */
2680 /* if have=0 then generate answer=0 */
2681 /* Return (Result is defined by Var1) */
2683 /* ------------------------------------------------------------------ */
2684 /* We need two working buffers during the long division; one (digits+ */
2685 /* 1) to accumulate the result, and the other (up to 2*digits+1) for */
2686 /* long subtractions. These are acc and var1 respectively. */
2687 /* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
2688 /* ------------------------------------------------------------------ */
2690 decDivideOp (decNumber * res,
2691 decNumber * lhs, decNumber * rhs,
2692 decContext * set, Flag op, uInt * status)
2694 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
2695 decNumber *allocrhs = NULL; /* .., rhs */
2696 Unit accbuff[D2U (DECBUFFER + DECDPUN)]; /* local buffer */
2697 Unit *acc = accbuff; /* -> accumulator array for result */
2698 Unit *allocacc = NULL; /* -> allocated buffer, iff allocated */
2699 Unit *accnext; /* -> where next digit will go */
2700 Int acclength; /* length of acc needed [Units] */
2701 Int accunits; /* count of units accumulated */
2702 Int accdigits; /* count of digits accumulated */
2704 Unit varbuff[D2U (DECBUFFER * 2 + DECDPUN) * sizeof (Unit)]; /* buffer for var1 */
2705 Unit *var1 = varbuff; /* -> var1 array for long subtraction */
2706 Unit *varalloc = NULL; /* -> allocated buffer, iff used */
2708 Unit *var2; /* -> var2 array */
2710 Int var1units, var2units; /* actual lengths */
2711 Int var2ulen; /* logical length (units) */
2712 Int var1initpad = 0; /* var1 initial padding (digits) */
2713 Unit *msu1, *msu2; /* -> msu of each var */
2714 Int msu2plus; /* msu2 plus one [does not vary] */
2715 eInt msu2pair; /* msu2 pair plus one [does not vary] */
2716 Int maxdigits; /* longest LHS or required acc length */
2717 Int mult; /* multiplier for subtraction */
2718 Unit thisunit; /* current unit being accumulated */
2719 Int residue; /* for rounding */
2720 Int reqdigits = set->digits; /* requested DIGITS */
2721 Int exponent; /* working exponent */
2722 Int maxexponent = 0; /* DIVIDE maximum exponent if unrounded */
2723 uByte bits; /* working sign */
2724 uByte merged; /* merged flags */
2725 Unit *target, *source; /* work */
2726 uInt const *pow; /* .. */
2727 Int shift, cut; /* .. */
2729 Int dropped; /* work */
2733 if (decCheckOperands (res, lhs, rhs, set))
2738 { /* protect allocated storage */
2742 /* reduce operands and set lostDigits status, as needed */
2743 if (lhs->digits > reqdigits)
2745 alloclhs = decRoundOperand (lhs, set, status);
2746 if (alloclhs == NULL)
2750 if (rhs->digits > reqdigits)
2752 allocrhs = decRoundOperand (rhs, set, status);
2753 if (allocrhs == NULL)
2759 /* [following code does not require input rounding] */
2761 bits = (lhs->bits ^ rhs->bits) & DECNEG; /* assumed sign for divisions */
2763 /* handle infinities and NaNs */
2764 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2766 { /* a special bit set */
2767 if (merged & (DECSNAN | DECNAN))
2768 { /* one or two NaNs */
2769 decNaNs (res, lhs, rhs, status);
2772 /* one or two infinities */
2773 if (decNumberIsInfinite (lhs))
2774 { /* LHS (dividend) is infinite */
2775 if (decNumberIsInfinite (rhs) || /* two infinities are invalid .. */
2776 op & (REMAINDER | REMNEAR))
2777 { /* as is remainder of infinity */
2778 *status |= DEC_Invalid_operation;
2781 /* [Note that infinity/0 raises no exceptions] */
2782 decNumberZero (res);
2783 res->bits = bits | DECINF; /* set +/- infinity */
2787 { /* RHS (divisor) is infinite */
2789 if (op & (REMAINDER | REMNEAR))
2791 /* result is [finished clone of] lhs */
2792 decCopyFit (res, lhs, set, &residue, status);
2796 decNumberZero (res);
2797 res->bits = bits; /* set +/- zero */
2798 /* for DIVIDEINT the exponent is always 0. For DIVIDE, result */
2799 /* is a 0 with infinitely negative exponent, clamped to minimum */
2802 res->exponent = set->emin - set->digits + 1;
2803 *status |= DEC_Clamped;
2806 decFinish (res, set, &residue, status);
2811 /* handle 0 rhs (x/0) */
2813 { /* x/0 is always exceptional */
2816 decNumberZero (res); /* [after lhs test] */
2817 *status |= DEC_Division_undefined; /* 0/0 will become NaN */
2821 decNumberZero (res);
2822 if (op & (REMAINDER | REMNEAR))
2823 *status |= DEC_Invalid_operation;
2826 *status |= DEC_Division_by_zero; /* x/0 */
2827 res->bits = bits | DECINF; /* .. is +/- Infinity */
2833 /* handle 0 lhs (0/x) */
2838 decNumberZero (res);
2845 exponent = lhs->exponent - rhs->exponent; /* ideal exponent */
2846 decNumberCopy (res, lhs); /* [zeros always fit] */
2847 res->bits = bits; /* sign as computed */
2848 res->exponent = exponent; /* exponent, too */
2849 decFinalize (res, set, &residue, status); /* check exponent */
2851 else if (op & DIVIDEINT)
2853 decNumberZero (res); /* integer 0 */
2854 res->bits = bits; /* sign as computed */
2858 exponent = rhs->exponent; /* [save in case overwrite] */
2859 decNumberCopy (res, lhs); /* [zeros always fit] */
2860 if (exponent < res->exponent)
2861 res->exponent = exponent; /* use lower */
2869 /* Precalculate exponent. This starts off adjusted (and hence fits */
2870 /* in 31 bits) and becomes the usual unadjusted exponent as the */
2871 /* division proceeds. The order of evaluation is important, here, */
2872 /* to avoid wrap. */
2874 (lhs->exponent + lhs->digits) - (rhs->exponent + rhs->digits);
2876 /* If the working exponent is -ve, then some quick exits are */
2877 /* possible because the quotient is known to be <1 */
2878 /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
2879 if (exponent < 0 && !(op == DIVIDE))
2883 decNumberZero (res); /* integer part is 0 */
2887 res->bits = bits; /* set +/- zero */
2890 /* we can fastpath remainders so long as the lhs has the */
2891 /* smaller (or equal) exponent */
2892 if (lhs->exponent <= rhs->exponent)
2894 if (op & REMAINDER || exponent < -1)
2896 /* It is REMAINDER or safe REMNEAR; result is [finished */
2897 /* clone of] lhs (r = x - 0*y) */
2899 decCopyFit (res, lhs, set, &residue, status);
2900 decFinish (res, set, &residue, status);
2903 /* [unsafe REMNEAR drops through] */
2907 /* We need long (slow) division; roll up the sleeves... */
2909 /* The accumulator will hold the quotient of the division. */
2910 /* If it needs to be too long for stack storage, then allocate. */
2911 acclength = D2U (reqdigits + DECDPUN); /* in Units */
2912 if (acclength * sizeof (Unit) > sizeof (accbuff))
2914 allocacc = (Unit *) malloc (acclength * sizeof (Unit));
2915 if (allocacc == NULL)
2916 { /* hopeless -- abandon */
2917 *status |= DEC_Insufficient_storage;
2920 acc = allocacc; /* use the allocated space */
2923 /* var1 is the padded LHS ready for subtractions. */
2924 /* If it needs to be too long for stack storage, then allocate. */
2925 /* The maximum units we need for var1 (long subtraction) is: */
2927 /* (rhs->digits+reqdigits-1) -- to allow full slide to right */
2928 /* or (lhs->digits) -- to allow for long lhs */
2929 /* whichever is larger */
2930 /* +1 -- for rounding of slide to right */
2931 /* +1 -- for leading 0s */
2932 /* +1 -- for pre-adjust if a remainder or DIVIDEINT */
2933 /* [Note: unused units do not participate in decUnitAddSub data] */
2934 maxdigits = rhs->digits + reqdigits - 1;
2935 if (lhs->digits > maxdigits)
2936 maxdigits = lhs->digits;
2937 var1units = D2U (maxdigits) + 2;
2938 /* allocate a guard unit above msu1 for REMAINDERNEAR */
2941 if ((var1units + 1) * sizeof (Unit) > sizeof (varbuff))
2943 varalloc = (Unit *) malloc ((var1units + 1) * sizeof (Unit));
2944 if (varalloc == NULL)
2945 { /* hopeless -- abandon */
2946 *status |= DEC_Insufficient_storage;
2949 var1 = varalloc; /* use the allocated space */
2952 /* Extend the lhs and rhs to full long subtraction length. The lhs */
2953 /* is truly extended into the var1 buffer, with 0 padding, so we can */
2954 /* subtract in place. The rhs (var2) has virtual padding */
2955 /* (implemented by decUnitAddSub). */
2956 /* We allocated one guard unit above msu1 for rem=rem+rem in REMAINDERNEAR */
2957 msu1 = var1 + var1units - 1; /* msu of var1 */
2958 source = lhs->lsu + D2U (lhs->digits) - 1; /* msu of input array */
2959 for (target = msu1; source >= lhs->lsu; source--, target--)
2961 for (; target >= var1; target--)
2964 /* rhs (var2) is left-aligned with var1 at the start */
2965 var2ulen = var1units; /* rhs logical length (units) */
2966 var2units = D2U (rhs->digits); /* rhs actual length (units) */
2967 var2 = rhs->lsu; /* -> rhs array */
2968 msu2 = var2 + var2units - 1; /* -> msu of var2 [never changes] */
2969 /* now set up the variables which we'll use for estimating the */
2970 /* multiplication factor. If these variables are not exact, we add */
2971 /* 1 to make sure that we never overestimate the multiplier. */
2972 msu2plus = *msu2; /* it's value .. */
2974 msu2plus++; /* .. +1 if any more */
2975 msu2pair = (eInt) * msu2 * (DECDPUNMAX + 1); /* top two pair .. */
2977 { /* .. [else treat 2nd as 0] */
2978 msu2pair += *(msu2 - 1); /* .. */
2980 msu2pair++; /* .. +1 if any more */
2983 /* Since we are working in units, the units may have leading zeros, */
2984 /* but we calculated the exponent on the assumption that they are */
2985 /* both left-aligned. Adjust the exponent to compensate: add the */
2986 /* number of leading zeros in var1 msu and subtract those in var2 msu. */
2987 /* [We actually do this by counting the digits and negating, as */
2988 /* lead1=DECDPUN-digits1, and similarly for lead2.] */
2989 for (pow = &powers[1]; *msu1 >= *pow; pow++)
2991 for (pow = &powers[1]; *msu2 >= *pow; pow++)
2994 /* Now, if doing an integer divide or remainder, we want to ensure */
2995 /* that the result will be Unit-aligned. To do this, we shift the */
2996 /* var1 accumulator towards least if need be. (It's much easier to */
2997 /* do this now than to reassemble the residue afterwards, if we are */
2998 /* doing a remainder.) Also ensure the exponent is not negative. */
3002 /* save the initial 'false' padding of var1, in digits */
3003 var1initpad = (var1units - D2U (lhs->digits)) * DECDPUN;
3004 /* Determine the shift to do. */
3008 cut = DECDPUN - exponent % DECDPUN;
3009 decShiftToLeast (var1, var1units, cut);
3010 exponent += cut; /* maintain numerical value */
3011 var1initpad -= cut; /* .. and reduce padding */
3012 /* clean any most-significant units we just emptied */
3013 for (u = msu1; cut >= DECDPUN; cut -= DECDPUN, u--)
3018 maxexponent = lhs->exponent - rhs->exponent; /* save */
3019 /* optimization: if the first iteration will just produce 0, */
3020 /* preadjust to skip it [valid for DIVIDE only] */
3023 var2ulen--; /* shift down */
3024 exponent -= DECDPUN; /* update the exponent */
3028 /* ---- start the long-division loops ------------------------------ */
3029 accunits = 0; /* no units accumulated yet */
3030 accdigits = 0; /* .. or digits */
3031 accnext = acc + acclength - 1; /* -> msu of acc [NB: allows digits+1] */
3033 { /* outer forever loop */
3034 thisunit = 0; /* current unit assumed 0 */
3035 /* find the next unit */
3037 { /* inner forever loop */
3038 /* strip leading zero units [from either pre-adjust or from */
3039 /* subtract last time around]. Leave at least one unit. */
3040 for (; *msu1 == 0 && msu1 > var1; msu1--)
3043 if (var1units < var2ulen)
3044 break; /* var1 too low for subtract */
3045 if (var1units == var2ulen)
3046 { /* unit-by-unit compare needed */
3047 /* compare the two numbers, from msu */
3048 Unit *pv1, *pv2, v2; /* units to compare */
3049 pv2 = msu2; /* -> msu */
3050 for (pv1 = msu1;; pv1--, pv2--)
3052 /* v1=*pv1 -- always OK */
3053 v2 = 0; /* assume in padding */
3055 v2 = *pv2; /* in range */
3057 break; /* no longer the same */
3059 break; /* done; leave pv1 as is */
3061 /* here when all inspected or a difference seen */
3063 break; /* var1 too low to subtract */
3065 { /* var1 == var2 */
3066 /* reach here if var1 and var2 are identical; subtraction */
3067 /* would increase digit by one, and the residue will be 0 so */
3068 /* we are done; leave the loop with residue set to 0. */
3069 thisunit++; /* as though subtracted */
3070 *var1 = 0; /* set var1 to 0 */
3071 var1units = 1; /* .. */
3072 break; /* from inner */
3073 } /* var1 == var2 */
3074 /* *pv1>v2. Prepare for real subtraction; the lengths are equal */
3075 /* Estimate the multiplier (there's always a msu1-1)... */
3076 /* Bring in two units of var2 to provide a good estimate. */
3078 (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3079 *(msu1 - 1)) / msu2pair);
3080 } /* lengths the same */
3082 { /* var1units > var2ulen, so subtraction is safe */
3083 /* The var2 msu is one unit towards the lsu of the var1 msu, */
3084 /* so we can only use one unit for var2. */
3086 (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3087 *(msu1 - 1)) / msu2plus);
3090 mult = 1; /* must always be at least 1 */
3091 /* subtraction needed; var1 is > var2 */
3092 thisunit = (Unit) (thisunit + mult); /* accumulate */
3093 /* subtract var1-var2, into var1; only the overlap needs */
3094 /* processing, as we are in place */
3095 shift = var2ulen - var2units;
3097 decDumpAr ('1', &var1[shift], var1units - shift);
3098 decDumpAr ('2', var2, var2units);
3099 printf ("m=%d\n", -mult);
3101 decUnitAddSub (&var1[shift], var1units - shift,
3102 var2, var2units, 0, &var1[shift], -mult);
3104 decDumpAr ('#', &var1[shift], var1units - shift);
3106 /* var1 now probably has leading zeros; these are removed at the */
3107 /* top of the inner loop. */
3110 /* We have the next unit; unless it's a leading zero, add to acc */
3111 if (accunits != 0 || thisunit != 0)
3112 { /* put the unit we got */
3113 *accnext = thisunit; /* store in accumulator */
3114 /* account exactly for the digits we got */
3117 accdigits++; /* at least one */
3118 for (pow = &powers[1]; thisunit >= *pow; pow++)
3122 accdigits += DECDPUN;
3123 accunits++; /* update count */
3124 accnext--; /* ready for next */
3125 if (accdigits > reqdigits)
3126 break; /* we have all we need */
3129 /* if the residue is zero, we're done (unless divide or */
3130 /* divideInteger and we haven't got enough digits yet) */
3131 if (*var1 == 0 && var1units == 1)
3132 { /* residue is 0 */
3133 if (op & (REMAINDER | REMNEAR))
3135 if ((op & DIVIDE) && (exponent <= maxexponent))
3137 /* [drop through if divideInteger] */
3139 /* we've also done enough if calculating remainder or integer */
3140 /* divide and we just did the last ('units') unit */
3141 if (exponent == 0 && !(op & DIVIDE))
3144 /* to get here, var1 is less than var2, so divide var2 by the per- */
3145 /* Unit power of ten and go for the next digit */
3146 var2ulen--; /* shift down */
3147 exponent -= DECDPUN; /* update the exponent */
3150 /* ---- division is complete --------------------------------------- */
3151 /* here: acc has at least reqdigits+1 of good results (or fewer */
3152 /* if early stop), starting at accnext+1 (its lsu) */
3153 /* var1 has any residue at the stopping point */
3154 /* accunits is the number of digits we collected in acc */
3157 accunits = 1; /* show we have one .. */
3158 accdigits = 1; /* .. */
3159 *accnext = 0; /* .. whose value is 0 */
3162 accnext++; /* back to last placed */
3163 /* accnext now -> lowest unit of result */
3165 residue = 0; /* assume no residue */
3168 /* record the presence of any residue, for rounding */
3169 if (*var1 != 0 || var1units > 1)
3173 /* We had an exact division; clean up spurious trailing 0s. */
3174 /* There will be at most DECDPUN-1, from the final multiply, */
3175 /* and then only if the result is non-0 (and even) and the */
3176 /* exponent is 'loose'. */
3178 Unit lsu = *accnext;
3179 if (!(lsu & 0x01) && (lsu != 0))
3181 /* count the trailing zeros */
3184 { /* [will terminate because lsu!=0] */
3185 if (exponent >= maxexponent)
3186 break; /* don't chop real 0s */
3188 if ((lsu - QUOT10 (lsu, drop + 1)
3189 * powers[drop + 1]) != 0)
3190 break; /* found non-0 digit */
3192 if (lsu % powers[drop + 1] != 0)
3193 break; /* found non-0 digit */
3199 accunits = decShiftToLeast (accnext, accunits, drop);
3200 accdigits = decGetDigits (accnext, accunits);
3201 accunits = D2U (accdigits);
3202 /* [exponent was adjusted in the loop] */
3204 } /* neither odd nor 0 */
3206 } /* exact divide */
3208 else /* op!=DIVIDE */
3210 /* check for coefficient overflow */
3211 if (accdigits + exponent > reqdigits)
3213 *status |= DEC_Division_impossible;
3216 if (op & (REMAINDER | REMNEAR))
3218 /* [Here, the exponent will be 0, because we adjusted var1 */
3219 /* appropriately.] */
3220 Int postshift; /* work */
3221 Flag wasodd = 0; /* integer was odd */
3222 Unit *quotlsu; /* for save */
3223 Int quotdigits; /* .. */
3225 /* Fastpath when residue is truly 0 is worthwhile [and */
3226 /* simplifies the code below] */
3227 if (*var1 == 0 && var1units == 1)
3228 { /* residue is 0 */
3229 Int exp = lhs->exponent; /* save min(exponents) */
3230 if (rhs->exponent < exp)
3231 exp = rhs->exponent;
3232 decNumberZero (res); /* 0 coefficient */
3236 res->exponent = exp; /* .. with proper exponent */
3239 /* note if the quotient was odd */
3240 if (*accnext & 0x01)
3241 wasodd = 1; /* acc is odd */
3242 quotlsu = accnext; /* save in case need to reinspect */
3243 quotdigits = accdigits; /* .. */
3245 /* treat the residue, in var1, as the value to return, via acc */
3246 /* calculate the unused zero digits. This is the smaller of: */
3247 /* var1 initial padding (saved above) */
3248 /* var2 residual padding, which happens to be given by: */
3250 var1initpad + exponent - lhs->exponent + rhs->exponent;
3251 /* [the 'exponent' term accounts for the shifts during divide] */
3252 if (var1initpad < postshift)
3253 postshift = var1initpad;
3255 /* shift var1 the requested amount, and adjust its digits */
3256 var1units = decShiftToLeast (var1, var1units, postshift);
3258 accdigits = decGetDigits (var1, var1units);
3259 accunits = D2U (accdigits);
3261 exponent = lhs->exponent; /* exponent is smaller of lhs & rhs */
3262 if (rhs->exponent < exponent)
3263 exponent = rhs->exponent;
3264 bits = lhs->bits; /* remainder sign is always as lhs */
3266 /* Now correct the result if we are doing remainderNear; if it */
3267 /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
3268 /* the integer was odd then the result should be rem-rhs. */
3271 Int compare, tarunits; /* work */
3275 /* calculate remainder*2 into the var1 buffer (which has */
3276 /* 'headroom' of an extra unit and hence enough space) */
3277 /* [a dedicated 'double' loop would be faster, here] */
3279 decUnitAddSub (accnext, accunits, accnext, accunits, 0,
3281 /* decDumpAr('r', accnext, tarunits); */
3283 /* Here, accnext (var1) holds tarunits Units with twice the */
3284 /* remainder's coefficient, which we must now compare to the */
3285 /* RHS. The remainder's exponent may be smaller than the RHS's. */
3287 decUnitCompare (accnext, tarunits, rhs->lsu,
3289 rhs->exponent - exponent);
3290 if (compare == BADINT)
3291 { /* deep trouble */
3292 *status |= DEC_Insufficient_storage;
3296 /* now restore the remainder by dividing by two; we know the */
3298 for (up = accnext; up < accnext + tarunits; up++)
3300 Int half; /* half to add to lower unit */
3302 *up /= 2; /* [shift] */
3305 *(up - 1) += (DECDPUNMAX + 1) / 2;
3307 /* [accunits still describes the original remainder length] */
3309 if (compare > 0 || (compare == 0 && wasodd))
3310 { /* adjustment needed */
3311 Int exp, expunits, exprem; /* work */
3312 /* This is effectively causing round-up of the quotient, */
3313 /* so if it was the rare case where it was full and all */
3314 /* nines, it would overflow and hence division-impossible */
3315 /* should be raised */
3316 Flag allnines = 0; /* 1 if quotient all nines */
3317 if (quotdigits == reqdigits)
3318 { /* could be borderline */
3319 for (up = quotlsu;; up++)
3321 if (quotdigits > DECDPUN)
3323 if (*up != DECDPUNMAX)
3324 break; /* non-nines */
3327 { /* this is the last Unit */
3328 if (*up == powers[quotdigits] - 1)
3332 quotdigits -= DECDPUN; /* checked those digits */
3334 } /* borderline check */
3337 *status |= DEC_Division_impossible;
3341 /* we need rem-rhs; the sign will invert. Again we can */
3342 /* safely use var1 for the working Units array. */
3343 exp = rhs->exponent - exponent; /* RHS padding needed */
3344 /* Calculate units and remainder from exponent. */
3345 expunits = exp / DECDPUN;
3346 exprem = exp % DECDPUN;
3347 /* subtract [A+B*(-m)]; the result will always be negative */
3348 accunits = -decUnitAddSub (accnext, accunits,
3349 rhs->lsu, D2U (rhs->digits),
3351 -(Int) powers[exprem]);
3352 accdigits = decGetDigits (accnext, accunits); /* count digits exactly */
3353 accunits = D2U (accdigits); /* and recalculate the units for copy */
3354 /* [exponent is as for original remainder] */
3355 bits ^= DECNEG; /* flip the sign */
3358 } /* REMAINDER or REMNEAR */
3361 /* Set exponent and bits */
3362 res->exponent = exponent;
3363 res->bits = (uByte) (bits & DECNEG); /* [cleaned] */
3365 /* Now the coefficient. */
3366 decSetCoeff (res, set, accnext, accdigits, &residue, status);
3368 decFinish (res, set, &residue, status); /* final cleanup */
3371 /* If a divide then strip trailing zeros if subset [after round] */
3372 if (!set->extended && (op == DIVIDE))
3373 decTrim (res, 0, &dropped);
3376 while (0); /* end protected */
3378 if (varalloc != NULL)
3379 free (varalloc); /* drop any storage we used */
3380 if (allocacc != NULL)
3381 free (allocacc); /* .. */
3382 if (allocrhs != NULL)
3383 free (allocrhs); /* .. */
3384 if (alloclhs != NULL)
3385 free (alloclhs); /* .. */
3389 /* ------------------------------------------------------------------ */
3390 /* decMultiplyOp -- multiplication operation */
3392 /* This routine performs the multiplication C=A x B. */
3394 /* res is C, the result. C may be A and/or B (e.g., X=X*X) */
3397 /* set is the context */
3398 /* status is the usual accumulator */
3400 /* C must have space for set->digits digits. */
3402 /* ------------------------------------------------------------------ */
3403 /* Note: We use 'long' multiplication rather than Karatsuba, as the */
3404 /* latter would give only a minor improvement for the short numbers */
3405 /* we expect to handle most (and uses much more memory). */
3407 /* We always have to use a buffer for the accumulator. */
3408 /* ------------------------------------------------------------------ */
3410 decMultiplyOp (decNumber * res, decNumber * lhs,
3411 decNumber * rhs, decContext * set, uInt * status)
3413 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3414 decNumber *allocrhs = NULL; /* .., rhs */
3415 Unit accbuff[D2U (DECBUFFER * 2 + 1)]; /* local buffer (+1 in case DECBUFFER==0) */
3416 Unit *acc = accbuff; /* -> accumulator array for exact result */
3417 Unit *allocacc = NULL; /* -> allocated buffer, iff allocated */
3418 Unit *mer, *mermsup; /* work */
3419 Int accunits; /* Units of accumulator in use */
3420 Int madlength; /* Units in multiplicand */
3421 Int shift; /* Units to shift multiplicand by */
3422 Int need; /* Accumulator units needed */
3423 Int exponent; /* work */
3424 Int residue = 0; /* rounding residue */
3425 uByte bits; /* result sign */
3426 uByte merged; /* merged flags */
3429 if (decCheckOperands (res, lhs, rhs, set))
3434 { /* protect allocated storage */
3438 /* reduce operands and set lostDigits status, as needed */
3439 if (lhs->digits > set->digits)
3441 alloclhs = decRoundOperand (lhs, set, status);
3442 if (alloclhs == NULL)
3446 if (rhs->digits > set->digits)
3448 allocrhs = decRoundOperand (rhs, set, status);
3449 if (allocrhs == NULL)
3455 /* [following code does not require input rounding] */
3457 /* precalculate result sign */
3458 bits = (uByte) ((lhs->bits ^ rhs->bits) & DECNEG);
3460 /* handle infinities and NaNs */
3461 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3463 { /* a special bit set */
3464 if (merged & (DECSNAN | DECNAN))
3465 { /* one or two NaNs */
3466 decNaNs (res, lhs, rhs, status);
3469 /* one or two infinities. Infinity * 0 is invalid */
3470 if (((lhs->bits & DECSPECIAL) == 0 && ISZERO (lhs))
3471 || ((rhs->bits & DECSPECIAL) == 0 && ISZERO (rhs)))
3473 *status |= DEC_Invalid_operation;
3476 decNumberZero (res);
3477 res->bits = bits | DECINF; /* infinity */
3481 /* For best speed, as in DMSRCN, we use the shorter number as the */
3482 /* multiplier (rhs) and the longer as the multiplicand (lhs) */
3483 if (lhs->digits < rhs->digits)
3485 decNumber *hold = lhs;
3490 /* if accumulator is too long for local storage, then allocate */
3491 need = D2U (lhs->digits) + D2U (rhs->digits); /* maximum units in result */
3492 if (need * sizeof (Unit) > sizeof (accbuff))
3494 allocacc = (Unit *) malloc (need * sizeof (Unit));
3495 if (allocacc == NULL)
3497 *status |= DEC_Insufficient_storage;
3500 acc = allocacc; /* use the allocated space */
3503 /* Now the main long multiplication loop */
3504 /* Unlike the equivalent in the IBM Java implementation, there */
3505 /* is no advantage in calculating from msu to lsu. So we do it */
3506 /* by the book, as it were. */
3507 /* Each iteration calculates ACC=ACC+MULTAND*MULT */
3508 accunits = 1; /* accumulator starts at '0' */
3509 *acc = 0; /* .. (lsu=0) */
3510 shift = 0; /* no multiplicand shift at first */
3511 madlength = D2U (lhs->digits); /* we know this won't change */
3512 mermsup = rhs->lsu + D2U (rhs->digits); /* -> msu+1 of multiplier */
3514 for (mer = rhs->lsu; mer < mermsup; mer++)
3516 /* Here, *mer is the next Unit in the multiplier to use */
3517 /* If non-zero [optimization] add it... */
3521 decUnitAddSub (&acc[shift], accunits - shift, lhs->lsu,
3522 madlength, 0, &acc[shift], *mer) + shift;
3525 { /* extend acc with a 0; we'll use it shortly */
3526 /* [this avoids length of <=0 later] */
3527 *(acc + accunits) = 0;
3530 /* multiply multiplicand by 10**DECDPUN for next Unit to left */
3531 shift++; /* add this for 'logical length' */
3534 /* Show exact result */
3535 decDumpAr ('*', acc, accunits);
3538 /* acc now contains the exact result of the multiplication */
3539 /* Build a decNumber from it, noting if any residue */
3540 res->bits = bits; /* set sign */
3541 res->digits = decGetDigits (acc, accunits); /* count digits exactly */
3543 /* We might have a 31-bit wrap in calculating the exponent. */
3544 /* This can only happen if both input exponents are negative and */
3545 /* both their magnitudes are large. If we did wrap, we set a safe */
3546 /* very negative exponent, from which decFinalize() will raise a */
3547 /* hard underflow. */
3548 exponent = lhs->exponent + rhs->exponent; /* calculate exponent */
3549 if (lhs->exponent < 0 && rhs->exponent < 0 && exponent > 0)
3550 exponent = -2 * DECNUMMAXE; /* force underflow */
3551 res->exponent = exponent; /* OK to overwrite now */
3553 /* Set the coefficient. If any rounding, residue records */
3554 decSetCoeff (res, set, acc, res->digits, &residue, status);
3556 decFinish (res, set, &residue, status); /* final cleanup */
3558 while (0); /* end protected */
3560 if (allocacc != NULL)
3561 free (allocacc); /* drop any storage we used */
3562 if (allocrhs != NULL)
3563 free (allocrhs); /* .. */
3564 if (alloclhs != NULL)
3565 free (alloclhs); /* .. */
3569 /* ------------------------------------------------------------------ */
3570 /* decQuantizeOp -- force exponent to requested value */
3572 /* This computes C = op(A, B), where op adjusts the coefficient */
3573 /* of C (by rounding or shifting) such that the exponent (-scale) */
3574 /* of C has the value B or matches the exponent of B. */
3575 /* The numerical value of C will equal A, except for the effects of */
3576 /* any rounding that occurred. */
3578 /* res is C, the result. C may be A or B */
3579 /* lhs is A, the number to adjust */
3580 /* rhs is B, the requested exponent */
3581 /* set is the context */
3582 /* quant is 1 for quantize or 0 for rescale */
3583 /* status is the status accumulator (this can be called without */
3584 /* risk of control loss) */
3586 /* C must have space for set->digits digits. */
3588 /* Unless there is an error or the result is infinite, the exponent */
3589 /* after the operation is guaranteed to be that requested. */
3590 /* ------------------------------------------------------------------ */
3592 decQuantizeOp (decNumber * res, decNumber * lhs,
3593 decNumber * rhs, decContext * set, Flag quant, uInt * status)
3595 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3596 decNumber *allocrhs = NULL; /* .., rhs */
3597 decNumber *inrhs = rhs; /* save original rhs */
3598 Int reqdigits = set->digits; /* requested DIGITS */
3599 Int reqexp; /* requested exponent [-scale] */
3600 Int residue = 0; /* rounding residue */
3601 uByte merged; /* merged flags */
3602 Int etiny = set->emin - (set->digits - 1);
3605 if (decCheckOperands (res, lhs, rhs, set))
3610 { /* protect allocated storage */
3614 /* reduce operands and set lostDigits status, as needed */
3615 if (lhs->digits > reqdigits)
3617 alloclhs = decRoundOperand (lhs, set, status);
3618 if (alloclhs == NULL)
3622 if (rhs->digits > reqdigits)
3623 { /* [this only checks lostDigits] */
3624 allocrhs = decRoundOperand (rhs, set, status);
3625 if (allocrhs == NULL)
3631 /* [following code does not require input rounding] */
3633 /* Handle special values */
3634 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3635 if ((lhs->bits | rhs->bits) & DECSPECIAL)
3637 /* NaNs get usual processing */
3638 if (merged & (DECSNAN | DECNAN))
3639 decNaNs (res, lhs, rhs, status);
3640 /* one infinity but not both is bad */
3641 else if ((lhs->bits ^ rhs->bits) & DECINF)
3642 *status |= DEC_Invalid_operation;
3643 /* both infinity: return lhs */
3645 decNumberCopy (res, lhs); /* [nop if in place] */
3649 /* set requested exponent */
3651 reqexp = inrhs->exponent; /* quantize -- match exponents */
3653 { /* rescale -- use value of rhs */
3654 /* Original rhs must be an integer that fits and is in range */
3656 reqexp = decGetInt (inrhs, set);
3658 reqexp = decGetInt (inrhs);
3664 etiny = set->emin; /* no subnormals */
3667 if (reqexp == BADINT /* bad (rescale only) or .. */
3668 || (reqexp < etiny) /* < lowest */
3669 || (reqexp > set->emax))
3671 *status |= DEC_Invalid_operation;
3675 /* we've processed the RHS, so we can overwrite it now if necessary */
3677 { /* zero coefficient unchanged */
3678 decNumberCopy (res, lhs); /* [nop if in place] */
3679 res->exponent = reqexp; /* .. just set exponent */
3682 res->bits = 0; /* subset specification; no -0 */
3686 { /* non-zero lhs */
3687 Int adjust = reqexp - lhs->exponent; /* digit adjustment needed */
3688 /* if adjusted coefficient will not fit, give up now */
3689 if ((lhs->digits - adjust) > reqdigits)
3691 *status |= DEC_Invalid_operation;
3696 { /* increasing exponent */
3697 /* this will decrease the length of the coefficient by adjust */
3698 /* digits, and must round as it does so */
3699 decContext workset; /* work */
3700 workset = *set; /* clone rounding, etc. */
3701 workset.digits = lhs->digits - adjust; /* set requested length */
3702 /* [note that the latter can be <1, here] */
3703 decCopyFit (res, lhs, &workset, &residue, status); /* fit to result */
3704 decApplyRound (res, &workset, residue, status); /* .. and round */
3705 residue = 0; /* [used] */
3706 /* If we rounded a 999s case, exponent will be off by one; */
3707 /* adjust back if so. */
3708 if (res->exponent > reqexp)
3710 res->digits = decShiftToMost (res->lsu, res->digits, 1); /* shift */
3711 res->exponent--; /* (re)adjust the exponent. */
3714 if (ISZERO (res) && !set->extended)
3715 res->bits = 0; /* subset; no -0 */
3718 else /* adjust<=0 */
3719 { /* decreasing or = exponent */
3720 /* this will increase the length of the coefficient by -adjust */
3721 /* digits, by adding trailing zeros. */
3722 decNumberCopy (res, lhs); /* [it will fit] */
3723 /* if padding needed (adjust<0), add it now... */
3727 decShiftToMost (res->lsu, res->digits, -adjust);
3728 res->exponent += adjust; /* adjust the exponent */
3733 /* Check for overflow [do not use Finalize in this case, as an */
3734 /* overflow here is a "don't fit" situation] */
3735 if (res->exponent > set->emax - res->digits + 1)
3737 *status |= DEC_Invalid_operation;
3742 decFinalize (res, set, &residue, status); /* set subnormal flags */
3743 *status &= ~DEC_Underflow; /* suppress Underflow [754r] */
3746 while (0); /* end protected */
3748 if (allocrhs != NULL)
3749 free (allocrhs); /* drop any storage we used */
3750 if (alloclhs != NULL)
3751 free (alloclhs); /* .. */
3755 /* ------------------------------------------------------------------ */
3756 /* decCompareOp -- compare, min, or max two Numbers */
3758 /* This computes C = A ? B and returns the signum (as a Number) */
3759 /* for COMPARE or the maximum or minimum (for COMPMAX and COMPMIN). */
3761 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
3764 /* set is the context */
3765 /* op is the operation flag */
3766 /* status is the usual accumulator */
3768 /* C must have space for one digit for COMPARE or set->digits for */
3769 /* COMPMAX and COMPMIN. */
3770 /* ------------------------------------------------------------------ */
3771 /* The emphasis here is on speed for common cases, and avoiding */
3772 /* coefficient comparison if possible. */
3773 /* ------------------------------------------------------------------ */
3775 decCompareOp (decNumber * res, decNumber * lhs, decNumber * rhs,
3776 decContext * set, Flag op, uInt * status)
3778 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3779 decNumber *allocrhs = NULL; /* .., rhs */
3780 Int result = 0; /* default result value */
3781 uByte merged; /* merged flags */
3782 uByte bits = 0; /* non-0 for NaN */
3785 if (decCheckOperands (res, lhs, rhs, set))
3790 { /* protect allocated storage */
3794 /* reduce operands and set lostDigits status, as needed */
3795 if (lhs->digits > set->digits)
3797 alloclhs = decRoundOperand (lhs, set, status);
3798 if (alloclhs == NULL)
3805 if (rhs->digits > set->digits)
3807 allocrhs = decRoundOperand (rhs, set, status);
3808 if (allocrhs == NULL)
3817 /* [following code does not require input rounding] */
3819 /* handle NaNs now; let infinities drop through */
3820 /* +++ review sNaN handling with 754r, for now assumes sNaN */
3821 /* (even just one) leads to NaN. */
3822 merged = (lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
3824 { /* a NaN bit set */
3826 else if (merged & DECSNAN);
3828 { /* 754r rules for MIN and MAX ignore single NaN */
3829 /* here if MIN or MAX, and one or two quiet NaNs */
3830 if (lhs->bits & rhs->bits & DECNAN);
3832 { /* just one quiet NaN */
3833 /* force choice to be the non-NaN operand */
3835 if (lhs->bits & DECNAN)
3836 result = -1; /* pick rhs */
3838 result = +1; /* pick lhs */
3842 op = COMPNAN; /* use special path */
3843 decNaNs (res, lhs, rhs, status);
3847 result = decCompare (lhs, rhs); /* we have numbers */
3849 while (0); /* end protected */
3851 if (result == BADINT)
3852 *status |= DEC_Insufficient_storage; /* rare */
3856 { /* return signum */
3857 decNumberZero (res); /* [always a valid result] */
3859 res->bits = bits; /* (maybe qNaN) */
3867 else if (op == COMPNAN); /* special, drop through */
3869 { /* MAX or MIN, non-NaN result */
3870 Int residue = 0; /* rounding accumulator */
3871 /* choose the operand for the result */
3874 { /* operands are numerically equal */
3875 /* choose according to sign then exponent (see 754r) */
3876 uByte slhs = (lhs->bits & DECNEG);
3877 uByte srhs = (rhs->bits & DECNEG);
3880 { /* subset: force left-hand */
3887 { /* signs differ */
3889 result = -1; /* rhs is max */
3891 result = +1; /* lhs is max */
3893 else if (slhs && srhs)
3894 { /* both negative */
3895 if (lhs->exponent < rhs->exponent)
3899 /* [if equal, we use lhs, technically identical] */
3902 { /* both positive */
3903 if (lhs->exponent > rhs->exponent)
3909 } /* numerically equal */
3910 /* here result will be non-0 */
3912 result = -result; /* reverse if looking for MIN */
3913 choice = (result > 0 ? lhs : rhs); /* choose */
3914 /* copy chosen to result, rounding if need be */
3915 decCopyFit (res, choice, set, &residue, status);
3916 decFinish (res, set, &residue, status);
3919 if (allocrhs != NULL)
3920 free (allocrhs); /* free any storage we used */
3921 if (alloclhs != NULL)
3922 free (alloclhs); /* .. */
3926 /* ------------------------------------------------------------------ */
3927 /* decCompare -- compare two decNumbers by numerical value */
3929 /* This routine compares A ? B without altering them. */
3931 /* Arg1 is A, a decNumber which is not a NaN */
3932 /* Arg2 is B, a decNumber which is not a NaN */
3934 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
3935 /* (the only possible failure is an allocation error) */
3936 /* ------------------------------------------------------------------ */
3937 /* This could be merged into decCompareOp */
3939 decCompare (decNumber * lhs, decNumber * rhs)
3941 Int result; /* result value */
3942 Int sigr; /* rhs signum */
3943 Int compare; /* work */
3944 result = 1; /* assume signum(lhs) */
3947 else if (decNumberIsNegative (lhs))
3949 sigr = 1; /* compute signum(rhs) */
3952 else if (decNumberIsNegative (rhs))
3955 return +1; /* L > R, return 1 */
3957 return -1; /* R < L, return -1 */
3959 /* signums are the same */
3961 return 0; /* both 0 */
3963 if ((lhs->bits | rhs->bits) & DECINF)
3964 { /* one or more infinities */
3965 if (lhs->bits == rhs->bits)
3966 result = 0; /* both the same */
3967 else if (decNumberIsInfinite (rhs))
3972 /* we must compare the coefficients, allowing for exponents */
3973 if (lhs->exponent > rhs->exponent)
3974 { /* LHS exponent larger */
3975 /* swap sides, and sign */
3976 decNumber *temp = lhs;
3982 compare = decUnitCompare (lhs->lsu, D2U (lhs->digits),
3983 rhs->lsu, D2U (rhs->digits),
3984 rhs->exponent - lhs->exponent);
3986 if (compare != BADINT)
3987 compare *= result; /* comparison succeeded */
3988 return compare; /* what we got */
3991 /* ------------------------------------------------------------------ */
3992 /* decUnitCompare -- compare two >=0 integers in Unit arrays */
3994 /* This routine compares A ? B*10**E where A and B are unit arrays */
3995 /* A is a plain integer */
3996 /* B has an exponent of E (which must be non-negative) */
3998 /* Arg1 is A first Unit (lsu) */
3999 /* Arg2 is A length in Units */
4000 /* Arg3 is B first Unit (lsu) */
4001 /* Arg4 is B length in Units */
4004 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
4005 /* (the only possible failure is an allocation error) */
4006 /* ------------------------------------------------------------------ */
4008 decUnitCompare (Unit * a, Int alength, Unit * b, Int blength, Int exp)
4010 Unit *acc; /* accumulator for result */
4011 Unit accbuff[D2U (DECBUFFER + 1)]; /* local buffer */
4012 Unit *allocacc = NULL; /* -> allocated acc buffer, iff allocated */
4013 Int accunits, need; /* units in use or needed for acc */
4014 Unit *l, *r, *u; /* work */
4015 Int expunits, exprem, result; /* .. */
4018 { /* aligned; fastpath */
4019 if (alength > blength)
4021 if (alength < blength)
4023 /* same number of units in both -- need unit-by-unit compare */
4024 l = a + alength - 1;
4025 r = b + alength - 1;
4026 for (; l >= a; l--, r--)
4033 return 0; /* all units match */
4036 /* Unaligned. If one is >1 unit longer than the other, padded */
4037 /* approximately, then we can return easily */
4038 if (alength > blength + (Int) D2U (exp))
4040 if (alength + 1 < blength + (Int) D2U (exp))
4043 /* We need to do a real subtract. For this, we need a result buffer */
4044 /* even though we only are interested in the sign. Its length needs */
4045 /* to be the larger of alength and padded blength, +2 */
4046 need = blength + D2U (exp); /* maximum real length of B */
4050 acc = accbuff; /* assume use local buffer */
4051 if (need * sizeof (Unit) > sizeof (accbuff))
4053 allocacc = (Unit *) malloc (need * sizeof (Unit));
4054 if (allocacc == NULL)
4055 return BADINT; /* hopeless -- abandon */
4058 /* Calculate units and remainder from exponent. */
4059 expunits = exp / DECDPUN;
4060 exprem = exp % DECDPUN;
4061 /* subtract [A+B*(-m)] */
4062 accunits = decUnitAddSub (a, alength, b, blength, expunits, acc,
4063 -(Int) powers[exprem]);
4064 /* [UnitAddSub result may have leading zeros, even on zero] */
4066 result = -1; /* negative result */
4068 { /* non-negative result */
4069 /* check units of the result before freeing any storage */
4070 for (u = acc; u < acc + accunits - 1 && *u == 0;)
4072 result = (*u == 0 ? 0 : +1);
4074 /* clean up and return the result */
4075 if (allocacc != NULL)
4076 free (allocacc); /* drop any storage we used */
4080 /* ------------------------------------------------------------------ */
4081 /* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays */
4083 /* This routine performs the calculation: */
4087 /* Where M is in the range -DECDPUNMAX through +DECDPUNMAX. */
4089 /* A may be shorter or longer than B. */
4091 /* Leading zeros are not removed after a calculation. The result is */
4092 /* either the same length as the longer of A and B (adding any */
4093 /* shift), or one Unit longer than that (if a Unit carry occurred). */
4095 /* A and B content are not altered unless C is also A or B. */
4096 /* C may be the same array as A or B, but only if no zero padding is */
4097 /* requested (that is, C may be B only if bshift==0). */
4098 /* C is filled from the lsu; only those units necessary to complete */
4099 /* the calculation are referenced. */
4101 /* Arg1 is A first Unit (lsu) */
4102 /* Arg2 is A length in Units */
4103 /* Arg3 is B first Unit (lsu) */
4104 /* Arg4 is B length in Units */
4105 /* Arg5 is B shift in Units (>=0; pads with 0 units if positive) */
4106 /* Arg6 is C first Unit (lsu) */
4107 /* Arg7 is M, the multiplier */
4109 /* returns the count of Units written to C, which will be non-zero */
4110 /* and negated if the result is negative. That is, the sign of the */
4111 /* returned Int is the sign of the result (positive for zero) and */
4112 /* the absolute value of the Int is the count of Units. */
4114 /* It is the caller's responsibility to make sure that C size is */
4115 /* safe, allowing space if necessary for a one-Unit carry. */
4117 /* This routine is severely performance-critical; *any* change here */
4118 /* must be measured (timed) to assure no performance degradation. */
4119 /* In particular, trickery here tends to be counter-productive, as */
4120 /* increased complexity of code hurts register optimizations on */
4121 /* register-poor architectures. Avoiding divisions is nearly */
4122 /* always a Good Idea, however. */
4124 /* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark */
4125 /* (IBM Warwick, UK) for some of the ideas used in this routine. */
4126 /* ------------------------------------------------------------------ */
4128 decUnitAddSub (Unit * a, Int alength,
4129 Unit * b, Int blength, Int bshift, Unit * c, Int m)
4131 Unit *alsu = a; /* A lsu [need to remember it] */
4132 Unit *clsu = c; /* C ditto */
4133 Unit *minC; /* low water mark for C */
4134 Unit *maxC; /* high water mark for C */
4135 eInt carry = 0; /* carry integer (could be Long) */
4137 #if DECDPUN==4 /* myriadal */
4138 Int est; /* estimated quotient */
4142 if (alength < 1 || blength < 1)
4143 printf ("decUnitAddSub: alen blen m %d %d [%d]\n", alength, blength, m);
4146 maxC = c + alength; /* A is usually the longer */
4147 minC = c + blength; /* .. and B the shorter */
4149 { /* B is shifted; low As copy across */
4151 /* if in place [common], skip copy unless there's a gap [rare] */
4152 if (a == c && bshift <= alength)
4158 for (; c < clsu + bshift; a++, c++)
4160 if (a < alsu + alength)
4173 /* For speed, we do the addition as two loops; the first where both A */
4174 /* and B contribute, and the second (if necessary) where only one or */
4175 /* other of the numbers contribute. */
4176 /* Carry handling is the same (i.e., duplicated) in each case. */
4177 for (; c < minC; c++)
4181 carry += ((eInt) * b) * m; /* [special-casing m=1/-1 */
4182 b++; /* here is not a win] */
4183 /* here carry is new Unit of digits; it could be +ve or -ve */
4184 if ((ueInt) carry <= DECDPUNMAX)
4185 { /* fastpath 0-DECDPUNMAX */
4190 /* remainder operator is undefined if negative, so we must test */
4191 #if DECDPUN==4 /* use divide-by-multiply */
4194 est = (((ueInt) carry >> 11) * 53687) >> 18;
4195 *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
4196 carry = est; /* likely quotient [89%] */
4197 if (*c < DECDPUNMAX + 1)
4198 continue; /* estimate was correct */
4200 *c -= DECDPUNMAX + 1;
4204 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4205 est = (((ueInt) carry >> 11) * 53687) >> 18;
4206 *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4207 carry = est - (DECDPUNMAX + 1); /* correctly negative */
4208 if (*c < DECDPUNMAX + 1)
4209 continue; /* was OK */
4211 *c -= DECDPUNMAX + 1;
4213 if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4214 { /* fastpath carry +1 */
4215 *c = (Unit) (carry - (DECDPUNMAX + 1)); /* [helps additions] */
4221 *c = (Unit) (carry % (DECDPUNMAX + 1));
4222 carry = carry / (DECDPUNMAX + 1);
4226 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4227 *c = (Unit) (carry % (DECDPUNMAX + 1));
4228 carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4232 /* we now may have one or other to complete */
4233 /* [pretest to avoid loop setup/shutdown] */
4235 for (; c < maxC; c++)
4237 if (a < alsu + alength)
4244 carry += ((eInt) * b) * m;
4247 /* here carry is new Unit of digits; it could be +ve or -ve and */
4248 /* magnitude up to DECDPUNMAX squared */
4249 if ((ueInt) carry <= DECDPUNMAX)
4250 { /* fastpath 0-DECDPUNMAX */
4255 /* result for this unit is negative or >DECDPUNMAX */
4256 #if DECDPUN==4 /* use divide-by-multiply */
4257 /* remainder is undefined if negative, so we must test */
4260 est = (((ueInt) carry >> 11) * 53687) >> 18;
4261 *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
4262 carry = est; /* likely quotient [79.7%] */
4263 if (*c < DECDPUNMAX + 1)
4264 continue; /* estimate was correct */
4266 *c -= DECDPUNMAX + 1;
4270 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4271 est = (((ueInt) carry >> 11) * 53687) >> 18;
4272 *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4273 carry = est - (DECDPUNMAX + 1); /* correctly negative */
4274 if (*c < DECDPUNMAX + 1)
4275 continue; /* was OK */
4277 *c -= DECDPUNMAX + 1;
4279 if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4280 { /* fastpath carry 1 */
4281 *c = (Unit) (carry - (DECDPUNMAX + 1));
4285 /* remainder is undefined if negative, so we must test */
4288 *c = (Unit) (carry % (DECDPUNMAX + 1));
4289 carry = carry / (DECDPUNMAX + 1);
4293 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4294 *c = (Unit) (carry % (DECDPUNMAX + 1));
4295 carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4299 /* OK, all A and B processed; might still have carry or borrow */
4300 /* return number of Units in the result, negated if a borrow */
4302 return c - clsu; /* no carry, we're done */
4304 { /* positive carry */
4305 *c = (Unit) carry; /* place as new unit */
4309 /* -ve carry: it's a borrow; complement needed */
4310 add = 1; /* temporary carry... */
4311 for (c = clsu; c < maxC; c++)
4313 add = DECDPUNMAX + add - *c;
4314 if (add <= DECDPUNMAX)
4325 /* add an extra unit iff it would be non-zero */
4327 printf ("UAS borrow: add %d, carry %d\n", add, carry);
4329 if ((add - carry - 1) != 0)
4331 *c = (Unit) (add - carry - 1);
4332 c++; /* interesting, include it */
4334 return clsu - c; /* -ve result indicates borrowed */
4337 /* ------------------------------------------------------------------ */
4338 /* decTrim -- trim trailing zeros or normalize */
4340 /* dn is the number to trim or normalize */
4341 /* all is 1 to remove all trailing zeros, 0 for just fraction ones */
4342 /* dropped returns the number of discarded trailing zeros */
4345 /* All fields are updated as required. This is a utility operation, */
4346 /* so special values are unchanged and no error is possible. */
4347 /* ------------------------------------------------------------------ */
4349 decTrim (decNumber * dn, Flag all, Int * dropped)
4351 Int d, exp; /* work */
4353 Unit *up; /* -> current Unit */
4356 if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
4360 *dropped = 0; /* assume no zeros dropped */
4361 if ((dn->bits & DECSPECIAL) /* fast exit if special .. */
4362 || (*dn->lsu & 0x01))
4363 return dn; /* .. or odd */
4366 dn->exponent = 0; /* (sign is preserved) */
4370 /* we have a finite number which is even */
4372 cut = 1; /* digit (1-DECDPUN) in Unit */
4373 up = dn->lsu; /* -> current Unit */
4374 for (d = 0; d < dn->digits - 1; d++)
4375 { /* [don't strip the final digit] */
4376 /* slice by powers */
4378 uInt quot = QUOT10 (*up, cut);
4379 if ((*up - quot * powers[cut]) != 0)
4380 break; /* found non-0 digit */
4382 if (*up % powers[cut] != 0)
4383 break; /* found non-0 digit */
4385 /* have a trailing 0 */
4388 /* [if exp>0 then all trailing 0s are significant for trim] */
4390 { /* if digit might be significant */
4392 break; /* then quit */
4393 exp++; /* next digit might be significant */
4396 cut++; /* next power */
4398 { /* need new Unit */
4404 return dn; /* none dropped */
4406 /* effect the drop */
4407 decShiftToLeast (dn->lsu, D2U (dn->digits), d);
4408 dn->exponent += d; /* maintain numerical value */
4409 dn->digits -= d; /* new length */
4410 *dropped = d; /* report the count */
4414 /* ------------------------------------------------------------------ */
4415 /* decShiftToMost -- shift digits in array towards most significant */
4417 /* uar is the array */
4418 /* digits is the count of digits in use in the array */
4419 /* shift is the number of zeros to pad with (least significant); */
4420 /* it must be zero or positive */
4422 /* returns the new length of the integer in the array, in digits */
4424 /* No overflow is permitted (that is, the uar array must be known to */
4425 /* be large enough to hold the result, after shifting). */
4426 /* ------------------------------------------------------------------ */
4428 decShiftToMost (Unit * uar, Int digits, Int shift)
4430 Unit *target, *source, *first; /* work */
4431 uInt rem; /* for division */
4432 Int cut; /* odd 0's to add */
4433 uInt next; /* work */
4436 return digits; /* [fastpath] nothing to do */
4437 if ((digits + shift) <= DECDPUN)
4438 { /* [fastpath] single-unit case */
4439 *uar = (Unit) (*uar * powers[shift]);
4440 return digits + shift;
4443 cut = (DECDPUN - shift % DECDPUN) % DECDPUN;
4444 source = uar + D2U (digits) - 1; /* where msu comes from */
4445 first = uar + D2U (digits + shift) - 1; /* where msu of source will end up */
4446 target = source + D2U (shift); /* where upper part of first cut goes */
4449 for (; source >= uar; source--, target--)
4451 /* split the source Unit and accumulate remainder for next */
4453 uInt quot = QUOT10 (*source, cut);
4454 rem = *source - quot * powers[cut];
4457 rem = *source % powers[cut];
4458 next += *source / powers[cut];
4460 if (target <= first)
4461 *target = (Unit) next; /* write to target iff valid */
4462 next = rem * powers[DECDPUN - cut]; /* save remainder for next Unit */
4464 /* propagate to one below and clear the rest */
4465 for (; target >= uar; target--)
4467 *target = (Unit) next;
4470 return digits + shift;
4473 /* ------------------------------------------------------------------ */
4474 /* decShiftToLeast -- shift digits in array towards least significant */
4476 /* uar is the array */
4477 /* units is length of the array, in units */
4478 /* shift is the number of digits to remove from the lsu end; it */
4479 /* must be zero or positive and less than units*DECDPUN. */
4481 /* returns the new length of the integer in the array, in units */
4483 /* Removed digits are discarded (lost). Units not required to hold */
4484 /* the final result are unchanged. */
4485 /* ------------------------------------------------------------------ */
4487 decShiftToLeast (Unit * uar, Int units, Int shift)
4489 Unit *target, *up; /* work */
4490 Int cut, count; /* work */
4491 Int quot, rem; /* for division */
4494 return units; /* [fastpath] nothing to do */
4496 up = uar + shift / DECDPUN; /* source; allow for whole Units */
4497 cut = shift % DECDPUN; /* odd 0's to drop */
4498 target = uar; /* both paths */
4500 { /* whole units shift */
4501 for (; up < uar + units; target++, up++)
4503 return target - uar;
4506 count = units * DECDPUN - shift; /* the maximum new length */
4508 quot = QUOT10 (*up, cut);
4510 quot = *up / powers[cut];
4514 *target = (Unit) quot;
4515 count -= (DECDPUN - cut);
4521 quot = QUOT10 (quot, cut);
4522 rem = *up - quot * powers[cut];
4524 rem = quot % powers[cut];
4525 quot = quot / powers[cut];
4527 *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4532 return target - uar + 1;
4536 /* ------------------------------------------------------------------ */
4537 /* decRoundOperand -- round an operand [used for subset only] */
4539 /* dn is the number to round (dn->digits is > set->digits) */
4540 /* set is the relevant context */
4541 /* status is the status accumulator */
4543 /* returns an allocated decNumber with the rounded result. */
4545 /* lostDigits and other status may be set by this. */
4547 /* Since the input is an operand, we are not permitted to modify it. */
4548 /* We therefore return an allocated decNumber, rounded as required. */
4549 /* It is the caller's responsibility to free the allocated storage. */
4551 /* If no storage is available then the result cannot be used, so NULL */
4553 /* ------------------------------------------------------------------ */
4555 decRoundOperand (decNumber * dn, decContext * set, uInt * status)
4557 decNumber *res; /* result structure */
4558 uInt newstatus = 0; /* status from round */
4559 Int residue = 0; /* rounding accumulator */
4561 /* Allocate storage for the returned decNumber, big enough for the */
4562 /* length specified by the context */
4563 res = (decNumber *) malloc (sizeof (decNumber)
4564 + (D2U (set->digits) - 1) * sizeof (Unit));
4567 *status |= DEC_Insufficient_storage;
4570 decCopyFit (res, dn, set, &residue, &newstatus);
4571 decApplyRound (res, set, residue, &newstatus);
4573 /* If that set Inexact then we "lost digits" */
4574 if (newstatus & DEC_Inexact)
4575 newstatus |= DEC_Lost_digits;
4576 *status |= newstatus;
4581 /* ------------------------------------------------------------------ */
4582 /* decCopyFit -- copy a number, shortening the coefficient if needed */
4584 /* dest is the target decNumber */
4585 /* src is the source decNumber */
4586 /* set is the context [used for length (digits) and rounding mode] */
4587 /* residue is the residue accumulator */
4588 /* status contains the current status to be updated */
4590 /* (dest==src is allowed and will be a no-op if fits) */
4591 /* All fields are updated as required. */
4592 /* ------------------------------------------------------------------ */
4594 decCopyFit (decNumber * dest, decNumber * src, decContext * set,
4595 Int * residue, uInt * status)
4597 dest->bits = src->bits;
4598 dest->exponent = src->exponent;
4599 decSetCoeff (dest, set, src->lsu, src->digits, residue, status);
4602 /* ------------------------------------------------------------------ */
4603 /* decSetCoeff -- set the coefficient of a number */
4605 /* dn is the number whose coefficient array is to be set. */
4606 /* It must have space for set->digits digits */
4607 /* set is the context [for size] */
4608 /* lsu -> lsu of the source coefficient [may be dn->lsu] */
4609 /* len is digits in the source coefficient [may be dn->digits] */
4610 /* residue is the residue accumulator. This has values as in */
4611 /* decApplyRound, and will be unchanged unless the */
4612 /* target size is less than len. In this case, the */
4613 /* coefficient is truncated and the residue is updated to */
4614 /* reflect the previous residue and the dropped digits. */
4615 /* status is the status accumulator, as usual */
4617 /* The coefficient may already be in the number, or it can be an */
4618 /* external intermediate array. If it is in the number, lsu must == */
4619 /* dn->lsu and len must == dn->digits. */
4621 /* Note that the coefficient length (len) may be < set->digits, and */
4622 /* in this case this merely copies the coefficient (or is a no-op */
4623 /* if dn->lsu==lsu). */
4625 /* Note also that (only internally, from decNumberRescale and */
4626 /* decSetSubnormal) the value of set->digits may be less than one, */
4627 /* indicating a round to left. */
4628 /* This routine handles that case correctly; caller ensures space. */
4630 /* dn->digits, dn->lsu (and as required), and dn->exponent are */
4631 /* updated as necessary. dn->bits (sign) is unchanged. */
4633 /* DEC_Rounded status is set if any digits are discarded. */
4634 /* DEC_Inexact status is set if any non-zero digits are discarded, or */
4635 /* incoming residue was non-0 (implies rounded) */
4636 /* ------------------------------------------------------------------ */
4637 /* mapping array: maps 0-9 to canonical residues, so that we can */
4638 /* adjust by a residue in range [-1, +1] and achieve correct rounding */
4639 /* 0 1 2 3 4 5 6 7 8 9 */
4640 static const uByte resmap[10] = { 0, 3, 3, 3, 3, 5, 7, 7, 7, 7 };
4642 decSetCoeff (decNumber * dn, decContext * set, Unit * lsu,
4643 Int len, Int * residue, uInt * status)
4645 Int discard; /* number of digits to discard */
4646 uInt discard1; /* first discarded digit */
4647 uInt cut; /* cut point in Unit */
4648 uInt quot, rem; /* for divisions */
4649 Unit *up, *target; /* work */
4655 discard = len - set->digits; /* digits to discard */
4657 { /* no digits are being discarded */
4660 /* copy the coefficient array to the result number; no shift needed */
4662 for (target = dn->lsu; target < dn->lsu + D2U (len); target++, up++)
4666 dn->digits = len; /* set the new length */
4668 /* dn->exponent and residue are unchanged */
4670 *status |= (DEC_Inexact | DEC_Rounded); /* record inexactitude */
4674 /* we have to discard some digits */
4675 *status |= DEC_Rounded; /* accumulate Rounded status */
4677 *residue = 1; /* previous residue now to right, so -1 to +1 */
4680 { /* everything, +1, is being discarded */
4681 /* guard digit is 0 */
4682 /* residue is all the number [NB could be all 0s] */
4684 for (up = lsu + D2U (len) - 1; up >= lsu; up--)
4687 { /* found a non-0 */
4689 break; /* no need to check any others */
4693 *status |= DEC_Inexact; /* record inexactitude */
4694 *dn->lsu = 0; /* coefficient will now be 0 */
4695 dn->digits = 1; /* .. */
4696 dn->exponent += discard; /* maintain numerical value */
4698 } /* total discard */
4700 /* partial discard [most common case] */
4701 /* here, at least the first (most significant) discarded digit exists */
4703 /* spin up the number, noting residue as we pass, until we get to */
4704 /* the Unit with the first discarded digit. When we get there, */
4705 /* extract it and remember where we're at */
4707 for (up = lsu;; up++)
4710 if (count >= discard)
4711 break; /* full ones all checked */
4716 /* here up -> Unit with discarded digit */
4717 cut = discard - (count - DECDPUN) - 1;
4718 if (cut == DECDPUN - 1)
4719 { /* discard digit is at top */
4721 discard1 = QUOT10 (*up, DECDPUN - 1);
4722 rem = *up - discard1 * powers[DECDPUN - 1];
4724 rem = *up % powers[DECDPUN - 1];
4725 discard1 = *up / powers[DECDPUN - 1];
4729 up++; /* move to next */
4730 cut = 0; /* bottom digit of result */
4731 quot = 0; /* keep a certain compiler happy */
4735 /* discard digit is in low digit(s), not top digit */
4739 { /* it's not at bottom of Unit */
4741 quot = QUOT10 (*up, cut);
4742 rem = *up - quot * powers[cut];
4744 rem = *up % powers[cut];
4745 quot = *up / powers[cut];
4750 /* discard digit is now at bottom of quot */
4752 temp = (quot * 6554) >> 16; /* fast /10 */
4753 /* Vowels algorithm here not a win (9 instructions) */
4754 discard1 = quot - X10 (temp);
4757 discard1 = quot % 10;
4760 cut++; /* update cut */
4763 /* here: up -> Unit of the array with discarded digit */
4764 /* cut is the division point for each Unit */
4765 /* quot holds the uncut high-order digits for the current */
4766 /* Unit, unless cut==0 in which case it's still in *up */
4767 /* copy the coefficient array to the result number, shifting as we go */
4768 count = set->digits; /* digits to end up with */
4770 { /* special for Rescale/Subnormal :-( */
4771 *dn->lsu = 0; /* .. result is 0 */
4772 dn->digits = 1; /* .. */
4775 { /* shift to least */
4776 /* [this is similar to decShiftToLeast code, with copy] */
4777 dn->digits = count; /* set the new length */
4780 /* on unit boundary, so simple shift down copy loop suffices */
4781 for (target = dn->lsu; target < dn->lsu + D2U (count);
4788 for (target = dn->lsu;; target++)
4790 *target = (Unit) quot;
4791 count -= (DECDPUN - cut);
4797 quot = QUOT10 (quot, cut);
4798 rem = *up - quot * powers[cut];
4800 rem = quot % powers[cut];
4801 quot = quot / powers[cut];
4803 *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4808 } /* shift to least needed */
4809 dn->exponent += discard; /* maintain numerical value */
4811 /* here, discard1 is the guard digit, and residue is everything else */
4812 /* [use mapping to accumulate residue safely] */
4813 *residue += resmap[discard1];
4816 *status |= DEC_Inexact; /* record inexactitude */
4820 /* ------------------------------------------------------------------ */
4821 /* decApplyRound -- apply pending rounding to a number */
4823 /* dn is the number, with space for set->digits digits */
4824 /* set is the context [for size and rounding mode] */
4825 /* residue indicates pending rounding, being any accumulated */
4826 /* guard and sticky information. It may be: */
4827 /* 6-9: rounding digit is >5 */
4828 /* 5: rounding digit is exactly half-way */
4829 /* 1-4: rounding digit is <5 and >0 */
4830 /* 0: the coefficient is exact */
4831 /* -1: as 1, but the hidden digits are subtractive, that */
4832 /* is, of the opposite sign to dn. In this case the */
4833 /* coefficient must be non-0. */
4834 /* status is the status accumulator, as usual */
4836 /* This routine applies rounding while keeping the length of the */
4837 /* coefficient constant. The exponent and status are unchanged */
4840 /* -- the coefficient was increased and is all nines (in which */
4841 /* case Overflow could occur, and is handled directly here so */
4842 /* the caller does not need to re-test for overflow) */
4844 /* -- the coefficient was decreased and becomes all nines (in which */
4845 /* case Underflow could occur, and is also handled directly). */
4847 /* All fields in dn are updated as required. */
4849 /* ------------------------------------------------------------------ */
4851 decApplyRound (decNumber * dn, decContext * set, Int residue, uInt * status)
4853 Int bump; /* 1 if coefficient needs to be incremented */
4854 /* -1 if coefficient needs to be decremented */
4857 return; /* nothing to apply */
4859 bump = 0; /* assume a smooth ride */
4861 /* now decide whether, and how, to round, depending on mode */
4864 case DEC_ROUND_DOWN:
4866 /* no change, except if negative residue */
4872 case DEC_ROUND_HALF_DOWN:
4879 case DEC_ROUND_HALF_EVEN:
4882 bump = 1; /* >0.5 goes up */
4883 else if (residue == 5)
4884 { /* exactly 0.5000... */
4885 /* 0.5 goes up iff [new] lsd is odd */
4886 if (*dn->lsu & 0x01)
4892 case DEC_ROUND_HALF_UP:
4906 case DEC_ROUND_CEILING:
4908 /* same as _UP for positive numbers, and as _DOWN for negatives */
4909 /* [negative residue cannot occur on 0] */
4910 if (decNumberIsNegative (dn))
4923 case DEC_ROUND_FLOOR:
4925 /* same as _UP for negative numbers, and as _DOWN for positive */
4926 /* [negative residue cannot occur on 0] */
4927 if (!decNumberIsNegative (dn))
4941 { /* e.g., DEC_ROUND_MAX */
4942 *status |= DEC_Invalid_context;
4944 printf ("Unknown rounding mode: %d\n", set->round);
4950 /* now bump the number, up or down, if need be */
4952 return; /* no action required */
4954 /* Simply use decUnitAddSub unless we are bumping up and the number */
4955 /* is all nines. In this special case we set to 1000... and adjust */
4956 /* the exponent by one (as otherwise we could overflow the array) */
4957 /* Similarly handle all-nines result if bumping down. */
4960 Unit *up; /* work */
4961 uInt count = dn->digits; /* digits to be checked */
4962 for (up = dn->lsu;; up++)
4964 if (count <= DECDPUN)
4966 /* this is the last Unit (the msu) */
4967 if (*up != powers[count] - 1)
4968 break; /* not still 9s */
4969 /* here if it, too, is all nines */
4970 *up = (Unit) powers[count - 1]; /* here 999 -> 100 etc. */
4971 for (up = up - 1; up >= dn->lsu; up--)
4972 *up = 0; /* others all to 0 */
4973 dn->exponent++; /* and bump exponent */
4974 /* [which, very rarely, could cause Overflow...] */
4975 if ((dn->exponent + dn->digits) > set->emax + 1)
4977 decSetOverflow (dn, set, status);
4981 /* a full unit to check, with more to come */
4982 if (*up != DECDPUNMAX)
4983 break; /* not still 9s */
4989 /* here we are lookng for a pre-bump of 1000... (leading 1, */
4990 /* all other digits zero) */
4991 Unit *up, *sup; /* work */
4992 uInt count = dn->digits; /* digits to be checked */
4993 for (up = dn->lsu;; up++)
4995 if (count <= DECDPUN)
4997 /* this is the last Unit (the msu) */
4998 if (*up != powers[count - 1])
4999 break; /* not 100.. */
5000 /* here if we have the 1000... case */
5001 sup = up; /* save msu pointer */
5002 *up = (Unit) powers[count] - 1; /* here 100 in msu -> 999 */
5003 /* others all to all-nines, too */
5004 for (up = up - 1; up >= dn->lsu; up--)
5005 *up = (Unit) powers[DECDPUN] - 1;
5006 dn->exponent--; /* and bump exponent */
5008 /* iff the number was at the subnormal boundary (exponent=etiny) */
5009 /* then the exponent is now out of range, so it will in fact get */
5010 /* clamped to etiny and the final 9 dropped. */
5011 /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
5012 /* dn->exponent, set->digits); */
5013 if (dn->exponent + 1 == set->emin - set->digits + 1)
5015 if (count == 1 && dn->digits == 1)
5016 *sup = 0; /* here 9 -> 0[.9] */
5019 *sup = (Unit) powers[count - 1] - 1; /* here 999.. in msu -> 99.. */
5024 DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5029 /* a full unit to check, with more to come */
5031 break; /* not still 0s */
5037 /* Actual bump needed. Do it. */
5038 decUnitAddSub (dn->lsu, D2U (dn->digits), one, 1, 0, dn->lsu, bump);
5042 /* ------------------------------------------------------------------ */
5043 /* decFinish -- finish processing a number */
5045 /* dn is the number */
5046 /* set is the context */
5047 /* residue is the rounding accumulator (as in decApplyRound) */
5048 /* status is the accumulator */
5050 /* This finishes off the current number by: */
5051 /* 1. If not extended: */
5052 /* a. Converting a zero result to clean '0' */
5053 /* b. Reducing positive exponents to 0, if would fit in digits */
5054 /* 2. Checking for overflow and subnormals (always) */
5055 /* Note this is just Finalize when no subset arithmetic. */
5056 /* All fields are updated as required. */
5057 /* ------------------------------------------------------------------ */
5059 decFinish (decNumber * dn, decContext * set, Int * residue, uInt * status)
5065 { /* value is zero */
5066 dn->exponent = 0; /* clean exponent .. */
5067 dn->bits = 0; /* .. and sign */
5068 return; /* no error possible */
5070 if (dn->exponent >= 0)
5071 { /* non-negative exponent */
5072 /* >0; reduce to integer if possible */
5073 if (set->digits >= (dn->exponent + dn->digits))
5075 dn->digits = decShiftToMost (dn->lsu, dn->digits, dn->exponent);
5081 decFinalize (dn, set, residue, status);
5085 /* ------------------------------------------------------------------ */
5086 /* decFinalize -- final check, clamp, and round of a number */
5088 /* dn is the number */
5089 /* set is the context */
5090 /* residue is the rounding accumulator (as in decApplyRound) */
5091 /* status is the status accumulator */
5093 /* This finishes off the current number by checking for subnormal */
5094 /* results, applying any pending rounding, checking for overflow, */
5095 /* and applying any clamping. */
5096 /* Underflow and overflow conditions are raised as appropriate. */
5097 /* All fields are updated as required. */
5098 /* ------------------------------------------------------------------ */
5100 decFinalize (decNumber * dn, decContext * set, Int * residue, uInt * status)
5102 Int shift; /* shift needed if clamping */
5104 /* We have to be careful when checking the exponent as the adjusted */
5105 /* exponent could overflow 31 bits [because it may already be up */
5106 /* to twice the expected]. */
5108 /* First test for subnormal. This must be done before any final */
5109 /* round as the result could be rounded to Nmin or 0. */
5110 if (dn->exponent < 0 /* negative exponent */
5111 && (dn->exponent < set->emin - dn->digits + 1))
5113 /* Go handle subnormals; this will apply round if needed. */
5114 decSetSubnormal (dn, set, residue, status);
5118 /* now apply any pending round (this could raise overflow). */
5120 decApplyRound (dn, set, *residue, status);
5122 /* Check for overflow [redundant in the 'rare' case] or clamp */
5123 if (dn->exponent <= set->emax - set->digits + 1)
5124 return; /* neither needed */
5126 /* here when we might have an overflow or clamp to do */
5127 if (dn->exponent > set->emax - dn->digits + 1)
5129 decSetOverflow (dn, set, status);
5132 /* here when the result is normal but in clamp range */
5136 /* here when we need to apply the IEEE exponent clamp (fold-down) */
5137 shift = dn->exponent - (set->emax - set->digits + 1);
5139 /* shift coefficient (if non-zero) */
5142 dn->digits = decShiftToMost (dn->lsu, dn->digits, shift);
5144 dn->exponent -= shift; /* adjust the exponent to match */
5145 *status |= DEC_Clamped; /* and record the dirty deed */
5149 /* ------------------------------------------------------------------ */
5150 /* decSetOverflow -- set number to proper overflow value */
5152 /* dn is the number (used for sign [only] and result) */
5153 /* set is the context [used for the rounding mode] */
5154 /* status contains the current status to be updated */
5156 /* This sets the sign of a number and sets its value to either */
5157 /* Infinity or the maximum finite value, depending on the sign of */
5158 /* dn and therounding mode, following IEEE 854 rules. */
5159 /* ------------------------------------------------------------------ */
5161 decSetOverflow (decNumber * dn, decContext * set, uInt * status)
5163 Flag needmax = 0; /* result is maximum finite value */
5164 uByte sign = dn->bits & DECNEG; /* clean and save sign bit */
5167 { /* zero does not overflow magnitude */
5168 Int emax = set->emax; /* limit value */
5170 emax -= set->digits - 1; /* lower if clamping */
5171 if (dn->exponent > emax)
5172 { /* clamp required */
5173 dn->exponent = emax;
5174 *status |= DEC_Clamped;
5182 case DEC_ROUND_DOWN:
5184 needmax = 1; /* never Infinity */
5187 case DEC_ROUND_CEILING:
5190 needmax = 1; /* Infinity if non-negative */
5193 case DEC_ROUND_FLOOR:
5196 needmax = 1; /* Infinity if negative */
5200 break; /* Infinity in all other cases */
5204 Unit *up; /* work */
5205 Int count = set->digits; /* nines to add */
5207 /* fill in all nines to set maximum value */
5208 for (up = dn->lsu;; up++)
5210 if (count > DECDPUN)
5211 *up = DECDPUNMAX; /* unit full o'nines */
5213 { /* this is the msu */
5214 *up = (Unit) (powers[count] - 1);
5217 count -= DECDPUN; /* we filled those digits */
5219 dn->bits = sign; /* sign */
5220 dn->exponent = set->emax - set->digits + 1;
5223 dn->bits = sign | DECINF; /* Value is +/-Infinity */
5224 *status |= DEC_Overflow | DEC_Inexact | DEC_Rounded;
5227 /* ------------------------------------------------------------------ */
5228 /* decSetSubnormal -- process value whose exponent is <Emin */
5230 /* dn is the number (used as input as well as output; it may have */
5231 /* an allowed subnormal value, which may need to be rounded) */
5232 /* set is the context [used for the rounding mode] */
5233 /* residue is any pending residue */
5234 /* status contains the current status to be updated */
5236 /* If subset mode, set result to zero and set Underflow flags. */
5238 /* Value may be zero with a low exponent; this does not set Subnormal */
5239 /* but the exponent will be clamped to Etiny. */
5241 /* Otherwise ensure exponent is not out of range, and round as */
5242 /* necessary. Underflow is set if the result is Inexact. */
5243 /* ------------------------------------------------------------------ */
5245 decSetSubnormal (decNumber * dn, decContext * set, Int * residue,
5248 decContext workset; /* work */
5249 Int etiny, adjust; /* .. */
5252 /* simple set to zero and 'hard underflow' for subset */
5256 /* always full overflow */
5257 *status |= DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5262 /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
5263 /* (Etiny) if needed */
5264 etiny = set->emin - (set->digits - 1); /* smallest allowed exponent */
5268 { /* value is zero */
5269 /* residue can never be non-zero here */
5273 printf ("++ Subnormal 0 residue %d\n", *residue);
5274 *status |= DEC_Invalid_operation;
5277 if (dn->exponent < etiny)
5278 { /* clamp required */
5279 dn->exponent = etiny;
5280 *status |= DEC_Clamped;
5285 *status |= DEC_Subnormal; /* we have a non-zero subnormal */
5287 adjust = etiny - dn->exponent; /* calculate digits to remove */
5289 { /* not out of range; unrounded */
5290 /* residue can never be non-zero here, so fast-path out */
5294 printf ("++ Subnormal no-adjust residue %d\n", *residue);
5295 *status |= DEC_Invalid_operation;
5298 /* it may already be inexact (from setting the coefficient) */
5299 if (*status & DEC_Inexact)
5300 *status |= DEC_Underflow;
5304 /* adjust>0. we need to rescale the result so exponent becomes Etiny */
5305 /* [this code is similar to that in rescale] */
5306 workset = *set; /* clone rounding, etc. */
5307 workset.digits = dn->digits - adjust; /* set requested length */
5308 workset.emin -= adjust; /* and adjust emin to match */
5309 /* [note that the latter can be <1, here, similar to Rescale case] */
5310 decSetCoeff (dn, &workset, dn->lsu, dn->digits, residue, status);
5311 decApplyRound (dn, &workset, *residue, status);
5313 /* Use 754R/854 default rule: Underflow is set iff Inexact */
5314 /* [independent of whether trapped] */
5315 if (*status & DEC_Inexact)
5316 *status |= DEC_Underflow;
5318 /* if we rounded up a 999s case, exponent will be off by one; adjust */
5319 /* back if so [it will fit, because we shortened] */
5320 if (dn->exponent > etiny)
5322 dn->digits = decShiftToMost (dn->lsu, dn->digits, 1);
5323 dn->exponent--; /* (re)adjust the exponent. */
5327 /* ------------------------------------------------------------------ */
5328 /* decGetInt -- get integer from a number */
5330 /* dn is the number [which will not be altered] */
5331 /* set is the context [requested digits], subset only */
5332 /* returns the converted integer, or BADINT if error */
5334 /* This checks and gets a whole number from the input decNumber. */
5335 /* The magnitude of the integer must be <2^31. */
5336 /* Any discarded fractional part must be 0. */
5337 /* If subset it must also fit in set->digits */
5338 /* ------------------------------------------------------------------ */
5341 decGetInt (decNumber * dn, decContext * set)
5345 decGetInt (decNumber * dn)
5348 Int theInt; /* result accumulator */
5349 Unit *up; /* work */
5350 Int got; /* digits (real or not) processed */
5351 Int ilength = dn->digits + dn->exponent; /* integral length */
5353 /* The number must be an integer that fits in 10 digits */
5354 /* Assert, here, that 10 is enough for any rescale Etiny */
5355 #if DEC_MAX_EMAX > 999999999
5356 #error GetInt may need updating [for Emax]
5358 #if DEC_MIN_EMIN < -999999999
5359 #error GetInt may need updating [for Emin]
5362 return 0; /* zeros are OK, with any exponent */
5364 return BADINT; /* always too big */
5366 if (!set->extended && ilength > set->digits)
5370 up = dn->lsu; /* ready for lsu */
5371 theInt = 0; /* ready to accumulate */
5372 if (dn->exponent >= 0)
5373 { /* relatively easy */
5374 /* no fractional part [usual]; allow for positive exponent */
5378 { /* -ve exponent; some fractional part to check and discard */
5379 Int count = -dn->exponent; /* digits to discard */
5380 /* spin up whole units until we get to the Unit with the unit digit */
5381 for (; count >= DECDPUN; up++)
5384 return BADINT; /* non-zero Unit to discard */
5388 got = 0; /* [a multiple of DECDPUN] */
5390 { /* [not multiple of DECDPUN] */
5392 /* slice off fraction digits and check for non-zero */
5394 theInt = QUOT10 (*up, count);
5395 rem = *up - theInt * powers[count];
5397 rem = *up % powers[count]; /* slice off discards */
5398 theInt = *up / powers[count];
5401 return BADINT; /* non-zero fraction */
5402 /* OK, we're good */
5403 got = DECDPUN - count; /* number of digits so far */
5404 up++; /* ready for next */
5407 /* collect the rest */
5408 for (; got < ilength; up++)
5410 theInt += *up * powers[got];
5413 if ((ilength == 10) /* check no wrap */
5414 && (theInt / (Int) powers[got - DECDPUN] != *(up - 1)))
5416 /* [that test also disallows the BADINT result case] */
5418 /* apply any sign and return */
5419 if (decNumberIsNegative (dn))
5424 /* ------------------------------------------------------------------ */
5425 /* decStrEq -- caseless comparison of strings */
5427 /* str1 is one of the strings to compare */
5428 /* str2 is the other */
5430 /* returns 1 if strings caseless-compare equal, 0 otherwise */
5432 /* Note that the strings must be the same length if they are to */
5433 /* compare equal; there is no padding. */
5434 /* ------------------------------------------------------------------ */
5435 /* [strcmpi is not in ANSI C] */
5437 decStrEq (const char *str1, const char *str2)
5439 for (;; str1++, str2++)
5448 if (tolower (*str1) != tolower (*str2))
5455 /* ------------------------------------------------------------------ */
5456 /* decNaNs -- handle NaN operand or operands */
5458 /* res is the result number */
5459 /* lhs is the first operand */
5460 /* rhs is the second operand, or NULL if none */
5461 /* status contains the current status */
5462 /* returns res in case convenient */
5464 /* Called when one or both operands is a NaN, and propagates the */
5465 /* appropriate result to res. When an sNaN is found, it is changed */
5466 /* to a qNaN and Invalid operation is set. */
5467 /* ------------------------------------------------------------------ */
5469 decNaNs (decNumber * res, decNumber * lhs, decNumber * rhs, uInt * status)
5471 /* This decision tree ends up with LHS being the source pointer, */
5472 /* and status updated if need be */
5473 if (lhs->bits & DECSNAN)
5474 *status |= DEC_Invalid_operation | DEC_sNaN;
5475 else if (rhs == NULL);
5476 else if (rhs->bits & DECSNAN)
5479 *status |= DEC_Invalid_operation | DEC_sNaN;
5481 else if (lhs->bits & DECNAN);
5485 decNumberCopy (res, lhs);
5486 res->bits &= ~DECSNAN; /* convert any sNaN to NaN, while */
5487 res->bits |= DECNAN; /* .. preserving sign */
5488 res->exponent = 0; /* clean exponent */
5489 /* [coefficient was copied] */
5493 /* ------------------------------------------------------------------ */
5494 /* decStatus -- apply non-zero status */
5496 /* dn is the number to set if error */
5497 /* status contains the current status (not yet in context) */
5498 /* set is the context */
5500 /* If the status is an error status, the number is set to a NaN, */
5501 /* unless the error was an overflow, divide-by-zero, or underflow, */
5502 /* in which case the number will have already been set. */
5504 /* The context status is then updated with the new status. Note that */
5505 /* this may raise a signal, so control may never return from this */
5506 /* routine (hence resources must be recovered before it is called). */
5507 /* ------------------------------------------------------------------ */
5509 decStatus (decNumber * dn, uInt status, decContext * set)
5511 if (status & DEC_NaNs)
5512 { /* error status -> NaN */
5513 /* if cause was an sNaN, clear and propagate [NaN is already set up] */
5514 if (status & DEC_sNaN)
5515 status &= ~DEC_sNaN;
5518 decNumberZero (dn); /* other error: clean throughout */
5519 dn->bits = DECNAN; /* and make a quiet NaN */
5522 decContextSetStatus (set, status);
5526 /* ------------------------------------------------------------------ */
5527 /* decGetDigits -- count digits in a Units array */
5529 /* uar is the Unit array holding the number [this is often an */
5530 /* accumulator of some sort] */
5531 /* len is the length of the array in units */
5533 /* returns the number of (significant) digits in the array */
5535 /* All leading zeros are excluded, except the last if the array has */
5536 /* only zero Units. */
5537 /* ------------------------------------------------------------------ */
5538 /* This may be called twice during some operations. */
5540 decGetDigits (Unit * uar, Int len)
5542 Unit *up = uar + len - 1; /* -> msu */
5543 Int digits = len * DECDPUN; /* maximum possible digits */
5544 uInt const *pow; /* work */
5546 for (; up >= uar; up--)
5552 continue; /* more to check */
5553 /* all units were 0 */
5554 digits++; /* .. so bump digits to 1 */
5557 /* found the first non-zero Unit */
5560 break; /* fastpath 1-9 */
5562 for (pow = &powers[2]; *up >= *pow; pow++)
5571 #if DECTRACE | DECCHECK
5572 /* ------------------------------------------------------------------ */
5573 /* decNumberShow -- display a number [debug aid] */
5574 /* dn is the number to show */
5576 /* Shows: sign, exponent, coefficient (msu first), digits */
5577 /* or: sign, special-value */
5578 /* ------------------------------------------------------------------ */
5579 /* this is public so other modules can use it */
5581 decNumberShow (decNumber * dn)
5583 Unit *up; /* work */
5586 char isign = '+'; /* main sign */
5592 if (decNumberIsNegative (dn))
5594 printf (" >> %c ", isign);
5595 if (dn->bits & DECSPECIAL)
5596 { /* Is a special value */
5597 if (decNumberIsInfinite (dn))
5598 printf ("Infinity");
5601 if (dn->bits & DECSNAN)
5602 printf ("sNaN"); /* signalling NaN */
5606 /* if coefficient and exponent are 0, we're done */
5607 if (dn->exponent == 0 && dn->digits == 1 && *dn->lsu == 0)
5612 /* drop through to report other information */
5616 /* now carefully display the coefficient */
5617 up = dn->lsu + D2U (dn->digits) - 1; /* msu */
5619 for (up = up - 1; up >= dn->lsu; up--)
5623 for (cut = DECDPUN - 1; cut >= 0; cut--)
5625 d = u / powers[cut];
5626 u -= d * powers[cut];
5630 if (dn->exponent != 0)
5633 if (dn->exponent < 0)
5635 printf (" E%c%d", esign, abs (dn->exponent));
5637 printf (" [%d]\n", dn->digits);
5641 #if DECTRACE || DECCHECK
5642 /* ------------------------------------------------------------------ */
5643 /* decDumpAr -- display a unit array [debug aid] */
5644 /* name is a single-character tag name */
5645 /* ar is the array to display */
5646 /* len is the length of the array in Units */
5647 /* ------------------------------------------------------------------ */
5649 decDumpAr (char name, Unit * ar, Int len)
5653 char *spec = "%04d ";
5657 printf (" :%c: ", name);
5658 for (i = len - 1; i >= 0; i--)
5661 printf ("%d ", ar[i]);
5663 printf (spec, ar[i]);
5671 /* ------------------------------------------------------------------ */
5672 /* decCheckOperands -- check operand(s) to a routine */
5673 /* res is the result structure (not checked; it will be set to */
5674 /* quiet NaN if error found (and it is not NULL)) */
5675 /* lhs is the first operand (may be DECUNUSED) */
5676 /* rhs is the second (may be DECUNUSED) */
5677 /* set is the context (may be DECUNUSED) */
5678 /* returns 0 if both operands, and the context are clean, or 1 */
5679 /* otherwise (in which case the context will show an error, */
5680 /* unless NULL). Note that res is not cleaned; caller should */
5681 /* handle this so res=NULL case is safe. */
5682 /* The caller is expected to abandon immediately if 1 is returned. */
5683 /* ------------------------------------------------------------------ */
5685 decCheckOperands (decNumber * res, decNumber * lhs,
5686 decNumber * rhs, decContext * set)
5690 { /* oops; hopeless */
5692 printf ("Context is NULL.\n");
5697 else if (set != DECUNUSED
5698 && (set->digits < 1 || set->round < 0
5699 || set->round >= DEC_ROUND_MAX))
5703 printf ("Bad context [digits=%d round=%d].\n", set->digits, set->round);
5712 printf ("Bad result [is NULL].\n");
5715 if (!bad && lhs != DECUNUSED)
5716 bad = (decCheckNumber (lhs, set));
5717 if (!bad && rhs != DECUNUSED)
5718 bad = (decCheckNumber (rhs, set));
5722 if (set != DECUNUSED)
5723 decContextSetStatus (set, DEC_Invalid_operation);
5724 if (res != DECUNUSED && res != NULL)
5726 decNumberZero (res);
5727 res->bits = DECNAN; /* qNaN */
5733 /* ------------------------------------------------------------------ */
5734 /* decCheckNumber -- check a number */
5735 /* dn is the number to check */
5736 /* set is the context (may be DECUNUSED) */
5737 /* returns 0 if the number is clean, or 1 otherwise */
5739 /* The number is considered valid if it could be a result from some */
5740 /* operation in some valid context (not necessarily the current one). */
5741 /* ------------------------------------------------------------------ */
5743 decCheckNumber (decNumber * dn, decContext * set)
5745 Unit *up; /* work */
5746 uInt maxuint; /* .. */
5747 Int ae, d, digits; /* .. */
5748 Int emin, emax; /* .. */
5753 printf ("Reference to decNumber is NULL.\n");
5758 /* check special values */
5759 if (dn->bits & DECSPECIAL)
5761 if (dn->exponent != 0)
5764 printf ("Exponent %d (not 0) for a special value.\n", dn->exponent);
5769 /* 2003.09.08: NaNs may now have coefficients, so next tests Inf only */
5770 if (decNumberIsInfinite (dn))
5772 if (dn->digits != 1)
5775 printf ("Digits %d (not 1) for an infinity.\n", dn->digits);
5782 printf ("LSU %d (not 0) for an infinity.\n", *dn->lsu);
5787 /* 2002.12.26: negative NaNs can now appear through proposed IEEE */
5788 /* concrete formats (decimal64, etc.), though they are */
5789 /* never visible in strings. */
5792 /* if ((dn->bits & DECINF) || (dn->bits & DECNEG)==0) return 0; */
5794 /* printf("Negative NaN in number.\n"); */
5799 /* check the coefficient */
5800 if (dn->digits < 1 || dn->digits > DECNUMMAXP)
5803 printf ("Digits %d in number.\n", dn->digits);
5810 for (up = dn->lsu; d > 0; up++)
5813 maxuint = DECDPUNMAX;
5815 { /* we are at the msu */
5816 maxuint = powers[d] - 1;
5817 if (dn->digits > 1 && *up < powers[d - 1])
5820 printf ("Leading 0 in number.\n");
5829 printf ("Bad Unit [%08x] in number at offset %d [maxuint %d].\n",
5830 *up, up - dn->lsu, maxuint);
5837 /* check the exponent. Note that input operands can have exponents */
5838 /* which are out of the set->emin/set->emax and set->digits range */
5839 /* (just as they can have more digits than set->digits). */
5840 ae = dn->exponent + dn->digits - 1; /* adjusted exponent */
5843 digits = DECNUMMAXP;
5844 if (ae < emin - (digits - 1))
5847 printf ("Adjusted exponent underflow [%d].\n", ae);
5855 printf ("Adjusted exponent overflow [%d].\n", ae);
5861 return 0; /* it's OK */
5868 /* ------------------------------------------------------------------ */
5869 /* decMalloc -- accountable allocation routine */
5870 /* n is the number of bytes to allocate */
5872 /* Semantics is the same as the stdlib malloc routine, but bytes */
5873 /* allocated are accounted for globally, and corruption fences are */
5874 /* added before and after the 'actual' storage. */
5875 /* ------------------------------------------------------------------ */
5876 /* This routine allocates storage with an extra twelve bytes; 8 are */
5877 /* at the start and hold: */
5878 /* 0-3 the original length requested */
5879 /* 4-7 buffer corruption detection fence (DECFENCE, x4) */
5880 /* The 4 bytes at the end also hold a corruption fence (DECFENCE, x4) */
5881 /* ------------------------------------------------------------------ */
5885 uInt size = n + 12; /* true size */
5886 void *alloc; /* -> allocated storage */
5888 uByte *b, *b0; /* .. */
5890 alloc = malloc (size); /* -> allocated storage */
5892 return NULL; /* out of strorage */
5893 b0 = (uByte *) alloc; /* as bytes */
5894 decAllocBytes += n; /* account for storage */
5895 j = (uInt *) alloc; /* -> first four bytes */
5896 *j = n; /* save n */
5897 /* printf("++ alloc(%d)\n", n); */
5898 for (b = b0 + 4; b < b0 + 8; b++)
5900 for (b = b0 + n + 8; b < b0 + n + 12; b++)
5902 return b0 + 8; /* -> play area */
5905 /* ------------------------------------------------------------------ */
5906 /* decFree -- accountable free routine */
5907 /* alloc is the storage to free */
5909 /* Semantics is the same as the stdlib malloc routine, except that */
5910 /* the global storage accounting is updated and the fences are */
5911 /* checked to ensure that no routine has written 'out of bounds'. */
5912 /* ------------------------------------------------------------------ */
5913 /* This routine first checks that the fences have not been corrupted. */
5914 /* It then frees the storage using the 'truw' storage address (that */
5915 /* is, offset by 8). */
5916 /* ------------------------------------------------------------------ */
5918 decFree (void *alloc)
5920 uInt *j, n; /* pointer, original length */
5921 uByte *b, *b0; /* work */
5924 return; /* allowed; it's a nop */
5925 b0 = (uByte *) alloc; /* as bytes */
5926 b0 -= 8; /* -> true start of storage */
5927 j = (uInt *) b0; /* -> first four bytes */
5929 for (b = b0 + 4; b < b0 + 8; b++)
5931 printf ("=== Corrupt byte [%02x] at offset %d from %d ===\n", *b,
5932 b - b0 - 8, (Int) b0);
5933 for (b = b0 + n + 8; b < b0 + n + 12; b++)
5935 printf ("=== Corrupt byte [%02x] at offset +%d from %d, n=%d ===\n", *b,
5936 b - b0 - 8, (Int) b0, n);
5937 free (b0); /* drop the storage */
5938 decAllocBytes -= n; /* account for storage */
OSDN Git Service
