* Makefile.in (TEXI_GCCINSTALL_FILES): Add gcc-vers.texi dependency.
[official-gcc.git] / libdecnumber / decNumber.c
blob039d19a459213a4cdd2a14e894d0c9e8ce7ba1bf
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
10 version.
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
15 for more details.
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, 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301, USA. */
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') */
27 /* arithmetic. */
28 /* */
29 /* Usage notes: */
30 /* */
31 /* 1. This code is ANSI C89 except: */
32 /* */
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.) */
36 /* */
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. */
40 /* */
41 /* c) If DECDPUN>4, non-ANSI 64-bit 'long long' types are used. */
42 /* To avoid these, set DECDPUN <= 4 (see documentation). */
43 /* */
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). */
51 /* */
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. */
55 /* */
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). */
60 /* */
61 /* It is the responsibility of the caller to clear the status */
62 /* flags as required. */
63 /* */
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). */
67 /* */
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. */
73 /* */
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.) */
80 /* */
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. */
84 /* */
85 /* 7. Subset arithmetic is available only if DECSUBSET is set to 1. */
86 /* ------------------------------------------------------------------ */
87 /* Implementation notes for maintenance of this module: */
88 /* */
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. */
95 /* */
96 /* Storage leak accounting can be enabled using DECALLOC. */
97 /* */
98 /* 2. All loops use the for(;;) construct. Any do construct is for */
99 /* protection as just described. */
100 /* */
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). */
107 /* */
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). */
113 /* */
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). */
121 /* */
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). */
127 /* */
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. */
134 /* */
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). */
140 /* */
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 /* Some of glibc's string inlines cause warnings. Plus we'd rather
153 rely on (and therefore test) GCC's string builtins. */
154 #define __NO_STRING_INLINES
156 #include <stdlib.h> /* for malloc, free, etc. */
157 #include <stdio.h> /* for printf [if needed] */
158 #include <string.h> /* for strcpy */
159 #include <ctype.h> /* for lower */
160 #include "config.h"
161 #include "decNumber.h" /* base number library */
162 #include "decNumberLocal.h" /* decNumber local types, etc. */
164 /* Constants */
165 /* Public constant array: powers of ten (powers[n]==10**n) */
166 const uInt powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
167 10000000, 100000000, 1000000000
170 /* Local constants */
171 #define DIVIDE 0x80 /* Divide operators */
172 #define REMAINDER 0x40 /* .. */
173 #define DIVIDEINT 0x20 /* .. */
174 #define REMNEAR 0x10 /* .. */
175 #define COMPARE 0x01 /* Compare operators */
176 #define COMPMAX 0x02 /* .. */
177 #define COMPMIN 0x03 /* .. */
178 #define COMPNAN 0x04 /* .. [NaN processing] */
180 #define DEC_sNaN 0x40000000 /* local status: sNaN signal */
181 #define BADINT (Int)0x80000000 /* most-negative Int; error indicator */
183 static Unit one[] = { 1 }; /* Unit array of 1, used for incrementing */
185 /* Granularity-dependent code */
186 #if DECDPUN<=4
187 #define eInt Int /* extended integer */
188 #define ueInt uInt /* unsigned extended integer */
189 /* Constant multipliers for divide-by-power-of five using reciprocal */
190 /* multiply, after removing powers of 2 by shifting, and final shift */
191 /* of 17 [we only need up to **4] */
192 static const uInt multies[] = { 131073, 26215, 5243, 1049, 210 };
194 /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
195 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
196 #else
197 /* For DECDPUN>4 we currently use non-ANSI 64-bit types. These could */
198 /* be replaced by subroutine calls later. */
199 #ifdef long
200 #undef long
201 #endif
202 typedef signed long long Long;
203 typedef unsigned long long uLong;
204 #define eInt Long /* extended integer */
205 #define ueInt uLong /* unsigned extended integer */
206 #endif
208 /* Local routines */
209 static decNumber *decAddOp (decNumber *, const decNumber *,
210 const decNumber *, decContext *,
211 uByte, uInt *);
212 static void decApplyRound (decNumber *, decContext *, Int, uInt *);
213 static Int decCompare (const decNumber * lhs, const decNumber * rhs);
214 static decNumber *decCompareOp (decNumber *, const decNumber *, const decNumber *,
215 decContext *, Flag, uInt *);
216 static void decCopyFit (decNumber *, const decNumber *, decContext *,
217 Int *, uInt *);
218 static decNumber *decDivideOp (decNumber *, const decNumber *, const decNumber *,
219 decContext *, Flag, uInt *);
220 static void decFinalize (decNumber *, decContext *, Int *, uInt *);
221 static Int decGetDigits (const Unit *, Int);
222 #if DECSUBSET
223 static Int decGetInt (const decNumber *, decContext *);
224 #else
225 static Int decGetInt (const decNumber *);
226 #endif
227 static decNumber *decMultiplyOp (decNumber *, const decNumber *,
228 const decNumber *, decContext *, uInt *);
229 static decNumber *decNaNs (decNumber *, const decNumber *, const decNumber *, uInt *);
230 static decNumber *decQuantizeOp (decNumber *, const decNumber *,
231 const decNumber *, decContext *, Flag, uInt *);
232 static void decSetCoeff (decNumber *, decContext *, const Unit *,
233 Int, Int *, uInt *);
234 static void decSetOverflow (decNumber *, decContext *, uInt *);
235 static void decSetSubnormal (decNumber *, decContext *, Int *, uInt *);
236 static Int decShiftToLeast (Unit *, Int, Int);
237 static Int decShiftToMost (Unit *, Int, Int);
238 static void decStatus (decNumber *, uInt, decContext *);
239 static Flag decStrEq (const char *, const char *);
240 static void decToString (const decNumber *, char[], Flag);
241 static decNumber *decTrim (decNumber *, Flag, Int *);
242 static Int decUnitAddSub (const Unit *, Int, const Unit *, Int, Int, Unit *, Int);
243 static Int decUnitCompare (const Unit *, Int, const Unit *, Int, Int);
245 #if !DECSUBSET
246 /* decFinish == decFinalize when no subset arithmetic needed */
247 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
248 #else
249 static void decFinish (decNumber *, decContext *, Int *, uInt *);
250 static decNumber *decRoundOperand (const decNumber *, decContext *, uInt *);
251 #endif
253 /* Diagnostic macros, etc. */
254 #if DECALLOC
255 /* Handle malloc/free accounting. If enabled, our accountable routines */
256 /* are used; otherwise the code just goes straight to the system malloc */
257 /* and free routines. */
258 #define malloc(a) decMalloc(a)
259 #define free(a) decFree(a)
260 #define DECFENCE 0x5a /* corruption detector */
261 /* 'Our' malloc and free: */
262 static void *decMalloc (size_t);
263 static void decFree (void *);
264 uInt decAllocBytes = 0; /* count of bytes allocated */
265 /* Note that DECALLOC code only checks for storage buffer overflow. */
266 /* To check for memory leaks, the decAllocBytes variable should be */
267 /* checked to be 0 at appropriate times (e.g., after the test */
268 /* harness completes a set of tests). This checking may be unreliable */
269 /* if the testing is done in a multi-thread environment. */
270 #endif
272 #if DECCHECK
273 /* Optional operand checking routines. Enabling these means that */
274 /* decNumber and decContext operands to operator routines are checked */
275 /* for correctness. This roughly doubles the execution time of the */
276 /* fastest routines (and adds 600+ bytes), so should not normally be */
277 /* used in 'production'. */
278 #define DECUNUSED (void *)(0xffffffff)
279 static Flag decCheckOperands (decNumber *, const decNumber *,
280 const decNumber *, decContext *);
281 static Flag decCheckNumber (const decNumber *, decContext *);
282 #endif
284 #if DECTRACE || DECCHECK
285 /* Optional trace/debugging routines. */
286 void decNumberShow (const decNumber *); /* displays the components of a number */
287 static void decDumpAr (char, const Unit *, Int);
288 #endif
290 /* ================================================================== */
291 /* Conversions */
292 /* ================================================================== */
294 /* ------------------------------------------------------------------ */
295 /* to-scientific-string -- conversion to numeric string */
296 /* to-engineering-string -- conversion to numeric string */
297 /* */
298 /* decNumberToString(dn, string); */
299 /* decNumberToEngString(dn, string); */
300 /* */
301 /* dn is the decNumber to convert */
302 /* string is the string where the result will be laid out */
303 /* */
304 /* string must be at least dn->digits+14 characters long */
305 /* */
306 /* No error is possible, and no status can be set. */
307 /* ------------------------------------------------------------------ */
308 char *
309 decNumberToString (const decNumber * dn, char *string)
311 decToString (dn, string, 0);
312 return string;
315 char *
316 decNumberToEngString (const decNumber * dn, char *string)
318 decToString (dn, string, 1);
319 return string;
322 /* ------------------------------------------------------------------ */
323 /* to-number -- conversion from numeric string */
324 /* */
325 /* decNumberFromString -- convert string to decNumber */
326 /* dn -- the number structure to fill */
327 /* chars[] -- the string to convert ('\0' terminated) */
328 /* set -- the context used for processing any error, */
329 /* determining the maximum precision available */
330 /* (set.digits), determining the maximum and minimum */
331 /* exponent (set.emax and set.emin), determining if */
332 /* extended values are allowed, and checking the */
333 /* rounding mode if overflow occurs or rounding is */
334 /* needed. */
335 /* */
336 /* The length of the coefficient and the size of the exponent are */
337 /* checked by this routine, so the correct error (Underflow or */
338 /* Overflow) can be reported or rounding applied, as necessary. */
339 /* */
340 /* If bad syntax is detected, the result will be a quiet NaN. */
341 /* ------------------------------------------------------------------ */
342 decNumber *
343 decNumberFromString (decNumber * dn, const char chars[], decContext * set)
345 Int exponent = 0; /* working exponent [assume 0] */
346 uByte bits = 0; /* working flags [assume +ve] */
347 Unit *res; /* where result will be built */
348 Unit resbuff[D2U (DECBUFFER + 1)]; /* local buffer in case need temporary */
349 Unit *allocres = NULL; /* -> allocated result, iff allocated */
350 Int need; /* units needed for result */
351 Int d = 0; /* count of digits found in decimal part */
352 const char *dotchar = NULL; /* where dot was found */
353 const char *cfirst; /* -> first character of decimal part */
354 const char *last = NULL; /* -> last digit of decimal part */
355 const char *firstexp; /* -> first significant exponent digit */
356 const char *c; /* work */
357 Unit *up; /* .. */
358 #if DECDPUN>1
359 Int i; /* .. */
360 #endif
361 Int residue = 0; /* rounding residue */
362 uInt status = 0; /* error code */
364 #if DECCHECK
365 if (decCheckOperands (DECUNUSED, DECUNUSED, DECUNUSED, set))
366 return decNumberZero (dn);
367 #endif
370 { /* status & malloc protection */
371 c = chars; /* -> input character */
372 if (*c == '-')
373 { /* handle leading '-' */
374 bits = DECNEG;
375 c++;
377 else if (*c == '+')
378 c++; /* step over leading '+' */
379 /* We're at the start of the number [we think] */
380 cfirst = c; /* save */
381 for (;; c++)
383 if (*c >= '0' && *c <= '9')
384 { /* test for Arabic digit */
385 last = c;
386 d++; /* count of real digits */
387 continue; /* still in decimal part */
389 if (*c != '.')
390 break; /* done with decimal part */
391 /* dot: record, check, and ignore */
392 if (dotchar != NULL)
393 { /* two dots */
394 last = NULL; /* indicate bad */
395 break;
396 } /* .. and go report */
397 dotchar = c; /* offset into decimal part */
398 } /* c */
400 if (last == NULL)
401 { /* no decimal digits, or >1 . */
402 #if DECSUBSET
403 /* If subset then infinities and NaNs are not allowed */
404 if (!set->extended)
406 status = DEC_Conversion_syntax;
407 break; /* all done */
409 else
411 #endif
412 /* Infinities and NaNs are possible, here */
413 decNumberZero (dn); /* be optimistic */
414 if (decStrEq (c, "Infinity") || decStrEq (c, "Inf"))
416 dn->bits = bits | DECINF;
417 break; /* all done */
419 else
420 { /* a NaN expected */
421 /* 2003.09.10 NaNs are now permitted to have a sign */
422 status = DEC_Conversion_syntax; /* assume the worst */
423 dn->bits = bits | DECNAN; /* assume simple NaN */
424 if (*c == 's' || *c == 'S')
425 { /* looks like an` sNaN */
426 c++;
427 dn->bits = bits | DECSNAN;
429 if (*c != 'n' && *c != 'N')
430 break; /* check caseless "NaN" */
431 c++;
432 if (*c != 'a' && *c != 'A')
433 break; /* .. */
434 c++;
435 if (*c != 'n' && *c != 'N')
436 break; /* .. */
437 c++;
438 /* now nothing, or nnnn, expected */
439 /* -> start of integer and skip leading 0s [including plain 0] */
440 for (cfirst = c; *cfirst == '0';)
441 cfirst++;
442 if (*cfirst == '\0')
443 { /* "NaN" or "sNaN", maybe with all 0s */
444 status = 0; /* it's good */
445 break; /* .. */
447 /* something other than 0s; setup last and d as usual [no dots] */
448 for (c = cfirst;; c++, d++)
450 if (*c < '0' || *c > '9')
451 break; /* test for Arabic digit */
452 last = c;
454 if (*c != '\0')
455 break; /* not all digits */
456 if (d > set->digits)
457 break; /* too many digits */
458 /* good; drop through and convert the integer */
459 status = 0;
460 bits = dn->bits; /* for copy-back */
461 } /* NaN expected */
462 #if DECSUBSET
464 #endif
465 } /* last==NULL */
467 if (*c != '\0')
468 { /* more there; exponent expected... */
469 Flag nege = 0; /* 1=negative exponent */
470 if (*c != 'e' && *c != 'E')
472 status = DEC_Conversion_syntax;
473 break;
476 /* Found 'e' or 'E' -- now process explicit exponent */
477 /* 1998.07.11: sign no longer required */
478 c++; /* to (expected) sign */
479 if (*c == '-')
481 nege = 1;
482 c++;
484 else if (*c == '+')
485 c++;
486 if (*c == '\0')
488 status = DEC_Conversion_syntax;
489 break;
492 for (; *c == '0' && *(c + 1) != '\0';)
493 c++; /* strip insignificant zeros */
494 firstexp = c; /* save exponent digit place */
495 for (;; c++)
497 if (*c < '0' || *c > '9')
498 break; /* not a digit */
499 exponent = X10 (exponent) + (Int) * c - (Int) '0';
500 } /* c */
501 /* if we didn't end on '\0' must not be a digit */
502 if (*c != '\0')
504 status = DEC_Conversion_syntax;
505 break;
508 /* (this next test must be after the syntax check) */
509 /* if it was too long the exponent may have wrapped, so check */
510 /* carefully and set it to a certain overflow if wrap possible */
511 if (c >= firstexp + 9 + 1)
513 if (c > firstexp + 9 + 1 || *firstexp > '1')
514 exponent = DECNUMMAXE * 2;
515 /* [up to 1999999999 is OK, for example 1E-1000000998] */
517 if (nege)
518 exponent = -exponent; /* was negative */
519 } /* had exponent */
520 /* Here when all inspected; syntax is good */
522 /* Handle decimal point... */
523 if (dotchar != NULL && dotchar < last) /* embedded . found, so */
524 exponent = exponent - (last - dotchar); /* .. adjust exponent */
525 /* [we can now ignore the .] */
527 /* strip leading zeros/dot (leave final if all 0's) */
528 for (c = cfirst; c < last; c++)
530 if (*c == '0')
531 d--; /* 0 stripped */
532 else if (*c != '.')
533 break;
534 cfirst++; /* step past leader */
535 } /* c */
537 #if DECSUBSET
538 /* We can now make a rapid exit for zeros if !extended */
539 if (*cfirst == '0' && !set->extended)
541 decNumberZero (dn); /* clean result */
542 break; /* [could be return] */
544 #endif
546 /* OK, the digits string is good. Copy to the decNumber, or to
547 a temporary decNumber if rounding is needed */
548 if (d <= set->digits)
549 res = dn->lsu; /* fits into given decNumber */
550 else
551 { /* rounding needed */
552 need = D2U (d); /* units needed */
553 res = resbuff; /* assume use local buffer */
554 if (need * sizeof (Unit) > sizeof (resbuff))
555 { /* too big for local */
556 allocres = (Unit *) malloc (need * sizeof (Unit));
557 if (allocres == NULL)
559 status |= DEC_Insufficient_storage;
560 break;
562 res = allocres;
565 /* res now -> number lsu, buffer, or allocated storage for Unit array */
567 /* Place the coefficient into the selected Unit array */
568 #if DECDPUN>1
569 i = d % DECDPUN; /* digits in top unit */
570 if (i == 0)
571 i = DECDPUN;
572 up = res + D2U (d) - 1; /* -> msu */
573 *up = 0;
574 for (c = cfirst;; c++)
575 { /* along the digits */
576 if (*c == '.')
577 { /* ignore . [don't decrement i] */
578 if (c != last)
579 continue;
580 break;
582 *up = (Unit) (X10 (*up) + (Int) * c - (Int) '0');
583 i--;
584 if (i > 0)
585 continue; /* more for this unit */
586 if (up == res)
587 break; /* just filled the last unit */
588 i = DECDPUN;
589 up--;
590 *up = 0;
591 } /* c */
592 #else
593 /* DECDPUN==1 */
594 up = res; /* -> lsu */
595 for (c = last; c >= cfirst; c--)
596 { /* over each character, from least */
597 if (*c == '.')
598 continue; /* ignore . [don't step b] */
599 *up = (Unit) ((Int) * c - (Int) '0');
600 up++;
601 } /* c */
602 #endif
604 dn->bits = bits;
605 dn->exponent = exponent;
606 dn->digits = d;
608 /* if not in number (too long) shorten into the number */
609 if (d > set->digits)
610 decSetCoeff (dn, set, res, d, &residue, &status);
612 /* Finally check for overflow or subnormal and round as needed */
613 decFinalize (dn, set, &residue, &status);
614 /* decNumberShow(dn); */
616 while (0); /* [for break] */
618 if (allocres != NULL)
619 free (allocres); /* drop any storage we used */
620 if (status != 0)
621 decStatus (dn, status, set);
622 return dn;
625 /* ================================================================== */
626 /* Operators */
627 /* ================================================================== */
629 /* ------------------------------------------------------------------ */
630 /* decNumberAbs -- absolute value operator */
631 /* */
632 /* This computes C = abs(A) */
633 /* */
634 /* res is C, the result. C may be A */
635 /* rhs is A */
636 /* set is the context */
637 /* */
638 /* C must have space for set->digits digits. */
639 /* ------------------------------------------------------------------ */
640 /* This has the same effect as decNumberPlus unless A is negative, */
641 /* in which case it has the same effect as decNumberMinus. */
642 /* ------------------------------------------------------------------ */
643 decNumber *
644 decNumberAbs (decNumber * res, const decNumber * rhs, decContext * set)
646 decNumber dzero; /* for 0 */
647 uInt status = 0; /* accumulator */
649 #if DECCHECK
650 if (decCheckOperands (res, DECUNUSED, rhs, set))
651 return res;
652 #endif
654 decNumberZero (&dzero); /* set 0 */
655 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
656 decAddOp (res, &dzero, rhs, set, (uByte) (rhs->bits & DECNEG), &status);
657 if (status != 0)
658 decStatus (res, status, set);
659 return res;
662 /* ------------------------------------------------------------------ */
663 /* decNumberAdd -- add two Numbers */
664 /* */
665 /* This computes C = A + B */
666 /* */
667 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
668 /* lhs is A */
669 /* rhs is B */
670 /* set is the context */
671 /* */
672 /* C must have space for set->digits digits. */
673 /* ------------------------------------------------------------------ */
674 /* This just calls the routine shared with Subtract */
675 decNumber *
676 decNumberAdd (decNumber * res, const decNumber * lhs,
677 const decNumber * rhs, decContext * set)
679 uInt status = 0; /* accumulator */
680 decAddOp (res, lhs, rhs, set, 0, &status);
681 if (status != 0)
682 decStatus (res, status, set);
683 return res;
686 /* ------------------------------------------------------------------ */
687 /* decNumberCompare -- compare two Numbers */
688 /* */
689 /* This computes C = A ? B */
690 /* */
691 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
692 /* lhs is A */
693 /* rhs is B */
694 /* set is the context */
695 /* */
696 /* C must have space for one digit. */
697 /* ------------------------------------------------------------------ */
698 decNumber *
699 decNumberCompare (decNumber * res, const decNumber * lhs,
700 const decNumber * rhs, decContext * set)
702 uInt status = 0; /* accumulator */
703 decCompareOp (res, lhs, rhs, set, COMPARE, &status);
704 if (status != 0)
705 decStatus (res, status, set);
706 return res;
709 /* ------------------------------------------------------------------ */
710 /* decNumberDivide -- divide one number by another */
711 /* */
712 /* This computes C = A / B */
713 /* */
714 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
715 /* lhs is A */
716 /* rhs is B */
717 /* set is the context */
718 /* */
719 /* C must have space for set->digits digits. */
720 /* ------------------------------------------------------------------ */
721 decNumber *
722 decNumberDivide (decNumber * res, const decNumber * lhs,
723 const decNumber * rhs, decContext * set)
725 uInt status = 0; /* accumulator */
726 decDivideOp (res, lhs, rhs, set, DIVIDE, &status);
727 if (status != 0)
728 decStatus (res, status, set);
729 return res;
732 /* ------------------------------------------------------------------ */
733 /* decNumberDivideInteger -- divide and return integer quotient */
734 /* */
735 /* This computes C = A # B, where # is the integer divide operator */
736 /* */
737 /* res is C, the result. C may be A and/or B (e.g., X=X#X) */
738 /* lhs is A */
739 /* rhs is B */
740 /* set is the context */
741 /* */
742 /* C must have space for set->digits digits. */
743 /* ------------------------------------------------------------------ */
744 decNumber *
745 decNumberDivideInteger (decNumber * res, const decNumber * lhs,
746 const decNumber * rhs, decContext * set)
748 uInt status = 0; /* accumulator */
749 decDivideOp (res, lhs, rhs, set, DIVIDEINT, &status);
750 if (status != 0)
751 decStatus (res, status, set);
752 return res;
755 /* ------------------------------------------------------------------ */
756 /* decNumberMax -- compare two Numbers and return the maximum */
757 /* */
758 /* This computes C = A ? B, returning the maximum or A if equal */
759 /* */
760 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
761 /* lhs is A */
762 /* rhs is B */
763 /* set is the context */
764 /* */
765 /* C must have space for set->digits digits. */
766 /* ------------------------------------------------------------------ */
767 decNumber *
768 decNumberMax (decNumber * res, const decNumber * lhs,
769 const decNumber * rhs, decContext * set)
771 uInt status = 0; /* accumulator */
772 decCompareOp (res, lhs, rhs, set, COMPMAX, &status);
773 if (status != 0)
774 decStatus (res, status, set);
775 return res;
778 /* ------------------------------------------------------------------ */
779 /* decNumberMin -- compare two Numbers and return the minimum */
780 /* */
781 /* This computes C = A ? B, returning the minimum or A if equal */
782 /* */
783 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
784 /* lhs is A */
785 /* rhs is B */
786 /* set is the context */
787 /* */
788 /* C must have space for set->digits digits. */
789 /* ------------------------------------------------------------------ */
790 decNumber *
791 decNumberMin (decNumber * res, const decNumber * lhs,
792 const decNumber * rhs, decContext * set)
794 uInt status = 0; /* accumulator */
795 decCompareOp (res, lhs, rhs, set, COMPMIN, &status);
796 if (status != 0)
797 decStatus (res, status, set);
798 return res;
801 /* ------------------------------------------------------------------ */
802 /* decNumberMinus -- prefix minus operator */
803 /* */
804 /* This computes C = 0 - A */
805 /* */
806 /* res is C, the result. C may be A */
807 /* rhs is A */
808 /* set is the context */
809 /* */
810 /* C must have space for set->digits digits. */
811 /* ------------------------------------------------------------------ */
812 /* We simply use AddOp for the subtract, which will do the necessary. */
813 /* ------------------------------------------------------------------ */
814 decNumber *
815 decNumberMinus (decNumber * res, const decNumber * rhs, decContext * set)
817 decNumber dzero;
818 uInt status = 0; /* accumulator */
820 #if DECCHECK
821 if (decCheckOperands (res, DECUNUSED, rhs, set))
822 return res;
823 #endif
825 decNumberZero (&dzero); /* make 0 */
826 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
827 decAddOp (res, &dzero, rhs, set, DECNEG, &status);
828 if (status != 0)
829 decStatus (res, status, set);
830 return res;
833 /* ------------------------------------------------------------------ */
834 /* decNumberPlus -- prefix plus operator */
835 /* */
836 /* This computes C = 0 + A */
837 /* */
838 /* res is C, the result. C may be A */
839 /* rhs is A */
840 /* set is the context */
841 /* */
842 /* C must have space for set->digits digits. */
843 /* ------------------------------------------------------------------ */
844 /* We simply use AddOp; Add will take fast path after preparing A. */
845 /* Performance is a concern here, as this routine is often used to */
846 /* check operands and apply rounding and overflow/underflow testing. */
847 /* ------------------------------------------------------------------ */
848 decNumber *
849 decNumberPlus (decNumber * res, const decNumber * rhs, decContext * set)
851 decNumber dzero;
852 uInt status = 0; /* accumulator */
854 #if DECCHECK
855 if (decCheckOperands (res, DECUNUSED, rhs, set))
856 return res;
857 #endif
859 decNumberZero (&dzero); /* make 0 */
860 dzero.exponent = rhs->exponent; /* [no coefficient expansion] */
861 decAddOp (res, &dzero, rhs, set, 0, &status);
862 if (status != 0)
863 decStatus (res, status, set);
864 return res;
867 /* ------------------------------------------------------------------ */
868 /* decNumberMultiply -- multiply two Numbers */
869 /* */
870 /* This computes C = A x B */
871 /* */
872 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
873 /* lhs is A */
874 /* rhs is B */
875 /* set is the context */
876 /* */
877 /* C must have space for set->digits digits. */
878 /* ------------------------------------------------------------------ */
879 decNumber *
880 decNumberMultiply (decNumber * res, const decNumber * lhs,
881 const decNumber * rhs, decContext * set)
883 uInt status = 0; /* accumulator */
884 decMultiplyOp (res, lhs, rhs, set, &status);
885 if (status != 0)
886 decStatus (res, status, set);
887 return res;
890 /* ------------------------------------------------------------------ */
891 /* decNumberNormalize -- remove trailing zeros */
892 /* */
893 /* This computes C = 0 + A, and normalizes the result */
894 /* */
895 /* res is C, the result. C may be A */
896 /* rhs is A */
897 /* set is the context */
898 /* */
899 /* C must have space for set->digits digits. */
900 /* ------------------------------------------------------------------ */
901 decNumber *
902 decNumberNormalize (decNumber * res, const decNumber * rhs, decContext * set)
904 decNumber *allocrhs = NULL; /* non-NULL if rounded rhs allocated */
905 uInt status = 0; /* as usual */
906 Int residue = 0; /* as usual */
907 Int dropped; /* work */
909 #if DECCHECK
910 if (decCheckOperands (res, DECUNUSED, rhs, set))
911 return res;
912 #endif
915 { /* protect allocated storage */
916 #if DECSUBSET
917 if (!set->extended)
919 /* reduce operand and set lostDigits status, as needed */
920 if (rhs->digits > set->digits)
922 allocrhs = decRoundOperand (rhs, set, &status);
923 if (allocrhs == NULL)
924 break;
925 rhs = allocrhs;
928 #endif
929 /* [following code does not require input rounding] */
931 /* specials copy through, except NaNs need care */
932 if (decNumberIsNaN (rhs))
934 decNaNs (res, rhs, NULL, &status);
935 break;
938 /* reduce result to the requested length and copy to result */
939 decCopyFit (res, rhs, set, &residue, &status); /* copy & round */
940 decFinish (res, set, &residue, &status); /* cleanup/set flags */
941 decTrim (res, 1, &dropped); /* normalize in place */
943 while (0); /* end protected */
945 if (allocrhs != NULL)
946 free (allocrhs); /* .. */
947 if (status != 0)
948 decStatus (res, status, set); /* then report status */
949 return res;
952 /* ------------------------------------------------------------------ */
953 /* decNumberPower -- raise a number to an integer power */
954 /* */
955 /* This computes C = A ** B */
956 /* */
957 /* res is C, the result. C may be A and/or B (e.g., X=X**X) */
958 /* lhs is A */
959 /* rhs is B */
960 /* set is the context */
961 /* */
962 /* C must have space for set->digits digits. */
963 /* */
964 /* Specification restriction: abs(n) must be <=999999999 */
965 /* ------------------------------------------------------------------ */
966 decNumber *
967 decNumberPower (decNumber * res, const decNumber * lhs,
968 const decNumber * rhs, decContext * set)
970 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
971 decNumber *allocrhs = NULL; /* .., rhs */
972 decNumber *allocdac = NULL; /* -> allocated acc buffer, iff used */
973 const decNumber *inrhs = rhs; /* save original rhs */
974 Int reqdigits = set->digits; /* requested DIGITS */
975 Int n; /* RHS in binary */
976 Int i; /* work */
977 #if DECSUBSET
978 Int dropped; /* .. */
979 #endif
980 uInt needbytes; /* buffer size needed */
981 Flag seenbit; /* seen a bit while powering */
982 Int residue = 0; /* rounding residue */
983 uInt status = 0; /* accumulator */
984 uByte bits = 0; /* result sign if errors */
985 decContext workset; /* working context */
986 decNumber dnOne; /* work value 1... */
987 /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
988 uByte dacbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
989 /* same again for possible 1/lhs calculation */
990 uByte lhsbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
991 decNumber *dac = (decNumber *) dacbuff; /* -> result accumulator */
993 #if DECCHECK
994 if (decCheckOperands (res, lhs, rhs, set))
995 return res;
996 #endif
999 { /* protect allocated storage */
1000 #if DECSUBSET
1001 if (!set->extended)
1003 /* reduce operands and set lostDigits status, as needed */
1004 if (lhs->digits > reqdigits)
1006 alloclhs = decRoundOperand (lhs, set, &status);
1007 if (alloclhs == NULL)
1008 break;
1009 lhs = alloclhs;
1011 /* rounding won't affect the result, but we might signal lostDigits */
1012 /* as well as the error for non-integer [x**y would need this too] */
1013 if (rhs->digits > reqdigits)
1015 allocrhs = decRoundOperand (rhs, set, &status);
1016 if (allocrhs == NULL)
1017 break;
1018 rhs = allocrhs;
1021 #endif
1022 /* [following code does not require input rounding] */
1024 /* handle rhs Infinity */
1025 if (decNumberIsInfinite (rhs))
1027 status |= DEC_Invalid_operation; /* bad */
1028 break;
1030 /* handle NaNs */
1031 if ((lhs->bits | rhs->bits) & (DECNAN | DECSNAN))
1033 decNaNs (res, lhs, rhs, &status);
1034 break;
1037 /* Original rhs must be an integer that fits and is in range */
1038 #if DECSUBSET
1039 n = decGetInt (inrhs, set);
1040 #else
1041 n = decGetInt (inrhs);
1042 #endif
1043 if (n == BADINT || n > 999999999 || n < -999999999)
1045 status |= DEC_Invalid_operation;
1046 break;
1048 if (n < 0)
1049 { /* negative */
1050 n = -n; /* use the absolute value */
1052 if (decNumberIsNegative (lhs) /* -x .. */
1053 && (n & 0x00000001))
1054 bits = DECNEG; /* .. to an odd power */
1056 /* handle LHS infinity */
1057 if (decNumberIsInfinite (lhs))
1058 { /* [NaNs already handled] */
1059 uByte rbits = rhs->bits; /* save */
1060 decNumberZero (res);
1061 if (n == 0)
1062 *res->lsu = 1; /* [-]Inf**0 => 1 */
1063 else
1065 if (!(rbits & DECNEG))
1066 bits |= DECINF; /* was not a **-n */
1067 /* [otherwise will be 0 or -0] */
1068 res->bits = bits;
1070 break;
1073 /* clone the context */
1074 workset = *set; /* copy all fields */
1075 /* calculate the working DIGITS */
1076 workset.digits = reqdigits + (inrhs->digits + inrhs->exponent) + 1;
1077 /* it's an error if this is more than we can handle */
1078 if (workset.digits > DECNUMMAXP)
1080 status |= DEC_Invalid_operation;
1081 break;
1084 /* workset.digits is the count of digits for the accumulator we need */
1085 /* if accumulator is too long for local storage, then allocate */
1086 needbytes =
1087 sizeof (decNumber) + (D2U (workset.digits) - 1) * sizeof (Unit);
1088 /* [needbytes also used below if 1/lhs needed] */
1089 if (needbytes > sizeof (dacbuff))
1091 allocdac = (decNumber *) malloc (needbytes);
1092 if (allocdac == NULL)
1093 { /* hopeless -- abandon */
1094 status |= DEC_Insufficient_storage;
1095 break;
1097 dac = allocdac; /* use the allocated space */
1099 decNumberZero (dac); /* acc=1 */
1100 *dac->lsu = 1; /* .. */
1102 if (n == 0)
1103 { /* x**0 is usually 1 */
1104 /* 0**0 is bad unless subset, when it becomes 1 */
1105 if (ISZERO (lhs)
1106 #if DECSUBSET
1107 && set->extended
1108 #endif
1110 status |= DEC_Invalid_operation;
1111 else
1112 decNumberCopy (res, dac); /* copy the 1 */
1113 break;
1116 /* if a negative power we'll need the constant 1, and if not subset */
1117 /* we'll invert the lhs now rather than inverting the result later */
1118 if (decNumberIsNegative (rhs))
1119 { /* was a **-n [hence digits>0] */
1120 decNumber * newlhs;
1121 decNumberCopy (&dnOne, dac); /* dnOne=1; [needed now or later] */
1122 #if DECSUBSET
1123 if (set->extended)
1124 { /* need to calculate 1/lhs */
1125 #endif
1126 /* divide lhs into 1, putting result in dac [dac=1/dac] */
1127 decDivideOp (dac, &dnOne, lhs, &workset, DIVIDE, &status);
1128 if (alloclhs != NULL)
1130 free (alloclhs); /* done with intermediate */
1131 alloclhs = NULL; /* indicate freed */
1133 /* now locate or allocate space for the inverted lhs */
1134 if (needbytes > sizeof (lhsbuff))
1136 alloclhs = (decNumber *) malloc (needbytes);
1137 if (alloclhs == NULL)
1138 { /* hopeless -- abandon */
1139 status |= DEC_Insufficient_storage;
1140 break;
1142 newlhs = alloclhs; /* use the allocated space */
1144 else
1145 newlhs = (decNumber *) lhsbuff; /* use stack storage */
1146 /* [lhs now points to buffer or allocated storage] */
1147 decNumberCopy (newlhs, dac); /* copy the 1/lhs */
1148 decNumberCopy (dac, &dnOne); /* restore acc=1 */
1149 lhs = newlhs;
1150 #if DECSUBSET
1152 #endif
1155 /* Raise-to-the-power loop... */
1156 seenbit = 0; /* set once we've seen a 1-bit */
1157 for (i = 1;; i++)
1158 { /* for each bit [top bit ignored] */
1159 /* abandon if we have had overflow or terminal underflow */
1160 if (status & (DEC_Overflow | DEC_Underflow))
1161 { /* interesting? */
1162 if (status & DEC_Overflow || ISZERO (dac))
1163 break;
1165 /* [the following two lines revealed an optimizer bug in a C++ */
1166 /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
1167 n = n << 1; /* move next bit to testable position */
1168 if (n < 0)
1169 { /* top bit is set */
1170 seenbit = 1; /* OK, we're off */
1171 decMultiplyOp (dac, dac, lhs, &workset, &status); /* dac=dac*x */
1173 if (i == 31)
1174 break; /* that was the last bit */
1175 if (!seenbit)
1176 continue; /* we don't have to square 1 */
1177 decMultiplyOp (dac, dac, dac, &workset, &status); /* dac=dac*dac [square] */
1178 } /*i *//* 32 bits */
1180 /* complete internal overflow or underflow processing */
1181 if (status & (DEC_Overflow | DEC_Subnormal))
1183 #if DECSUBSET
1184 /* If subset, and power was negative, reverse the kind of -erflow */
1185 /* [1/x not yet done] */
1186 if (!set->extended && decNumberIsNegative (rhs))
1188 if (status & DEC_Overflow)
1189 status ^= DEC_Overflow | DEC_Underflow | DEC_Subnormal;
1190 else
1191 { /* trickier -- Underflow may or may not be set */
1192 status &= ~(DEC_Underflow | DEC_Subnormal); /* [one or both] */
1193 status |= DEC_Overflow;
1196 #endif
1197 dac->bits = (dac->bits & ~DECNEG) | bits; /* force correct sign */
1198 /* round subnormals [to set.digits rather than workset.digits] */
1199 /* or set overflow result similarly as required */
1200 decFinalize (dac, set, &residue, &status);
1201 decNumberCopy (res, dac); /* copy to result (is now OK length) */
1202 break;
1205 #if DECSUBSET
1206 if (!set->extended && /* subset math */
1207 decNumberIsNegative (rhs))
1208 { /* was a **-n [hence digits>0] */
1209 /* so divide result into 1 [dac=1/dac] */
1210 decDivideOp (dac, &dnOne, dac, &workset, DIVIDE, &status);
1212 #endif
1214 /* reduce result to the requested length and copy to result */
1215 decCopyFit (res, dac, set, &residue, &status);
1216 decFinish (res, set, &residue, &status); /* final cleanup */
1217 #if DECSUBSET
1218 if (!set->extended)
1219 decTrim (res, 0, &dropped); /* trailing zeros */
1220 #endif
1222 while (0); /* end protected */
1224 if (allocdac != NULL)
1225 free (allocdac); /* drop any storage we used */
1226 if (allocrhs != NULL)
1227 free (allocrhs); /* .. */
1228 if (alloclhs != NULL)
1229 free (alloclhs); /* .. */
1230 if (status != 0)
1231 decStatus (res, status, set);
1232 return res;
1235 /* ------------------------------------------------------------------ */
1236 /* decNumberQuantize -- force exponent to requested value */
1237 /* */
1238 /* This computes C = op(A, B), where op adjusts the coefficient */
1239 /* of C (by rounding or shifting) such that the exponent (-scale) */
1240 /* of C has exponent of B. The numerical value of C will equal A, */
1241 /* except for the effects of any rounding that occurred. */
1242 /* */
1243 /* res is C, the result. C may be A or B */
1244 /* lhs is A, the number to adjust */
1245 /* rhs is B, the number with exponent to match */
1246 /* set is the context */
1247 /* */
1248 /* C must have space for set->digits digits. */
1249 /* */
1250 /* Unless there is an error or the result is infinite, the exponent */
1251 /* after the operation is guaranteed to be equal to that of B. */
1252 /* ------------------------------------------------------------------ */
1253 decNumber *
1254 decNumberQuantize (decNumber * res, const decNumber * lhs,
1255 const decNumber * rhs, decContext * set)
1257 uInt status = 0; /* accumulator */
1258 decQuantizeOp (res, lhs, rhs, set, 1, &status);
1259 if (status != 0)
1260 decStatus (res, status, set);
1261 return res;
1264 /* ------------------------------------------------------------------ */
1265 /* decNumberRescale -- force exponent to requested value */
1266 /* */
1267 /* This computes C = op(A, B), where op adjusts the coefficient */
1268 /* of C (by rounding or shifting) such that the exponent (-scale) */
1269 /* of C has the value B. The numerical value of C will equal A, */
1270 /* except for the effects of any rounding that occurred. */
1271 /* */
1272 /* res is C, the result. C may be A or B */
1273 /* lhs is A, the number to adjust */
1274 /* rhs is B, the requested exponent */
1275 /* set is the context */
1276 /* */
1277 /* C must have space for set->digits digits. */
1278 /* */
1279 /* Unless there is an error or the result is infinite, the exponent */
1280 /* after the operation is guaranteed to be equal to B. */
1281 /* ------------------------------------------------------------------ */
1282 decNumber *
1283 decNumberRescale (decNumber * res, const decNumber * lhs,
1284 const decNumber * rhs, decContext * set)
1286 uInt status = 0; /* accumulator */
1287 decQuantizeOp (res, lhs, rhs, set, 0, &status);
1288 if (status != 0)
1289 decStatus (res, status, set);
1290 return res;
1293 /* ------------------------------------------------------------------ */
1294 /* decNumberRemainder -- divide and return remainder */
1295 /* */
1296 /* This computes C = A % B */
1297 /* */
1298 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1299 /* lhs is A */
1300 /* rhs is B */
1301 /* set is the context */
1302 /* */
1303 /* C must have space for set->digits digits. */
1304 /* ------------------------------------------------------------------ */
1305 decNumber *
1306 decNumberRemainder (decNumber * res, const decNumber * lhs,
1307 const decNumber * rhs, decContext * set)
1309 uInt status = 0; /* accumulator */
1310 decDivideOp (res, lhs, rhs, set, REMAINDER, &status);
1311 if (status != 0)
1312 decStatus (res, status, set);
1313 return res;
1316 /* ------------------------------------------------------------------ */
1317 /* decNumberRemainderNear -- divide and return remainder from nearest */
1318 /* */
1319 /* This computes C = A % B, where % is the IEEE remainder operator */
1320 /* */
1321 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1322 /* lhs is A */
1323 /* rhs is B */
1324 /* set is the context */
1325 /* */
1326 /* C must have space for set->digits digits. */
1327 /* ------------------------------------------------------------------ */
1328 decNumber *
1329 decNumberRemainderNear (decNumber * res, const decNumber * lhs,
1330 const decNumber * rhs, decContext * set)
1332 uInt status = 0; /* accumulator */
1333 decDivideOp (res, lhs, rhs, set, REMNEAR, &status);
1334 if (status != 0)
1335 decStatus (res, status, set);
1336 return res;
1339 /* ------------------------------------------------------------------ */
1340 /* decNumberSameQuantum -- test for equal exponents */
1341 /* */
1342 /* res is the result number, which will contain either 0 or 1 */
1343 /* lhs is a number to test */
1344 /* rhs is the second (usually a pattern) */
1345 /* */
1346 /* No errors are possible and no context is needed. */
1347 /* ------------------------------------------------------------------ */
1348 decNumber *
1349 decNumberSameQuantum (decNumber * res, const decNumber * lhs, const decNumber * rhs)
1351 uByte merged; /* merged flags */
1352 Unit ret = 0; /* return value */
1354 #if DECCHECK
1355 if (decCheckOperands (res, lhs, rhs, DECUNUSED))
1356 return res;
1357 #endif
1359 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
1360 if (merged)
1362 if (decNumberIsNaN (lhs) && decNumberIsNaN (rhs))
1363 ret = 1;
1364 else if (decNumberIsInfinite (lhs) && decNumberIsInfinite (rhs))
1365 ret = 1;
1366 /* [anything else with a special gives 0] */
1368 else if (lhs->exponent == rhs->exponent)
1369 ret = 1;
1371 decNumberZero (res); /* OK to overwrite an operand */
1372 *res->lsu = ret;
1373 return res;
1376 /* ------------------------------------------------------------------ */
1377 /* decNumberSquareRoot -- square root operator */
1378 /* */
1379 /* This computes C = squareroot(A) */
1380 /* */
1381 /* res is C, the result. C may be A */
1382 /* rhs is A */
1383 /* set is the context; note that rounding mode has no effect */
1384 /* */
1385 /* C must have space for set->digits digits. */
1386 /* ------------------------------------------------------------------ */
1387 /* This uses the following varying-precision algorithm in: */
1388 /* */
1389 /* Properly Rounded Variable Precision Square Root, T. E. Hull and */
1390 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
1391 /* pp229-237, ACM, September 1985. */
1392 /* */
1393 /* % [Reformatted original Numerical Turing source code follows.] */
1394 /* function sqrt(x : real) : real */
1395 /* % sqrt(x) returns the properly rounded approximation to the square */
1396 /* % root of x, in the precision of the calling environment, or it */
1397 /* % fails if x < 0. */
1398 /* % t e hull and a abrham, august, 1984 */
1399 /* if x <= 0 then */
1400 /* if x < 0 then */
1401 /* assert false */
1402 /* else */
1403 /* result 0 */
1404 /* end if */
1405 /* end if */
1406 /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
1407 /* var e := getexp(x) % exponent part of x */
1408 /* var approx : real */
1409 /* if e mod 2 = 0 then */
1410 /* approx := .259 + .819 * f % approx to root of f */
1411 /* else */
1412 /* f := f/l0 % adjustments */
1413 /* e := e + 1 % for odd */
1414 /* approx := .0819 + 2.59 * f % exponent */
1415 /* end if */
1416 /* */
1417 /* var p:= 3 */
1418 /* const maxp := currentprecision + 2 */
1419 /* loop */
1420 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
1421 /* precision p */
1422 /* approx := .5 * (approx + f/approx) */
1423 /* exit when p = maxp */
1424 /* end loop */
1425 /* */
1426 /* % approx is now within 1 ulp of the properly rounded square root */
1427 /* % of f; to ensure proper rounding, compare squares of (approx - */
1428 /* % l/2 ulp) and (approx + l/2 ulp) with f. */
1429 /* p := currentprecision */
1430 /* begin */
1431 /* precision p + 2 */
1432 /* const approxsubhalf := approx - setexp(.5, -p) */
1433 /* if mulru(approxsubhalf, approxsubhalf) > f then */
1434 /* approx := approx - setexp(.l, -p + 1) */
1435 /* else */
1436 /* const approxaddhalf := approx + setexp(.5, -p) */
1437 /* if mulrd(approxaddhalf, approxaddhalf) < f then */
1438 /* approx := approx + setexp(.l, -p + 1) */
1439 /* end if */
1440 /* end if */
1441 /* end */
1442 /* result setexp(approx, e div 2) % fix exponent */
1443 /* end sqrt */
1444 /* ------------------------------------------------------------------ */
1445 decNumber *
1446 decNumberSquareRoot (decNumber * res, const decNumber * rhs, decContext * set)
1448 decContext workset, approxset; /* work contexts */
1449 decNumber dzero; /* used for constant zero */
1450 Int maxp = set->digits + 2; /* largest working precision */
1451 Int residue = 0; /* rounding residue */
1452 uInt status = 0, ignore = 0; /* status accumulators */
1453 Int exp; /* working exponent */
1454 Int ideal; /* ideal (preferred) exponent */
1455 uInt needbytes; /* work */
1456 Int dropped; /* .. */
1458 decNumber *allocrhs = NULL; /* non-NULL if rounded rhs allocated */
1459 /* buffer for f [needs +1 in case DECBUFFER 0] */
1460 uByte buff[sizeof (decNumber) + (D2U (DECBUFFER + 1) - 1) * sizeof (Unit)];
1461 /* buffer for a [needs +2 to match maxp] */
1462 uByte bufa[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1463 /* buffer for temporary, b [must be same size as a] */
1464 uByte bufb[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1465 decNumber *allocbuff = NULL; /* -> allocated buff, iff allocated */
1466 decNumber *allocbufa = NULL; /* -> allocated bufa, iff allocated */
1467 decNumber *allocbufb = NULL; /* -> allocated bufb, iff allocated */
1468 decNumber *f = (decNumber *) buff; /* reduced fraction */
1469 decNumber *a = (decNumber *) bufa; /* approximation to result */
1470 decNumber *b = (decNumber *) bufb; /* intermediate result */
1471 /* buffer for temporary variable, up to 3 digits */
1472 uByte buft[sizeof (decNumber) + (D2U (3) - 1) * sizeof (Unit)];
1473 decNumber *t = (decNumber *) buft; /* up-to-3-digit constant or work */
1475 #if DECCHECK
1476 if (decCheckOperands (res, DECUNUSED, rhs, set))
1477 return res;
1478 #endif
1481 { /* protect allocated storage */
1482 #if DECSUBSET
1483 if (!set->extended)
1485 /* reduce operand and set lostDigits status, as needed */
1486 if (rhs->digits > set->digits)
1488 allocrhs = decRoundOperand (rhs, set, &status);
1489 if (allocrhs == NULL)
1490 break;
1491 /* [Note: 'f' allocation below could reuse this buffer if */
1492 /* used, but as this is rare we keep them separate for clarity.] */
1493 rhs = allocrhs;
1496 #endif
1497 /* [following code does not require input rounding] */
1499 /* handle infinities and NaNs */
1500 if (rhs->bits & DECSPECIAL)
1502 if (decNumberIsInfinite (rhs))
1503 { /* an infinity */
1504 if (decNumberIsNegative (rhs))
1505 status |= DEC_Invalid_operation;
1506 else
1507 decNumberCopy (res, rhs); /* +Infinity */
1509 else
1510 decNaNs (res, rhs, NULL, &status); /* a NaN */
1511 break;
1514 /* calculate the ideal (preferred) exponent [floor(exp/2)] */
1515 /* [We would like to write: ideal=rhs->exponent>>1, but this */
1516 /* generates a compiler warning. Generated code is the same.] */
1517 ideal = (rhs->exponent & ~1) / 2; /* target */
1519 /* handle zeros */
1520 if (ISZERO (rhs))
1522 decNumberCopy (res, rhs); /* could be 0 or -0 */
1523 res->exponent = ideal; /* use the ideal [safe] */
1524 break;
1527 /* any other -x is an oops */
1528 if (decNumberIsNegative (rhs))
1530 status |= DEC_Invalid_operation;
1531 break;
1534 /* we need space for three working variables */
1535 /* f -- the same precision as the RHS, reduced to 0.01->0.99... */
1536 /* a -- Hull's approx -- precision, when assigned, is */
1537 /* currentprecision (we allow +2 for use as temporary) */
1538 /* b -- intermediate temporary result */
1539 /* if any is too long for local storage, then allocate */
1540 needbytes =
1541 sizeof (decNumber) + (D2U (rhs->digits) - 1) * sizeof (Unit);
1542 if (needbytes > sizeof (buff))
1544 allocbuff = (decNumber *) malloc (needbytes);
1545 if (allocbuff == NULL)
1546 { /* hopeless -- abandon */
1547 status |= DEC_Insufficient_storage;
1548 break;
1550 f = allocbuff; /* use the allocated space */
1552 /* a and b both need to be able to hold a maxp-length number */
1553 needbytes = sizeof (decNumber) + (D2U (maxp) - 1) * sizeof (Unit);
1554 if (needbytes > sizeof (bufa))
1555 { /* [same applies to b] */
1556 allocbufa = (decNumber *) malloc (needbytes);
1557 allocbufb = (decNumber *) malloc (needbytes);
1558 if (allocbufa == NULL || allocbufb == NULL)
1559 { /* hopeless */
1560 status |= DEC_Insufficient_storage;
1561 break;
1563 a = allocbufa; /* use the allocated space */
1564 b = allocbufb; /* .. */
1567 /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
1568 decNumberCopy (f, rhs);
1569 exp = f->exponent + f->digits; /* adjusted to Hull rules */
1570 f->exponent = -(f->digits); /* to range */
1572 /* set up working contexts (the second is used for Numerical */
1573 /* Turing assignment) */
1574 decContextDefault (&workset, DEC_INIT_DECIMAL64);
1575 decContextDefault (&approxset, DEC_INIT_DECIMAL64);
1576 approxset.digits = set->digits; /* approx's length */
1578 /* [Until further notice, no error is possible and status bits */
1579 /* (Rounded, etc.) should be ignored, not accumulated.] */
1581 /* Calculate initial approximation, and allow for odd exponent */
1582 workset.digits = set->digits; /* p for initial calculation */
1583 t->bits = 0;
1584 t->digits = 3;
1585 a->bits = 0;
1586 a->digits = 3;
1587 if ((exp & 1) == 0)
1588 { /* even exponent */
1589 /* Set t=0.259, a=0.819 */
1590 t->exponent = -3;
1591 a->exponent = -3;
1592 #if DECDPUN>=3
1593 t->lsu[0] = 259;
1594 a->lsu[0] = 819;
1595 #elif DECDPUN==2
1596 t->lsu[0] = 59;
1597 t->lsu[1] = 2;
1598 a->lsu[0] = 19;
1599 a->lsu[1] = 8;
1600 #else
1601 t->lsu[0] = 9;
1602 t->lsu[1] = 5;
1603 t->lsu[2] = 2;
1604 a->lsu[0] = 9;
1605 a->lsu[1] = 1;
1606 a->lsu[2] = 8;
1607 #endif
1609 else
1610 { /* odd exponent */
1611 /* Set t=0.0819, a=2.59 */
1612 f->exponent--; /* f=f/10 */
1613 exp++; /* e=e+1 */
1614 t->exponent = -4;
1615 a->exponent = -2;
1616 #if DECDPUN>=3
1617 t->lsu[0] = 819;
1618 a->lsu[0] = 259;
1619 #elif DECDPUN==2
1620 t->lsu[0] = 19;
1621 t->lsu[1] = 8;
1622 a->lsu[0] = 59;
1623 a->lsu[1] = 2;
1624 #else
1625 t->lsu[0] = 9;
1626 t->lsu[1] = 1;
1627 t->lsu[2] = 8;
1628 a->lsu[0] = 9;
1629 a->lsu[1] = 5;
1630 a->lsu[2] = 2;
1631 #endif
1633 decMultiplyOp (a, a, f, &workset, &ignore); /* a=a*f */
1634 decAddOp (a, a, t, &workset, 0, &ignore); /* ..+t */
1635 /* [a is now the initial approximation for sqrt(f), calculated with */
1636 /* currentprecision, which is also a's precision.] */
1638 /* the main calculation loop */
1639 decNumberZero (&dzero); /* make 0 */
1640 decNumberZero (t); /* set t = 0.5 */
1641 t->lsu[0] = 5; /* .. */
1642 t->exponent = -1; /* .. */
1643 workset.digits = 3; /* initial p */
1644 for (;;)
1646 /* set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp] */
1647 workset.digits = workset.digits * 2 - 2;
1648 if (workset.digits > maxp)
1649 workset.digits = maxp;
1650 /* a = 0.5 * (a + f/a) */
1651 /* [calculated at p then rounded to currentprecision] */
1652 decDivideOp (b, f, a, &workset, DIVIDE, &ignore); /* b=f/a */
1653 decAddOp (b, b, a, &workset, 0, &ignore); /* b=b+a */
1654 decMultiplyOp (a, b, t, &workset, &ignore); /* a=b*0.5 */
1655 /* assign to approx [round to length] */
1656 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1657 if (workset.digits == maxp)
1658 break; /* just did final */
1659 } /* loop */
1661 /* a is now at currentprecision and within 1 ulp of the properly */
1662 /* rounded square root of f; to ensure proper rounding, compare */
1663 /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
1664 /* Here workset.digits=maxp and t=0.5 */
1665 workset.digits--; /* maxp-1 is OK now */
1666 t->exponent = -set->digits - 1; /* make 0.5 ulp */
1667 decNumberCopy (b, a);
1668 decAddOp (b, b, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */
1669 workset.round = DEC_ROUND_UP;
1670 decMultiplyOp (b, b, b, &workset, &ignore); /* b = mulru(b, b) */
1671 decCompareOp (b, f, b, &workset, COMPARE, &ignore); /* b ? f, reversed */
1672 if (decNumberIsNegative (b))
1673 { /* f < b [i.e., b > f] */
1674 /* this is the more common adjustment, though both are rare */
1675 t->exponent++; /* make 1.0 ulp */
1676 t->lsu[0] = 1; /* .. */
1677 decAddOp (a, a, t, &workset, DECNEG, &ignore); /* a = a - 1 ulp */
1678 /* assign to approx [round to length] */
1679 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1681 else
1683 decNumberCopy (b, a);
1684 decAddOp (b, b, t, &workset, 0, &ignore); /* b = a + 0.5 ulp */
1685 workset.round = DEC_ROUND_DOWN;
1686 decMultiplyOp (b, b, b, &workset, &ignore); /* b = mulrd(b, b) */
1687 decCompareOp (b, b, f, &workset, COMPARE, &ignore); /* b ? f */
1688 if (decNumberIsNegative (b))
1689 { /* b < f */
1690 t->exponent++; /* make 1.0 ulp */
1691 t->lsu[0] = 1; /* .. */
1692 decAddOp (a, a, t, &workset, 0, &ignore); /* a = a + 1 ulp */
1693 /* assign to approx [round to length] */
1694 decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1697 /* [no errors are possible in the above, and rounding/inexact during */
1698 /* estimation are irrelevant, so status was not accumulated] */
1700 /* Here, 0.1 <= a < 1 [Hull] */
1701 a->exponent += exp / 2; /* set correct exponent */
1703 /* Process Subnormals */
1704 decFinalize (a, set, &residue, &status);
1706 /* count dropable zeros [after any subnormal rounding] */
1707 decNumberCopy (b, a);
1708 decTrim (b, 1, &dropped); /* [drops trailing zeros] */
1710 /* Finally set Inexact and Rounded. The answer can only be exact if */
1711 /* it is short enough so that squaring it could fit in set->digits, */
1712 /* so this is the only (relatively rare) time we have to check */
1713 /* carefully */
1714 if (b->digits * 2 - 1 > set->digits)
1715 { /* cannot fit */
1716 status |= DEC_Inexact | DEC_Rounded;
1718 else
1719 { /* could be exact/unrounded */
1720 uInt mstatus = 0; /* local status */
1721 decMultiplyOp (b, b, b, &workset, &mstatus); /* try the multiply */
1722 if (mstatus != 0)
1723 { /* result won't fit */
1724 status |= DEC_Inexact | DEC_Rounded;
1726 else
1727 { /* plausible */
1728 decCompareOp (t, b, rhs, &workset, COMPARE, &mstatus); /* b ? rhs */
1729 if (!ISZERO (t))
1731 status |= DEC_Inexact | DEC_Rounded;
1733 else
1734 { /* is Exact */
1735 /* here, dropped is the count of trailing zeros in 'a' */
1736 /* use closest exponent to ideal... */
1737 Int todrop = ideal - a->exponent; /* most we can drop */
1739 if (todrop < 0)
1740 { /* ideally would add 0s */
1741 status |= DEC_Rounded;
1743 else
1744 { /* unrounded */
1745 if (dropped < todrop)
1746 todrop = dropped; /* clamp to those available */
1747 if (todrop > 0)
1748 { /* OK, some to drop */
1749 decShiftToLeast (a->lsu, D2U (a->digits), todrop);
1750 a->exponent += todrop; /* maintain numerical value */
1751 a->digits -= todrop; /* new length */
1757 decNumberCopy (res, a); /* assume this is the result */
1759 while (0); /* end protected */
1761 if (allocbuff != NULL)
1762 free (allocbuff); /* drop any storage we used */
1763 if (allocbufa != NULL)
1764 free (allocbufa); /* .. */
1765 if (allocbufb != NULL)
1766 free (allocbufb); /* .. */
1767 if (allocrhs != NULL)
1768 free (allocrhs); /* .. */
1769 if (status != 0)
1770 decStatus (res, status, set); /* then report status */
1771 return res;
1774 /* ------------------------------------------------------------------ */
1775 /* decNumberSubtract -- subtract two Numbers */
1776 /* */
1777 /* This computes C = A - B */
1778 /* */
1779 /* res is C, the result. C may be A and/or B (e.g., X=X-X) */
1780 /* lhs is A */
1781 /* rhs is B */
1782 /* set is the context */
1783 /* */
1784 /* C must have space for set->digits digits. */
1785 /* ------------------------------------------------------------------ */
1786 decNumber *
1787 decNumberSubtract (decNumber * res, const decNumber * lhs,
1788 const decNumber * rhs, decContext * set)
1790 uInt status = 0; /* accumulator */
1792 decAddOp (res, lhs, rhs, set, DECNEG, &status);
1793 if (status != 0)
1794 decStatus (res, status, set);
1795 return res;
1798 /* ------------------------------------------------------------------ */
1799 /* decNumberToIntegralValue -- round-to-integral-value */
1800 /* */
1801 /* res is the result */
1802 /* rhs is input number */
1803 /* set is the context */
1804 /* */
1805 /* res must have space for any value of rhs. */
1806 /* */
1807 /* This implements the IEEE special operator and therefore treats */
1808 /* special values as valid, and also never sets Inexact. For finite */
1809 /* numbers it returns rescale(rhs, 0) if rhs->exponent is <0. */
1810 /* Otherwise the result is rhs (so no error is possible). */
1811 /* */
1812 /* The context is used for rounding mode and status after sNaN, but */
1813 /* the digits setting is ignored. */
1814 /* ------------------------------------------------------------------ */
1815 decNumber *
1816 decNumberToIntegralValue (decNumber * res, const decNumber * rhs, decContext * set)
1818 decNumber dn;
1819 decContext workset; /* working context */
1821 #if DECCHECK
1822 if (decCheckOperands (res, DECUNUSED, rhs, set))
1823 return res;
1824 #endif
1826 /* handle infinities and NaNs */
1827 if (rhs->bits & DECSPECIAL)
1829 uInt status = 0;
1830 if (decNumberIsInfinite (rhs))
1831 decNumberCopy (res, rhs); /* an Infinity */
1832 else
1833 decNaNs (res, rhs, NULL, &status); /* a NaN */
1834 if (status != 0)
1835 decStatus (res, status, set);
1836 return res;
1839 /* we have a finite number; no error possible */
1840 if (rhs->exponent >= 0)
1841 return decNumberCopy (res, rhs);
1842 /* that was easy, but if negative exponent we have work to do... */
1843 workset = *set; /* clone rounding, etc. */
1844 workset.digits = rhs->digits; /* no length rounding */
1845 workset.traps = 0; /* no traps */
1846 decNumberZero (&dn); /* make a number with exponent 0 */
1847 return decNumberQuantize (res, rhs, &dn, &workset);
1850 /* ================================================================== */
1851 /* Utility routines */
1852 /* ================================================================== */
1854 /* ------------------------------------------------------------------ */
1855 /* decNumberCopy -- copy a number */
1856 /* */
1857 /* dest is the target decNumber */
1858 /* src is the source decNumber */
1859 /* returns dest */
1860 /* */
1861 /* (dest==src is allowed and is a no-op) */
1862 /* All fields are updated as required. This is a utility operation, */
1863 /* so special values are unchanged and no error is possible. */
1864 /* ------------------------------------------------------------------ */
1865 decNumber *
1866 decNumberCopy (decNumber * dest, const decNumber * src)
1869 #if DECCHECK
1870 if (src == NULL)
1871 return decNumberZero (dest);
1872 #endif
1874 if (dest == src)
1875 return dest; /* no copy required */
1877 /* We use explicit assignments here as structure assignment can copy */
1878 /* more than just the lsu (for small DECDPUN). This would not affect */
1879 /* the value of the results, but would disturb test harness spill */
1880 /* checking. */
1881 dest->bits = src->bits;
1882 dest->exponent = src->exponent;
1883 dest->digits = src->digits;
1884 dest->lsu[0] = src->lsu[0];
1885 if (src->digits > DECDPUN)
1886 { /* more Units to come */
1887 Unit *d; /* work */
1888 const Unit *s, *smsup; /* work */
1889 /* memcpy for the remaining Units would be safe as they cannot */
1890 /* overlap. However, this explicit loop is faster in short cases. */
1891 d = dest->lsu + 1; /* -> first destination */
1892 smsup = src->lsu + D2U (src->digits); /* -> source msu+1 */
1893 for (s = src->lsu + 1; s < smsup; s++, d++)
1894 *d = *s;
1896 return dest;
1899 /* ------------------------------------------------------------------ */
1900 /* decNumberTrim -- remove insignificant zeros */
1901 /* */
1902 /* dn is the number to trim */
1903 /* returns dn */
1904 /* */
1905 /* All fields are updated as required. This is a utility operation, */
1906 /* so special values are unchanged and no error is possible. */
1907 /* ------------------------------------------------------------------ */
1908 decNumber *
1909 decNumberTrim (decNumber * dn)
1911 Int dropped; /* work */
1912 return decTrim (dn, 0, &dropped);
1915 /* ------------------------------------------------------------------ */
1916 /* decNumberVersion -- return the name and version of this module */
1917 /* */
1918 /* No error is possible. */
1919 /* ------------------------------------------------------------------ */
1920 const char *
1921 decNumberVersion (void)
1923 return DECVERSION;
1926 /* ------------------------------------------------------------------ */
1927 /* decNumberZero -- set a number to 0 */
1928 /* */
1929 /* dn is the number to set, with space for one digit */
1930 /* returns dn */
1931 /* */
1932 /* No error is possible. */
1933 /* ------------------------------------------------------------------ */
1934 /* Memset is not used as it is much slower in some environments. */
1935 decNumber *
1936 decNumberZero (decNumber * dn)
1939 #if DECCHECK
1940 if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
1941 return dn;
1942 #endif
1944 dn->bits = 0;
1945 dn->exponent = 0;
1946 dn->digits = 1;
1947 dn->lsu[0] = 0;
1948 return dn;
1951 /* ================================================================== */
1952 /* Local routines */
1953 /* ================================================================== */
1955 /* ------------------------------------------------------------------ */
1956 /* decToString -- lay out a number into a string */
1957 /* */
1958 /* dn is the number to lay out */
1959 /* string is where to lay out the number */
1960 /* eng is 1 if Engineering, 0 if Scientific */
1961 /* */
1962 /* str must be at least dn->digits+14 characters long */
1963 /* No error is possible. */
1964 /* */
1965 /* Note that this routine can generate a -0 or 0.000. These are */
1966 /* never generated in subset to-number or arithmetic, but can occur */
1967 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234). */
1968 /* ------------------------------------------------------------------ */
1969 /* If DECCHECK is enabled the string "?" is returned if a number is */
1970 /* invalid. */
1972 /* TODIGIT -- macro to remove the leading digit from the unsigned */
1973 /* integer u at column cut (counting from the right, LSD=0) and place */
1974 /* it as an ASCII character into the character pointed to by c. Note */
1975 /* that cut must be <= 9, and the maximum value for u is 2,000,000,000 */
1976 /* (as is needed for negative exponents of subnormals). The unsigned */
1977 /* integer pow is used as a temporary variable. */
1978 #define TODIGIT(u, cut, c) { \
1979 *(c)='0'; \
1980 pow=powers[cut]*2; \
1981 if ((u)>pow) { \
1982 pow*=4; \
1983 if ((u)>=pow) {(u)-=pow; *(c)+=8;} \
1984 pow/=2; \
1985 if ((u)>=pow) {(u)-=pow; *(c)+=4;} \
1986 pow/=2; \
1988 if ((u)>=pow) {(u)-=pow; *(c)+=2;} \
1989 pow/=2; \
1990 if ((u)>=pow) {(u)-=pow; *(c)+=1;} \
1993 static void
1994 decToString (const decNumber * dn, char *string, Flag eng)
1996 Int exp = dn->exponent; /* local copy */
1997 Int e; /* E-part value */
1998 Int pre; /* digits before the '.' */
1999 Int cut; /* for counting digits in a Unit */
2000 char *c = string; /* work [output pointer] */
2001 const Unit *up = dn->lsu + D2U (dn->digits) - 1; /* -> msu [input pointer] */
2002 uInt u, pow; /* work */
2004 #if DECCHECK
2005 if (decCheckOperands (DECUNUSED, dn, DECUNUSED, DECUNUSED))
2007 strcpy (string, "?");
2008 return;
2010 #endif
2012 if (decNumberIsNegative (dn))
2013 { /* Negatives get a minus (except */
2014 *c = '-'; /* NaNs, which remove the '-' below) */
2015 c++;
2017 if (dn->bits & DECSPECIAL)
2018 { /* Is a special value */
2019 if (decNumberIsInfinite (dn))
2021 strcpy (c, "Infinity");
2022 return;
2024 /* a NaN */
2025 if (dn->bits & DECSNAN)
2026 { /* signalling NaN */
2027 *c = 's';
2028 c++;
2030 strcpy (c, "NaN");
2031 c += 3; /* step past */
2032 /* if not a clean non-zero coefficient, that's all we have in a */
2033 /* NaN string */
2034 if (exp != 0 || (*dn->lsu == 0 && dn->digits == 1))
2035 return;
2036 /* [drop through to add integer] */
2039 /* calculate how many digits in msu, and hence first cut */
2040 cut = dn->digits % DECDPUN;
2041 if (cut == 0)
2042 cut = DECDPUN; /* msu is full */
2043 cut--; /* power of ten for digit */
2045 if (exp == 0)
2046 { /* simple integer [common fastpath, */
2047 /* used for NaNs, too] */
2048 for (; up >= dn->lsu; up--)
2049 { /* each Unit from msu */
2050 u = *up; /* contains DECDPUN digits to lay out */
2051 for (; cut >= 0; c++, cut--)
2052 TODIGIT (u, cut, c);
2053 cut = DECDPUN - 1; /* next Unit has all digits */
2055 *c = '\0'; /* terminate the string */
2056 return;
2059 /* non-0 exponent -- assume plain form */
2060 pre = dn->digits + exp; /* digits before '.' */
2061 e = 0; /* no E */
2062 if ((exp > 0) || (pre < -5))
2063 { /* need exponential form */
2064 e = exp + dn->digits - 1; /* calculate E value */
2065 pre = 1; /* assume one digit before '.' */
2066 if (eng && (e != 0))
2067 { /* may need to adjust */
2068 Int adj; /* adjustment */
2069 /* The C remainder operator is undefined for negative numbers, so */
2070 /* we must use positive remainder calculation here */
2071 if (e < 0)
2073 adj = (-e) % 3;
2074 if (adj != 0)
2075 adj = 3 - adj;
2077 else
2078 { /* e>0 */
2079 adj = e % 3;
2081 e = e - adj;
2082 /* if we are dealing with zero we will use exponent which is a */
2083 /* multiple of three, as expected, but there will only be the */
2084 /* one zero before the E, still. Otherwise note the padding. */
2085 if (!ISZERO (dn))
2086 pre += adj;
2087 else
2088 { /* is zero */
2089 if (adj != 0)
2090 { /* 0.00Esnn needed */
2091 e = e + 3;
2092 pre = -(2 - adj);
2094 } /* zero */
2095 } /* eng */
2098 /* lay out the digits of the coefficient, adding 0s and . as needed */
2099 u = *up;
2100 if (pre > 0)
2101 { /* xxx.xxx or xx00 (engineering) form */
2102 for (; pre > 0; pre--, c++, cut--)
2104 if (cut < 0)
2105 { /* need new Unit */
2106 if (up == dn->lsu)
2107 break; /* out of input digits (pre>digits) */
2108 up--;
2109 cut = DECDPUN - 1;
2110 u = *up;
2112 TODIGIT (u, cut, c);
2114 if (up > dn->lsu || (up == dn->lsu && cut >= 0))
2115 { /* more to come, after '.' */
2116 *c = '.';
2117 c++;
2118 for (;; c++, cut--)
2120 if (cut < 0)
2121 { /* need new Unit */
2122 if (up == dn->lsu)
2123 break; /* out of input digits */
2124 up--;
2125 cut = DECDPUN - 1;
2126 u = *up;
2128 TODIGIT (u, cut, c);
2131 else
2132 for (; pre > 0; pre--, c++)
2133 *c = '0'; /* 0 padding (for engineering) needed */
2135 else
2136 { /* 0.xxx or 0.000xxx form */
2137 *c = '0';
2138 c++;
2139 *c = '.';
2140 c++;
2141 for (; pre < 0; pre++, c++)
2142 *c = '0'; /* add any 0's after '.' */
2143 for (;; c++, cut--)
2145 if (cut < 0)
2146 { /* need new Unit */
2147 if (up == dn->lsu)
2148 break; /* out of input digits */
2149 up--;
2150 cut = DECDPUN - 1;
2151 u = *up;
2153 TODIGIT (u, cut, c);
2157 /* Finally add the E-part, if needed. It will never be 0, has a
2158 base maximum and minimum of +999999999 through -999999999, but
2159 could range down to -1999999998 for subnormal numbers */
2160 if (e != 0)
2162 Flag had = 0; /* 1=had non-zero */
2163 *c = 'E';
2164 c++;
2165 *c = '+';
2166 c++; /* assume positive */
2167 u = e; /* .. */
2168 if (e < 0)
2170 *(c - 1) = '-'; /* oops, need - */
2171 u = -e; /* uInt, please */
2173 /* layout the exponent (_itoa is not ANSI C) */
2174 for (cut = 9; cut >= 0; cut--)
2176 TODIGIT (u, cut, c);
2177 if (*c == '0' && !had)
2178 continue; /* skip leading zeros */
2179 had = 1; /* had non-0 */
2180 c++; /* step for next */
2181 } /* cut */
2183 *c = '\0'; /* terminate the string (all paths) */
2184 return;
2187 /* ------------------------------------------------------------------ */
2188 /* decAddOp -- add/subtract operation */
2189 /* */
2190 /* This computes C = A + B */
2191 /* */
2192 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
2193 /* lhs is A */
2194 /* rhs is B */
2195 /* set is the context */
2196 /* negate is DECNEG if rhs should be negated, or 0 otherwise */
2197 /* status accumulates status for the caller */
2198 /* */
2199 /* C must have space for set->digits digits. */
2200 /* ------------------------------------------------------------------ */
2201 /* If possible, we calculate the coefficient directly into C. */
2202 /* However, if: */
2203 /* -- we need a digits+1 calculation because numbers are unaligned */
2204 /* and span more than set->digits digits */
2205 /* -- a carry to digits+1 digits looks possible */
2206 /* -- C is the same as A or B, and the result would destructively */
2207 /* overlap the A or B coefficient */
2208 /* then we must calculate into a temporary buffer. In this latter */
2209 /* case we use the local (stack) buffer if possible, and only if too */
2210 /* long for that do we resort to malloc. */
2211 /* */
2212 /* Misalignment is handled as follows: */
2213 /* Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp. */
2214 /* BPad: Apply the padding by a combination of shifting (whole */
2215 /* units) and multiplication (part units). */
2216 /* */
2217 /* Addition, especially x=x+1, is speed-critical, so we take pains */
2218 /* to make returning as fast as possible, by flagging any allocation. */
2219 /* ------------------------------------------------------------------ */
2220 static decNumber *
2221 decAddOp (decNumber * res, const decNumber * lhs,
2222 const decNumber * rhs, decContext * set, uByte negate, uInt * status)
2224 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
2225 decNumber *allocrhs = NULL; /* .., rhs */
2226 Int rhsshift; /* working shift (in Units) */
2227 Int maxdigits; /* longest logical length */
2228 Int mult; /* multiplier */
2229 Int residue; /* rounding accumulator */
2230 uByte bits; /* result bits */
2231 Flag diffsign; /* non-0 if arguments have different sign */
2232 Unit *acc; /* accumulator for result */
2233 Unit accbuff[D2U (DECBUFFER + 1)]; /* local buffer [+1 is for possible */
2234 /* final carry digit or DECBUFFER=0] */
2235 Unit *allocacc = NULL; /* -> allocated acc buffer, iff allocated */
2236 Flag alloced = 0; /* set non-0 if any allocations */
2237 Int reqdigits = set->digits; /* local copy; requested DIGITS */
2238 uByte merged; /* merged flags */
2239 Int padding; /* work */
2241 #if DECCHECK
2242 if (decCheckOperands (res, lhs, rhs, set))
2243 return res;
2244 #endif
2247 { /* protect allocated storage */
2248 #if DECSUBSET
2249 if (!set->extended)
2251 /* reduce operands and set lostDigits status, as needed */
2252 if (lhs->digits > reqdigits)
2254 alloclhs = decRoundOperand (lhs, set, status);
2255 if (alloclhs == NULL)
2256 break;
2257 lhs = alloclhs;
2258 alloced = 1;
2260 if (rhs->digits > reqdigits)
2262 allocrhs = decRoundOperand (rhs, set, status);
2263 if (allocrhs == NULL)
2264 break;
2265 rhs = allocrhs;
2266 alloced = 1;
2269 #endif
2270 /* [following code does not require input rounding] */
2272 /* note whether signs differ */
2273 diffsign = (Flag) ((lhs->bits ^ rhs->bits ^ negate) & DECNEG);
2275 /* handle infinities and NaNs */
2276 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2277 if (merged)
2278 { /* a special bit set */
2279 if (merged & (DECSNAN | DECNAN)) /* a NaN */
2280 decNaNs (res, lhs, rhs, status);
2281 else
2282 { /* one or two infinities */
2283 if (decNumberIsInfinite (lhs))
2284 { /* LHS is infinity */
2285 /* two infinities with different signs is invalid */
2286 if (decNumberIsInfinite (rhs) && diffsign)
2288 *status |= DEC_Invalid_operation;
2289 break;
2291 bits = lhs->bits & DECNEG; /* get sign from LHS */
2293 else
2294 bits = (rhs->bits ^ negate) & DECNEG; /* RHS must be Infinity */
2295 bits |= DECINF;
2296 decNumberZero (res);
2297 res->bits = bits; /* set +/- infinity */
2298 } /* an infinity */
2299 break;
2302 /* Quick exit for add 0s; return the non-0, modified as need be */
2303 if (ISZERO (lhs))
2305 Int adjust; /* work */
2306 Int lexp = lhs->exponent; /* save in case LHS==RES */
2307 bits = lhs->bits; /* .. */
2308 residue = 0; /* clear accumulator */
2309 decCopyFit (res, rhs, set, &residue, status); /* copy (as needed) */
2310 res->bits ^= negate; /* flip if rhs was negated */
2311 #if DECSUBSET
2312 if (set->extended)
2313 { /* exponents on zeros count */
2314 #endif
2315 /* exponent will be the lower of the two */
2316 adjust = lexp - res->exponent; /* adjustment needed [if -ve] */
2317 if (ISZERO (res))
2318 { /* both 0: special IEEE 854 rules */
2319 if (adjust < 0)
2320 res->exponent = lexp; /* set exponent */
2321 /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
2322 if (diffsign)
2324 if (set->round != DEC_ROUND_FLOOR)
2325 res->bits = 0;
2326 else
2327 res->bits = DECNEG; /* preserve 0 sign */
2330 else
2331 { /* non-0 res */
2332 if (adjust < 0)
2333 { /* 0-padding needed */
2334 if ((res->digits - adjust) > set->digits)
2336 adjust = res->digits - set->digits; /* to fit exactly */
2337 *status |= DEC_Rounded; /* [but exact] */
2339 res->digits =
2340 decShiftToMost (res->lsu, res->digits, -adjust);
2341 res->exponent += adjust; /* set the exponent. */
2343 } /* non-0 res */
2344 #if DECSUBSET
2345 } /* extended */
2346 #endif
2347 decFinish (res, set, &residue, status); /* clean and finalize */
2348 break;
2351 if (ISZERO (rhs))
2352 { /* [lhs is non-zero] */
2353 Int adjust; /* work */
2354 Int rexp = rhs->exponent; /* save in case RHS==RES */
2355 bits = rhs->bits; /* be clean */
2356 residue = 0; /* clear accumulator */
2357 decCopyFit (res, lhs, set, &residue, status); /* copy (as needed) */
2358 #if DECSUBSET
2359 if (set->extended)
2360 { /* exponents on zeros count */
2361 #endif
2362 /* exponent will be the lower of the two */
2363 /* [0-0 case handled above] */
2364 adjust = rexp - res->exponent; /* adjustment needed [if -ve] */
2365 if (adjust < 0)
2366 { /* 0-padding needed */
2367 if ((res->digits - adjust) > set->digits)
2369 adjust = res->digits - set->digits; /* to fit exactly */
2370 *status |= DEC_Rounded; /* [but exact] */
2372 res->digits =
2373 decShiftToMost (res->lsu, res->digits, -adjust);
2374 res->exponent += adjust; /* set the exponent. */
2376 #if DECSUBSET
2377 } /* extended */
2378 #endif
2379 decFinish (res, set, &residue, status); /* clean and finalize */
2380 break;
2382 /* [both fastpath and mainpath code below assume these cases */
2383 /* (notably 0-0) have already been handled] */
2385 /* calculate the padding needed to align the operands */
2386 padding = rhs->exponent - lhs->exponent;
2388 /* Fastpath cases where the numbers are aligned and normal, the RHS */
2389 /* is all in one unit, no operand rounding is needed, and no carry, */
2390 /* lengthening, or borrow is needed */
2391 if (rhs->digits <= DECDPUN && padding == 0 && rhs->exponent >= set->emin /* [some normals drop through] */
2392 && rhs->digits <= reqdigits && lhs->digits <= reqdigits)
2394 Int partial = *lhs->lsu;
2395 if (!diffsign)
2396 { /* adding */
2397 Int maxv = DECDPUNMAX; /* highest no-overflow */
2398 if (lhs->digits < DECDPUN)
2399 maxv = powers[lhs->digits] - 1;
2400 partial += *rhs->lsu;
2401 if (partial <= maxv)
2402 { /* no carry */
2403 if (res != lhs)
2404 decNumberCopy (res, lhs); /* not in place */
2405 *res->lsu = (Unit) partial; /* [copy could have overwritten RHS] */
2406 break;
2408 /* else drop out for careful add */
2410 else
2411 { /* signs differ */
2412 partial -= *rhs->lsu;
2413 if (partial > 0)
2414 { /* no borrow needed, and non-0 result */
2415 if (res != lhs)
2416 decNumberCopy (res, lhs); /* not in place */
2417 *res->lsu = (Unit) partial;
2418 /* this could have reduced digits [but result>0] */
2419 res->digits = decGetDigits (res->lsu, D2U (res->digits));
2420 break;
2422 /* else drop out for careful subtract */
2426 /* Now align (pad) the lhs or rhs so we can add or subtract them, as
2427 necessary. If one number is much larger than the other (that is,
2428 if in plain form there is a least one digit between the lowest
2429 digit or one and the highest of the other) we need to pad with up
2430 to DIGITS-1 trailing zeros, and then apply rounding (as exotic
2431 rounding modes may be affected by the residue).
2433 rhsshift = 0; /* rhs shift to left (padding) in Units */
2434 bits = lhs->bits; /* assume sign is that of LHS */
2435 mult = 1; /* likely multiplier */
2437 /* if padding==0 the operands are aligned; no padding needed */
2438 if (padding != 0)
2440 /* some padding needed */
2441 /* We always pad the RHS, as we can then effect any required */
2442 /* padding by a combination of shifts and a multiply */
2443 Flag swapped = 0;
2444 if (padding < 0)
2445 { /* LHS needs the padding */
2446 const decNumber *t;
2447 padding = -padding; /* will be +ve */
2448 bits = (uByte) (rhs->bits ^ negate); /* assumed sign is now that of RHS */
2449 t = lhs;
2450 lhs = rhs;
2451 rhs = t;
2452 swapped = 1;
2455 /* If, after pad, rhs would be longer than lhs by digits+1 or */
2456 /* more then lhs cannot affect the answer, except as a residue, */
2457 /* so we only need to pad up to a length of DIGITS+1. */
2458 if (rhs->digits + padding > lhs->digits + reqdigits + 1)
2460 /* The RHS is sufficient */
2461 /* for residue we use the relative sign indication... */
2462 Int shift = reqdigits - rhs->digits; /* left shift needed */
2463 residue = 1; /* residue for rounding */
2464 if (diffsign)
2465 residue = -residue; /* signs differ */
2466 /* copy, shortening if necessary */
2467 decCopyFit (res, rhs, set, &residue, status);
2468 /* if it was already shorter, then need to pad with zeros */
2469 if (shift > 0)
2471 res->digits = decShiftToMost (res->lsu, res->digits, shift);
2472 res->exponent -= shift; /* adjust the exponent. */
2474 /* flip the result sign if unswapped and rhs was negated */
2475 if (!swapped)
2476 res->bits ^= negate;
2477 decFinish (res, set, &residue, status); /* done */
2478 break;
2481 /* LHS digits may affect result */
2482 rhsshift = D2U (padding + 1) - 1; /* this much by Unit shift .. */
2483 mult = powers[padding - (rhsshift * DECDPUN)]; /* .. this by multiplication */
2484 } /* padding needed */
2486 if (diffsign)
2487 mult = -mult; /* signs differ */
2489 /* determine the longer operand */
2490 maxdigits = rhs->digits + padding; /* virtual length of RHS */
2491 if (lhs->digits > maxdigits)
2492 maxdigits = lhs->digits;
2494 /* Decide on the result buffer to use; if possible place directly */
2495 /* into result. */
2496 acc = res->lsu; /* assume build direct */
2497 /* If destructive overlap, or the number is too long, or a carry or */
2498 /* borrow to DIGITS+1 might be possible we must use a buffer. */
2499 /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
2500 if ((maxdigits >= reqdigits) /* is, or could be, too large */
2501 || (res == rhs && rhsshift > 0))
2502 { /* destructive overlap */
2503 /* buffer needed; choose it */
2504 /* we'll need units for maxdigits digits, +1 Unit for carry or borrow */
2505 Int need = D2U (maxdigits) + 1;
2506 acc = accbuff; /* assume use local buffer */
2507 if (need * sizeof (Unit) > sizeof (accbuff))
2509 allocacc = (Unit *) malloc (need * sizeof (Unit));
2510 if (allocacc == NULL)
2511 { /* hopeless -- abandon */
2512 *status |= DEC_Insufficient_storage;
2513 break;
2515 acc = allocacc;
2516 alloced = 1;
2520 res->bits = (uByte) (bits & DECNEG); /* it's now safe to overwrite.. */
2521 res->exponent = lhs->exponent; /* .. operands (even if aliased) */
2523 #if DECTRACE
2524 decDumpAr ('A', lhs->lsu, D2U (lhs->digits));
2525 decDumpAr ('B', rhs->lsu, D2U (rhs->digits));
2526 printf (" :h: %d %d\n", rhsshift, mult);
2527 #endif
2529 /* add [A+B*m] or subtract [A+B*(-m)] */
2530 res->digits = decUnitAddSub (lhs->lsu, D2U (lhs->digits), rhs->lsu, D2U (rhs->digits), rhsshift, acc, mult) * DECDPUN; /* [units -> digits] */
2531 if (res->digits < 0)
2532 { /* we borrowed */
2533 res->digits = -res->digits;
2534 res->bits ^= DECNEG; /* flip the sign */
2536 #if DECTRACE
2537 decDumpAr ('+', acc, D2U (res->digits));
2538 #endif
2540 /* If we used a buffer we need to copy back, possibly shortening */
2541 /* (If we didn't use buffer it must have fit, so can't need rounding */
2542 /* and residue must be 0.) */
2543 residue = 0; /* clear accumulator */
2544 if (acc != res->lsu)
2546 #if DECSUBSET
2547 if (set->extended)
2548 { /* round from first significant digit */
2549 #endif
2550 /* remove leading zeros that we added due to rounding up to */
2551 /* integral Units -- before the test for rounding. */
2552 if (res->digits > reqdigits)
2553 res->digits = decGetDigits (acc, D2U (res->digits));
2554 decSetCoeff (res, set, acc, res->digits, &residue, status);
2555 #if DECSUBSET
2557 else
2558 { /* subset arithmetic rounds from original significant digit */
2559 /* We may have an underestimate. This only occurs when both */
2560 /* numbers fit in DECDPUN digits and we are padding with a */
2561 /* negative multiple (-10, -100...) and the top digit(s) become */
2562 /* 0. (This only matters if we are using X3.274 rules where the */
2563 /* leading zero could be included in the rounding.) */
2564 if (res->digits < maxdigits)
2566 *(acc + D2U (res->digits)) = 0; /* ensure leading 0 is there */
2567 res->digits = maxdigits;
2569 else
2571 /* remove leading zeros that we added due to rounding up to */
2572 /* integral Units (but only those in excess of the original */
2573 /* maxdigits length, unless extended) before test for rounding. */
2574 if (res->digits > reqdigits)
2576 res->digits = decGetDigits (acc, D2U (res->digits));
2577 if (res->digits < maxdigits)
2578 res->digits = maxdigits;
2581 decSetCoeff (res, set, acc, res->digits, &residue, status);
2582 /* Now apply rounding if needed before removing leading zeros. */
2583 /* This is safe because subnormals are not a possibility */
2584 if (residue != 0)
2586 decApplyRound (res, set, residue, status);
2587 residue = 0; /* we did what we had to do */
2589 } /* subset */
2590 #endif
2591 } /* used buffer */
2593 /* strip leading zeros [these were left on in case of subset subtract] */
2594 res->digits = decGetDigits (res->lsu, D2U (res->digits));
2596 /* apply checks and rounding */
2597 decFinish (res, set, &residue, status);
2599 /* "When the sum of two operands with opposite signs is exactly */
2600 /* zero, the sign of that sum shall be '+' in all rounding modes */
2601 /* except round toward -Infinity, in which mode that sign shall be */
2602 /* '-'." [Subset zeros also never have '-', set by decFinish.] */
2603 if (ISZERO (res) && diffsign
2604 #if DECSUBSET
2605 && set->extended
2606 #endif
2607 && (*status & DEC_Inexact) == 0)
2609 if (set->round == DEC_ROUND_FLOOR)
2610 res->bits |= DECNEG; /* sign - */
2611 else
2612 res->bits &= ~DECNEG; /* sign + */
2615 while (0); /* end protected */
2617 if (alloced)
2619 if (allocacc != NULL)
2620 free (allocacc); /* drop any storage we used */
2621 if (allocrhs != NULL)
2622 free (allocrhs); /* .. */
2623 if (alloclhs != NULL)
2624 free (alloclhs); /* .. */
2626 return res;
2629 /* ------------------------------------------------------------------ */
2630 /* decDivideOp -- division operation */
2631 /* */
2632 /* This routine performs the calculations for all four division */
2633 /* operators (divide, divideInteger, remainder, remainderNear). */
2634 /* */
2635 /* C=A op B */
2636 /* */
2637 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
2638 /* lhs is A */
2639 /* rhs is B */
2640 /* set is the context */
2641 /* op is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively. */
2642 /* status is the usual accumulator */
2643 /* */
2644 /* C must have space for set->digits digits. */
2645 /* */
2646 /* ------------------------------------------------------------------ */
2647 /* The underlying algorithm of this routine is the same as in the */
2648 /* 1981 S/370 implementation, that is, non-restoring long division */
2649 /* with bi-unit (rather than bi-digit) estimation for each unit */
2650 /* multiplier. In this pseudocode overview, complications for the */
2651 /* Remainder operators and division residues for exact rounding are */
2652 /* omitted for clarity. */
2653 /* */
2654 /* Prepare operands and handle special values */
2655 /* Test for x/0 and then 0/x */
2656 /* Exp =Exp1 - Exp2 */
2657 /* Exp =Exp +len(var1) -len(var2) */
2658 /* Sign=Sign1 * Sign2 */
2659 /* Pad accumulator (Var1) to double-length with 0's (pad1) */
2660 /* Pad Var2 to same length as Var1 */
2661 /* msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round */
2662 /* have=0 */
2663 /* Do until (have=digits+1 OR residue=0) */
2664 /* if exp<0 then if integer divide/residue then leave */
2665 /* this_unit=0 */
2666 /* Do forever */
2667 /* compare numbers */
2668 /* if <0 then leave inner_loop */
2669 /* if =0 then (* quick exit without subtract *) do */
2670 /* this_unit=this_unit+1; output this_unit */
2671 /* leave outer_loop; end */
2672 /* Compare lengths of numbers (mantissae): */
2673 /* If same then tops2=msu2pair -- {units 1&2 of var2} */
2674 /* else tops2=msu2plus -- {0, unit 1 of var2} */
2675 /* tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
2676 /* mult=tops1/tops2 -- Good and safe guess at divisor */
2677 /* if mult=0 then mult=1 */
2678 /* this_unit=this_unit+mult */
2679 /* subtract */
2680 /* end inner_loop */
2681 /* if have\=0 | this_unit\=0 then do */
2682 /* output this_unit */
2683 /* have=have+1; end */
2684 /* var2=var2/10 */
2685 /* exp=exp-1 */
2686 /* end outer_loop */
2687 /* exp=exp+1 -- set the proper exponent */
2688 /* if have=0 then generate answer=0 */
2689 /* Return (Result is defined by Var1) */
2690 /* */
2691 /* ------------------------------------------------------------------ */
2692 /* We need two working buffers during the long division; one (digits+ */
2693 /* 1) to accumulate the result, and the other (up to 2*digits+1) for */
2694 /* long subtractions. These are acc and var1 respectively. */
2695 /* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
2696 /* ------------------------------------------------------------------ */
2697 static decNumber *
2698 decDivideOp (decNumber * res,
2699 const decNumber * lhs, const decNumber * rhs,
2700 decContext * set, Flag op, uInt * status)
2702 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
2703 decNumber *allocrhs = NULL; /* .., rhs */
2704 Unit accbuff[D2U (DECBUFFER + DECDPUN)]; /* local buffer */
2705 Unit *acc = accbuff; /* -> accumulator array for result */
2706 Unit *allocacc = NULL; /* -> allocated buffer, iff allocated */
2707 Unit *accnext; /* -> where next digit will go */
2708 Int acclength; /* length of acc needed [Units] */
2709 Int accunits; /* count of units accumulated */
2710 Int accdigits; /* count of digits accumulated */
2712 Unit varbuff[D2U (DECBUFFER * 2 + DECDPUN) * sizeof (Unit)]; /* buffer for var1 */
2713 Unit *var1 = varbuff; /* -> var1 array for long subtraction */
2714 Unit *varalloc = NULL; /* -> allocated buffer, iff used */
2716 const Unit *var2; /* -> var2 array */
2718 Int var1units, var2units; /* actual lengths */
2719 Int var2ulen; /* logical length (units) */
2720 Int var1initpad = 0; /* var1 initial padding (digits) */
2721 Unit *msu1; /* -> msu of each var */
2722 const Unit *msu2; /* -> msu of each var */
2723 Int msu2plus; /* msu2 plus one [does not vary] */
2724 eInt msu2pair; /* msu2 pair plus one [does not vary] */
2725 Int maxdigits; /* longest LHS or required acc length */
2726 Int mult; /* multiplier for subtraction */
2727 Unit thisunit; /* current unit being accumulated */
2728 Int residue; /* for rounding */
2729 Int reqdigits = set->digits; /* requested DIGITS */
2730 Int exponent; /* working exponent */
2731 Int maxexponent = 0; /* DIVIDE maximum exponent if unrounded */
2732 uByte bits; /* working sign */
2733 uByte merged; /* merged flags */
2734 Unit *target; /* work */
2735 const Unit *source; /* work */
2736 uInt const *pow; /* .. */
2737 Int shift, cut; /* .. */
2738 #if DECSUBSET
2739 Int dropped; /* work */
2740 #endif
2742 #if DECCHECK
2743 if (decCheckOperands (res, lhs, rhs, set))
2744 return res;
2745 #endif
2748 { /* protect allocated storage */
2749 #if DECSUBSET
2750 if (!set->extended)
2752 /* reduce operands and set lostDigits status, as needed */
2753 if (lhs->digits > reqdigits)
2755 alloclhs = decRoundOperand (lhs, set, status);
2756 if (alloclhs == NULL)
2757 break;
2758 lhs = alloclhs;
2760 if (rhs->digits > reqdigits)
2762 allocrhs = decRoundOperand (rhs, set, status);
2763 if (allocrhs == NULL)
2764 break;
2765 rhs = allocrhs;
2768 #endif
2769 /* [following code does not require input rounding] */
2771 bits = (lhs->bits ^ rhs->bits) & DECNEG; /* assumed sign for divisions */
2773 /* handle infinities and NaNs */
2774 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2775 if (merged)
2776 { /* a special bit set */
2777 if (merged & (DECSNAN | DECNAN))
2778 { /* one or two NaNs */
2779 decNaNs (res, lhs, rhs, status);
2780 break;
2782 /* one or two infinities */
2783 if (decNumberIsInfinite (lhs))
2784 { /* LHS (dividend) is infinite */
2785 if (decNumberIsInfinite (rhs) || /* two infinities are invalid .. */
2786 op & (REMAINDER | REMNEAR))
2787 { /* as is remainder of infinity */
2788 *status |= DEC_Invalid_operation;
2789 break;
2791 /* [Note that infinity/0 raises no exceptions] */
2792 decNumberZero (res);
2793 res->bits = bits | DECINF; /* set +/- infinity */
2794 break;
2796 else
2797 { /* RHS (divisor) is infinite */
2798 residue = 0;
2799 if (op & (REMAINDER | REMNEAR))
2801 /* result is [finished clone of] lhs */
2802 decCopyFit (res, lhs, set, &residue, status);
2804 else
2805 { /* a division */
2806 decNumberZero (res);
2807 res->bits = bits; /* set +/- zero */
2808 /* for DIVIDEINT the exponent is always 0. For DIVIDE, result */
2809 /* is a 0 with infinitely negative exponent, clamped to minimum */
2810 if (op & DIVIDE)
2812 res->exponent = set->emin - set->digits + 1;
2813 *status |= DEC_Clamped;
2816 decFinish (res, set, &residue, status);
2817 break;
2821 /* handle 0 rhs (x/0) */
2822 if (ISZERO (rhs))
2823 { /* x/0 is always exceptional */
2824 if (ISZERO (lhs))
2826 decNumberZero (res); /* [after lhs test] */
2827 *status |= DEC_Division_undefined; /* 0/0 will become NaN */
2829 else
2831 decNumberZero (res);
2832 if (op & (REMAINDER | REMNEAR))
2833 *status |= DEC_Invalid_operation;
2834 else
2836 *status |= DEC_Division_by_zero; /* x/0 */
2837 res->bits = bits | DECINF; /* .. is +/- Infinity */
2840 break;
2843 /* handle 0 lhs (0/x) */
2844 if (ISZERO (lhs))
2845 { /* 0/x [x!=0] */
2846 #if DECSUBSET
2847 if (!set->extended)
2848 decNumberZero (res);
2849 else
2851 #endif
2852 if (op & DIVIDE)
2854 residue = 0;
2855 exponent = lhs->exponent - rhs->exponent; /* ideal exponent */
2856 decNumberCopy (res, lhs); /* [zeros always fit] */
2857 res->bits = bits; /* sign as computed */
2858 res->exponent = exponent; /* exponent, too */
2859 decFinalize (res, set, &residue, status); /* check exponent */
2861 else if (op & DIVIDEINT)
2863 decNumberZero (res); /* integer 0 */
2864 res->bits = bits; /* sign as computed */
2866 else
2867 { /* a remainder */
2868 exponent = rhs->exponent; /* [save in case overwrite] */
2869 decNumberCopy (res, lhs); /* [zeros always fit] */
2870 if (exponent < res->exponent)
2871 res->exponent = exponent; /* use lower */
2873 #if DECSUBSET
2875 #endif
2876 break;
2879 /* Precalculate exponent. This starts off adjusted (and hence fits */
2880 /* in 31 bits) and becomes the usual unadjusted exponent as the */
2881 /* division proceeds. The order of evaluation is important, here, */
2882 /* to avoid wrap. */
2883 exponent =
2884 (lhs->exponent + lhs->digits) - (rhs->exponent + rhs->digits);
2886 /* If the working exponent is -ve, then some quick exits are */
2887 /* possible because the quotient is known to be <1 */
2888 /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
2889 if (exponent < 0 && !(op == DIVIDE))
2891 if (op & DIVIDEINT)
2893 decNumberZero (res); /* integer part is 0 */
2894 #if DECSUBSET
2895 if (set->extended)
2896 #endif
2897 res->bits = bits; /* set +/- zero */
2898 break;
2900 /* we can fastpath remainders so long as the lhs has the */
2901 /* smaller (or equal) exponent */
2902 if (lhs->exponent <= rhs->exponent)
2904 if (op & REMAINDER || exponent < -1)
2906 /* It is REMAINDER or safe REMNEAR; result is [finished */
2907 /* clone of] lhs (r = x - 0*y) */
2908 residue = 0;
2909 decCopyFit (res, lhs, set, &residue, status);
2910 decFinish (res, set, &residue, status);
2911 break;
2913 /* [unsafe REMNEAR drops through] */
2915 } /* fastpaths */
2917 /* We need long (slow) division; roll up the sleeves... */
2919 /* The accumulator will hold the quotient of the division. */
2920 /* If it needs to be too long for stack storage, then allocate. */
2921 acclength = D2U (reqdigits + DECDPUN); /* in Units */
2922 if (acclength * sizeof (Unit) > sizeof (accbuff))
2924 allocacc = (Unit *) malloc (acclength * sizeof (Unit));
2925 if (allocacc == NULL)
2926 { /* hopeless -- abandon */
2927 *status |= DEC_Insufficient_storage;
2928 break;
2930 acc = allocacc; /* use the allocated space */
2933 /* var1 is the padded LHS ready for subtractions. */
2934 /* If it needs to be too long for stack storage, then allocate. */
2935 /* The maximum units we need for var1 (long subtraction) is: */
2936 /* Enough for */
2937 /* (rhs->digits+reqdigits-1) -- to allow full slide to right */
2938 /* or (lhs->digits) -- to allow for long lhs */
2939 /* whichever is larger */
2940 /* +1 -- for rounding of slide to right */
2941 /* +1 -- for leading 0s */
2942 /* +1 -- for pre-adjust if a remainder or DIVIDEINT */
2943 /* [Note: unused units do not participate in decUnitAddSub data] */
2944 maxdigits = rhs->digits + reqdigits - 1;
2945 if (lhs->digits > maxdigits)
2946 maxdigits = lhs->digits;
2947 var1units = D2U (maxdigits) + 2;
2948 /* allocate a guard unit above msu1 for REMAINDERNEAR */
2949 if (!(op & DIVIDE))
2950 var1units++;
2951 if ((var1units + 1) * sizeof (Unit) > sizeof (varbuff))
2953 varalloc = (Unit *) malloc ((var1units + 1) * sizeof (Unit));
2954 if (varalloc == NULL)
2955 { /* hopeless -- abandon */
2956 *status |= DEC_Insufficient_storage;
2957 break;
2959 var1 = varalloc; /* use the allocated space */
2962 /* Extend the lhs and rhs to full long subtraction length. The lhs */
2963 /* is truly extended into the var1 buffer, with 0 padding, so we can */
2964 /* subtract in place. The rhs (var2) has virtual padding */
2965 /* (implemented by decUnitAddSub). */
2966 /* We allocated one guard unit above msu1 for rem=rem+rem in REMAINDERNEAR */
2967 msu1 = var1 + var1units - 1; /* msu of var1 */
2968 source = lhs->lsu + D2U (lhs->digits) - 1; /* msu of input array */
2969 for (target = msu1; source >= lhs->lsu; source--, target--)
2970 *target = *source;
2971 for (; target >= var1; target--)
2972 *target = 0;
2974 /* rhs (var2) is left-aligned with var1 at the start */
2975 var2ulen = var1units; /* rhs logical length (units) */
2976 var2units = D2U (rhs->digits); /* rhs actual length (units) */
2977 var2 = rhs->lsu; /* -> rhs array */
2978 msu2 = var2 + var2units - 1; /* -> msu of var2 [never changes] */
2979 /* now set up the variables which we'll use for estimating the */
2980 /* multiplication factor. If these variables are not exact, we add */
2981 /* 1 to make sure that we never overestimate the multiplier. */
2982 msu2plus = *msu2; /* it's value .. */
2983 if (var2units > 1)
2984 msu2plus++; /* .. +1 if any more */
2985 msu2pair = (eInt) * msu2 * (DECDPUNMAX + 1); /* top two pair .. */
2986 if (var2units > 1)
2987 { /* .. [else treat 2nd as 0] */
2988 msu2pair += *(msu2 - 1); /* .. */
2989 if (var2units > 2)
2990 msu2pair++; /* .. +1 if any more */
2993 /* Since we are working in units, the units may have leading zeros, */
2994 /* but we calculated the exponent on the assumption that they are */
2995 /* both left-aligned. Adjust the exponent to compensate: add the */
2996 /* number of leading zeros in var1 msu and subtract those in var2 msu. */
2997 /* [We actually do this by counting the digits and negating, as */
2998 /* lead1=DECDPUN-digits1, and similarly for lead2.] */
2999 for (pow = &powers[1]; *msu1 >= *pow; pow++)
3000 exponent--;
3001 for (pow = &powers[1]; *msu2 >= *pow; pow++)
3002 exponent++;
3004 /* Now, if doing an integer divide or remainder, we want to ensure */
3005 /* that the result will be Unit-aligned. To do this, we shift the */
3006 /* var1 accumulator towards least if need be. (It's much easier to */
3007 /* do this now than to reassemble the residue afterwards, if we are */
3008 /* doing a remainder.) Also ensure the exponent is not negative. */
3009 if (!(op & DIVIDE))
3011 Unit *u;
3012 /* save the initial 'false' padding of var1, in digits */
3013 var1initpad = (var1units - D2U (lhs->digits)) * DECDPUN;
3014 /* Determine the shift to do. */
3015 if (exponent < 0)
3016 cut = -exponent;
3017 else
3018 cut = DECDPUN - exponent % DECDPUN;
3019 decShiftToLeast (var1, var1units, cut);
3020 exponent += cut; /* maintain numerical value */
3021 var1initpad -= cut; /* .. and reduce padding */
3022 /* clean any most-significant units we just emptied */
3023 for (u = msu1; cut >= DECDPUN; cut -= DECDPUN, u--)
3024 *u = 0;
3025 } /* align */
3026 else
3027 { /* is DIVIDE */
3028 maxexponent = lhs->exponent - rhs->exponent; /* save */
3029 /* optimization: if the first iteration will just produce 0, */
3030 /* preadjust to skip it [valid for DIVIDE only] */
3031 if (*msu1 < *msu2)
3033 var2ulen--; /* shift down */
3034 exponent -= DECDPUN; /* update the exponent */
3038 /* ---- start the long-division loops ------------------------------ */
3039 accunits = 0; /* no units accumulated yet */
3040 accdigits = 0; /* .. or digits */
3041 accnext = acc + acclength - 1; /* -> msu of acc [NB: allows digits+1] */
3042 for (;;)
3043 { /* outer forever loop */
3044 thisunit = 0; /* current unit assumed 0 */
3045 /* find the next unit */
3046 for (;;)
3047 { /* inner forever loop */
3048 /* strip leading zero units [from either pre-adjust or from */
3049 /* subtract last time around]. Leave at least one unit. */
3050 for (; *msu1 == 0 && msu1 > var1; msu1--)
3051 var1units--;
3053 if (var1units < var2ulen)
3054 break; /* var1 too low for subtract */
3055 if (var1units == var2ulen)
3056 { /* unit-by-unit compare needed */
3057 /* compare the two numbers, from msu */
3058 Unit *pv1, v2; /* units to compare */
3059 const Unit *pv2; /* units to compare */
3060 pv2 = msu2; /* -> msu */
3061 for (pv1 = msu1;; pv1--, pv2--)
3063 /* v1=*pv1 -- always OK */
3064 v2 = 0; /* assume in padding */
3065 if (pv2 >= var2)
3066 v2 = *pv2; /* in range */
3067 if (*pv1 != v2)
3068 break; /* no longer the same */
3069 if (pv1 == var1)
3070 break; /* done; leave pv1 as is */
3072 /* here when all inspected or a difference seen */
3073 if (*pv1 < v2)
3074 break; /* var1 too low to subtract */
3075 if (*pv1 == v2)
3076 { /* var1 == var2 */
3077 /* reach here if var1 and var2 are identical; subtraction */
3078 /* would increase digit by one, and the residue will be 0 so */
3079 /* we are done; leave the loop with residue set to 0. */
3080 thisunit++; /* as though subtracted */
3081 *var1 = 0; /* set var1 to 0 */
3082 var1units = 1; /* .. */
3083 break; /* from inner */
3084 } /* var1 == var2 */
3085 /* *pv1>v2. Prepare for real subtraction; the lengths are equal */
3086 /* Estimate the multiplier (there's always a msu1-1)... */
3087 /* Bring in two units of var2 to provide a good estimate. */
3088 mult =
3089 (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3090 *(msu1 - 1)) / msu2pair);
3091 } /* lengths the same */
3092 else
3093 { /* var1units > var2ulen, so subtraction is safe */
3094 /* The var2 msu is one unit towards the lsu of the var1 msu, */
3095 /* so we can only use one unit for var2. */
3096 mult =
3097 (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3098 *(msu1 - 1)) / msu2plus);
3100 if (mult == 0)
3101 mult = 1; /* must always be at least 1 */
3102 /* subtraction needed; var1 is > var2 */
3103 thisunit = (Unit) (thisunit + mult); /* accumulate */
3104 /* subtract var1-var2, into var1; only the overlap needs */
3105 /* processing, as we are in place */
3106 shift = var2ulen - var2units;
3107 #if DECTRACE
3108 decDumpAr ('1', &var1[shift], var1units - shift);
3109 decDumpAr ('2', var2, var2units);
3110 printf ("m=%d\n", -mult);
3111 #endif
3112 decUnitAddSub (&var1[shift], var1units - shift,
3113 var2, var2units, 0, &var1[shift], -mult);
3114 #if DECTRACE
3115 decDumpAr ('#', &var1[shift], var1units - shift);
3116 #endif
3117 /* var1 now probably has leading zeros; these are removed at the */
3118 /* top of the inner loop. */
3119 } /* inner loop */
3121 /* We have the next unit; unless it's a leading zero, add to acc */
3122 if (accunits != 0 || thisunit != 0)
3123 { /* put the unit we got */
3124 *accnext = thisunit; /* store in accumulator */
3125 /* account exactly for the digits we got */
3126 if (accunits == 0)
3128 accdigits++; /* at least one */
3129 for (pow = &powers[1]; thisunit >= *pow; pow++)
3130 accdigits++;
3132 else
3133 accdigits += DECDPUN;
3134 accunits++; /* update count */
3135 accnext--; /* ready for next */
3136 if (accdigits > reqdigits)
3137 break; /* we have all we need */
3140 /* if the residue is zero, we're done (unless divide or */
3141 /* divideInteger and we haven't got enough digits yet) */
3142 if (*var1 == 0 && var1units == 1)
3143 { /* residue is 0 */
3144 if (op & (REMAINDER | REMNEAR))
3145 break;
3146 if ((op & DIVIDE) && (exponent <= maxexponent))
3147 break;
3148 /* [drop through if divideInteger] */
3150 /* we've also done enough if calculating remainder or integer */
3151 /* divide and we just did the last ('units') unit */
3152 if (exponent == 0 && !(op & DIVIDE))
3153 break;
3155 /* to get here, var1 is less than var2, so divide var2 by the per- */
3156 /* Unit power of ten and go for the next digit */
3157 var2ulen--; /* shift down */
3158 exponent -= DECDPUN; /* update the exponent */
3159 } /* outer loop */
3161 /* ---- division is complete --------------------------------------- */
3162 /* here: acc has at least reqdigits+1 of good results (or fewer */
3163 /* if early stop), starting at accnext+1 (its lsu) */
3164 /* var1 has any residue at the stopping point */
3165 /* accunits is the number of digits we collected in acc */
3166 if (accunits == 0)
3167 { /* acc is 0 */
3168 accunits = 1; /* show we have one .. */
3169 accdigits = 1; /* .. */
3170 *accnext = 0; /* .. whose value is 0 */
3172 else
3173 accnext++; /* back to last placed */
3174 /* accnext now -> lowest unit of result */
3176 residue = 0; /* assume no residue */
3177 if (op & DIVIDE)
3179 /* record the presence of any residue, for rounding */
3180 if (*var1 != 0 || var1units > 1)
3181 residue = 1;
3182 else
3183 { /* no residue */
3184 /* We had an exact division; clean up spurious trailing 0s. */
3185 /* There will be at most DECDPUN-1, from the final multiply, */
3186 /* and then only if the result is non-0 (and even) and the */
3187 /* exponent is 'loose'. */
3188 #if DECDPUN>1
3189 Unit lsu = *accnext;
3190 if (!(lsu & 0x01) && (lsu != 0))
3192 /* count the trailing zeros */
3193 Int drop = 0;
3194 for (;; drop++)
3195 { /* [will terminate because lsu!=0] */
3196 if (exponent >= maxexponent)
3197 break; /* don't chop real 0s */
3198 #if DECDPUN<=4
3199 if ((lsu - QUOT10 (lsu, drop + 1)
3200 * powers[drop + 1]) != 0)
3201 break; /* found non-0 digit */
3202 #else
3203 if (lsu % powers[drop + 1] != 0)
3204 break; /* found non-0 digit */
3205 #endif
3206 exponent++;
3208 if (drop > 0)
3210 accunits = decShiftToLeast (accnext, accunits, drop);
3211 accdigits = decGetDigits (accnext, accunits);
3212 accunits = D2U (accdigits);
3213 /* [exponent was adjusted in the loop] */
3215 } /* neither odd nor 0 */
3216 #endif
3217 } /* exact divide */
3218 } /* divide */
3219 else /* op!=DIVIDE */
3221 /* check for coefficient overflow */
3222 if (accdigits + exponent > reqdigits)
3224 *status |= DEC_Division_impossible;
3225 break;
3227 if (op & (REMAINDER | REMNEAR))
3229 /* [Here, the exponent will be 0, because we adjusted var1 */
3230 /* appropriately.] */
3231 Int postshift; /* work */
3232 Flag wasodd = 0; /* integer was odd */
3233 Unit *quotlsu; /* for save */
3234 Int quotdigits; /* .. */
3236 /* Fastpath when residue is truly 0 is worthwhile [and */
3237 /* simplifies the code below] */
3238 if (*var1 == 0 && var1units == 1)
3239 { /* residue is 0 */
3240 Int exp = lhs->exponent; /* save min(exponents) */
3241 if (rhs->exponent < exp)
3242 exp = rhs->exponent;
3243 decNumberZero (res); /* 0 coefficient */
3244 #if DECSUBSET
3245 if (set->extended)
3246 #endif
3247 res->exponent = exp; /* .. with proper exponent */
3248 break;
3250 /* note if the quotient was odd */
3251 if (*accnext & 0x01)
3252 wasodd = 1; /* acc is odd */
3253 quotlsu = accnext; /* save in case need to reinspect */
3254 quotdigits = accdigits; /* .. */
3256 /* treat the residue, in var1, as the value to return, via acc */
3257 /* calculate the unused zero digits. This is the smaller of: */
3258 /* var1 initial padding (saved above) */
3259 /* var2 residual padding, which happens to be given by: */
3260 postshift =
3261 var1initpad + exponent - lhs->exponent + rhs->exponent;
3262 /* [the 'exponent' term accounts for the shifts during divide] */
3263 if (var1initpad < postshift)
3264 postshift = var1initpad;
3266 /* shift var1 the requested amount, and adjust its digits */
3267 var1units = decShiftToLeast (var1, var1units, postshift);
3268 accnext = var1;
3269 accdigits = decGetDigits (var1, var1units);
3270 accunits = D2U (accdigits);
3272 exponent = lhs->exponent; /* exponent is smaller of lhs & rhs */
3273 if (rhs->exponent < exponent)
3274 exponent = rhs->exponent;
3275 bits = lhs->bits; /* remainder sign is always as lhs */
3277 /* Now correct the result if we are doing remainderNear; if it */
3278 /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
3279 /* the integer was odd then the result should be rem-rhs. */
3280 if (op & REMNEAR)
3282 Int compare, tarunits; /* work */
3283 Unit *up; /* .. */
3286 /* calculate remainder*2 into the var1 buffer (which has */
3287 /* 'headroom' of an extra unit and hence enough space) */
3288 /* [a dedicated 'double' loop would be faster, here] */
3289 tarunits =
3290 decUnitAddSub (accnext, accunits, accnext, accunits, 0,
3291 accnext, 1);
3292 /* decDumpAr('r', accnext, tarunits); */
3294 /* Here, accnext (var1) holds tarunits Units with twice the */
3295 /* remainder's coefficient, which we must now compare to the */
3296 /* RHS. The remainder's exponent may be smaller than the RHS's. */
3297 compare =
3298 decUnitCompare (accnext, tarunits, rhs->lsu,
3299 D2U (rhs->digits),
3300 rhs->exponent - exponent);
3301 if (compare == BADINT)
3302 { /* deep trouble */
3303 *status |= DEC_Insufficient_storage;
3304 break;
3307 /* now restore the remainder by dividing by two; we know the */
3308 /* lsu is even. */
3309 for (up = accnext; up < accnext + tarunits; up++)
3311 Int half; /* half to add to lower unit */
3312 half = *up & 0x01;
3313 *up /= 2; /* [shift] */
3314 if (!half)
3315 continue;
3316 *(up - 1) += (DECDPUNMAX + 1) / 2;
3318 /* [accunits still describes the original remainder length] */
3320 if (compare > 0 || (compare == 0 && wasodd))
3321 { /* adjustment needed */
3322 Int exp, expunits, exprem; /* work */
3323 /* This is effectively causing round-up of the quotient, */
3324 /* so if it was the rare case where it was full and all */
3325 /* nines, it would overflow and hence division-impossible */
3326 /* should be raised */
3327 Flag allnines = 0; /* 1 if quotient all nines */
3328 if (quotdigits == reqdigits)
3329 { /* could be borderline */
3330 for (up = quotlsu;; up++)
3332 if (quotdigits > DECDPUN)
3334 if (*up != DECDPUNMAX)
3335 break; /* non-nines */
3337 else
3338 { /* this is the last Unit */
3339 if (*up == powers[quotdigits] - 1)
3340 allnines = 1;
3341 break;
3343 quotdigits -= DECDPUN; /* checked those digits */
3344 } /* up */
3345 } /* borderline check */
3346 if (allnines)
3348 *status |= DEC_Division_impossible;
3349 break;
3352 /* we need rem-rhs; the sign will invert. Again we can */
3353 /* safely use var1 for the working Units array. */
3354 exp = rhs->exponent - exponent; /* RHS padding needed */
3355 /* Calculate units and remainder from exponent. */
3356 expunits = exp / DECDPUN;
3357 exprem = exp % DECDPUN;
3358 /* subtract [A+B*(-m)]; the result will always be negative */
3359 accunits = -decUnitAddSub (accnext, accunits,
3360 rhs->lsu, D2U (rhs->digits),
3361 expunits, accnext,
3362 -(Int) powers[exprem]);
3363 accdigits = decGetDigits (accnext, accunits); /* count digits exactly */
3364 accunits = D2U (accdigits); /* and recalculate the units for copy */
3365 /* [exponent is as for original remainder] */
3366 bits ^= DECNEG; /* flip the sign */
3368 } /* REMNEAR */
3369 } /* REMAINDER or REMNEAR */
3370 } /* not DIVIDE */
3372 /* Set exponent and bits */
3373 res->exponent = exponent;
3374 res->bits = (uByte) (bits & DECNEG); /* [cleaned] */
3376 /* Now the coefficient. */
3377 decSetCoeff (res, set, accnext, accdigits, &residue, status);
3379 decFinish (res, set, &residue, status); /* final cleanup */
3381 #if DECSUBSET
3382 /* If a divide then strip trailing zeros if subset [after round] */
3383 if (!set->extended && (op == DIVIDE))
3384 decTrim (res, 0, &dropped);
3385 #endif
3387 while (0); /* end protected */
3389 if (varalloc != NULL)
3390 free (varalloc); /* drop any storage we used */
3391 if (allocacc != NULL)
3392 free (allocacc); /* .. */
3393 if (allocrhs != NULL)
3394 free (allocrhs); /* .. */
3395 if (alloclhs != NULL)
3396 free (alloclhs); /* .. */
3397 return res;
3400 /* ------------------------------------------------------------------ */
3401 /* decMultiplyOp -- multiplication operation */
3402 /* */
3403 /* This routine performs the multiplication C=A x B. */
3404 /* */
3405 /* res is C, the result. C may be A and/or B (e.g., X=X*X) */
3406 /* lhs is A */
3407 /* rhs is B */
3408 /* set is the context */
3409 /* status is the usual accumulator */
3410 /* */
3411 /* C must have space for set->digits digits. */
3412 /* */
3413 /* ------------------------------------------------------------------ */
3414 /* Note: We use 'long' multiplication rather than Karatsuba, as the */
3415 /* latter would give only a minor improvement for the short numbers */
3416 /* we expect to handle most (and uses much more memory). */
3417 /* */
3418 /* We always have to use a buffer for the accumulator. */
3419 /* ------------------------------------------------------------------ */
3420 static decNumber *
3421 decMultiplyOp (decNumber * res, const decNumber * lhs,
3422 const decNumber * rhs, decContext * set, uInt * status)
3424 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3425 decNumber *allocrhs = NULL; /* .., rhs */
3426 Unit accbuff[D2U (DECBUFFER * 2 + 1)]; /* local buffer (+1 in case DECBUFFER==0) */
3427 Unit *acc = accbuff; /* -> accumulator array for exact result */
3428 Unit *allocacc = NULL; /* -> allocated buffer, iff allocated */
3429 const Unit *mer, *mermsup; /* work */
3430 Int accunits; /* Units of accumulator in use */
3431 Int madlength; /* Units in multiplicand */
3432 Int shift; /* Units to shift multiplicand by */
3433 Int need; /* Accumulator units needed */
3434 Int exponent; /* work */
3435 Int residue = 0; /* rounding residue */
3436 uByte bits; /* result sign */
3437 uByte merged; /* merged flags */
3439 #if DECCHECK
3440 if (decCheckOperands (res, lhs, rhs, set))
3441 return res;
3442 #endif
3445 { /* protect allocated storage */
3446 #if DECSUBSET
3447 if (!set->extended)
3449 /* reduce operands and set lostDigits status, as needed */
3450 if (lhs->digits > set->digits)
3452 alloclhs = decRoundOperand (lhs, set, status);
3453 if (alloclhs == NULL)
3454 break;
3455 lhs = alloclhs;
3457 if (rhs->digits > set->digits)
3459 allocrhs = decRoundOperand (rhs, set, status);
3460 if (allocrhs == NULL)
3461 break;
3462 rhs = allocrhs;
3465 #endif
3466 /* [following code does not require input rounding] */
3468 /* precalculate result sign */
3469 bits = (uByte) ((lhs->bits ^ rhs->bits) & DECNEG);
3471 /* handle infinities and NaNs */
3472 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3473 if (merged)
3474 { /* a special bit set */
3475 if (merged & (DECSNAN | DECNAN))
3476 { /* one or two NaNs */
3477 decNaNs (res, lhs, rhs, status);
3478 break;
3480 /* one or two infinities. Infinity * 0 is invalid */
3481 if (((lhs->bits & DECSPECIAL) == 0 && ISZERO (lhs))
3482 || ((rhs->bits & DECSPECIAL) == 0 && ISZERO (rhs)))
3484 *status |= DEC_Invalid_operation;
3485 break;
3487 decNumberZero (res);
3488 res->bits = bits | DECINF; /* infinity */
3489 break;
3492 /* For best speed, as in DMSRCN, we use the shorter number as the */
3493 /* multiplier (rhs) and the longer as the multiplicand (lhs) */
3494 if (lhs->digits < rhs->digits)
3495 { /* swap... */
3496 const decNumber *hold = lhs;
3497 lhs = rhs;
3498 rhs = hold;
3501 /* if accumulator is too long for local storage, then allocate */
3502 need = D2U (lhs->digits) + D2U (rhs->digits); /* maximum units in result */
3503 if (need * sizeof (Unit) > sizeof (accbuff))
3505 allocacc = (Unit *) malloc (need * sizeof (Unit));
3506 if (allocacc == NULL)
3508 *status |= DEC_Insufficient_storage;
3509 break;
3511 acc = allocacc; /* use the allocated space */
3514 /* Now the main long multiplication loop */
3515 /* Unlike the equivalent in the IBM Java implementation, there */
3516 /* is no advantage in calculating from msu to lsu. So we do it */
3517 /* by the book, as it were. */
3518 /* Each iteration calculates ACC=ACC+MULTAND*MULT */
3519 accunits = 1; /* accumulator starts at '0' */
3520 *acc = 0; /* .. (lsu=0) */
3521 shift = 0; /* no multiplicand shift at first */
3522 madlength = D2U (lhs->digits); /* we know this won't change */
3523 mermsup = rhs->lsu + D2U (rhs->digits); /* -> msu+1 of multiplier */
3525 for (mer = rhs->lsu; mer < mermsup; mer++)
3527 /* Here, *mer is the next Unit in the multiplier to use */
3528 /* If non-zero [optimization] add it... */
3529 if (*mer != 0)
3531 accunits =
3532 decUnitAddSub (&acc[shift], accunits - shift, lhs->lsu,
3533 madlength, 0, &acc[shift], *mer) + shift;
3535 else
3536 { /* extend acc with a 0; we'll use it shortly */
3537 /* [this avoids length of <=0 later] */
3538 *(acc + accunits) = 0;
3539 accunits++;
3541 /* multiply multiplicand by 10**DECDPUN for next Unit to left */
3542 shift++; /* add this for 'logical length' */
3543 } /* n */
3544 #if DECTRACE
3545 /* Show exact result */
3546 decDumpAr ('*', acc, accunits);
3547 #endif
3549 /* acc now contains the exact result of the multiplication */
3550 /* Build a decNumber from it, noting if any residue */
3551 res->bits = bits; /* set sign */
3552 res->digits = decGetDigits (acc, accunits); /* count digits exactly */
3554 /* We might have a 31-bit wrap in calculating the exponent. */
3555 /* This can only happen if both input exponents are negative and */
3556 /* both their magnitudes are large. If we did wrap, we set a safe */
3557 /* very negative exponent, from which decFinalize() will raise a */
3558 /* hard underflow. */
3559 exponent = lhs->exponent + rhs->exponent; /* calculate exponent */
3560 if (lhs->exponent < 0 && rhs->exponent < 0 && exponent > 0)
3561 exponent = -2 * DECNUMMAXE; /* force underflow */
3562 res->exponent = exponent; /* OK to overwrite now */
3564 /* Set the coefficient. If any rounding, residue records */
3565 decSetCoeff (res, set, acc, res->digits, &residue, status);
3567 decFinish (res, set, &residue, status); /* final cleanup */
3569 while (0); /* end protected */
3571 if (allocacc != NULL)
3572 free (allocacc); /* drop any storage we used */
3573 if (allocrhs != NULL)
3574 free (allocrhs); /* .. */
3575 if (alloclhs != NULL)
3576 free (alloclhs); /* .. */
3577 return res;
3580 /* ------------------------------------------------------------------ */
3581 /* decQuantizeOp -- force exponent to requested value */
3582 /* */
3583 /* This computes C = op(A, B), where op adjusts the coefficient */
3584 /* of C (by rounding or shifting) such that the exponent (-scale) */
3585 /* of C has the value B or matches the exponent of B. */
3586 /* The numerical value of C will equal A, except for the effects of */
3587 /* any rounding that occurred. */
3588 /* */
3589 /* res is C, the result. C may be A or B */
3590 /* lhs is A, the number to adjust */
3591 /* rhs is B, the requested exponent */
3592 /* set is the context */
3593 /* quant is 1 for quantize or 0 for rescale */
3594 /* status is the status accumulator (this can be called without */
3595 /* risk of control loss) */
3596 /* */
3597 /* C must have space for set->digits digits. */
3598 /* */
3599 /* Unless there is an error or the result is infinite, the exponent */
3600 /* after the operation is guaranteed to be that requested. */
3601 /* ------------------------------------------------------------------ */
3602 static decNumber *
3603 decQuantizeOp (decNumber * res, const decNumber * lhs,
3604 const decNumber * rhs, decContext * set, Flag quant, uInt * status)
3606 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3607 decNumber *allocrhs = NULL; /* .., rhs */
3608 const decNumber *inrhs = rhs; /* save original rhs */
3609 Int reqdigits = set->digits; /* requested DIGITS */
3610 Int reqexp; /* requested exponent [-scale] */
3611 Int residue = 0; /* rounding residue */
3612 uByte merged; /* merged flags */
3613 Int etiny = set->emin - (set->digits - 1);
3615 #if DECCHECK
3616 if (decCheckOperands (res, lhs, rhs, set))
3617 return res;
3618 #endif
3621 { /* protect allocated storage */
3622 #if DECSUBSET
3623 if (!set->extended)
3625 /* reduce operands and set lostDigits status, as needed */
3626 if (lhs->digits > reqdigits)
3628 alloclhs = decRoundOperand (lhs, set, status);
3629 if (alloclhs == NULL)
3630 break;
3631 lhs = alloclhs;
3633 if (rhs->digits > reqdigits)
3634 { /* [this only checks lostDigits] */
3635 allocrhs = decRoundOperand (rhs, set, status);
3636 if (allocrhs == NULL)
3637 break;
3638 rhs = allocrhs;
3641 #endif
3642 /* [following code does not require input rounding] */
3644 /* Handle special values */
3645 merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3646 if ((lhs->bits | rhs->bits) & DECSPECIAL)
3648 /* NaNs get usual processing */
3649 if (merged & (DECSNAN | DECNAN))
3650 decNaNs (res, lhs, rhs, status);
3651 /* one infinity but not both is bad */
3652 else if ((lhs->bits ^ rhs->bits) & DECINF)
3653 *status |= DEC_Invalid_operation;
3654 /* both infinity: return lhs */
3655 else
3656 decNumberCopy (res, lhs); /* [nop if in place] */
3657 break;
3660 /* set requested exponent */
3661 if (quant)
3662 reqexp = inrhs->exponent; /* quantize -- match exponents */
3663 else
3664 { /* rescale -- use value of rhs */
3665 /* Original rhs must be an integer that fits and is in range */
3666 #if DECSUBSET
3667 reqexp = decGetInt (inrhs, set);
3668 #else
3669 reqexp = decGetInt (inrhs);
3670 #endif
3673 #if DECSUBSET
3674 if (!set->extended)
3675 etiny = set->emin; /* no subnormals */
3676 #endif
3678 if (reqexp == BADINT /* bad (rescale only) or .. */
3679 || (reqexp < etiny) /* < lowest */
3680 || (reqexp > set->emax))
3681 { /* > Emax */
3682 *status |= DEC_Invalid_operation;
3683 break;
3686 /* we've processed the RHS, so we can overwrite it now if necessary */
3687 if (ISZERO (lhs))
3688 { /* zero coefficient unchanged */
3689 decNumberCopy (res, lhs); /* [nop if in place] */
3690 res->exponent = reqexp; /* .. just set exponent */
3691 #if DECSUBSET
3692 if (!set->extended)
3693 res->bits = 0; /* subset specification; no -0 */
3694 #endif
3696 else
3697 { /* non-zero lhs */
3698 Int adjust = reqexp - lhs->exponent; /* digit adjustment needed */
3699 /* if adjusted coefficient will not fit, give up now */
3700 if ((lhs->digits - adjust) > reqdigits)
3702 *status |= DEC_Invalid_operation;
3703 break;
3706 if (adjust > 0)
3707 { /* increasing exponent */
3708 /* this will decrease the length of the coefficient by adjust */
3709 /* digits, and must round as it does so */
3710 decContext workset; /* work */
3711 workset = *set; /* clone rounding, etc. */
3712 workset.digits = lhs->digits - adjust; /* set requested length */
3713 /* [note that the latter can be <1, here] */
3714 decCopyFit (res, lhs, &workset, &residue, status); /* fit to result */
3715 decApplyRound (res, &workset, residue, status); /* .. and round */
3716 residue = 0; /* [used] */
3717 /* If we rounded a 999s case, exponent will be off by one; */
3718 /* adjust back if so. */
3719 if (res->exponent > reqexp)
3721 res->digits = decShiftToMost (res->lsu, res->digits, 1); /* shift */
3722 res->exponent--; /* (re)adjust the exponent. */
3724 #if DECSUBSET
3725 if (ISZERO (res) && !set->extended)
3726 res->bits = 0; /* subset; no -0 */
3727 #endif
3728 } /* increase */
3729 else /* adjust<=0 */
3730 { /* decreasing or = exponent */
3731 /* this will increase the length of the coefficient by -adjust */
3732 /* digits, by adding trailing zeros. */
3733 decNumberCopy (res, lhs); /* [it will fit] */
3734 /* if padding needed (adjust<0), add it now... */
3735 if (adjust < 0)
3737 res->digits =
3738 decShiftToMost (res->lsu, res->digits, -adjust);
3739 res->exponent += adjust; /* adjust the exponent */
3741 } /* decrease */
3742 } /* non-zero */
3744 /* Check for overflow [do not use Finalize in this case, as an */
3745 /* overflow here is a "don't fit" situation] */
3746 if (res->exponent > set->emax - res->digits + 1)
3747 { /* too big */
3748 *status |= DEC_Invalid_operation;
3749 break;
3751 else
3753 decFinalize (res, set, &residue, status); /* set subnormal flags */
3754 *status &= ~DEC_Underflow; /* suppress Underflow [754r] */
3757 while (0); /* end protected */
3759 if (allocrhs != NULL)
3760 free (allocrhs); /* drop any storage we used */
3761 if (alloclhs != NULL)
3762 free (alloclhs); /* .. */
3763 return res;
3766 /* ------------------------------------------------------------------ */
3767 /* decCompareOp -- compare, min, or max two Numbers */
3768 /* */
3769 /* This computes C = A ? B and returns the signum (as a Number) */
3770 /* for COMPARE or the maximum or minimum (for COMPMAX and COMPMIN). */
3771 /* */
3772 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
3773 /* lhs is A */
3774 /* rhs is B */
3775 /* set is the context */
3776 /* op is the operation flag */
3777 /* status is the usual accumulator */
3778 /* */
3779 /* C must have space for one digit for COMPARE or set->digits for */
3780 /* COMPMAX and COMPMIN. */
3781 /* ------------------------------------------------------------------ */
3782 /* The emphasis here is on speed for common cases, and avoiding */
3783 /* coefficient comparison if possible. */
3784 /* ------------------------------------------------------------------ */
3785 decNumber *
3786 decCompareOp (decNumber * res, const decNumber * lhs, const decNumber * rhs,
3787 decContext * set, Flag op, uInt * status)
3789 decNumber *alloclhs = NULL; /* non-NULL if rounded lhs allocated */
3790 decNumber *allocrhs = NULL; /* .., rhs */
3791 Int result = 0; /* default result value */
3792 uByte merged; /* merged flags */
3793 uByte bits = 0; /* non-0 for NaN */
3795 #if DECCHECK
3796 if (decCheckOperands (res, lhs, rhs, set))
3797 return res;
3798 #endif
3801 { /* protect allocated storage */
3802 #if DECSUBSET
3803 if (!set->extended)
3805 /* reduce operands and set lostDigits status, as needed */
3806 if (lhs->digits > set->digits)
3808 alloclhs = decRoundOperand (lhs, set, status);
3809 if (alloclhs == NULL)
3811 result = BADINT;
3812 break;
3814 lhs = alloclhs;
3816 if (rhs->digits > set->digits)
3818 allocrhs = decRoundOperand (rhs, set, status);
3819 if (allocrhs == NULL)
3821 result = BADINT;
3822 break;
3824 rhs = allocrhs;
3827 #endif
3828 /* [following code does not require input rounding] */
3830 /* handle NaNs now; let infinities drop through */
3831 /* +++ review sNaN handling with 754r, for now assumes sNaN */
3832 /* (even just one) leads to NaN. */
3833 merged = (lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
3834 if (merged)
3835 { /* a NaN bit set */
3836 if (op == COMPARE);
3837 else if (merged & DECSNAN);
3838 else
3839 { /* 754r rules for MIN and MAX ignore single NaN */
3840 /* here if MIN or MAX, and one or two quiet NaNs */
3841 if (lhs->bits & rhs->bits & DECNAN);
3842 else
3843 { /* just one quiet NaN */
3844 /* force choice to be the non-NaN operand */
3845 op = COMPMAX;
3846 if (lhs->bits & DECNAN)
3847 result = -1; /* pick rhs */
3848 else
3849 result = +1; /* pick lhs */
3850 break;
3853 op = COMPNAN; /* use special path */
3854 decNaNs (res, lhs, rhs, status);
3855 break;
3858 result = decCompare (lhs, rhs); /* we have numbers */
3860 while (0); /* end protected */
3862 if (result == BADINT)
3863 *status |= DEC_Insufficient_storage; /* rare */
3864 else
3866 if (op == COMPARE)
3867 { /* return signum */
3868 decNumberZero (res); /* [always a valid result] */
3869 if (result == 0)
3870 res->bits = bits; /* (maybe qNaN) */
3871 else
3873 *res->lsu = 1;
3874 if (result < 0)
3875 res->bits = DECNEG;
3878 else if (op == COMPNAN); /* special, drop through */
3879 else
3880 { /* MAX or MIN, non-NaN result */
3881 Int residue = 0; /* rounding accumulator */
3882 /* choose the operand for the result */
3883 const decNumber *choice;
3884 if (result == 0)
3885 { /* operands are numerically equal */
3886 /* choose according to sign then exponent (see 754r) */
3887 uByte slhs = (lhs->bits & DECNEG);
3888 uByte srhs = (rhs->bits & DECNEG);
3889 #if DECSUBSET
3890 if (!set->extended)
3891 { /* subset: force left-hand */
3892 op = COMPMAX;
3893 result = +1;
3895 else
3896 #endif
3897 if (slhs != srhs)
3898 { /* signs differ */
3899 if (slhs)
3900 result = -1; /* rhs is max */
3901 else
3902 result = +1; /* lhs is max */
3904 else if (slhs && srhs)
3905 { /* both negative */
3906 if (lhs->exponent < rhs->exponent)
3907 result = +1;
3908 else
3909 result = -1;
3910 /* [if equal, we use lhs, technically identical] */
3912 else
3913 { /* both positive */
3914 if (lhs->exponent > rhs->exponent)
3915 result = +1;
3916 else
3917 result = -1;
3918 /* [ditto] */
3920 } /* numerically equal */
3921 /* here result will be non-0 */
3922 if (op == COMPMIN)
3923 result = -result; /* reverse if looking for MIN */
3924 choice = (result > 0 ? lhs : rhs); /* choose */
3925 /* copy chosen to result, rounding if need be */
3926 decCopyFit (res, choice, set, &residue, status);
3927 decFinish (res, set, &residue, status);
3930 if (allocrhs != NULL)
3931 free (allocrhs); /* free any storage we used */
3932 if (alloclhs != NULL)
3933 free (alloclhs); /* .. */
3934 return res;
3937 /* ------------------------------------------------------------------ */
3938 /* decCompare -- compare two decNumbers by numerical value */
3939 /* */
3940 /* This routine compares A ? B without altering them. */
3941 /* */
3942 /* Arg1 is A, a decNumber which is not a NaN */
3943 /* Arg2 is B, a decNumber which is not a NaN */
3944 /* */
3945 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
3946 /* (the only possible failure is an allocation error) */
3947 /* ------------------------------------------------------------------ */
3948 /* This could be merged into decCompareOp */
3949 static Int
3950 decCompare (const decNumber * lhs, const decNumber * rhs)
3952 Int result; /* result value */
3953 Int sigr; /* rhs signum */
3954 Int compare; /* work */
3955 result = 1; /* assume signum(lhs) */
3956 if (ISZERO (lhs))
3957 result = 0;
3958 else if (decNumberIsNegative (lhs))
3959 result = -1;
3960 sigr = 1; /* compute signum(rhs) */
3961 if (ISZERO (rhs))
3962 sigr = 0;
3963 else if (decNumberIsNegative (rhs))
3964 sigr = -1;
3965 if (result > sigr)
3966 return +1; /* L > R, return 1 */
3967 if (result < sigr)
3968 return -1; /* R < L, return -1 */
3970 /* signums are the same */
3971 if (result == 0)
3972 return 0; /* both 0 */
3973 /* Both non-zero */
3974 if ((lhs->bits | rhs->bits) & DECINF)
3975 { /* one or more infinities */
3976 if (lhs->bits == rhs->bits)
3977 result = 0; /* both the same */
3978 else if (decNumberIsInfinite (rhs))
3979 result = -result;
3980 return result;
3983 /* we must compare the coefficients, allowing for exponents */
3984 if (lhs->exponent > rhs->exponent)
3985 { /* LHS exponent larger */
3986 /* swap sides, and sign */
3987 const decNumber *temp = lhs;
3988 lhs = rhs;
3989 rhs = temp;
3990 result = -result;
3993 compare = decUnitCompare (lhs->lsu, D2U (lhs->digits),
3994 rhs->lsu, D2U (rhs->digits),
3995 rhs->exponent - lhs->exponent);
3997 if (compare != BADINT)
3998 compare *= result; /* comparison succeeded */
3999 return compare; /* what we got */
4002 /* ------------------------------------------------------------------ */
4003 /* decUnitCompare -- compare two >=0 integers in Unit arrays */
4004 /* */
4005 /* This routine compares A ? B*10**E where A and B are unit arrays */
4006 /* A is a plain integer */
4007 /* B has an exponent of E (which must be non-negative) */
4008 /* */
4009 /* Arg1 is A first Unit (lsu) */
4010 /* Arg2 is A length in Units */
4011 /* Arg3 is B first Unit (lsu) */
4012 /* Arg4 is B length in Units */
4013 /* Arg5 is E */
4014 /* */
4015 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
4016 /* (the only possible failure is an allocation error) */
4017 /* ------------------------------------------------------------------ */
4018 static Int
4019 decUnitCompare (const Unit * a, Int alength, const Unit * b, Int blength, Int exp)
4021 Unit *acc; /* accumulator for result */
4022 Unit accbuff[D2U (DECBUFFER + 1)]; /* local buffer */
4023 Unit *allocacc = NULL; /* -> allocated acc buffer, iff allocated */
4024 Int accunits, need; /* units in use or needed for acc */
4025 const Unit *l, *r, *u; /* work */
4026 Int expunits, exprem, result; /* .. */
4028 if (exp == 0)
4029 { /* aligned; fastpath */
4030 if (alength > blength)
4031 return 1;
4032 if (alength < blength)
4033 return -1;
4034 /* same number of units in both -- need unit-by-unit compare */
4035 l = a + alength - 1;
4036 r = b + alength - 1;
4037 for (; l >= a; l--, r--)
4039 if (*l > *r)
4040 return 1;
4041 if (*l < *r)
4042 return -1;
4044 return 0; /* all units match */
4045 } /* aligned */
4047 /* Unaligned. If one is >1 unit longer than the other, padded */
4048 /* approximately, then we can return easily */
4049 if (alength > blength + (Int) D2U (exp))
4050 return 1;
4051 if (alength + 1 < blength + (Int) D2U (exp))
4052 return -1;
4054 /* We need to do a real subtract. For this, we need a result buffer */
4055 /* even though we only are interested in the sign. Its length needs */
4056 /* to be the larger of alength and padded blength, +2 */
4057 need = blength + D2U (exp); /* maximum real length of B */
4058 if (need < alength)
4059 need = alength;
4060 need += 2;
4061 acc = accbuff; /* assume use local buffer */
4062 if (need * sizeof (Unit) > sizeof (accbuff))
4064 allocacc = (Unit *) malloc (need * sizeof (Unit));
4065 if (allocacc == NULL)
4066 return BADINT; /* hopeless -- abandon */
4067 acc = allocacc;
4069 /* Calculate units and remainder from exponent. */
4070 expunits = exp / DECDPUN;
4071 exprem = exp % DECDPUN;
4072 /* subtract [A+B*(-m)] */
4073 accunits = decUnitAddSub (a, alength, b, blength, expunits, acc,
4074 -(Int) powers[exprem]);
4075 /* [UnitAddSub result may have leading zeros, even on zero] */
4076 if (accunits < 0)
4077 result = -1; /* negative result */
4078 else
4079 { /* non-negative result */
4080 /* check units of the result before freeing any storage */
4081 for (u = acc; u < acc + accunits - 1 && *u == 0;)
4082 u++;
4083 result = (*u == 0 ? 0 : +1);
4085 /* clean up and return the result */
4086 if (allocacc != NULL)
4087 free (allocacc); /* drop any storage we used */
4088 return result;
4091 /* ------------------------------------------------------------------ */
4092 /* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays */
4093 /* */
4094 /* This routine performs the calculation: */
4095 /* */
4096 /* C=A+(B*M) */
4097 /* */
4098 /* Where M is in the range -DECDPUNMAX through +DECDPUNMAX. */
4099 /* */
4100 /* A may be shorter or longer than B. */
4101 /* */
4102 /* Leading zeros are not removed after a calculation. The result is */
4103 /* either the same length as the longer of A and B (adding any */
4104 /* shift), or one Unit longer than that (if a Unit carry occurred). */
4105 /* */
4106 /* A and B content are not altered unless C is also A or B. */
4107 /* C may be the same array as A or B, but only if no zero padding is */
4108 /* requested (that is, C may be B only if bshift==0). */
4109 /* C is filled from the lsu; only those units necessary to complete */
4110 /* the calculation are referenced. */
4111 /* */
4112 /* Arg1 is A first Unit (lsu) */
4113 /* Arg2 is A length in Units */
4114 /* Arg3 is B first Unit (lsu) */
4115 /* Arg4 is B length in Units */
4116 /* Arg5 is B shift in Units (>=0; pads with 0 units if positive) */
4117 /* Arg6 is C first Unit (lsu) */
4118 /* Arg7 is M, the multiplier */
4119 /* */
4120 /* returns the count of Units written to C, which will be non-zero */
4121 /* and negated if the result is negative. That is, the sign of the */
4122 /* returned Int is the sign of the result (positive for zero) and */
4123 /* the absolute value of the Int is the count of Units. */
4124 /* */
4125 /* It is the caller's responsibility to make sure that C size is */
4126 /* safe, allowing space if necessary for a one-Unit carry. */
4127 /* */
4128 /* This routine is severely performance-critical; *any* change here */
4129 /* must be measured (timed) to assure no performance degradation. */
4130 /* In particular, trickery here tends to be counter-productive, as */
4131 /* increased complexity of code hurts register optimizations on */
4132 /* register-poor architectures. Avoiding divisions is nearly */
4133 /* always a Good Idea, however. */
4134 /* */
4135 /* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark */
4136 /* (IBM Warwick, UK) for some of the ideas used in this routine. */
4137 /* ------------------------------------------------------------------ */
4138 static Int
4139 decUnitAddSub (const Unit * a, Int alength,
4140 const Unit * b, Int blength, Int bshift, Unit * c, Int m)
4142 const Unit *alsu = a; /* A lsu [need to remember it] */
4143 Unit *clsu = c; /* C ditto */
4144 Unit *minC; /* low water mark for C */
4145 Unit *maxC; /* high water mark for C */
4146 eInt carry = 0; /* carry integer (could be Long) */
4147 Int add; /* work */
4148 #if DECDPUN==4 /* myriadal */
4149 Int est; /* estimated quotient */
4150 #endif
4152 #if DECTRACE
4153 if (alength < 1 || blength < 1)
4154 printf ("decUnitAddSub: alen blen m %d %d [%d]\n", alength, blength, m);
4155 #endif
4157 maxC = c + alength; /* A is usually the longer */
4158 minC = c + blength; /* .. and B the shorter */
4159 if (bshift != 0)
4160 { /* B is shifted; low As copy across */
4161 minC += bshift;
4162 /* if in place [common], skip copy unless there's a gap [rare] */
4163 if (a == c && bshift <= alength)
4165 c += bshift;
4166 a += bshift;
4168 else
4169 for (; c < clsu + bshift; a++, c++)
4170 { /* copy needed */
4171 if (a < alsu + alength)
4172 *c = *a;
4173 else
4174 *c = 0;
4177 if (minC > maxC)
4178 { /* swap */
4179 Unit *hold = minC;
4180 minC = maxC;
4181 maxC = hold;
4184 /* For speed, we do the addition as two loops; the first where both A */
4185 /* and B contribute, and the second (if necessary) where only one or */
4186 /* other of the numbers contribute. */
4187 /* Carry handling is the same (i.e., duplicated) in each case. */
4188 for (; c < minC; c++)
4190 carry += *a;
4191 a++;
4192 carry += ((eInt) * b) * m; /* [special-casing m=1/-1 */
4193 b++; /* here is not a win] */
4194 /* here carry is new Unit of digits; it could be +ve or -ve */
4195 if ((ueInt) carry <= DECDPUNMAX)
4196 { /* fastpath 0-DECDPUNMAX */
4197 *c = (Unit) carry;
4198 carry = 0;
4199 continue;
4201 /* remainder operator is undefined if negative, so we must test */
4202 #if DECDPUN==4 /* use divide-by-multiply */
4203 if (carry >= 0)
4205 est = (((ueInt) carry >> 11) * 53687) >> 18;
4206 *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
4207 carry = est; /* likely quotient [89%] */
4208 if (*c < DECDPUNMAX + 1)
4209 continue; /* estimate was correct */
4210 carry++;
4211 *c -= DECDPUNMAX + 1;
4212 continue;
4214 /* negative case */
4215 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4216 est = (((ueInt) carry >> 11) * 53687) >> 18;
4217 *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4218 carry = est - (DECDPUNMAX + 1); /* correctly negative */
4219 if (*c < DECDPUNMAX + 1)
4220 continue; /* was OK */
4221 carry++;
4222 *c -= DECDPUNMAX + 1;
4223 #else
4224 if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4225 { /* fastpath carry +1 */
4226 *c = (Unit) (carry - (DECDPUNMAX + 1)); /* [helps additions] */
4227 carry = 1;
4228 continue;
4230 if (carry >= 0)
4232 *c = (Unit) (carry % (DECDPUNMAX + 1));
4233 carry = carry / (DECDPUNMAX + 1);
4234 continue;
4236 /* negative case */
4237 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4238 *c = (Unit) (carry % (DECDPUNMAX + 1));
4239 carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4240 #endif
4241 } /* c */
4243 /* we now may have one or other to complete */
4244 /* [pretest to avoid loop setup/shutdown] */
4245 if (c < maxC)
4246 for (; c < maxC; c++)
4248 if (a < alsu + alength)
4249 { /* still in A */
4250 carry += *a;
4251 a++;
4253 else
4254 { /* inside B */
4255 carry += ((eInt) * b) * m;
4256 b++;
4258 /* here carry is new Unit of digits; it could be +ve or -ve and */
4259 /* magnitude up to DECDPUNMAX squared */
4260 if ((ueInt) carry <= DECDPUNMAX)
4261 { /* fastpath 0-DECDPUNMAX */
4262 *c = (Unit) carry;
4263 carry = 0;
4264 continue;
4266 /* result for this unit is negative or >DECDPUNMAX */
4267 #if DECDPUN==4 /* use divide-by-multiply */
4268 /* remainder is undefined if negative, so we must test */
4269 if (carry >= 0)
4271 est = (((ueInt) carry >> 11) * 53687) >> 18;
4272 *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
4273 carry = est; /* likely quotient [79.7%] */
4274 if (*c < DECDPUNMAX + 1)
4275 continue; /* estimate was correct */
4276 carry++;
4277 *c -= DECDPUNMAX + 1;
4278 continue;
4280 /* negative case */
4281 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4282 est = (((ueInt) carry >> 11) * 53687) >> 18;
4283 *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4284 carry = est - (DECDPUNMAX + 1); /* correctly negative */
4285 if (*c < DECDPUNMAX + 1)
4286 continue; /* was OK */
4287 carry++;
4288 *c -= DECDPUNMAX + 1;
4289 #else
4290 if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4291 { /* fastpath carry 1 */
4292 *c = (Unit) (carry - (DECDPUNMAX + 1));
4293 carry = 1;
4294 continue;
4296 /* remainder is undefined if negative, so we must test */
4297 if (carry >= 0)
4299 *c = (Unit) (carry % (DECDPUNMAX + 1));
4300 carry = carry / (DECDPUNMAX + 1);
4301 continue;
4303 /* negative case */
4304 carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1); /* make positive */
4305 *c = (Unit) (carry % (DECDPUNMAX + 1));
4306 carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4307 #endif
4308 } /* c */
4310 /* OK, all A and B processed; might still have carry or borrow */
4311 /* return number of Units in the result, negated if a borrow */
4312 if (carry == 0)
4313 return c - clsu; /* no carry, we're done */
4314 if (carry > 0)
4315 { /* positive carry */
4316 *c = (Unit) carry; /* place as new unit */
4317 c++; /* .. */
4318 return c - clsu;
4320 /* -ve carry: it's a borrow; complement needed */
4321 add = 1; /* temporary carry... */
4322 for (c = clsu; c < maxC; c++)
4324 add = DECDPUNMAX + add - *c;
4325 if (add <= DECDPUNMAX)
4327 *c = (Unit) add;
4328 add = 0;
4330 else
4332 *c = 0;
4333 add = 1;
4336 /* add an extra unit iff it would be non-zero */
4337 #if DECTRACE
4338 printf ("UAS borrow: add %d, carry %d\n", add, carry);
4339 #endif
4340 if ((add - carry - 1) != 0)
4342 *c = (Unit) (add - carry - 1);
4343 c++; /* interesting, include it */
4345 return clsu - c; /* -ve result indicates borrowed */
4348 /* ------------------------------------------------------------------ */
4349 /* decTrim -- trim trailing zeros or normalize */
4350 /* */
4351 /* dn is the number to trim or normalize */
4352 /* all is 1 to remove all trailing zeros, 0 for just fraction ones */
4353 /* dropped returns the number of discarded trailing zeros */
4354 /* returns dn */
4355 /* */
4356 /* All fields are updated as required. This is a utility operation, */
4357 /* so special values are unchanged and no error is possible. */
4358 /* ------------------------------------------------------------------ */
4359 static decNumber *
4360 decTrim (decNumber * dn, Flag all, Int * dropped)
4362 Int d, exp; /* work */
4363 uInt cut; /* .. */
4364 Unit *up; /* -> current Unit */
4366 #if DECCHECK
4367 if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
4368 return dn;
4369 #endif
4371 *dropped = 0; /* assume no zeros dropped */
4372 if ((dn->bits & DECSPECIAL) /* fast exit if special .. */
4373 || (*dn->lsu & 0x01))
4374 return dn; /* .. or odd */
4375 if (ISZERO (dn))
4376 { /* .. or 0 */
4377 dn->exponent = 0; /* (sign is preserved) */
4378 return dn;
4381 /* we have a finite number which is even */
4382 exp = dn->exponent;
4383 cut = 1; /* digit (1-DECDPUN) in Unit */
4384 up = dn->lsu; /* -> current Unit */
4385 for (d = 0; d < dn->digits - 1; d++)
4386 { /* [don't strip the final digit] */
4387 /* slice by powers */
4388 #if DECDPUN<=4
4389 uInt quot = QUOT10 (*up, cut);
4390 if ((*up - quot * powers[cut]) != 0)
4391 break; /* found non-0 digit */
4392 #else
4393 if (*up % powers[cut] != 0)
4394 break; /* found non-0 digit */
4395 #endif
4396 /* have a trailing 0 */
4397 if (!all)
4398 { /* trimming */
4399 /* [if exp>0 then all trailing 0s are significant for trim] */
4400 if (exp <= 0)
4401 { /* if digit might be significant */
4402 if (exp == 0)
4403 break; /* then quit */
4404 exp++; /* next digit might be significant */
4407 cut++; /* next power */
4408 if (cut > DECDPUN)
4409 { /* need new Unit */
4410 up++;
4411 cut = 1;
4413 } /* d */
4414 if (d == 0)
4415 return dn; /* none dropped */
4417 /* effect the drop */
4418 decShiftToLeast (dn->lsu, D2U (dn->digits), d);
4419 dn->exponent += d; /* maintain numerical value */
4420 dn->digits -= d; /* new length */
4421 *dropped = d; /* report the count */
4422 return dn;
4425 /* ------------------------------------------------------------------ */
4426 /* decShiftToMost -- shift digits in array towards most significant */
4427 /* */
4428 /* uar is the array */
4429 /* digits is the count of digits in use in the array */
4430 /* shift is the number of zeros to pad with (least significant); */
4431 /* it must be zero or positive */
4432 /* */
4433 /* returns the new length of the integer in the array, in digits */
4434 /* */
4435 /* No overflow is permitted (that is, the uar array must be known to */
4436 /* be large enough to hold the result, after shifting). */
4437 /* ------------------------------------------------------------------ */
4438 static Int
4439 decShiftToMost (Unit * uar, Int digits, Int shift)
4441 Unit *target, *source, *first; /* work */
4442 uInt rem; /* for division */
4443 Int cut; /* odd 0's to add */
4444 uInt next; /* work */
4446 if (shift == 0)
4447 return digits; /* [fastpath] nothing to do */
4448 if ((digits + shift) <= DECDPUN)
4449 { /* [fastpath] single-unit case */
4450 *uar = (Unit) (*uar * powers[shift]);
4451 return digits + shift;
4454 cut = (DECDPUN - shift % DECDPUN) % DECDPUN;
4455 source = uar + D2U (digits) - 1; /* where msu comes from */
4456 first = uar + D2U (digits + shift) - 1; /* where msu of source will end up */
4457 target = source + D2U (shift); /* where upper part of first cut goes */
4458 next = 0;
4460 for (; source >= uar; source--, target--)
4462 /* split the source Unit and accumulate remainder for next */
4463 #if DECDPUN<=4
4464 uInt quot = QUOT10 (*source, cut);
4465 rem = *source - quot * powers[cut];
4466 next += quot;
4467 #else
4468 rem = *source % powers[cut];
4469 next += *source / powers[cut];
4470 #endif
4471 if (target <= first)
4472 *target = (Unit) next; /* write to target iff valid */
4473 next = rem * powers[DECDPUN - cut]; /* save remainder for next Unit */
4475 /* propagate to one below and clear the rest */
4476 for (; target >= uar; target--)
4478 *target = (Unit) next;
4479 next = 0;
4481 return digits + shift;
4484 /* ------------------------------------------------------------------ */
4485 /* decShiftToLeast -- shift digits in array towards least significant */
4486 /* */
4487 /* uar is the array */
4488 /* units is length of the array, in units */
4489 /* shift is the number of digits to remove from the lsu end; it */
4490 /* must be zero or positive and less than units*DECDPUN. */
4491 /* */
4492 /* returns the new length of the integer in the array, in units */
4493 /* */
4494 /* Removed digits are discarded (lost). Units not required to hold */
4495 /* the final result are unchanged. */
4496 /* ------------------------------------------------------------------ */
4497 static Int
4498 decShiftToLeast (Unit * uar, Int units, Int shift)
4500 Unit *target, *up; /* work */
4501 Int cut, count; /* work */
4502 Int quot, rem; /* for division */
4504 if (shift == 0)
4505 return units; /* [fastpath] nothing to do */
4507 up = uar + shift / DECDPUN; /* source; allow for whole Units */
4508 cut = shift % DECDPUN; /* odd 0's to drop */
4509 target = uar; /* both paths */
4510 if (cut == 0)
4511 { /* whole units shift */
4512 for (; up < uar + units; target++, up++)
4513 *target = *up;
4514 return target - uar;
4516 /* messier */
4517 count = units * DECDPUN - shift; /* the maximum new length */
4518 #if DECDPUN<=4
4519 quot = QUOT10 (*up, cut);
4520 #else
4521 quot = *up / powers[cut];
4522 #endif
4523 for (;; target++)
4525 *target = (Unit) quot;
4526 count -= (DECDPUN - cut);
4527 if (count <= 0)
4528 break;
4529 up++;
4530 quot = *up;
4531 #if DECDPUN<=4
4532 quot = QUOT10 (quot, cut);
4533 rem = *up - quot * powers[cut];
4534 #else
4535 rem = quot % powers[cut];
4536 quot = quot / powers[cut];
4537 #endif
4538 *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4539 count -= cut;
4540 if (count <= 0)
4541 break;
4543 return target - uar + 1;
4546 #if DECSUBSET
4547 /* ------------------------------------------------------------------ */
4548 /* decRoundOperand -- round an operand [used for subset only] */
4549 /* */
4550 /* dn is the number to round (dn->digits is > set->digits) */
4551 /* set is the relevant context */
4552 /* status is the status accumulator */
4553 /* */
4554 /* returns an allocated decNumber with the rounded result. */
4555 /* */
4556 /* lostDigits and other status may be set by this. */
4557 /* */
4558 /* Since the input is an operand, we are not permitted to modify it. */
4559 /* We therefore return an allocated decNumber, rounded as required. */
4560 /* It is the caller's responsibility to free the allocated storage. */
4561 /* */
4562 /* If no storage is available then the result cannot be used, so NULL */
4563 /* is returned. */
4564 /* ------------------------------------------------------------------ */
4565 static decNumber *
4566 decRoundOperand (const decNumber * dn, decContext * set, uInt * status)
4568 decNumber *res; /* result structure */
4569 uInt newstatus = 0; /* status from round */
4570 Int residue = 0; /* rounding accumulator */
4572 /* Allocate storage for the returned decNumber, big enough for the */
4573 /* length specified by the context */
4574 res = (decNumber *) malloc (sizeof (decNumber)
4575 + (D2U (set->digits) - 1) * sizeof (Unit));
4576 if (res == NULL)
4578 *status |= DEC_Insufficient_storage;
4579 return NULL;
4581 decCopyFit (res, dn, set, &residue, &newstatus);
4582 decApplyRound (res, set, residue, &newstatus);
4584 /* If that set Inexact then we "lost digits" */
4585 if (newstatus & DEC_Inexact)
4586 newstatus |= DEC_Lost_digits;
4587 *status |= newstatus;
4588 return res;
4590 #endif
4592 /* ------------------------------------------------------------------ */
4593 /* decCopyFit -- copy a number, shortening the coefficient if needed */
4594 /* */
4595 /* dest is the target decNumber */
4596 /* src is the source decNumber */
4597 /* set is the context [used for length (digits) and rounding mode] */
4598 /* residue is the residue accumulator */
4599 /* status contains the current status to be updated */
4600 /* */
4601 /* (dest==src is allowed and will be a no-op if fits) */
4602 /* All fields are updated as required. */
4603 /* ------------------------------------------------------------------ */
4604 static void
4605 decCopyFit (decNumber * dest, const decNumber * src, decContext * set,
4606 Int * residue, uInt * status)
4608 dest->bits = src->bits;
4609 dest->exponent = src->exponent;
4610 decSetCoeff (dest, set, src->lsu, src->digits, residue, status);
4613 /* ------------------------------------------------------------------ */
4614 /* decSetCoeff -- set the coefficient of a number */
4615 /* */
4616 /* dn is the number whose coefficient array is to be set. */
4617 /* It must have space for set->digits digits */
4618 /* set is the context [for size] */
4619 /* lsu -> lsu of the source coefficient [may be dn->lsu] */
4620 /* len is digits in the source coefficient [may be dn->digits] */
4621 /* residue is the residue accumulator. This has values as in */
4622 /* decApplyRound, and will be unchanged unless the */
4623 /* target size is less than len. In this case, the */
4624 /* coefficient is truncated and the residue is updated to */
4625 /* reflect the previous residue and the dropped digits. */
4626 /* status is the status accumulator, as usual */
4627 /* */
4628 /* The coefficient may already be in the number, or it can be an */
4629 /* external intermediate array. If it is in the number, lsu must == */
4630 /* dn->lsu and len must == dn->digits. */
4631 /* */
4632 /* Note that the coefficient length (len) may be < set->digits, and */
4633 /* in this case this merely copies the coefficient (or is a no-op */
4634 /* if dn->lsu==lsu). */
4635 /* */
4636 /* Note also that (only internally, from decNumberRescale and */
4637 /* decSetSubnormal) the value of set->digits may be less than one, */
4638 /* indicating a round to left. */
4639 /* This routine handles that case correctly; caller ensures space. */
4640 /* */
4641 /* dn->digits, dn->lsu (and as required), and dn->exponent are */
4642 /* updated as necessary. dn->bits (sign) is unchanged. */
4643 /* */
4644 /* DEC_Rounded status is set if any digits are discarded. */
4645 /* DEC_Inexact status is set if any non-zero digits are discarded, or */
4646 /* incoming residue was non-0 (implies rounded) */
4647 /* ------------------------------------------------------------------ */
4648 /* mapping array: maps 0-9 to canonical residues, so that we can */
4649 /* adjust by a residue in range [-1, +1] and achieve correct rounding */
4650 /* 0 1 2 3 4 5 6 7 8 9 */
4651 static const uByte resmap[10] = { 0, 3, 3, 3, 3, 5, 7, 7, 7, 7 };
4652 static void
4653 decSetCoeff (decNumber * dn, decContext * set, const Unit * lsu,
4654 Int len, Int * residue, uInt * status)
4656 Int discard; /* number of digits to discard */
4657 uInt discard1; /* first discarded digit */
4658 uInt cut; /* cut point in Unit */
4659 uInt quot, rem; /* for divisions */
4660 Unit *target; /* work */
4661 const Unit *up; /* work */
4662 Int count; /* .. */
4663 #if DECDPUN<=4
4664 uInt temp; /* .. */
4665 #endif
4667 discard = len - set->digits; /* digits to discard */
4668 if (discard <= 0)
4669 { /* no digits are being discarded */
4670 if (dn->lsu != lsu)
4671 { /* copy needed */
4672 /* copy the coefficient array to the result number; no shift needed */
4673 up = lsu;
4674 for (target = dn->lsu; target < dn->lsu + D2U (len); target++, up++)
4676 *target = *up;
4678 dn->digits = len; /* set the new length */
4680 /* dn->exponent and residue are unchanged */
4681 if (*residue != 0)
4682 *status |= (DEC_Inexact | DEC_Rounded); /* record inexactitude */
4683 return;
4686 /* we have to discard some digits */
4687 *status |= DEC_Rounded; /* accumulate Rounded status */
4688 if (*residue > 1)
4689 *residue = 1; /* previous residue now to right, so -1 to +1 */
4691 if (discard > len)
4692 { /* everything, +1, is being discarded */
4693 /* guard digit is 0 */
4694 /* residue is all the number [NB could be all 0s] */
4695 if (*residue <= 0)
4696 for (up = lsu + D2U (len) - 1; up >= lsu; up--)
4698 if (*up != 0)
4699 { /* found a non-0 */
4700 *residue = 1;
4701 break; /* no need to check any others */
4704 if (*residue != 0)
4705 *status |= DEC_Inexact; /* record inexactitude */
4706 *dn->lsu = 0; /* coefficient will now be 0 */
4707 dn->digits = 1; /* .. */
4708 dn->exponent += discard; /* maintain numerical value */
4709 return;
4710 } /* total discard */
4712 /* partial discard [most common case] */
4713 /* here, at least the first (most significant) discarded digit exists */
4715 /* spin up the number, noting residue as we pass, until we get to */
4716 /* the Unit with the first discarded digit. When we get there, */
4717 /* extract it and remember where we're at */
4718 count = 0;
4719 for (up = lsu;; up++)
4721 count += DECDPUN;
4722 if (count >= discard)
4723 break; /* full ones all checked */
4724 if (*up != 0)
4725 *residue = 1;
4726 } /* up */
4728 /* here up -> Unit with discarded digit */
4729 cut = discard - (count - DECDPUN) - 1;
4730 if (cut == DECDPUN - 1)
4731 { /* discard digit is at top */
4732 #if DECDPUN<=4
4733 discard1 = QUOT10 (*up, DECDPUN - 1);
4734 rem = *up - discard1 * powers[DECDPUN - 1];
4735 #else
4736 rem = *up % powers[DECDPUN - 1];
4737 discard1 = *up / powers[DECDPUN - 1];
4738 #endif
4739 if (rem != 0)
4740 *residue = 1;
4741 up++; /* move to next */
4742 cut = 0; /* bottom digit of result */
4743 quot = 0; /* keep a certain compiler happy */
4745 else
4747 /* discard digit is in low digit(s), not top digit */
4748 if (cut == 0)
4749 quot = *up;
4750 else /* cut>0 */
4751 { /* it's not at bottom of Unit */
4752 #if DECDPUN<=4
4753 quot = QUOT10 (*up, cut);
4754 rem = *up - quot * powers[cut];
4755 #else
4756 rem = *up % powers[cut];
4757 quot = *up / powers[cut];
4758 #endif
4759 if (rem != 0)
4760 *residue = 1;
4762 /* discard digit is now at bottom of quot */
4763 #if DECDPUN<=4
4764 temp = (quot * 6554) >> 16; /* fast /10 */
4765 /* Vowels algorithm here not a win (9 instructions) */
4766 discard1 = quot - X10 (temp);
4767 quot = temp;
4768 #else
4769 discard1 = quot % 10;
4770 quot = quot / 10;
4771 #endif
4772 cut++; /* update cut */
4775 /* here: up -> Unit of the array with discarded digit */
4776 /* cut is the division point for each Unit */
4777 /* quot holds the uncut high-order digits for the current */
4778 /* Unit, unless cut==0 in which case it's still in *up */
4779 /* copy the coefficient array to the result number, shifting as we go */
4780 count = set->digits; /* digits to end up with */
4781 if (count <= 0)
4782 { /* special for Rescale/Subnormal :-( */
4783 *dn->lsu = 0; /* .. result is 0 */
4784 dn->digits = 1; /* .. */
4786 else
4787 { /* shift to least */
4788 /* [this is similar to decShiftToLeast code, with copy] */
4789 dn->digits = count; /* set the new length */
4790 if (cut == 0)
4792 /* on unit boundary, so simple shift down copy loop suffices */
4793 for (target = dn->lsu; target < dn->lsu + D2U (count);
4794 target++, up++)
4796 *target = *up;
4799 else
4800 for (target = dn->lsu;; target++)
4802 *target = (Unit) quot;
4803 count -= (DECDPUN - cut);
4804 if (count <= 0)
4805 break;
4806 up++;
4807 quot = *up;
4808 #if DECDPUN<=4
4809 quot = QUOT10 (quot, cut);
4810 rem = *up - quot * powers[cut];
4811 #else
4812 rem = quot % powers[cut];
4813 quot = quot / powers[cut];
4814 #endif
4815 *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4816 count -= cut;
4817 if (count <= 0)
4818 break;
4820 } /* shift to least needed */
4821 dn->exponent += discard; /* maintain numerical value */
4823 /* here, discard1 is the guard digit, and residue is everything else */
4824 /* [use mapping to accumulate residue safely] */
4825 *residue += resmap[discard1];
4827 if (*residue != 0)
4828 *status |= DEC_Inexact; /* record inexactitude */
4829 return;
4832 /* ------------------------------------------------------------------ */
4833 /* decApplyRound -- apply pending rounding to a number */
4834 /* */
4835 /* dn is the number, with space for set->digits digits */
4836 /* set is the context [for size and rounding mode] */
4837 /* residue indicates pending rounding, being any accumulated */
4838 /* guard and sticky information. It may be: */
4839 /* 6-9: rounding digit is >5 */
4840 /* 5: rounding digit is exactly half-way */
4841 /* 1-4: rounding digit is <5 and >0 */
4842 /* 0: the coefficient is exact */
4843 /* -1: as 1, but the hidden digits are subtractive, that */
4844 /* is, of the opposite sign to dn. In this case the */
4845 /* coefficient must be non-0. */
4846 /* status is the status accumulator, as usual */
4847 /* */
4848 /* This routine applies rounding while keeping the length of the */
4849 /* coefficient constant. The exponent and status are unchanged */
4850 /* except if: */
4851 /* */
4852 /* -- the coefficient was increased and is all nines (in which */
4853 /* case Overflow could occur, and is handled directly here so */
4854 /* the caller does not need to re-test for overflow) */
4855 /* */
4856 /* -- the coefficient was decreased and becomes all nines (in which */
4857 /* case Underflow could occur, and is also handled directly). */
4858 /* */
4859 /* All fields in dn are updated as required. */
4860 /* */
4861 /* ------------------------------------------------------------------ */
4862 static void
4863 decApplyRound (decNumber * dn, decContext * set, Int residue, uInt * status)
4865 Int bump; /* 1 if coefficient needs to be incremented */
4866 /* -1 if coefficient needs to be decremented */
4868 if (residue == 0)
4869 return; /* nothing to apply */
4871 bump = 0; /* assume a smooth ride */
4873 /* now decide whether, and how, to round, depending on mode */
4874 switch (set->round)
4876 case DEC_ROUND_DOWN:
4878 /* no change, except if negative residue */
4879 if (residue < 0)
4880 bump = -1;
4881 break;
4882 } /* r-d */
4884 case DEC_ROUND_HALF_DOWN:
4886 if (residue > 5)
4887 bump = 1;
4888 break;
4889 } /* r-h-d */
4891 case DEC_ROUND_HALF_EVEN:
4893 if (residue > 5)
4894 bump = 1; /* >0.5 goes up */
4895 else if (residue == 5)
4896 { /* exactly 0.5000... */
4897 /* 0.5 goes up iff [new] lsd is odd */
4898 if (*dn->lsu & 0x01)
4899 bump = 1;
4901 break;
4902 } /* r-h-e */
4904 case DEC_ROUND_HALF_UP:
4906 if (residue >= 5)
4907 bump = 1;
4908 break;
4909 } /* r-h-u */
4911 case DEC_ROUND_UP:
4913 if (residue > 0)
4914 bump = 1;
4915 break;
4916 } /* r-u */
4918 case DEC_ROUND_CEILING:
4920 /* same as _UP for positive numbers, and as _DOWN for negatives */
4921 /* [negative residue cannot occur on 0] */
4922 if (decNumberIsNegative (dn))
4924 if (residue < 0)
4925 bump = -1;
4927 else
4929 if (residue > 0)
4930 bump = 1;
4932 break;
4933 } /* r-c */
4935 case DEC_ROUND_FLOOR:
4937 /* same as _UP for negative numbers, and as _DOWN for positive */
4938 /* [negative residue cannot occur on 0] */
4939 if (!decNumberIsNegative (dn))
4941 if (residue < 0)
4942 bump = -1;
4944 else
4946 if (residue > 0)
4947 bump = 1;
4949 break;
4950 } /* r-f */
4952 default:
4953 { /* e.g., DEC_ROUND_MAX */
4954 *status |= DEC_Invalid_context;
4955 #if DECTRACE
4956 printf ("Unknown rounding mode: %d\n", set->round);
4957 #endif
4958 break;
4960 } /* switch */
4962 /* now bump the number, up or down, if need be */
4963 if (bump == 0)
4964 return; /* no action required */
4966 /* Simply use decUnitAddSub unless we are bumping up and the number */
4967 /* is all nines. In this special case we set to 1000... and adjust */
4968 /* the exponent by one (as otherwise we could overflow the array) */
4969 /* Similarly handle all-nines result if bumping down. */
4970 if (bump > 0)
4972 Unit *up; /* work */
4973 uInt count = dn->digits; /* digits to be checked */
4974 for (up = dn->lsu;; up++)
4976 if (count <= DECDPUN)
4978 /* this is the last Unit (the msu) */
4979 if (*up != powers[count] - 1)
4980 break; /* not still 9s */
4981 /* here if it, too, is all nines */
4982 *up = (Unit) powers[count - 1]; /* here 999 -> 100 etc. */
4983 for (up = up - 1; up >= dn->lsu; up--)
4984 *up = 0; /* others all to 0 */
4985 dn->exponent++; /* and bump exponent */
4986 /* [which, very rarely, could cause Overflow...] */
4987 if ((dn->exponent + dn->digits) > set->emax + 1)
4989 decSetOverflow (dn, set, status);
4991 return; /* done */
4993 /* a full unit to check, with more to come */
4994 if (*up != DECDPUNMAX)
4995 break; /* not still 9s */
4996 count -= DECDPUN;
4997 } /* up */
4998 } /* bump>0 */
4999 else
5000 { /* -1 */
5001 /* here we are lookng for a pre-bump of 1000... (leading 1, */
5002 /* all other digits zero) */
5003 Unit *up, *sup; /* work */
5004 uInt count = dn->digits; /* digits to be checked */
5005 for (up = dn->lsu;; up++)
5007 if (count <= DECDPUN)
5009 /* this is the last Unit (the msu) */
5010 if (*up != powers[count - 1])
5011 break; /* not 100.. */
5012 /* here if we have the 1000... case */
5013 sup = up; /* save msu pointer */
5014 *up = (Unit) powers[count] - 1; /* here 100 in msu -> 999 */
5015 /* others all to all-nines, too */
5016 for (up = up - 1; up >= dn->lsu; up--)
5017 *up = (Unit) powers[DECDPUN] - 1;
5018 dn->exponent--; /* and bump exponent */
5020 /* iff the number was at the subnormal boundary (exponent=etiny) */
5021 /* then the exponent is now out of range, so it will in fact get */
5022 /* clamped to etiny and the final 9 dropped. */
5023 /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
5024 /* dn->exponent, set->digits); */
5025 if (dn->exponent + 1 == set->emin - set->digits + 1)
5027 if (count == 1 && dn->digits == 1)
5028 *sup = 0; /* here 9 -> 0[.9] */
5029 else
5031 *sup = (Unit) powers[count - 1] - 1; /* here 999.. in msu -> 99.. */
5032 dn->digits--;
5034 dn->exponent++;
5035 *status |=
5036 DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5038 return; /* done */
5041 /* a full unit to check, with more to come */
5042 if (*up != 0)
5043 break; /* not still 0s */
5044 count -= DECDPUN;
5045 } /* up */
5047 } /* bump<0 */
5049 /* Actual bump needed. Do it. */
5050 decUnitAddSub (dn->lsu, D2U (dn->digits), one, 1, 0, dn->lsu, bump);
5053 #if DECSUBSET
5054 /* ------------------------------------------------------------------ */
5055 /* decFinish -- finish processing a number */
5056 /* */
5057 /* dn is the number */
5058 /* set is the context */
5059 /* residue is the rounding accumulator (as in decApplyRound) */
5060 /* status is the accumulator */
5061 /* */
5062 /* This finishes off the current number by: */
5063 /* 1. If not extended: */
5064 /* a. Converting a zero result to clean '0' */
5065 /* b. Reducing positive exponents to 0, if would fit in digits */
5066 /* 2. Checking for overflow and subnormals (always) */
5067 /* Note this is just Finalize when no subset arithmetic. */
5068 /* All fields are updated as required. */
5069 /* ------------------------------------------------------------------ */
5070 static void
5071 decFinish (decNumber * dn, decContext * set, Int * residue, uInt * status)
5073 if (!set->extended)
5075 if ISZERO
5076 (dn)
5077 { /* value is zero */
5078 dn->exponent = 0; /* clean exponent .. */
5079 dn->bits = 0; /* .. and sign */
5080 return; /* no error possible */
5082 if (dn->exponent >= 0)
5083 { /* non-negative exponent */
5084 /* >0; reduce to integer if possible */
5085 if (set->digits >= (dn->exponent + dn->digits))
5087 dn->digits = decShiftToMost (dn->lsu, dn->digits, dn->exponent);
5088 dn->exponent = 0;
5091 } /* !extended */
5093 decFinalize (dn, set, residue, status);
5095 #endif
5097 /* ------------------------------------------------------------------ */
5098 /* decFinalize -- final check, clamp, and round of a number */
5099 /* */
5100 /* dn is the number */
5101 /* set is the context */
5102 /* residue is the rounding accumulator (as in decApplyRound) */
5103 /* status is the status accumulator */
5104 /* */
5105 /* This finishes off the current number by checking for subnormal */
5106 /* results, applying any pending rounding, checking for overflow, */
5107 /* and applying any clamping. */
5108 /* Underflow and overflow conditions are raised as appropriate. */
5109 /* All fields are updated as required. */
5110 /* ------------------------------------------------------------------ */
5111 static void
5112 decFinalize (decNumber * dn, decContext * set, Int * residue, uInt * status)
5114 Int shift; /* shift needed if clamping */
5116 /* We have to be careful when checking the exponent as the adjusted */
5117 /* exponent could overflow 31 bits [because it may already be up */
5118 /* to twice the expected]. */
5120 /* First test for subnormal. This must be done before any final */
5121 /* round as the result could be rounded to Nmin or 0. */
5122 if (dn->exponent < 0 /* negative exponent */
5123 && (dn->exponent < set->emin - dn->digits + 1))
5125 /* Go handle subnormals; this will apply round if needed. */
5126 decSetSubnormal (dn, set, residue, status);
5127 return;
5130 /* now apply any pending round (this could raise overflow). */
5131 if (*residue != 0)
5132 decApplyRound (dn, set, *residue, status);
5134 /* Check for overflow [redundant in the 'rare' case] or clamp */
5135 if (dn->exponent <= set->emax - set->digits + 1)
5136 return; /* neither needed */
5138 /* here when we might have an overflow or clamp to do */
5139 if (dn->exponent > set->emax - dn->digits + 1)
5140 { /* too big */
5141 decSetOverflow (dn, set, status);
5142 return;
5144 /* here when the result is normal but in clamp range */
5145 if (!set->clamp)
5146 return;
5148 /* here when we need to apply the IEEE exponent clamp (fold-down) */
5149 shift = dn->exponent - (set->emax - set->digits + 1);
5151 /* shift coefficient (if non-zero) */
5152 if (!ISZERO (dn))
5154 dn->digits = decShiftToMost (dn->lsu, dn->digits, shift);
5156 dn->exponent -= shift; /* adjust the exponent to match */
5157 *status |= DEC_Clamped; /* and record the dirty deed */
5158 return;
5161 /* ------------------------------------------------------------------ */
5162 /* decSetOverflow -- set number to proper overflow value */
5163 /* */
5164 /* dn is the number (used for sign [only] and result) */
5165 /* set is the context [used for the rounding mode] */
5166 /* status contains the current status to be updated */
5167 /* */
5168 /* This sets the sign of a number and sets its value to either */
5169 /* Infinity or the maximum finite value, depending on the sign of */
5170 /* dn and therounding mode, following IEEE 854 rules. */
5171 /* ------------------------------------------------------------------ */
5172 static void
5173 decSetOverflow (decNumber * dn, decContext * set, uInt * status)
5175 Flag needmax = 0; /* result is maximum finite value */
5176 uByte sign = dn->bits & DECNEG; /* clean and save sign bit */
5178 if (ISZERO (dn))
5179 { /* zero does not overflow magnitude */
5180 Int emax = set->emax; /* limit value */
5181 if (set->clamp)
5182 emax -= set->digits - 1; /* lower if clamping */
5183 if (dn->exponent > emax)
5184 { /* clamp required */
5185 dn->exponent = emax;
5186 *status |= DEC_Clamped;
5188 return;
5191 decNumberZero (dn);
5192 switch (set->round)
5194 case DEC_ROUND_DOWN:
5196 needmax = 1; /* never Infinity */
5197 break;
5198 } /* r-d */
5199 case DEC_ROUND_CEILING:
5201 if (sign)
5202 needmax = 1; /* Infinity if non-negative */
5203 break;
5204 } /* r-c */
5205 case DEC_ROUND_FLOOR:
5207 if (!sign)
5208 needmax = 1; /* Infinity if negative */
5209 break;
5210 } /* r-f */
5211 default:
5212 break; /* Infinity in all other cases */
5214 if (needmax)
5216 Unit *up; /* work */
5217 Int count = set->digits; /* nines to add */
5218 dn->digits = count;
5219 /* fill in all nines to set maximum value */
5220 for (up = dn->lsu;; up++)
5222 if (count > DECDPUN)
5223 *up = DECDPUNMAX; /* unit full o'nines */
5224 else
5225 { /* this is the msu */
5226 *up = (Unit) (powers[count] - 1);
5227 break;
5229 count -= DECDPUN; /* we filled those digits */
5230 } /* up */
5231 dn->bits = sign; /* sign */
5232 dn->exponent = set->emax - set->digits + 1;
5234 else
5235 dn->bits = sign | DECINF; /* Value is +/-Infinity */
5236 *status |= DEC_Overflow | DEC_Inexact | DEC_Rounded;
5239 /* ------------------------------------------------------------------ */
5240 /* decSetSubnormal -- process value whose exponent is <Emin */
5241 /* */
5242 /* dn is the number (used as input as well as output; it may have */
5243 /* an allowed subnormal value, which may need to be rounded) */
5244 /* set is the context [used for the rounding mode] */
5245 /* residue is any pending residue */
5246 /* status contains the current status to be updated */
5247 /* */
5248 /* If subset mode, set result to zero and set Underflow flags. */
5249 /* */
5250 /* Value may be zero with a low exponent; this does not set Subnormal */
5251 /* but the exponent will be clamped to Etiny. */
5252 /* */
5253 /* Otherwise ensure exponent is not out of range, and round as */
5254 /* necessary. Underflow is set if the result is Inexact. */
5255 /* ------------------------------------------------------------------ */
5256 static void
5257 decSetSubnormal (decNumber * dn, decContext * set,
5258 Int * residue, uInt * status)
5260 decContext workset; /* work */
5261 Int etiny, adjust; /* .. */
5263 #if DECSUBSET
5264 /* simple set to zero and 'hard underflow' for subset */
5265 if (!set->extended)
5267 decNumberZero (dn);
5268 /* always full overflow */
5269 *status |= DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5270 return;
5272 #endif
5274 /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
5275 /* (Etiny) if needed */
5276 etiny = set->emin - (set->digits - 1); /* smallest allowed exponent */
5278 if ISZERO
5279 (dn)
5280 { /* value is zero */
5281 /* residue can never be non-zero here */
5282 #if DECCHECK
5283 if (*residue != 0)
5285 printf ("++ Subnormal 0 residue %d\n", *residue);
5286 *status |= DEC_Invalid_operation;
5288 #endif
5289 if (dn->exponent < etiny)
5290 { /* clamp required */
5291 dn->exponent = etiny;
5292 *status |= DEC_Clamped;
5294 return;
5297 *status |= DEC_Subnormal; /* we have a non-zero subnormal */
5299 adjust = etiny - dn->exponent; /* calculate digits to remove */
5300 if (adjust <= 0)
5301 { /* not out of range; unrounded */
5302 /* residue can never be non-zero here, so fast-path out */
5303 #if DECCHECK
5304 if (*residue != 0)
5306 printf ("++ Subnormal no-adjust residue %d\n", *residue);
5307 *status |= DEC_Invalid_operation;
5309 #endif
5310 /* it may already be inexact (from setting the coefficient) */
5311 if (*status & DEC_Inexact)
5312 *status |= DEC_Underflow;
5313 return;
5316 /* adjust>0. we need to rescale the result so exponent becomes Etiny */
5317 /* [this code is similar to that in rescale] */
5318 workset = *set; /* clone rounding, etc. */
5319 workset.digits = dn->digits - adjust; /* set requested length */
5320 workset.emin -= adjust; /* and adjust emin to match */
5321 /* [note that the latter can be <1, here, similar to Rescale case] */
5322 decSetCoeff (dn, &workset, dn->lsu, dn->digits, residue, status);
5323 decApplyRound (dn, &workset, *residue, status);
5325 /* Use 754R/854 default rule: Underflow is set iff Inexact */
5326 /* [independent of whether trapped] */
5327 if (*status & DEC_Inexact)
5328 *status |= DEC_Underflow;
5330 /* if we rounded up a 999s case, exponent will be off by one; adjust */
5331 /* back if so [it will fit, because we shortened] */
5332 if (dn->exponent > etiny)
5334 dn->digits = decShiftToMost (dn->lsu, dn->digits, 1);
5335 dn->exponent--; /* (re)adjust the exponent. */
5339 /* ------------------------------------------------------------------ */
5340 /* decGetInt -- get integer from a number */
5341 /* */
5342 /* dn is the number [which will not be altered] */
5343 /* set is the context [requested digits], subset only */
5344 /* returns the converted integer, or BADINT if error */
5345 /* */
5346 /* This checks and gets a whole number from the input decNumber. */
5347 /* The magnitude of the integer must be <2^31. */
5348 /* Any discarded fractional part must be 0. */
5349 /* If subset it must also fit in set->digits */
5350 /* ------------------------------------------------------------------ */
5351 #if DECSUBSET
5352 static Int
5353 decGetInt (const decNumber * dn, decContext * set)
5355 #else
5356 static Int
5357 decGetInt (const decNumber * dn)
5359 #endif
5360 Int theInt; /* result accumulator */
5361 const Unit *up; /* work */
5362 Int got; /* digits (real or not) processed */
5363 Int ilength = dn->digits + dn->exponent; /* integral length */
5365 /* The number must be an integer that fits in 10 digits */
5366 /* Assert, here, that 10 is enough for any rescale Etiny */
5367 #if DEC_MAX_EMAX > 999999999
5368 #error GetInt may need updating [for Emax]
5369 #endif
5370 #if DEC_MIN_EMIN < -999999999
5371 #error GetInt may need updating [for Emin]
5372 #endif
5373 if (ISZERO (dn))
5374 return 0; /* zeros are OK, with any exponent */
5375 if (ilength > 10)
5376 return BADINT; /* always too big */
5377 #if DECSUBSET
5378 if (!set->extended && ilength > set->digits)
5379 return BADINT;
5380 #endif
5382 up = dn->lsu; /* ready for lsu */
5383 theInt = 0; /* ready to accumulate */
5384 if (dn->exponent >= 0)
5385 { /* relatively easy */
5386 /* no fractional part [usual]; allow for positive exponent */
5387 got = dn->exponent;
5389 else
5390 { /* -ve exponent; some fractional part to check and discard */
5391 Int count = -dn->exponent; /* digits to discard */
5392 /* spin up whole units until we get to the Unit with the unit digit */
5393 for (; count >= DECDPUN; up++)
5395 if (*up != 0)
5396 return BADINT; /* non-zero Unit to discard */
5397 count -= DECDPUN;
5399 if (count == 0)
5400 got = 0; /* [a multiple of DECDPUN] */
5401 else
5402 { /* [not multiple of DECDPUN] */
5403 Int rem; /* work */
5404 /* slice off fraction digits and check for non-zero */
5405 #if DECDPUN<=4
5406 theInt = QUOT10 (*up, count);
5407 rem = *up - theInt * powers[count];
5408 #else
5409 rem = *up % powers[count]; /* slice off discards */
5410 theInt = *up / powers[count];
5411 #endif
5412 if (rem != 0)
5413 return BADINT; /* non-zero fraction */
5414 /* OK, we're good */
5415 got = DECDPUN - count; /* number of digits so far */
5416 up++; /* ready for next */
5419 /* collect the rest */
5420 for (; got < ilength; up++)
5422 theInt += *up * powers[got];
5423 got += DECDPUN;
5425 if ((ilength == 10) /* check no wrap */
5426 && (theInt / (Int) powers[got - DECDPUN] != *(up - 1)))
5427 return BADINT;
5428 /* [that test also disallows the BADINT result case] */
5430 /* apply any sign and return */
5431 if (decNumberIsNegative (dn))
5432 theInt = -theInt;
5433 return theInt;
5436 /* ------------------------------------------------------------------ */
5437 /* decStrEq -- caseless comparison of strings */
5438 /* */
5439 /* str1 is one of the strings to compare */
5440 /* str2 is the other */
5441 /* */
5442 /* returns 1 if strings caseless-compare equal, 0 otherwise */
5443 /* */
5444 /* Note that the strings must be the same length if they are to */
5445 /* compare equal; there is no padding. */
5446 /* ------------------------------------------------------------------ */
5447 /* [strcmpi is not in ANSI C] */
5448 static Flag
5449 decStrEq (const char *str1, const char *str2)
5451 for (;; str1++, str2++)
5453 unsigned char u1 = (unsigned char) *str1;
5454 unsigned char u2 = (unsigned char) *str2;
5455 if (u1 == u2)
5457 if (u1 == '\0')
5458 break;
5460 else
5462 if (tolower (u1) != tolower (u2))
5463 return 0;
5465 } /* stepping */
5466 return 1;
5469 /* ------------------------------------------------------------------ */
5470 /* decNaNs -- handle NaN operand or operands */
5471 /* */
5472 /* res is the result number */
5473 /* lhs is the first operand */
5474 /* rhs is the second operand, or NULL if none */
5475 /* status contains the current status */
5476 /* returns res in case convenient */
5477 /* */
5478 /* Called when one or both operands is a NaN, and propagates the */
5479 /* appropriate result to res. When an sNaN is found, it is changed */
5480 /* to a qNaN and Invalid operation is set. */
5481 /* ------------------------------------------------------------------ */
5482 static decNumber *
5483 decNaNs (decNumber * res, const decNumber * lhs, const decNumber * rhs, uInt * status)
5485 /* This decision tree ends up with LHS being the source pointer, */
5486 /* and status updated if need be */
5487 if (lhs->bits & DECSNAN)
5488 *status |= DEC_Invalid_operation | DEC_sNaN;
5489 else if (rhs == NULL);
5490 else if (rhs->bits & DECSNAN)
5492 lhs = rhs;
5493 *status |= DEC_Invalid_operation | DEC_sNaN;
5495 else if (lhs->bits & DECNAN);
5496 else
5497 lhs = rhs;
5499 decNumberCopy (res, lhs);
5500 res->bits &= ~DECSNAN; /* convert any sNaN to NaN, while */
5501 res->bits |= DECNAN; /* .. preserving sign */
5502 res->exponent = 0; /* clean exponent */
5503 /* [coefficient was copied] */
5504 return res;
5507 /* ------------------------------------------------------------------ */
5508 /* decStatus -- apply non-zero status */
5509 /* */
5510 /* dn is the number to set if error */
5511 /* status contains the current status (not yet in context) */
5512 /* set is the context */
5513 /* */
5514 /* If the status is an error status, the number is set to a NaN, */
5515 /* unless the error was an overflow, divide-by-zero, or underflow, */
5516 /* in which case the number will have already been set. */
5517 /* */
5518 /* The context status is then updated with the new status. Note that */
5519 /* this may raise a signal, so control may never return from this */
5520 /* routine (hence resources must be recovered before it is called). */
5521 /* ------------------------------------------------------------------ */
5522 static void
5523 decStatus (decNumber * dn, uInt status, decContext * set)
5525 if (status & DEC_NaNs)
5526 { /* error status -> NaN */
5527 /* if cause was an sNaN, clear and propagate [NaN is already set up] */
5528 if (status & DEC_sNaN)
5529 status &= ~DEC_sNaN;
5530 else
5532 decNumberZero (dn); /* other error: clean throughout */
5533 dn->bits = DECNAN; /* and make a quiet NaN */
5536 decContextSetStatus (set, status);
5537 return;
5540 /* ------------------------------------------------------------------ */
5541 /* decGetDigits -- count digits in a Units array */
5542 /* */
5543 /* uar is the Unit array holding the number [this is often an */
5544 /* accumulator of some sort] */
5545 /* len is the length of the array in units */
5546 /* */
5547 /* returns the number of (significant) digits in the array */
5548 /* */
5549 /* All leading zeros are excluded, except the last if the array has */
5550 /* only zero Units. */
5551 /* ------------------------------------------------------------------ */
5552 /* This may be called twice during some operations. */
5553 static Int
5554 decGetDigits (const Unit * uar, Int len)
5556 const Unit *up = uar + len - 1; /* -> msu */
5557 Int digits = len * DECDPUN; /* maximum possible digits */
5558 uInt const *pow; /* work */
5560 for (; up >= uar; up--)
5562 digits -= DECDPUN;
5563 if (*up == 0)
5564 { /* unit is 0 */
5565 if (digits != 0)
5566 continue; /* more to check */
5567 /* all units were 0 */
5568 digits++; /* .. so bump digits to 1 */
5569 break;
5571 /* found the first non-zero Unit */
5572 digits++;
5573 if (*up < 10)
5574 break; /* fastpath 1-9 */
5575 digits++;
5576 for (pow = &powers[2]; *up >= *pow; pow++)
5577 digits++;
5578 break;
5579 } /* up */
5581 return digits;
5585 #if DECTRACE | DECCHECK
5586 /* ------------------------------------------------------------------ */
5587 /* decNumberShow -- display a number [debug aid] */
5588 /* dn is the number to show */
5589 /* */
5590 /* Shows: sign, exponent, coefficient (msu first), digits */
5591 /* or: sign, special-value */
5592 /* ------------------------------------------------------------------ */
5593 /* this is public so other modules can use it */
5594 void
5595 decNumberShow (const decNumber * dn)
5597 const Unit *up; /* work */
5598 uInt u, d; /* .. */
5599 Int cut; /* .. */
5600 char isign = '+'; /* main sign */
5601 if (dn == NULL)
5603 printf ("NULL\n");
5604 return;
5606 if (decNumberIsNegative (dn))
5607 isign = '-';
5608 printf (" >> %c ", isign);
5609 if (dn->bits & DECSPECIAL)
5610 { /* Is a special value */
5611 if (decNumberIsInfinite (dn))
5612 printf ("Infinity");
5613 else
5614 { /* a NaN */
5615 if (dn->bits & DECSNAN)
5616 printf ("sNaN"); /* signalling NaN */
5617 else
5618 printf ("NaN");
5620 /* if coefficient and exponent are 0, we're done */
5621 if (dn->exponent == 0 && dn->digits == 1 && *dn->lsu == 0)
5623 printf ("\n");
5624 return;
5626 /* drop through to report other information */
5627 printf (" ");
5630 /* now carefully display the coefficient */
5631 up = dn->lsu + D2U (dn->digits) - 1; /* msu */
5632 printf ("%d", *up);
5633 for (up = up - 1; up >= dn->lsu; up--)
5635 u = *up;
5636 printf (":");
5637 for (cut = DECDPUN - 1; cut >= 0; cut--)
5639 d = u / powers[cut];
5640 u -= d * powers[cut];
5641 printf ("%d", d);
5642 } /* cut */
5643 } /* up */
5644 if (dn->exponent != 0)
5646 char esign = '+';
5647 if (dn->exponent < 0)
5648 esign = '-';
5649 printf (" E%c%d", esign, abs (dn->exponent));
5651 printf (" [%d]\n", dn->digits);
5653 #endif
5655 #if DECTRACE || DECCHECK
5656 /* ------------------------------------------------------------------ */
5657 /* decDumpAr -- display a unit array [debug aid] */
5658 /* name is a single-character tag name */
5659 /* ar is the array to display */
5660 /* len is the length of the array in Units */
5661 /* ------------------------------------------------------------------ */
5662 static void
5663 decDumpAr (char name, const Unit * ar, Int len)
5665 Int i;
5666 #if DECDPUN==4
5667 const char *spec = "%04d ";
5668 #else
5669 const char *spec = "%d ";
5670 #endif
5671 printf (" :%c: ", name);
5672 for (i = len - 1; i >= 0; i--)
5674 if (i == len - 1)
5675 printf ("%d ", ar[i]);
5676 else
5677 printf (spec, ar[i]);
5679 printf ("\n");
5680 return;
5682 #endif
5684 #if DECCHECK
5685 /* ------------------------------------------------------------------ */
5686 /* decCheckOperands -- check operand(s) to a routine */
5687 /* res is the result structure (not checked; it will be set to */
5688 /* quiet NaN if error found (and it is not NULL)) */
5689 /* lhs is the first operand (may be DECUNUSED) */
5690 /* rhs is the second (may be DECUNUSED) */
5691 /* set is the context (may be DECUNUSED) */
5692 /* returns 0 if both operands, and the context are clean, or 1 */
5693 /* otherwise (in which case the context will show an error, */
5694 /* unless NULL). Note that res is not cleaned; caller should */
5695 /* handle this so res=NULL case is safe. */
5696 /* The caller is expected to abandon immediately if 1 is returned. */
5697 /* ------------------------------------------------------------------ */
5698 static Flag
5699 decCheckOperands (decNumber * res, const decNumber * lhs,
5700 const decNumber * rhs, decContext * set)
5702 Flag bad = 0;
5703 if (set == NULL)
5704 { /* oops; hopeless */
5705 #if DECTRACE
5706 printf ("Context is NULL.\n");
5707 #endif
5708 bad = 1;
5709 return 1;
5711 else if (set != DECUNUSED
5712 && (set->digits < 1 || set->round < 0
5713 || set->round >= DEC_ROUND_MAX))
5715 bad = 1;
5716 #if DECTRACE
5717 printf ("Bad context [digits=%d round=%d].\n", set->digits, set->round);
5718 #endif
5720 else
5722 if (res == NULL)
5724 bad = 1;
5725 #if DECTRACE
5726 printf ("Bad result [is NULL].\n");
5727 #endif
5729 if (!bad && lhs != DECUNUSED)
5730 bad = (decCheckNumber (lhs, set));
5731 if (!bad && rhs != DECUNUSED)
5732 bad = (decCheckNumber (rhs, set));
5734 if (bad)
5736 if (set != DECUNUSED)
5737 decContextSetStatus (set, DEC_Invalid_operation);
5738 if (res != DECUNUSED && res != NULL)
5740 decNumberZero (res);
5741 res->bits = DECNAN; /* qNaN */
5744 return bad;
5747 /* ------------------------------------------------------------------ */
5748 /* decCheckNumber -- check a number */
5749 /* dn is the number to check */
5750 /* set is the context (may be DECUNUSED) */
5751 /* returns 0 if the number is clean, or 1 otherwise */
5752 /* */
5753 /* The number is considered valid if it could be a result from some */
5754 /* operation in some valid context (not necessarily the current one). */
5755 /* ------------------------------------------------------------------ */
5756 Flag
5757 decCheckNumber (const decNumber * dn, decContext * set)
5759 const Unit *up; /* work */
5760 uInt maxuint; /* .. */
5761 Int ae, d, digits; /* .. */
5762 Int emin, emax; /* .. */
5764 if (dn == NULL)
5765 { /* hopeless */
5766 #if DECTRACE
5767 printf ("Reference to decNumber is NULL.\n");
5768 #endif
5769 return 1;
5772 /* check special values */
5773 if (dn->bits & DECSPECIAL)
5775 if (dn->exponent != 0)
5777 #if DECTRACE
5778 printf ("Exponent %d (not 0) for a special value.\n", dn->exponent);
5779 #endif
5780 return 1;
5783 /* 2003.09.08: NaNs may now have coefficients, so next tests Inf only */
5784 if (decNumberIsInfinite (dn))
5786 if (dn->digits != 1)
5788 #if DECTRACE
5789 printf ("Digits %d (not 1) for an infinity.\n", dn->digits);
5790 #endif
5791 return 1;
5793 if (*dn->lsu != 0)
5795 #if DECTRACE
5796 printf ("LSU %d (not 0) for an infinity.\n", *dn->lsu);
5797 #endif
5798 return 1;
5800 } /* Inf */
5801 /* 2002.12.26: negative NaNs can now appear through proposed IEEE */
5802 /* concrete formats (decimal64, etc.), though they are */
5803 /* never visible in strings. */
5804 return 0;
5806 /* if ((dn->bits & DECINF) || (dn->bits & DECNEG)==0) return 0; */
5807 /* #if DECTRACE */
5808 /* printf("Negative NaN in number.\n"); */
5809 /* #endif */
5810 /* return 1; */
5813 /* check the coefficient */
5814 if (dn->digits < 1 || dn->digits > DECNUMMAXP)
5816 #if DECTRACE
5817 printf ("Digits %d in number.\n", dn->digits);
5818 #endif
5819 return 1;
5822 d = dn->digits;
5824 for (up = dn->lsu; d > 0; up++)
5826 if (d > DECDPUN)
5827 maxuint = DECDPUNMAX;
5828 else
5829 { /* we are at the msu */
5830 maxuint = powers[d] - 1;
5831 if (dn->digits > 1 && *up < powers[d - 1])
5833 #if DECTRACE
5834 printf ("Leading 0 in number.\n");
5835 decNumberShow (dn);
5836 #endif
5837 return 1;
5840 if (*up > maxuint)
5842 #if DECTRACE
5843 printf ("Bad Unit [%08x] in number at offset %d [maxuint %d].\n",
5844 *up, up - dn->lsu, maxuint);
5845 #endif
5846 return 1;
5848 d -= DECDPUN;
5851 /* check the exponent. Note that input operands can have exponents */
5852 /* which are out of the set->emin/set->emax and set->digits range */
5853 /* (just as they can have more digits than set->digits). */
5854 ae = dn->exponent + dn->digits - 1; /* adjusted exponent */
5855 emax = DECNUMMAXE;
5856 emin = DECNUMMINE;
5857 digits = DECNUMMAXP;
5858 if (ae < emin - (digits - 1))
5860 #if DECTRACE
5861 printf ("Adjusted exponent underflow [%d].\n", ae);
5862 decNumberShow (dn);
5863 #endif
5864 return 1;
5866 if (ae > +emax)
5868 #if DECTRACE
5869 printf ("Adjusted exponent overflow [%d].\n", ae);
5870 decNumberShow (dn);
5871 #endif
5872 return 1;
5875 return 0; /* it's OK */
5877 #endif
5879 #if DECALLOC
5880 #undef malloc
5881 #undef free
5882 /* ------------------------------------------------------------------ */
5883 /* decMalloc -- accountable allocation routine */
5884 /* n is the number of bytes to allocate */
5885 /* */
5886 /* Semantics is the same as the stdlib malloc routine, but bytes */
5887 /* allocated are accounted for globally, and corruption fences are */
5888 /* added before and after the 'actual' storage. */
5889 /* ------------------------------------------------------------------ */
5890 /* This routine allocates storage with an extra twelve bytes; 8 are */
5891 /* at the start and hold: */
5892 /* 0-3 the original length requested */
5893 /* 4-7 buffer corruption detection fence (DECFENCE, x4) */
5894 /* The 4 bytes at the end also hold a corruption fence (DECFENCE, x4) */
5895 /* ------------------------------------------------------------------ */
5896 static void *
5897 decMalloc (uInt n)
5899 uInt size = n + 12; /* true size */
5900 void *alloc; /* -> allocated storage */
5901 uInt *j; /* work */
5902 uByte *b, *b0; /* .. */
5904 alloc = malloc (size); /* -> allocated storage */
5905 if (alloc == NULL)
5906 return NULL; /* out of strorage */
5907 b0 = (uByte *) alloc; /* as bytes */
5908 decAllocBytes += n; /* account for storage */
5909 j = (uInt *) alloc; /* -> first four bytes */
5910 *j = n; /* save n */
5911 /* printf("++ alloc(%d)\n", n); */
5912 for (b = b0 + 4; b < b0 + 8; b++)
5913 *b = DECFENCE;
5914 for (b = b0 + n + 8; b < b0 + n + 12; b++)
5915 *b = DECFENCE;
5916 return b0 + 8; /* -> play area */
5919 /* ------------------------------------------------------------------ */
5920 /* decFree -- accountable free routine */
5921 /* alloc is the storage to free */
5922 /* */
5923 /* Semantics is the same as the stdlib malloc routine, except that */
5924 /* the global storage accounting is updated and the fences are */
5925 /* checked to ensure that no routine has written 'out of bounds'. */
5926 /* ------------------------------------------------------------------ */
5927 /* This routine first checks that the fences have not been corrupted. */
5928 /* It then frees the storage using the 'truw' storage address (that */
5929 /* is, offset by 8). */
5930 /* ------------------------------------------------------------------ */
5931 static void
5932 decFree (void *alloc)
5934 uInt *j, n; /* pointer, original length */
5935 uByte *b, *b0; /* work */
5937 if (alloc == NULL)
5938 return; /* allowed; it's a nop */
5939 b0 = (uByte *) alloc; /* as bytes */
5940 b0 -= 8; /* -> true start of storage */
5941 j = (uInt *) b0; /* -> first four bytes */
5942 n = *j; /* lift */
5943 for (b = b0 + 4; b < b0 + 8; b++)
5944 if (*b != DECFENCE)
5945 printf ("=== Corrupt byte [%02x] at offset %d from %d ===\n", *b,
5946 b - b0 - 8, (Int) b0);
5947 for (b = b0 + n + 8; b < b0 + n + 12; b++)
5948 if (*b != DECFENCE)
5949 printf ("=== Corrupt byte [%02x] at offset +%d from %d, n=%d ===\n", *b,
5950 b - b0 - 8, (Int) b0, n);
5951 free (b0); /* drop the storage */
5952 decAllocBytes -= n; /* account for storage */
5954 #endif