OSDN Git Service

libgfortran:
[pf3gnuchains/gcc-fork.git] / libdecnumber / decNumber.c
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.
4
5    This file is part of GCC.
6
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.
11
12    In addition to the permissions in the GNU General Public License,
13    the Free Software Foundation gives you unlimited permission to link
14    the compiled version of this file into combinations with other
15    programs, and to distribute those combinations without any
16    restriction coming from the use of this file.  (The General Public
17    License restrictions do apply in other respects; for example, they
18    cover modification of the file, and distribution when not linked
19    into a combine executable.)
20
21    GCC is distributed in the hope that it will be useful, but WITHOUT ANY
22    WARRANTY; without even the implied warranty of MERCHANTABILITY or
23    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
24    for more details.
25
26    You should have received a copy of the GNU General Public License
27    along with GCC; see the file COPYING.  If not, write to the Free
28    Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
29    02110-1301, USA.  */
30
31 /* ------------------------------------------------------------------ */
32 /* This module comprises the routines for Standard Decimal Arithmetic */
33 /* as defined in the specification which may be found on the          */
34 /* http://www2.hursley.ibm.com/decimal web pages.  It implements both */
35 /* the full ('extended') arithmetic and the simpler ('subset')        */
36 /* arithmetic.                                                        */
37 /*                                                                    */
38 /* Usage notes:                                                       */
39 /*                                                                    */
40 /* 1. This code is ANSI C89 except:                                   */
41 /*                                                                    */
42 /*    a) Line comments (double forward slash) are used.  (Most C      */
43 /*       compilers accept these.  If yours does not, a simple script  */
44 /*       can be used to convert them to ANSI C comments.)             */
45 /*                                                                    */
46 /*    b) Types from C99 stdint.h are used.  If you do not have this   */
47 /*       header file, see the User's Guide section of the decNumber   */
48 /*       documentation; this lists the necessary definitions.         */
49 /*                                                                    */
50 /*    c) If DECDPUN>4, non-ANSI 64-bit 'long long' types are used.    */
51 /*       To avoid these, set DECDPUN <= 4 (see documentation).        */
52 /*                                                                    */
53 /* 2. The decNumber format which this library uses is optimized for   */
54 /*    efficient processing of relatively short numbers; in particular */
55 /*    it allows the use of fixed sized structures and minimizes copy  */
56 /*    and move operations.  It does, however, support arbitrary       */
57 /*    precision (up to 999,999,999 digits) and arbitrary exponent     */
58 /*    range (Emax in the range 0 through 999,999,999 and Emin in the  */
59 /*    range -999,999,999 through 0).                                  */
60 /*                                                                    */
61 /* 3. Operands to operator functions are never modified unless they   */
62 /*    are also specified to be the result number (which is always     */
63 /*    permitted).  Other than that case, operands may not overlap.    */
64 /*                                                                    */
65 /* 4. Error handling: the type of the error is ORed into the status   */
66 /*    flags in the current context (decContext structure).  The       */
67 /*    SIGFPE signal is then raised if the corresponding trap-enabler  */
68 /*    flag in the decContext is set (is 1).                           */
69 /*                                                                    */
70 /*    It is the responsibility of the caller to clear the status      */
71 /*    flags as required.                                              */
72 /*                                                                    */
73 /*    The result of any routine which returns a number will always    */
74 /*    be a valid number (which may be a special value, such as an     */
75 /*    Infinity or NaN).                                               */
76 /*                                                                    */
77 /* 5. The decNumber format is not an exchangeable concrete            */
78 /*    representation as it comprises fields which may be machine-     */
79 /*    dependent (big-endian or little-endian, for example).           */
80 /*    Canonical conversions to and from strings are provided; other   */
81 /*    conversions are available in separate modules.                  */
82 /*                                                                    */
83 /* 6. Normally, input operands are assumed to be valid.  Set DECCHECK */
84 /*    to 1 for extended operand checking (including NULL operands).   */
85 /*    Results are undefined if a badly-formed structure (or a NULL    */
86 /*    NULL pointer to a structure) is provided, though with DECCHECK  */
87 /*    enabled the operator routines are protected against exceptions. */
88 /*    (Except if the result pointer is NULL, which is unrecoverable.) */
89 /*                                                                    */
90 /*    However, the routines will never cause exceptions if they are   */
91 /*    given well-formed operands, even if the value of the operands   */
92 /*    is inappropriate for the operation and DECCHECK is not set.     */
93 /*                                                                    */
94 /* 7. Subset arithmetic is available only if DECSUBSET is set to 1.   */
95 /* ------------------------------------------------------------------ */
96 /* Implementation notes for maintenance of this module:               */
97 /*                                                                    */
98 /* 1. Storage leak protection:  Routines which use malloc are not     */
99 /*    permitted to use return for fastpath or error exits (i.e.,      */
100 /*    they follow strict structured programming conventions).         */
101 /*    Instead they have a do{}while(0); construct surrounding the     */
102 /*    code which is protected -- break may be used from this.         */
103 /*    Other routines are allowed to use the return statement inline.  */
104 /*                                                                    */
105 /*    Storage leak accounting can be enabled using DECALLOC.          */
106 /*                                                                    */
107 /* 2. All loops use the for(;;) construct.  Any do construct is for   */
108 /*    protection as just described.                                   */
109 /*                                                                    */
110 /* 3. Setting status in the context must always be the very last      */
111 /*    action in a routine, as non-0 status may raise a trap and hence */
112 /*    the call to set status may not return (if the handler uses long */
113 /*    jump).  Therefore all cleanup must be done first.  In general,  */
114 /*    to achieve this we accumulate status and only finally apply it  */
115 /*    by calling decContextSetStatus (via decStatus).                 */
116 /*                                                                    */
117 /*    Routines which allocate storage cannot, therefore, use the      */
118 /*    'top level' routines which could cause a non-returning          */
119 /*    transfer of control.  The decXxxxOp routines are safe (do not   */
120 /*    call decStatus even if traps are set in the context) and should */
121 /*    be used instead (they are also a little faster).                */
122 /*                                                                    */
123 /* 4. Exponent checking is minimized by allowing the exponent to      */
124 /*    grow outside its limits during calculations, provided that      */
125 /*    the decFinalize function is called later.  Multiplication and   */
126 /*    division, and intermediate calculations in exponentiation,      */
127 /*    require more careful checks because of the risk of 31-bit       */
128 /*    overflow (the most negative valid exponent is -1999999997, for  */
129 /*    a 999999999-digit number with adjusted exponent of -999999999). */
130 /*                                                                    */
131 /* 5. Rounding is deferred until finalization of results, with any    */
132 /*    'off to the right' data being represented as a single digit     */
133 /*    residue (in the range -1 through 9).  This avoids any double-   */
134 /*    rounding when more than one shortening takes place (for         */
135 /*    example, when a result is subnormal).                           */
136 /*                                                                    */
137 /* 6. The digits count is allowed to rise to a multiple of DECDPUN    */
138 /*    during many operations, so whole Units are handled and exact    */
139 /*    accounting of digits is not needed.  The correct digits value   */
140 /*    is found by decGetDigits, which accounts for leading zeros.     */
141 /*    This must be called before any rounding if the number of digits */
142 /*    is not known exactly.                                           */
143 /*                                                                    */
144 /* 7. We use the multiply-by-reciprocal 'trick' for partitioning      */
145 /*    numbers up to four digits, using appropriate constants.  This   */
146 /*    is not useful for longer numbers because overflow of 32 bits    */
147 /*    would lead to 4 multiplies, which is almost as expensive as     */
148 /*    a divide (unless we assumed floating-point multiply available). */
149 /*                                                                    */
150 /* 8. Unusual abbreviations possibly used in the commentary:          */
151 /*      lhs -- left hand side (operand, of an operation)              */
152 /*      lsd -- least significant digit (of coefficient)               */
153 /*      lsu -- least significant Unit (of coefficient)                */
154 /*      msd -- most significant digit (of coefficient)                */
155 /*      msu -- most significant Unit (of coefficient)                 */
156 /*      rhs -- right hand side (operand, of an operation)             */
157 /*      +ve -- positive                                               */
158 /*      -ve -- negative                                               */
159 /* ------------------------------------------------------------------ */
160
161 /* Some of glibc's string inlines cause warnings.  Plus we'd rather
162    rely on (and therefore test) GCC's string builtins.  */
163 #define __NO_STRING_INLINES
164
165 #include <stdlib.h>             /* for malloc, free, etc. */
166 #include <stdio.h>              /* for printf [if needed] */
167 #include <string.h>             /* for strcpy */
168 #include <ctype.h>              /* for lower */
169 #include "config.h"
170 #include "decNumber.h"          /* base number library */
171 #include "decNumberLocal.h"     /* decNumber local types, etc. */
172
173 /* Constants */
174 /* Public constant array: powers of ten (powers[n]==10**n) */
175 const uInt powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
176   10000000, 100000000, 1000000000
177 };
178
179 /* Local constants */
180 #define DIVIDE    0x80          /* Divide operators */
181 #define REMAINDER 0x40          /* .. */
182 #define DIVIDEINT 0x20          /* .. */
183 #define REMNEAR   0x10          /* .. */
184 #define COMPARE   0x01          /* Compare operators */
185 #define COMPMAX   0x02          /* .. */
186 #define COMPMIN   0x03          /* .. */
187 #define COMPNAN   0x04          /* .. [NaN processing] */
188
189 #define DEC_sNaN 0x40000000     /* local status: sNaN signal */
190 #define BADINT (Int)0x80000000  /* most-negative Int; error indicator */
191
192 static Unit one[] = { 1 };      /* Unit array of 1, used for incrementing */
193
194 /* Granularity-dependent code */
195 #if DECDPUN<=4
196 #define eInt  Int               /* extended integer */
197 #define ueInt uInt              /* unsigned extended integer */
198   /* Constant multipliers for divide-by-power-of five using reciprocal */
199   /* multiply, after removing powers of 2 by shifting, and final shift */
200   /* of 17 [we only need up to **4] */
201 static const uInt multies[] = { 131073, 26215, 5243, 1049, 210 };
202
203   /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
204 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
205 #else
206   /* For DECDPUN>4 we currently use non-ANSI 64-bit types.  These could */
207   /* be replaced by subroutine calls later. */
208 #ifdef long
209 #undef long
210 #endif
211 typedef signed long long Long;
212 typedef unsigned long long uLong;
213 #define eInt  Long              /* extended integer */
214 #define ueInt uLong             /* unsigned extended integer */
215 #endif
216
217 /* Local routines */
218 static decNumber *decAddOp (decNumber *, const decNumber *,
219                             const decNumber *, decContext *,
220                             uByte, uInt *);
221 static void decApplyRound (decNumber *, decContext *, Int, uInt *);
222 static Int decCompare (const decNumber * lhs, const decNumber * rhs);
223 static decNumber *decCompareOp (decNumber *, const decNumber *, const decNumber *,
224                                 decContext *, Flag, uInt *);
225 static void decCopyFit (decNumber *, const decNumber *, decContext *,
226                         Int *, uInt *);
227 static decNumber *decDivideOp (decNumber *, const decNumber *, const decNumber *,
228                                decContext *, Flag, uInt *);
229 static void decFinalize (decNumber *, decContext *, Int *, uInt *);
230 static Int decGetDigits (const Unit *, Int);
231 #if DECSUBSET
232 static Int decGetInt (const decNumber *, decContext *);
233 #else
234 static Int decGetInt (const decNumber *);
235 #endif
236 static decNumber *decMultiplyOp (decNumber *, const decNumber *,
237                                  const decNumber *, decContext *, uInt *);
238 static decNumber *decNaNs (decNumber *, const decNumber *, const decNumber *, uInt *);
239 static decNumber *decQuantizeOp (decNumber *, const decNumber *,
240                                  const decNumber *, decContext *, Flag, uInt *);
241 static void decSetCoeff (decNumber *, decContext *, const Unit *,
242                          Int, Int *, uInt *);
243 static void decSetOverflow (decNumber *, decContext *, uInt *);
244 static void decSetSubnormal (decNumber *, decContext *, Int *, uInt *);
245 static Int decShiftToLeast (Unit *, Int, Int);
246 static Int decShiftToMost (Unit *, Int, Int);
247 static void decStatus (decNumber *, uInt, decContext *);
248 static Flag decStrEq (const char *, const char *);
249 static void decToString (const decNumber *, char[], Flag);
250 static decNumber *decTrim (decNumber *, Flag, Int *);
251 static Int decUnitAddSub (const Unit *, Int, const Unit *, Int, Int, Unit *, Int);
252 static Int decUnitCompare (const Unit *, Int, const Unit *, Int, Int);
253
254 #if !DECSUBSET
255 /* decFinish == decFinalize when no subset arithmetic needed */
256 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
257 #else
258 static void decFinish (decNumber *, decContext *, Int *, uInt *);
259 static decNumber *decRoundOperand (const decNumber *, decContext *, uInt *);
260 #endif
261
262 /* Diagnostic macros, etc. */
263 #if DECALLOC
264 /* Handle malloc/free accounting.  If enabled, our accountable routines */
265 /* are used; otherwise the code just goes straight to the system malloc */
266 /* and free routines. */
267 #define malloc(a) decMalloc(a)
268 #define free(a) decFree(a)
269 #define DECFENCE 0x5a           /* corruption detector */
270 /* 'Our' malloc and free: */
271 static void *decMalloc (size_t);
272 static void decFree (void *);
273 uInt decAllocBytes = 0;         /* count of bytes allocated */
274 /* Note that DECALLOC code only checks for storage buffer overflow. */
275 /* To check for memory leaks, the decAllocBytes variable should be */
276 /* checked to be 0 at appropriate times (e.g., after the test */
277 /* harness completes a set of tests).  This checking may be unreliable */
278 /* if the testing is done in a multi-thread environment. */
279 #endif
280
281 #if DECCHECK
282 /* Optional operand checking routines.  Enabling these means that */
283 /* decNumber and decContext operands to operator routines are checked */
284 /* for correctness.  This roughly doubles the execution time of the */
285 /* fastest routines (and adds 600+ bytes), so should not normally be */
286 /* used in 'production'. */
287 #define DECUNUSED (void *)(0xffffffff)
288 static Flag decCheckOperands (decNumber *, const decNumber *,
289                               const decNumber *, decContext *);
290 static Flag decCheckNumber (const decNumber *, decContext *);
291 #endif
292
293 #if DECTRACE || DECCHECK
294 /* Optional trace/debugging routines. */
295 void decNumberShow (const decNumber *); /* displays the components of a number */
296 static void decDumpAr (char, const Unit *, Int);
297 #endif
298
299 /* ================================================================== */
300 /* Conversions                                                        */
301 /* ================================================================== */
302
303 /* ------------------------------------------------------------------ */
304 /* to-scientific-string -- conversion to numeric string               */
305 /* to-engineering-string -- conversion to numeric string              */
306 /*                                                                    */
307 /*   decNumberToString(dn, string);                                   */
308 /*   decNumberToEngString(dn, string);                                */
309 /*                                                                    */
310 /*  dn is the decNumber to convert                                    */
311 /*  string is the string where the result will be laid out            */
312 /*                                                                    */
313 /*  string must be at least dn->digits+14 characters long             */
314 /*                                                                    */
315 /*  No error is possible, and no status can be set.                   */
316 /* ------------------------------------------------------------------ */
317 char *
318 decNumberToString (const decNumber * dn, char *string)
319 {
320   decToString (dn, string, 0);
321   return string;
322 }
323
324 char *
325 decNumberToEngString (const decNumber * dn, char *string)
326 {
327   decToString (dn, string, 1);
328   return string;
329 }
330
331 /* ------------------------------------------------------------------ */
332 /* to-number -- conversion from numeric string                        */
333 /*                                                                    */
334 /* decNumberFromString -- convert string to decNumber                 */
335 /*   dn        -- the number structure to fill                        */
336 /*   chars[]   -- the string to convert ('\0' terminated)             */
337 /*   set       -- the context used for processing any error,          */
338 /*                determining the maximum precision available         */
339 /*                (set.digits), determining the maximum and minimum   */
340 /*                exponent (set.emax and set.emin), determining if    */
341 /*                extended values are allowed, and checking the       */
342 /*                rounding mode if overflow occurs or rounding is     */
343 /*                needed.                                             */
344 /*                                                                    */
345 /* The length of the coefficient and the size of the exponent are     */
346 /* checked by this routine, so the correct error (Underflow or        */
347 /* Overflow) can be reported or rounding applied, as necessary.       */
348 /*                                                                    */
349 /* If bad syntax is detected, the result will be a quiet NaN.         */
350 /* ------------------------------------------------------------------ */
351 decNumber *
352 decNumberFromString (decNumber * dn, const char chars[], decContext * set)
353 {
354   Int exponent = 0;             /* working exponent [assume 0] */
355   uByte bits = 0;               /* working flags [assume +ve] */
356   Unit *res;                    /* where result will be built */
357   Unit resbuff[D2U (DECBUFFER + 1)];    /* local buffer in case need temporary */
358   Unit *allocres = NULL;        /* -> allocated result, iff allocated */
359   Int need;                     /* units needed for result */
360   Int d = 0;                    /* count of digits found in decimal part */
361   const char *dotchar = NULL;   /* where dot was found */
362   const char *cfirst;           /* -> first character of decimal part */
363   const char *last = NULL;      /* -> last digit of decimal part */
364   const char *firstexp;         /* -> first significant exponent digit */
365   const char *c;                /* work */
366   Unit *up;                     /* .. */
367 #if DECDPUN>1
368   Int i;                        /* .. */
369 #endif
370   Int residue = 0;              /* rounding residue */
371   uInt status = 0;              /* error code */
372
373 #if DECCHECK
374   if (decCheckOperands (DECUNUSED, DECUNUSED, DECUNUSED, set))
375     return decNumberZero (dn);
376 #endif
377
378   do
379     {                           /* status & malloc protection */
380       c = chars;                /* -> input character */
381       if (*c == '-')
382         {                       /* handle leading '-' */
383           bits = DECNEG;
384           c++;
385         }
386       else if (*c == '+')
387         c++;                    /* step over leading '+' */
388       /* We're at the start of the number [we think] */
389       cfirst = c;               /* save */
390       for (;; c++)
391         {
392           if (*c >= '0' && *c <= '9')
393             {                   /* test for Arabic digit */
394               last = c;
395               d++;              /* count of real digits */
396               continue;         /* still in decimal part */
397             }
398           if (*c != '.')
399             break;              /* done with decimal part */
400           /* dot: record, check, and ignore */
401           if (dotchar != NULL)
402             {                   /* two dots */
403               last = NULL;      /* indicate bad */
404               break;
405             }                   /* .. and go report */
406           dotchar = c;          /* offset into decimal part */
407         }                       /* c */
408
409       if (last == NULL)
410         {                       /* no decimal digits, or >1 . */
411 #if DECSUBSET
412           /* If subset then infinities and NaNs are not allowed */
413           if (!set->extended)
414             {
415               status = DEC_Conversion_syntax;
416               break;            /* all done */
417             }
418           else
419             {
420 #endif
421               /* Infinities and NaNs are possible, here */
422               decNumberZero (dn);       /* be optimistic */
423               if (decStrEq (c, "Infinity") || decStrEq (c, "Inf"))
424                 {
425                   dn->bits = bits | DECINF;
426                   break;        /* all done */
427                 }
428               else
429                 {               /* a NaN expected */
430                   /* 2003.09.10 NaNs are now permitted to have a sign */
431                   status = DEC_Conversion_syntax;       /* assume the worst */
432                   dn->bits = bits | DECNAN;     /* assume simple NaN */
433                   if (*c == 's' || *c == 'S')
434                     {           /* looks like an` sNaN */
435                       c++;
436                       dn->bits = bits | DECSNAN;
437                     }
438                   if (*c != 'n' && *c != 'N')
439                     break;      /* check caseless "NaN" */
440                   c++;
441                   if (*c != 'a' && *c != 'A')
442                     break;      /* .. */
443                   c++;
444                   if (*c != 'n' && *c != 'N')
445                     break;      /* .. */
446                   c++;
447                   /* now nothing, or nnnn, expected */
448                   /* -> start of integer and skip leading 0s [including plain 0] */
449                   for (cfirst = c; *cfirst == '0';)
450                     cfirst++;
451                   if (*cfirst == '\0')
452                     {           /* "NaN" or "sNaN", maybe with all 0s */
453                       status = 0;       /* it's good */
454                       break;    /* .. */
455                     }
456                   /* something other than 0s; setup last and d as usual [no dots] */
457                   for (c = cfirst;; c++, d++)
458                     {
459                       if (*c < '0' || *c > '9')
460                         break;  /* test for Arabic digit */
461                       last = c;
462                     }
463                   if (*c != '\0')
464                     break;      /* not all digits */
465                   if (d > set->digits)
466                     break;      /* too many digits */
467                   /* good; drop through and convert the integer */
468                   status = 0;
469                   bits = dn->bits;      /* for copy-back */
470                 }               /* NaN expected */
471 #if DECSUBSET
472             }
473 #endif
474         }                       /* last==NULL */
475
476       if (*c != '\0')
477         {                       /* more there; exponent expected... */
478           Flag nege = 0;        /* 1=negative exponent */
479           if (*c != 'e' && *c != 'E')
480             {
481               status = DEC_Conversion_syntax;
482               break;
483             }
484
485           /* Found 'e' or 'E' -- now process explicit exponent */
486           /* 1998.07.11: sign no longer required */
487           c++;                  /* to (expected) sign */
488           if (*c == '-')
489             {
490               nege = 1;
491               c++;
492             }
493           else if (*c == '+')
494             c++;
495           if (*c == '\0')
496             {
497               status = DEC_Conversion_syntax;
498               break;
499             }
500
501           for (; *c == '0' && *(c + 1) != '\0';)
502             c++;                /* strip insignificant zeros */
503           firstexp = c;         /* save exponent digit place */
504           for (;; c++)
505             {
506               if (*c < '0' || *c > '9')
507                 break;          /* not a digit */
508               exponent = X10 (exponent) + (Int) * c - (Int) '0';
509             }                   /* c */
510           /* if we didn't end on '\0' must not be a digit */
511           if (*c != '\0')
512             {
513               status = DEC_Conversion_syntax;
514               break;
515             }
516
517           /* (this next test must be after the syntax check) */
518           /* if it was too long the exponent may have wrapped, so check */
519           /* carefully and set it to a certain overflow if wrap possible */
520           if (c >= firstexp + 9 + 1)
521             {
522               if (c > firstexp + 9 + 1 || *firstexp > '1')
523                 exponent = DECNUMMAXE * 2;
524               /* [up to 1999999999 is OK, for example 1E-1000000998] */
525             }
526           if (nege)
527             exponent = -exponent;       /* was negative */
528         }                       /* had exponent */
529       /* Here when all inspected; syntax is good */
530
531       /* Handle decimal point... */
532       if (dotchar != NULL && dotchar < last)    /* embedded . found, so */
533         exponent = exponent - (last - dotchar); /* .. adjust exponent */
534       /* [we can now ignore the .] */
535
536       /* strip leading zeros/dot (leave final if all 0's) */
537       for (c = cfirst; c < last; c++)
538         {
539           if (*c == '0')
540             d--;                /* 0 stripped */
541           else if (*c != '.')
542             break;
543           cfirst++;             /* step past leader */
544         }                       /* c */
545
546 #if DECSUBSET
547       /* We can now make a rapid exit for zeros if !extended */
548       if (*cfirst == '0' && !set->extended)
549         {
550           decNumberZero (dn);   /* clean result */
551           break;                /* [could be return] */
552         }
553 #endif
554
555       /* OK, the digits string is good.  Copy to the decNumber, or to
556          a temporary decNumber if rounding is needed */
557       if (d <= set->digits)
558         res = dn->lsu;          /* fits into given decNumber */
559       else
560         {                       /* rounding needed */
561           need = D2U (d);       /* units needed */
562           res = resbuff;        /* assume use local buffer */
563           if (need * sizeof (Unit) > sizeof (resbuff))
564             {                   /* too big for local */
565               allocres = (Unit *) malloc (need * sizeof (Unit));
566               if (allocres == NULL)
567                 {
568                   status |= DEC_Insufficient_storage;
569                   break;
570                 }
571               res = allocres;
572             }
573         }
574       /* res now -> number lsu, buffer, or allocated storage for Unit array */
575
576       /* Place the coefficient into the selected Unit array */
577 #if DECDPUN>1
578       i = d % DECDPUN;          /* digits in top unit */
579       if (i == 0)
580         i = DECDPUN;
581       up = res + D2U (d) - 1;   /* -> msu */
582       *up = 0;
583       for (c = cfirst;; c++)
584         {                       /* along the digits */
585           if (*c == '.')
586             {                   /* ignore . [don't decrement i] */
587               if (c != last)
588                 continue;
589               break;
590             }
591           *up = (Unit) (X10 (*up) + (Int) * c - (Int) '0');
592           i--;
593           if (i > 0)
594             continue;           /* more for this unit */
595           if (up == res)
596             break;              /* just filled the last unit */
597           i = DECDPUN;
598           up--;
599           *up = 0;
600         }                       /* c */
601 #else
602       /* DECDPUN==1 */
603       up = res;                 /* -> lsu */
604       for (c = last; c >= cfirst; c--)
605         {                       /* over each character, from least */
606           if (*c == '.')
607             continue;           /* ignore . [don't step b] */
608           *up = (Unit) ((Int) * c - (Int) '0');
609           up++;
610         }                       /* c */
611 #endif
612
613       dn->bits = bits;
614       dn->exponent = exponent;
615       dn->digits = d;
616
617       /* if not in number (too long) shorten into the number */
618       if (d > set->digits)
619         decSetCoeff (dn, set, res, d, &residue, &status);
620
621       /* Finally check for overflow or subnormal and round as needed */
622       decFinalize (dn, set, &residue, &status);
623       /* decNumberShow(dn); */
624     }
625   while (0);                    /* [for break] */
626
627   if (allocres != NULL)
628     free (allocres);            /* drop any storage we used */
629   if (status != 0)
630     decStatus (dn, status, set);
631   return dn;
632 }
633
634 /* ================================================================== */
635 /* Operators                                                          */
636 /* ================================================================== */
637
638 /* ------------------------------------------------------------------ */
639 /* decNumberAbs -- absolute value operator                            */
640 /*                                                                    */
641 /*   This computes C = abs(A)                                         */
642 /*                                                                    */
643 /*   res is C, the result.  C may be A                                */
644 /*   rhs is A                                                         */
645 /*   set is the context                                               */
646 /*                                                                    */
647 /* C must have space for set->digits digits.                          */
648 /* ------------------------------------------------------------------ */
649 /* This has the same effect as decNumberPlus unless A is negative,    */
650 /* in which case it has the same effect as decNumberMinus.            */
651 /* ------------------------------------------------------------------ */
652 decNumber *
653 decNumberAbs (decNumber * res, const decNumber * rhs, decContext * set)
654 {
655   decNumber dzero;              /* for 0 */
656   uInt status = 0;              /* accumulator */
657
658 #if DECCHECK
659   if (decCheckOperands (res, DECUNUSED, rhs, set))
660     return res;
661 #endif
662
663   decNumberZero (&dzero);       /* set 0 */
664   dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
665   decAddOp (res, &dzero, rhs, set, (uByte) (rhs->bits & DECNEG), &status);
666   if (status != 0)
667     decStatus (res, status, set);
668   return res;
669 }
670
671 /* ------------------------------------------------------------------ */
672 /* decNumberAdd -- add two Numbers                                    */
673 /*                                                                    */
674 /*   This computes C = A + B                                          */
675 /*                                                                    */
676 /*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
677 /*   lhs is A                                                         */
678 /*   rhs is B                                                         */
679 /*   set is the context                                               */
680 /*                                                                    */
681 /* C must have space for set->digits digits.                          */
682 /* ------------------------------------------------------------------ */
683 /* This just calls the routine shared with Subtract                   */
684 decNumber *
685 decNumberAdd (decNumber * res, const decNumber * lhs,
686               const decNumber * rhs, decContext * set)
687 {
688   uInt status = 0;              /* accumulator */
689   decAddOp (res, lhs, rhs, set, 0, &status);
690   if (status != 0)
691     decStatus (res, status, set);
692   return res;
693 }
694
695 /* ------------------------------------------------------------------ */
696 /* decNumberCompare -- compare two Numbers                            */
697 /*                                                                    */
698 /*   This computes C = A ? B                                          */
699 /*                                                                    */
700 /*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
701 /*   lhs is A                                                         */
702 /*   rhs is B                                                         */
703 /*   set is the context                                               */
704 /*                                                                    */
705 /* C must have space for one digit.                                   */
706 /* ------------------------------------------------------------------ */
707 decNumber *
708 decNumberCompare (decNumber * res, const decNumber * lhs,
709                   const decNumber * rhs, decContext * set)
710 {
711   uInt status = 0;              /* accumulator */
712   decCompareOp (res, lhs, rhs, set, COMPARE, &status);
713   if (status != 0)
714     decStatus (res, status, set);
715   return res;
716 }
717
718 /* ------------------------------------------------------------------ */
719 /* decNumberDivide -- divide one number by another                    */
720 /*                                                                    */
721 /*   This computes C = A / B                                          */
722 /*                                                                    */
723 /*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
724 /*   lhs is A                                                         */
725 /*   rhs is B                                                         */
726 /*   set is the context                                               */
727 /*                                                                    */
728 /* C must have space for set->digits digits.                          */
729 /* ------------------------------------------------------------------ */
730 decNumber *
731 decNumberDivide (decNumber * res, const decNumber * lhs,
732                  const decNumber * rhs, decContext * set)
733 {
734   uInt status = 0;              /* accumulator */
735   decDivideOp (res, lhs, rhs, set, DIVIDE, &status);
736   if (status != 0)
737     decStatus (res, status, set);
738   return res;
739 }
740
741 /* ------------------------------------------------------------------ */
742 /* decNumberDivideInteger -- divide and return integer quotient       */
743 /*                                                                    */
744 /*   This computes C = A # B, where # is the integer divide operator  */
745 /*                                                                    */
746 /*   res is C, the result.  C may be A and/or B (e.g., X=X#X)         */
747 /*   lhs is A                                                         */
748 /*   rhs is B                                                         */
749 /*   set is the context                                               */
750 /*                                                                    */
751 /* C must have space for set->digits digits.                          */
752 /* ------------------------------------------------------------------ */
753 decNumber *
754 decNumberDivideInteger (decNumber * res, const decNumber * lhs,
755                         const decNumber * rhs, decContext * set)
756 {
757   uInt status = 0;              /* accumulator */
758   decDivideOp (res, lhs, rhs, set, DIVIDEINT, &status);
759   if (status != 0)
760     decStatus (res, status, set);
761   return res;
762 }
763
764 /* ------------------------------------------------------------------ */
765 /* decNumberMax -- compare two Numbers and return the maximum         */
766 /*                                                                    */
767 /*   This computes C = A ? B, returning the maximum or A if equal     */
768 /*                                                                    */
769 /*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
770 /*   lhs is A                                                         */
771 /*   rhs is B                                                         */
772 /*   set is the context                                               */
773 /*                                                                    */
774 /* C must have space for set->digits digits.                          */
775 /* ------------------------------------------------------------------ */
776 decNumber *
777 decNumberMax (decNumber * res, const decNumber * lhs,
778               const decNumber * rhs, decContext * set)
779 {
780   uInt status = 0;              /* accumulator */
781   decCompareOp (res, lhs, rhs, set, COMPMAX, &status);
782   if (status != 0)
783     decStatus (res, status, set);
784   return res;
785 }
786
787 /* ------------------------------------------------------------------ */
788 /* decNumberMin -- compare two Numbers and return the minimum         */
789 /*                                                                    */
790 /*   This computes C = A ? B, returning the minimum or A if equal     */
791 /*                                                                    */
792 /*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
793 /*   lhs is A                                                         */
794 /*   rhs is B                                                         */
795 /*   set is the context                                               */
796 /*                                                                    */
797 /* C must have space for set->digits digits.                          */
798 /* ------------------------------------------------------------------ */
799 decNumber *
800 decNumberMin (decNumber * res, const decNumber * lhs,
801               const decNumber * rhs, decContext * set)
802 {
803   uInt status = 0;              /* accumulator */
804   decCompareOp (res, lhs, rhs, set, COMPMIN, &status);
805   if (status != 0)
806     decStatus (res, status, set);
807   return res;
808 }
809
810 /* ------------------------------------------------------------------ */
811 /* decNumberMinus -- prefix minus operator                            */
812 /*                                                                    */
813 /*   This computes C = 0 - A                                          */
814 /*                                                                    */
815 /*   res is C, the result.  C may be A                                */
816 /*   rhs is A                                                         */
817 /*   set is the context                                               */
818 /*                                                                    */
819 /* C must have space for set->digits digits.                          */
820 /* ------------------------------------------------------------------ */
821 /* We simply use AddOp for the subtract, which will do the necessary. */
822 /* ------------------------------------------------------------------ */
823 decNumber *
824 decNumberMinus (decNumber * res, const decNumber * rhs, decContext * set)
825 {
826   decNumber dzero;
827   uInt status = 0;              /* accumulator */
828
829 #if DECCHECK
830   if (decCheckOperands (res, DECUNUSED, rhs, set))
831     return res;
832 #endif
833
834   decNumberZero (&dzero);       /* make 0 */
835   dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
836   decAddOp (res, &dzero, rhs, set, DECNEG, &status);
837   if (status != 0)
838     decStatus (res, status, set);
839   return res;
840 }
841
842 /* ------------------------------------------------------------------ */
843 /* decNumberPlus -- prefix plus operator                              */
844 /*                                                                    */
845 /*   This computes C = 0 + A                                          */
846 /*                                                                    */
847 /*   res is C, the result.  C may be A                                */
848 /*   rhs is A                                                         */
849 /*   set is the context                                               */
850 /*                                                                    */
851 /* C must have space for set->digits digits.                          */
852 /* ------------------------------------------------------------------ */
853 /* We simply use AddOp; Add will take fast path after preparing A.    */
854 /* Performance is a concern here, as this routine is often used to    */
855 /* check operands and apply rounding and overflow/underflow testing.  */
856 /* ------------------------------------------------------------------ */
857 decNumber *
858 decNumberPlus (decNumber * res, const decNumber * rhs, decContext * set)
859 {
860   decNumber dzero;
861   uInt status = 0;              /* accumulator */
862
863 #if DECCHECK
864   if (decCheckOperands (res, DECUNUSED, rhs, set))
865     return res;
866 #endif
867
868   decNumberZero (&dzero);       /* make 0 */
869   dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
870   decAddOp (res, &dzero, rhs, set, 0, &status);
871   if (status != 0)
872     decStatus (res, status, set);
873   return res;
874 }
875
876 /* ------------------------------------------------------------------ */
877 /* decNumberMultiply -- multiply two Numbers                          */
878 /*                                                                    */
879 /*   This computes C = A x B                                          */
880 /*                                                                    */
881 /*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
882 /*   lhs is A                                                         */
883 /*   rhs is B                                                         */
884 /*   set is the context                                               */
885 /*                                                                    */
886 /* C must have space for set->digits digits.                          */
887 /* ------------------------------------------------------------------ */
888 decNumber *
889 decNumberMultiply (decNumber * res, const decNumber * lhs,
890                    const decNumber * rhs, decContext * set)
891 {
892   uInt status = 0;              /* accumulator */
893   decMultiplyOp (res, lhs, rhs, set, &status);
894   if (status != 0)
895     decStatus (res, status, set);
896   return res;
897 }
898
899 /* ------------------------------------------------------------------ */
900 /* decNumberNormalize -- remove trailing zeros                        */
901 /*                                                                    */
902 /*   This computes C = 0 + A, and normalizes the result               */
903 /*                                                                    */
904 /*   res is C, the result.  C may be A                                */
905 /*   rhs is A                                                         */
906 /*   set is the context                                               */
907 /*                                                                    */
908 /* C must have space for set->digits digits.                          */
909 /* ------------------------------------------------------------------ */
910 decNumber *
911 decNumberNormalize (decNumber * res, const decNumber * rhs, decContext * set)
912 {
913   decNumber *allocrhs = NULL;   /* non-NULL if rounded rhs allocated */
914   uInt status = 0;              /* as usual */
915   Int residue = 0;              /* as usual */
916   Int dropped;                  /* work */
917
918 #if DECCHECK
919   if (decCheckOperands (res, DECUNUSED, rhs, set))
920     return res;
921 #endif
922
923   do
924     {                           /* protect allocated storage */
925 #if DECSUBSET
926       if (!set->extended)
927         {
928           /* reduce operand and set lostDigits status, as needed */
929           if (rhs->digits > set->digits)
930             {
931               allocrhs = decRoundOperand (rhs, set, &status);
932               if (allocrhs == NULL)
933                 break;
934               rhs = allocrhs;
935             }
936         }
937 #endif
938       /* [following code does not require input rounding] */
939
940       /* specials copy through, except NaNs need care */
941       if (decNumberIsNaN (rhs))
942         {
943           decNaNs (res, rhs, NULL, &status);
944           break;
945         }
946
947       /* reduce result to the requested length and copy to result */
948       decCopyFit (res, rhs, set, &residue, &status);    /* copy & round */
949       decFinish (res, set, &residue, &status);  /* cleanup/set flags */
950       decTrim (res, 1, &dropped);       /* normalize in place */
951     }
952   while (0);                    /* end protected */
953
954   if (allocrhs != NULL)
955     free (allocrhs);            /* .. */
956   if (status != 0)
957     decStatus (res, status, set);       /* then report status */
958   return res;
959 }
960
961 /* ------------------------------------------------------------------ */
962 /* decNumberPower -- raise a number to an integer power               */
963 /*                                                                    */
964 /*   This computes C = A ** B                                         */
965 /*                                                                    */
966 /*   res is C, the result.  C may be A and/or B (e.g., X=X**X)        */
967 /*   lhs is A                                                         */
968 /*   rhs is B                                                         */
969 /*   set is the context                                               */
970 /*                                                                    */
971 /* C must have space for set->digits digits.                          */
972 /*                                                                    */
973 /* Specification restriction: abs(n) must be <=999999999              */
974 /* ------------------------------------------------------------------ */
975 decNumber *
976 decNumberPower (decNumber * res, const decNumber * lhs,
977                 const decNumber * rhs, decContext * set)
978 {
979   decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
980   decNumber *allocrhs = NULL;   /* .., rhs */
981   decNumber *allocdac = NULL;   /* -> allocated acc buffer, iff used */
982   const decNumber *inrhs = rhs; /* save original rhs */
983   Int reqdigits = set->digits;  /* requested DIGITS */
984   Int n;                        /* RHS in binary */
985   Int i;                        /* work */
986 #if DECSUBSET
987   Int dropped;                  /* .. */
988 #endif
989   uInt needbytes;               /* buffer size needed */
990   Flag seenbit;                 /* seen a bit while powering */
991   Int residue = 0;              /* rounding residue */
992   uInt status = 0;              /* accumulator */
993   uByte bits = 0;               /* result sign if errors */
994   decContext workset;           /* working context */
995   decNumber dnOne;              /* work value 1... */
996   /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
997   uByte dacbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
998   /* same again for possible 1/lhs calculation */
999   uByte lhsbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
1000   decNumber *dac = (decNumber *) dacbuff;       /* -> result accumulator */
1001
1002 #if DECCHECK
1003   if (decCheckOperands (res, lhs, rhs, set))
1004     return res;
1005 #endif
1006
1007   do
1008     {                           /* protect allocated storage */
1009 #if DECSUBSET
1010       if (!set->extended)
1011         {
1012           /* reduce operands and set lostDigits status, as needed */
1013           if (lhs->digits > reqdigits)
1014             {
1015               alloclhs = decRoundOperand (lhs, set, &status);
1016               if (alloclhs == NULL)
1017                 break;
1018               lhs = alloclhs;
1019             }
1020           /* rounding won't affect the result, but we might signal lostDigits */
1021           /* as well as the error for non-integer [x**y would need this too] */
1022           if (rhs->digits > reqdigits)
1023             {
1024               allocrhs = decRoundOperand (rhs, set, &status);
1025               if (allocrhs == NULL)
1026                 break;
1027               rhs = allocrhs;
1028             }
1029         }
1030 #endif
1031       /* [following code does not require input rounding] */
1032
1033       /* handle rhs Infinity */
1034       if (decNumberIsInfinite (rhs))
1035         {
1036           status |= DEC_Invalid_operation;      /* bad */
1037           break;
1038         }
1039       /* handle NaNs */
1040       if ((lhs->bits | rhs->bits) & (DECNAN | DECSNAN))
1041         {
1042           decNaNs (res, lhs, rhs, &status);
1043           break;
1044         }
1045
1046       /* Original rhs must be an integer that fits and is in range */
1047 #if DECSUBSET
1048       n = decGetInt (inrhs, set);
1049 #else
1050       n = decGetInt (inrhs);
1051 #endif
1052       if (n == BADINT || n > 999999999 || n < -999999999)
1053         {
1054           status |= DEC_Invalid_operation;
1055           break;
1056         }
1057       if (n < 0)
1058         {                       /* negative */
1059           n = -n;               /* use the absolute value */
1060         }
1061       if (decNumberIsNegative (lhs)     /* -x .. */
1062           && (n & 0x00000001))
1063         bits = DECNEG;          /* .. to an odd power */
1064
1065       /* handle LHS infinity */
1066       if (decNumberIsInfinite (lhs))
1067         {                       /* [NaNs already handled] */
1068           uByte rbits = rhs->bits;      /* save */
1069           decNumberZero (res);
1070           if (n == 0)
1071             *res->lsu = 1;      /* [-]Inf**0 => 1 */
1072           else
1073             {
1074               if (!(rbits & DECNEG))
1075                 bits |= DECINF; /* was not a **-n */
1076               /* [otherwise will be 0 or -0] */
1077               res->bits = bits;
1078             }
1079           break;
1080         }
1081
1082       /* clone the context */
1083       workset = *set;           /* copy all fields */
1084       /* calculate the working DIGITS */
1085       workset.digits = reqdigits + (inrhs->digits + inrhs->exponent) + 1;
1086       /* it's an error if this is more than we can handle */
1087       if (workset.digits > DECNUMMAXP)
1088         {
1089           status |= DEC_Invalid_operation;
1090           break;
1091         }
1092
1093       /* workset.digits is the count of digits for the accumulator we need */
1094       /* if accumulator is too long for local storage, then allocate */
1095       needbytes =
1096         sizeof (decNumber) + (D2U (workset.digits) - 1) * sizeof (Unit);
1097       /* [needbytes also used below if 1/lhs needed] */
1098       if (needbytes > sizeof (dacbuff))
1099         {
1100           allocdac = (decNumber *) malloc (needbytes);
1101           if (allocdac == NULL)
1102             {                   /* hopeless -- abandon */
1103               status |= DEC_Insufficient_storage;
1104               break;
1105             }
1106           dac = allocdac;       /* use the allocated space */
1107         }
1108       decNumberZero (dac);      /* acc=1 */
1109       *dac->lsu = 1;            /* .. */
1110
1111       if (n == 0)
1112         {                       /* x**0 is usually 1 */
1113           /* 0**0 is bad unless subset, when it becomes 1 */
1114           if (ISZERO (lhs)
1115 #if DECSUBSET
1116               && set->extended
1117 #endif
1118             )
1119             status |= DEC_Invalid_operation;
1120           else
1121             decNumberCopy (res, dac);   /* copy the 1 */
1122           break;
1123         }
1124
1125       /* if a negative power we'll need the constant 1, and if not subset */
1126       /* we'll invert the lhs now rather than inverting the result later */
1127       if (decNumberIsNegative (rhs))
1128         {                       /* was a **-n [hence digits>0] */
1129           decNumber * newlhs;
1130           decNumberCopy (&dnOne, dac);  /* dnOne=1;  [needed now or later] */
1131 #if DECSUBSET
1132           if (set->extended)
1133             {                   /* need to calculate 1/lhs */
1134 #endif
1135               /* divide lhs into 1, putting result in dac [dac=1/dac] */
1136               decDivideOp (dac, &dnOne, lhs, &workset, DIVIDE, &status);
1137               if (alloclhs != NULL)
1138                 {
1139                   free (alloclhs);      /* done with intermediate */
1140                   alloclhs = NULL;      /* indicate freed */
1141                 }
1142               /* now locate or allocate space for the inverted lhs */
1143               if (needbytes > sizeof (lhsbuff))
1144                 {
1145                   alloclhs = (decNumber *) malloc (needbytes);
1146                   if (alloclhs == NULL)
1147                     {           /* hopeless -- abandon */
1148                       status |= DEC_Insufficient_storage;
1149                       break;
1150                     }
1151                   newlhs = alloclhs;    /* use the allocated space */
1152                 }
1153               else
1154                 newlhs = (decNumber *) lhsbuff; /* use stack storage */
1155               /* [lhs now points to buffer or allocated storage] */
1156               decNumberCopy (newlhs, dac);      /* copy the 1/lhs */
1157               decNumberCopy (dac, &dnOne);      /* restore acc=1 */
1158               lhs = newlhs;
1159 #if DECSUBSET
1160             }
1161 #endif
1162         }
1163
1164       /* Raise-to-the-power loop... */
1165       seenbit = 0;              /* set once we've seen a 1-bit */
1166       for (i = 1;; i++)
1167         {                       /* for each bit [top bit ignored] */
1168           /* abandon if we have had overflow or terminal underflow */
1169           if (status & (DEC_Overflow | DEC_Underflow))
1170             {                   /* interesting? */
1171               if (status & DEC_Overflow || ISZERO (dac))
1172                 break;
1173             }
1174           /* [the following two lines revealed an optimizer bug in a C++ */
1175           /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
1176           n = n << 1;           /* move next bit to testable position */
1177           if (n < 0)
1178             {                   /* top bit is set */
1179               seenbit = 1;      /* OK, we're off */
1180               decMultiplyOp (dac, dac, lhs, &workset, &status); /* dac=dac*x */
1181             }
1182           if (i == 31)
1183             break;              /* that was the last bit */
1184           if (!seenbit)
1185             continue;           /* we don't have to square 1 */
1186           decMultiplyOp (dac, dac, dac, &workset, &status);     /* dac=dac*dac [square] */
1187         }                       /*i *//* 32 bits */
1188
1189       /* complete internal overflow or underflow processing */
1190       if (status & (DEC_Overflow | DEC_Subnormal))
1191         {
1192 #if DECSUBSET
1193           /* If subset, and power was negative, reverse the kind of -erflow */
1194           /* [1/x not yet done] */
1195           if (!set->extended && decNumberIsNegative (rhs))
1196             {
1197               if (status & DEC_Overflow)
1198                 status ^= DEC_Overflow | DEC_Underflow | DEC_Subnormal;
1199               else
1200                 {               /* trickier -- Underflow may or may not be set */
1201                   status &= ~(DEC_Underflow | DEC_Subnormal);   /* [one or both] */
1202                   status |= DEC_Overflow;
1203                 }
1204             }
1205 #endif
1206           dac->bits = (dac->bits & ~DECNEG) | bits;     /* force correct sign */
1207           /* round subnormals [to set.digits rather than workset.digits] */
1208           /* or set overflow result similarly as required */
1209           decFinalize (dac, set, &residue, &status);
1210           decNumberCopy (res, dac);     /* copy to result (is now OK length) */
1211           break;
1212         }
1213
1214 #if DECSUBSET
1215       if (!set->extended &&     /* subset math */
1216           decNumberIsNegative (rhs))
1217         {                       /* was a **-n [hence digits>0] */
1218           /* so divide result into 1 [dac=1/dac] */
1219           decDivideOp (dac, &dnOne, dac, &workset, DIVIDE, &status);
1220         }
1221 #endif
1222
1223       /* reduce result to the requested length and copy to result */
1224       decCopyFit (res, dac, set, &residue, &status);
1225       decFinish (res, set, &residue, &status);  /* final cleanup */
1226 #if DECSUBSET
1227       if (!set->extended)
1228         decTrim (res, 0, &dropped);     /* trailing zeros */
1229 #endif
1230     }
1231   while (0);                    /* end protected */
1232
1233   if (allocdac != NULL)
1234     free (allocdac);            /* drop any storage we used */
1235   if (allocrhs != NULL)
1236     free (allocrhs);            /* .. */
1237   if (alloclhs != NULL)
1238     free (alloclhs);            /* .. */
1239   if (status != 0)
1240     decStatus (res, status, set);
1241   return res;
1242 }
1243
1244 /* ------------------------------------------------------------------ */
1245 /* decNumberQuantize -- force exponent to requested value             */
1246 /*                                                                    */
1247 /*   This computes C = op(A, B), where op adjusts the coefficient     */
1248 /*   of C (by rounding or shifting) such that the exponent (-scale)   */
1249 /*   of C has exponent of B.  The numerical value of C will equal A,  */
1250 /*   except for the effects of any rounding that occurred.            */
1251 /*                                                                    */
1252 /*   res is C, the result.  C may be A or B                           */
1253 /*   lhs is A, the number to adjust                                   */
1254 /*   rhs is B, the number with exponent to match                      */
1255 /*   set is the context                                               */
1256 /*                                                                    */
1257 /* C must have space for set->digits digits.                          */
1258 /*                                                                    */
1259 /* Unless there is an error or the result is infinite, the exponent   */
1260 /* after the operation is guaranteed to be equal to that of B.        */
1261 /* ------------------------------------------------------------------ */
1262 decNumber *
1263 decNumberQuantize (decNumber * res, const decNumber * lhs,
1264                    const decNumber * rhs, decContext * set)
1265 {
1266   uInt status = 0;              /* accumulator */
1267   decQuantizeOp (res, lhs, rhs, set, 1, &status);
1268   if (status != 0)
1269     decStatus (res, status, set);
1270   return res;
1271 }
1272
1273 /* ------------------------------------------------------------------ */
1274 /* decNumberRescale -- force exponent to requested value              */
1275 /*                                                                    */
1276 /*   This computes C = op(A, B), where op adjusts the coefficient     */
1277 /*   of C (by rounding or shifting) such that the exponent (-scale)   */
1278 /*   of C has the value B.  The numerical value of C will equal A,    */
1279 /*   except for the effects of any rounding that occurred.            */
1280 /*                                                                    */
1281 /*   res is C, the result.  C may be A or B                           */
1282 /*   lhs is A, the number to adjust                                   */
1283 /*   rhs is B, the requested exponent                                 */
1284 /*   set is the context                                               */
1285 /*                                                                    */
1286 /* C must have space for set->digits digits.                          */
1287 /*                                                                    */
1288 /* Unless there is an error or the result is infinite, the exponent   */
1289 /* after the operation is guaranteed to be equal to B.                */
1290 /* ------------------------------------------------------------------ */
1291 decNumber *
1292 decNumberRescale (decNumber * res, const decNumber * lhs,
1293                   const decNumber * rhs, decContext * set)
1294 {
1295   uInt status = 0;              /* accumulator */
1296   decQuantizeOp (res, lhs, rhs, set, 0, &status);
1297   if (status != 0)
1298     decStatus (res, status, set);
1299   return res;
1300 }
1301
1302 /* ------------------------------------------------------------------ */
1303 /* decNumberRemainder -- divide and return remainder                  */
1304 /*                                                                    */
1305 /*   This computes C = A % B                                          */
1306 /*                                                                    */
1307 /*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
1308 /*   lhs is A                                                         */
1309 /*   rhs is B                                                         */
1310 /*   set is the context                                               */
1311 /*                                                                    */
1312 /* C must have space for set->digits digits.                          */
1313 /* ------------------------------------------------------------------ */
1314 decNumber *
1315 decNumberRemainder (decNumber * res, const decNumber * lhs,
1316                     const decNumber * rhs, decContext * set)
1317 {
1318   uInt status = 0;              /* accumulator */
1319   decDivideOp (res, lhs, rhs, set, REMAINDER, &status);
1320   if (status != 0)
1321     decStatus (res, status, set);
1322   return res;
1323 }
1324
1325 /* ------------------------------------------------------------------ */
1326 /* decNumberRemainderNear -- divide and return remainder from nearest */
1327 /*                                                                    */
1328 /*   This computes C = A % B, where % is the IEEE remainder operator  */
1329 /*                                                                    */
1330 /*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
1331 /*   lhs is A                                                         */
1332 /*   rhs is B                                                         */
1333 /*   set is the context                                               */
1334 /*                                                                    */
1335 /* C must have space for set->digits digits.                          */
1336 /* ------------------------------------------------------------------ */
1337 decNumber *
1338 decNumberRemainderNear (decNumber * res, const decNumber * lhs,
1339                         const decNumber * rhs, decContext * set)
1340 {
1341   uInt status = 0;              /* accumulator */
1342   decDivideOp (res, lhs, rhs, set, REMNEAR, &status);
1343   if (status != 0)
1344     decStatus (res, status, set);
1345   return res;
1346 }
1347
1348 /* ------------------------------------------------------------------ */
1349 /* decNumberSameQuantum -- test for equal exponents                   */
1350 /*                                                                    */
1351 /*   res is the result number, which will contain either 0 or 1       */
1352 /*   lhs is a number to test                                          */
1353 /*   rhs is the second (usually a pattern)                            */
1354 /*                                                                    */
1355 /* No errors are possible and no context is needed.                   */
1356 /* ------------------------------------------------------------------ */
1357 decNumber *
1358 decNumberSameQuantum (decNumber * res, const decNumber * lhs, const decNumber * rhs)
1359 {
1360   uByte merged;                 /* merged flags */
1361   Unit ret = 0;                 /* return value */
1362
1363 #if DECCHECK
1364   if (decCheckOperands (res, lhs, rhs, DECUNUSED))
1365     return res;
1366 #endif
1367
1368   merged = (lhs->bits | rhs->bits) & DECSPECIAL;
1369   if (merged)
1370     {
1371       if (decNumberIsNaN (lhs) && decNumberIsNaN (rhs))
1372         ret = 1;
1373       else if (decNumberIsInfinite (lhs) && decNumberIsInfinite (rhs))
1374         ret = 1;
1375       /* [anything else with a special gives 0] */
1376     }
1377   else if (lhs->exponent == rhs->exponent)
1378     ret = 1;
1379
1380   decNumberZero (res);          /* OK to overwrite an operand */
1381   *res->lsu = ret;
1382   return res;
1383 }
1384
1385 /* ------------------------------------------------------------------ */
1386 /* decNumberSquareRoot -- square root operator                        */
1387 /*                                                                    */
1388 /*   This computes C = squareroot(A)                                  */
1389 /*                                                                    */
1390 /*   res is C, the result.  C may be A                                */
1391 /*   rhs is A                                                         */
1392 /*   set is the context; note that rounding mode has no effect        */
1393 /*                                                                    */
1394 /* C must have space for set->digits digits.                          */
1395 /* ------------------------------------------------------------------ */
1396 /* This uses the following varying-precision algorithm in:            */
1397 /*                                                                    */
1398 /*   Properly Rounded Variable Precision Square Root, T. E. Hull and  */
1399 /*   A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
1400 /*   pp229-237, ACM, September 1985.                                  */
1401 /*                                                                    */
1402 /* % [Reformatted original Numerical Turing source code follows.]     */
1403 /* function sqrt(x : real) : real                                     */
1404 /* % sqrt(x) returns the properly rounded approximation to the square */
1405 /* % root of x, in the precision of the calling environment, or it    */
1406 /* % fails if x < 0.                                                  */
1407 /* % t e hull and a abrham, august, 1984                              */
1408 /* if x <= 0 then                                                     */
1409 /*   if x < 0 then                                                    */
1410 /*     assert false                                                   */
1411 /*   else                                                             */
1412 /*     result 0                                                       */
1413 /*   end if                                                           */
1414 /* end if                                                             */
1415 /* var f := setexp(x, 0)  % fraction part of x   [0.1 <= x < 1]       */
1416 /* var e := getexp(x)     % exponent part of x                        */
1417 /* var approx : real                                                  */
1418 /* if e mod 2 = 0  then                                               */
1419 /*   approx := .259 + .819 * f   % approx to root of f                */
1420 /* else                                                               */
1421 /*   f := f/l0                   % adjustments                        */
1422 /*   e := e + 1                  %   for odd                          */
1423 /*   approx := .0819 + 2.59 * f  %   exponent                         */
1424 /* end if                                                             */
1425 /*                                                                    */
1426 /* var p:= 3                                                          */
1427 /* const maxp := currentprecision + 2                                 */
1428 /* loop                                                               */
1429 /*   p := min(2*p - 2, maxp)     % p = 4,6,10, . . . , maxp           */
1430 /*   precision p                                                      */
1431 /*   approx := .5 * (approx + f/approx)                               */
1432 /*   exit when p = maxp                                               */
1433 /* end loop                                                           */
1434 /*                                                                    */
1435 /* % approx is now within 1 ulp of the properly rounded square root   */
1436 /* % of f; to ensure proper rounding, compare squares of (approx -    */
1437 /* % l/2 ulp) and (approx + l/2 ulp) with f.                          */
1438 /* p := currentprecision                                              */
1439 /* begin                                                              */
1440 /*   precision p + 2                                                  */
1441 /*   const approxsubhalf := approx - setexp(.5, -p)                   */
1442 /*   if mulru(approxsubhalf, approxsubhalf) > f then                  */
1443 /*     approx := approx - setexp(.l, -p + 1)                          */
1444 /*   else                                                             */
1445 /*     const approxaddhalf := approx + setexp(.5, -p)                 */
1446 /*     if mulrd(approxaddhalf, approxaddhalf) < f then                */
1447 /*       approx := approx + setexp(.l, -p + 1)                        */
1448 /*     end if                                                         */
1449 /*   end if                                                           */
1450 /* end                                                                */
1451 /* result setexp(approx, e div 2)  % fix exponent                     */
1452 /* end sqrt                                                           */
1453 /* ------------------------------------------------------------------ */
1454 decNumber *
1455 decNumberSquareRoot (decNumber * res, const decNumber * rhs, decContext * set)
1456 {
1457   decContext workset, approxset;        /* work contexts */
1458   decNumber dzero;              /* used for constant zero */
1459   Int maxp = set->digits + 2;   /* largest working precision */
1460   Int residue = 0;              /* rounding residue */
1461   uInt status = 0, ignore = 0;  /* status accumulators */
1462   Int exp;                      /* working exponent */
1463   Int ideal;                    /* ideal (preferred) exponent */
1464   uInt needbytes;               /* work */
1465   Int dropped;                  /* .. */
1466
1467   decNumber *allocrhs = NULL;   /* non-NULL if rounded rhs allocated */
1468   /* buffer for f [needs +1 in case DECBUFFER 0] */
1469   uByte buff[sizeof (decNumber) + (D2U (DECBUFFER + 1) - 1) * sizeof (Unit)];
1470   /* buffer for a [needs +2 to match maxp] */
1471   uByte bufa[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1472   /* buffer for temporary, b [must be same size as a] */
1473   uByte bufb[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1474   decNumber *allocbuff = NULL;  /* -> allocated buff, iff allocated */
1475   decNumber *allocbufa = NULL;  /* -> allocated bufa, iff allocated */
1476   decNumber *allocbufb = NULL;  /* -> allocated bufb, iff allocated */
1477   decNumber *f = (decNumber *) buff;    /* reduced fraction */
1478   decNumber *a = (decNumber *) bufa;    /* approximation to result */
1479   decNumber *b = (decNumber *) bufb;    /* intermediate result */
1480   /* buffer for temporary variable, up to 3 digits */
1481   uByte buft[sizeof (decNumber) + (D2U (3) - 1) * sizeof (Unit)];
1482   decNumber *t = (decNumber *) buft;    /* up-to-3-digit constant or work */
1483
1484 #if DECCHECK
1485   if (decCheckOperands (res, DECUNUSED, rhs, set))
1486     return res;
1487 #endif
1488
1489   do
1490     {                           /* protect allocated storage */
1491 #if DECSUBSET
1492       if (!set->extended)
1493         {
1494           /* reduce operand and set lostDigits status, as needed */
1495           if (rhs->digits > set->digits)
1496             {
1497               allocrhs = decRoundOperand (rhs, set, &status);
1498               if (allocrhs == NULL)
1499                 break;
1500               /* [Note: 'f' allocation below could reuse this buffer if */
1501               /* used, but as this is rare we keep them separate for clarity.] */
1502               rhs = allocrhs;
1503             }
1504         }
1505 #endif
1506       /* [following code does not require input rounding] */
1507
1508       /* handle infinities and NaNs */
1509       if (rhs->bits & DECSPECIAL)
1510         {
1511           if (decNumberIsInfinite (rhs))
1512             {                   /* an infinity */
1513               if (decNumberIsNegative (rhs))
1514                 status |= DEC_Invalid_operation;
1515               else
1516                 decNumberCopy (res, rhs);       /* +Infinity */
1517             }
1518           else
1519             decNaNs (res, rhs, NULL, &status);  /* a NaN */
1520           break;
1521         }
1522
1523       /* calculate the ideal (preferred) exponent [floor(exp/2)] */
1524       /* [We would like to write: ideal=rhs->exponent>>1, but this */
1525       /* generates a compiler warning.  Generated code is the same.] */
1526       ideal = (rhs->exponent & ~1) / 2; /* target */
1527
1528       /* handle zeros */
1529       if (ISZERO (rhs))
1530         {
1531           decNumberCopy (res, rhs);     /* could be 0 or -0 */
1532           res->exponent = ideal;        /* use the ideal [safe] */
1533           break;
1534         }
1535
1536       /* any other -x is an oops */
1537       if (decNumberIsNegative (rhs))
1538         {
1539           status |= DEC_Invalid_operation;
1540           break;
1541         }
1542
1543       /* we need space for three working variables */
1544       /*   f -- the same precision as the RHS, reduced to 0.01->0.99... */
1545       /*   a -- Hull's approx -- precision, when assigned, is */
1546       /*        currentprecision (we allow +2 for use as temporary) */
1547       /*   b -- intermediate temporary result */
1548       /* if any is too long for local storage, then allocate */
1549       needbytes =
1550         sizeof (decNumber) + (D2U (rhs->digits) - 1) * sizeof (Unit);
1551       if (needbytes > sizeof (buff))
1552         {
1553           allocbuff = (decNumber *) malloc (needbytes);
1554           if (allocbuff == NULL)
1555             {                   /* hopeless -- abandon */
1556               status |= DEC_Insufficient_storage;
1557               break;
1558             }
1559           f = allocbuff;        /* use the allocated space */
1560         }
1561       /* a and b both need to be able to hold a maxp-length number */
1562       needbytes = sizeof (decNumber) + (D2U (maxp) - 1) * sizeof (Unit);
1563       if (needbytes > sizeof (bufa))
1564         {                       /* [same applies to b] */
1565           allocbufa = (decNumber *) malloc (needbytes);
1566           allocbufb = (decNumber *) malloc (needbytes);
1567           if (allocbufa == NULL || allocbufb == NULL)
1568             {                   /* hopeless */
1569               status |= DEC_Insufficient_storage;
1570               break;
1571             }
1572           a = allocbufa;        /* use the allocated space */
1573           b = allocbufb;        /* .. */
1574         }
1575
1576       /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
1577       decNumberCopy (f, rhs);
1578       exp = f->exponent + f->digits;    /* adjusted to Hull rules */
1579       f->exponent = -(f->digits);       /* to range */
1580
1581       /* set up working contexts (the second is used for Numerical */
1582       /* Turing assignment) */
1583       decContextDefault (&workset, DEC_INIT_DECIMAL64);
1584       decContextDefault (&approxset, DEC_INIT_DECIMAL64);
1585       approxset.digits = set->digits;   /* approx's length */
1586
1587       /* [Until further notice, no error is possible and status bits */
1588       /* (Rounded, etc.) should be ignored, not accumulated.] */
1589
1590       /* Calculate initial approximation, and allow for odd exponent */
1591       workset.digits = set->digits;     /* p for initial calculation */
1592       t->bits = 0;
1593       t->digits = 3;
1594       a->bits = 0;
1595       a->digits = 3;
1596       if ((exp & 1) == 0)
1597         {                       /* even exponent */
1598           /* Set t=0.259, a=0.819 */
1599           t->exponent = -3;
1600           a->exponent = -3;
1601 #if DECDPUN>=3
1602           t->lsu[0] = 259;
1603           a->lsu[0] = 819;
1604 #elif DECDPUN==2
1605           t->lsu[0] = 59;
1606           t->lsu[1] = 2;
1607           a->lsu[0] = 19;
1608           a->lsu[1] = 8;
1609 #else
1610           t->lsu[0] = 9;
1611           t->lsu[1] = 5;
1612           t->lsu[2] = 2;
1613           a->lsu[0] = 9;
1614           a->lsu[1] = 1;
1615           a->lsu[2] = 8;
1616 #endif
1617         }
1618       else
1619         {                       /* odd exponent */
1620           /* Set t=0.0819, a=2.59 */
1621           f->exponent--;        /* f=f/10 */
1622           exp++;                /* e=e+1 */
1623           t->exponent = -4;
1624           a->exponent = -2;
1625 #if DECDPUN>=3
1626           t->lsu[0] = 819;
1627           a->lsu[0] = 259;
1628 #elif DECDPUN==2
1629           t->lsu[0] = 19;
1630           t->lsu[1] = 8;
1631           a->lsu[0] = 59;
1632           a->lsu[1] = 2;
1633 #else
1634           t->lsu[0] = 9;
1635           t->lsu[1] = 1;
1636           t->lsu[2] = 8;
1637           a->lsu[0] = 9;
1638           a->lsu[1] = 5;
1639           a->lsu[2] = 2;
1640 #endif
1641         }
1642       decMultiplyOp (a, a, f, &workset, &ignore);       /* a=a*f */
1643       decAddOp (a, a, t, &workset, 0, &ignore); /* ..+t */
1644       /* [a is now the initial approximation for sqrt(f), calculated with */
1645       /* currentprecision, which is also a's precision.] */
1646
1647       /* the main calculation loop */
1648       decNumberZero (&dzero);   /* make 0 */
1649       decNumberZero (t);        /* set t = 0.5 */
1650       t->lsu[0] = 5;            /* .. */
1651       t->exponent = -1;         /* .. */
1652       workset.digits = 3;       /* initial p */
1653       for (;;)
1654         {
1655           /* set p to min(2*p - 2, maxp)  [hence 3; or: 4, 6, 10, ... , maxp] */
1656           workset.digits = workset.digits * 2 - 2;
1657           if (workset.digits > maxp)
1658             workset.digits = maxp;
1659           /* a = 0.5 * (a + f/a) */
1660           /* [calculated at p then rounded to currentprecision] */
1661           decDivideOp (b, f, a, &workset, DIVIDE, &ignore);     /* b=f/a */
1662           decAddOp (b, b, a, &workset, 0, &ignore);     /* b=b+a */
1663           decMultiplyOp (a, b, t, &workset, &ignore);   /* a=b*0.5 */
1664           /* assign to approx [round to length] */
1665           decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1666           if (workset.digits == maxp)
1667             break;              /* just did final */
1668         }                       /* loop */
1669
1670       /* a is now at currentprecision and within 1 ulp of the properly */
1671       /* rounded square root of f; to ensure proper rounding, compare */
1672       /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
1673       /* Here workset.digits=maxp and t=0.5 */
1674       workset.digits--;         /* maxp-1 is OK now */
1675       t->exponent = -set->digits - 1;   /* make 0.5 ulp */
1676       decNumberCopy (b, a);
1677       decAddOp (b, b, t, &workset, DECNEG, &ignore);    /* b = a - 0.5 ulp */
1678       workset.round = DEC_ROUND_UP;
1679       decMultiplyOp (b, b, b, &workset, &ignore);       /* b = mulru(b, b) */
1680       decCompareOp (b, f, b, &workset, COMPARE, &ignore);       /* b ? f, reversed */
1681       if (decNumberIsNegative (b))
1682         {                       /* f < b [i.e., b > f] */
1683           /* this is the more common adjustment, though both are rare */
1684           t->exponent++;        /* make 1.0 ulp */
1685           t->lsu[0] = 1;        /* .. */
1686           decAddOp (a, a, t, &workset, DECNEG, &ignore);        /* a = a - 1 ulp */
1687           /* assign to approx [round to length] */
1688           decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1689         }
1690       else
1691         {
1692           decNumberCopy (b, a);
1693           decAddOp (b, b, t, &workset, 0, &ignore);     /* b = a + 0.5 ulp */
1694           workset.round = DEC_ROUND_DOWN;
1695           decMultiplyOp (b, b, b, &workset, &ignore);   /* b = mulrd(b, b) */
1696           decCompareOp (b, b, f, &workset, COMPARE, &ignore);   /* b ? f */
1697           if (decNumberIsNegative (b))
1698             {                   /* b < f */
1699               t->exponent++;    /* make 1.0 ulp */
1700               t->lsu[0] = 1;    /* .. */
1701               decAddOp (a, a, t, &workset, 0, &ignore); /* a = a + 1 ulp */
1702               /* assign to approx [round to length] */
1703               decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1704             }
1705         }
1706       /* [no errors are possible in the above, and rounding/inexact during */
1707       /* estimation are irrelevant, so status was not accumulated] */
1708
1709       /* Here, 0.1 <= a < 1  [Hull] */
1710       a->exponent += exp / 2;   /* set correct exponent */
1711
1712       /* Process Subnormals */
1713       decFinalize (a, set, &residue, &status);
1714
1715       /* count dropable zeros [after any subnormal rounding] */
1716       decNumberCopy (b, a);
1717       decTrim (b, 1, &dropped); /* [drops trailing zeros] */
1718
1719       /* Finally set Inexact and Rounded.  The answer can only be exact if */
1720       /* it is short enough so that squaring it could fit in set->digits, */
1721       /* so this is the only (relatively rare) time we have to check */
1722       /* carefully */
1723       if (b->digits * 2 - 1 > set->digits)
1724         {                       /* cannot fit */
1725           status |= DEC_Inexact | DEC_Rounded;
1726         }
1727       else
1728         {                       /* could be exact/unrounded */
1729           uInt mstatus = 0;     /* local status */
1730           decMultiplyOp (b, b, b, &workset, &mstatus);  /* try the multiply */
1731           if (mstatus != 0)
1732             {                   /* result won't fit */
1733               status |= DEC_Inexact | DEC_Rounded;
1734             }
1735           else
1736             {                   /* plausible */
1737               decCompareOp (t, b, rhs, &workset, COMPARE, &mstatus);    /* b ? rhs */
1738               if (!ISZERO (t))
1739                 {
1740                   status |= DEC_Inexact | DEC_Rounded;
1741                 }
1742               else
1743                 {               /* is Exact */
1744                   /* here, dropped is the count of trailing zeros in 'a' */
1745                   /* use closest exponent to ideal... */
1746                   Int todrop = ideal - a->exponent;     /* most we can drop */
1747
1748                   if (todrop < 0)
1749                     {           /* ideally would add 0s */
1750                       status |= DEC_Rounded;
1751                     }
1752                   else
1753                     {           /* unrounded */
1754                       if (dropped < todrop)
1755                         todrop = dropped;       /* clamp to those available */
1756                       if (todrop > 0)
1757                         {       /* OK, some to drop */
1758                           decShiftToLeast (a->lsu, D2U (a->digits), todrop);
1759                           a->exponent += todrop;        /* maintain numerical value */
1760                           a->digits -= todrop;  /* new length */
1761                         }
1762                     }
1763                 }
1764             }
1765         }
1766       decNumberCopy (res, a);   /* assume this is the result */
1767     }
1768   while (0);                    /* end protected */
1769
1770   if (allocbuff != NULL)
1771     free (allocbuff);           /* drop any storage we used */
1772   if (allocbufa != NULL)
1773     free (allocbufa);           /* .. */
1774   if (allocbufb != NULL)
1775     free (allocbufb);           /* .. */
1776   if (allocrhs != NULL)
1777     free (allocrhs);            /* .. */
1778   if (status != 0)
1779     decStatus (res, status, set);       /* then report status */
1780   return res;
1781 }
1782
1783 /* ------------------------------------------------------------------ */
1784 /* decNumberSubtract -- subtract two Numbers                          */
1785 /*                                                                    */
1786 /*   This computes C = A - B                                          */
1787 /*                                                                    */
1788 /*   res is C, the result.  C may be A and/or B (e.g., X=X-X)         */
1789 /*   lhs is A                                                         */
1790 /*   rhs is B                                                         */
1791 /*   set is the context                                               */
1792 /*                                                                    */
1793 /* C must have space for set->digits digits.                          */
1794 /* ------------------------------------------------------------------ */
1795 decNumber *
1796 decNumberSubtract (decNumber * res, const decNumber * lhs,
1797                    const decNumber * rhs, decContext * set)
1798 {
1799   uInt status = 0;              /* accumulator */
1800
1801   decAddOp (res, lhs, rhs, set, DECNEG, &status);
1802   if (status != 0)
1803     decStatus (res, status, set);
1804   return res;
1805 }
1806
1807 /* ------------------------------------------------------------------ */
1808 /* decNumberToIntegralValue -- round-to-integral-value                */
1809 /*                                                                    */
1810 /*   res is the result                                                */
1811 /*   rhs is input number                                              */
1812 /*   set is the context                                               */
1813 /*                                                                    */
1814 /* res must have space for any value of rhs.                          */
1815 /*                                                                    */
1816 /* This implements the IEEE special operator and therefore treats     */
1817 /* special values as valid, and also never sets Inexact.  For finite  */
1818 /* numbers it returns rescale(rhs, 0) if rhs->exponent is <0.         */
1819 /* Otherwise the result is rhs (so no error is possible).             */
1820 /*                                                                    */
1821 /* The context is used for rounding mode and status after sNaN, but   */
1822 /* the digits setting is ignored.                                     */
1823 /* ------------------------------------------------------------------ */
1824 decNumber *
1825 decNumberToIntegralValue (decNumber * res, const decNumber * rhs, decContext * set)
1826 {
1827   decNumber dn;
1828   decContext workset;           /* working context */
1829
1830 #if DECCHECK
1831   if (decCheckOperands (res, DECUNUSED, rhs, set))
1832     return res;
1833 #endif
1834
1835   /* handle infinities and NaNs */
1836   if (rhs->bits & DECSPECIAL)
1837     {
1838       uInt status = 0;
1839       if (decNumberIsInfinite (rhs))
1840         decNumberCopy (res, rhs);       /* an Infinity */
1841       else
1842         decNaNs (res, rhs, NULL, &status);      /* a NaN */
1843       if (status != 0)
1844         decStatus (res, status, set);
1845       return res;
1846     }
1847
1848   /* we have a finite number; no error possible */
1849   if (rhs->exponent >= 0)
1850     return decNumberCopy (res, rhs);
1851   /* that was easy, but if negative exponent we have work to do... */
1852   workset = *set;               /* clone rounding, etc. */
1853   workset.digits = rhs->digits; /* no length rounding */
1854   workset.traps = 0;            /* no traps */
1855   decNumberZero (&dn);          /* make a number with exponent 0 */
1856   return decNumberQuantize (res, rhs, &dn, &workset);
1857 }
1858
1859 /* ================================================================== */
1860 /* Utility routines                                                   */
1861 /* ================================================================== */
1862
1863 /* ------------------------------------------------------------------ */
1864 /* decNumberCopy -- copy a number                                     */
1865 /*                                                                    */
1866 /*   dest is the target decNumber                                     */
1867 /*   src  is the source decNumber                                     */
1868 /*   returns dest                                                     */
1869 /*                                                                    */
1870 /* (dest==src is allowed and is a no-op)                              */
1871 /* All fields are updated as required.  This is a utility operation,  */
1872 /* so special values are unchanged and no error is possible.          */
1873 /* ------------------------------------------------------------------ */
1874 decNumber *
1875 decNumberCopy (decNumber * dest, const decNumber * src)
1876 {
1877
1878 #if DECCHECK
1879   if (src == NULL)
1880     return decNumberZero (dest);
1881 #endif
1882
1883   if (dest == src)
1884     return dest;                /* no copy required */
1885
1886   /* We use explicit assignments here as structure assignment can copy */
1887   /* more than just the lsu (for small DECDPUN).  This would not affect */
1888   /* the value of the results, but would disturb test harness spill */
1889   /* checking. */
1890   dest->bits = src->bits;
1891   dest->exponent = src->exponent;
1892   dest->digits = src->digits;
1893   dest->lsu[0] = src->lsu[0];
1894   if (src->digits > DECDPUN)
1895     {                           /* more Units to come */
1896       Unit *d;                  /* work */
1897       const Unit *s, *smsup;    /* work */
1898       /* memcpy for the remaining Units would be safe as they cannot */
1899       /* overlap.  However, this explicit loop is faster in short cases. */
1900       d = dest->lsu + 1;        /* -> first destination */
1901       smsup = src->lsu + D2U (src->digits);     /* -> source msu+1 */
1902       for (s = src->lsu + 1; s < smsup; s++, d++)
1903         *d = *s;
1904     }
1905   return dest;
1906 }
1907
1908 /* ------------------------------------------------------------------ */
1909 /* decNumberTrim -- remove insignificant zeros                        */
1910 /*                                                                    */
1911 /*   dn is the number to trim                                         */
1912 /*   returns dn                                                       */
1913 /*                                                                    */
1914 /* All fields are updated as required.  This is a utility operation,  */
1915 /* so special values are unchanged and no error is possible.          */
1916 /* ------------------------------------------------------------------ */
1917 decNumber *
1918 decNumberTrim (decNumber * dn)
1919 {
1920   Int dropped;                  /* work */
1921   return decTrim (dn, 0, &dropped);
1922 }
1923
1924 /* ------------------------------------------------------------------ */
1925 /* decNumberVersion -- return the name and version of this module     */
1926 /*                                                                    */
1927 /* No error is possible.                                              */
1928 /* ------------------------------------------------------------------ */
1929 const char *
1930 decNumberVersion (void)
1931 {
1932   return DECVERSION;
1933 }
1934
1935 /* ------------------------------------------------------------------ */
1936 /* decNumberZero -- set a number to 0                                 */
1937 /*                                                                    */
1938 /*   dn is the number to set, with space for one digit                */
1939 /*   returns dn                                                       */
1940 /*                                                                    */
1941 /* No error is possible.                                              */
1942 /* ------------------------------------------------------------------ */
1943 /* Memset is not used as it is much slower in some environments. */
1944 decNumber *
1945 decNumberZero (decNumber * dn)
1946 {
1947
1948 #if DECCHECK
1949   if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
1950     return dn;
1951 #endif
1952
1953   dn->bits = 0;
1954   dn->exponent = 0;
1955   dn->digits = 1;
1956   dn->lsu[0] = 0;
1957   return dn;
1958 }
1959
1960 /* ================================================================== */
1961 /* Local routines                                                     */
1962 /* ================================================================== */
1963
1964 /* ------------------------------------------------------------------ */
1965 /* decToString -- lay out a number into a string                      */
1966 /*                                                                    */
1967 /*   dn     is the number to lay out                                  */
1968 /*   string is where to lay out the number                            */
1969 /*   eng    is 1 if Engineering, 0 if Scientific                      */
1970 /*                                                                    */
1971 /* str must be at least dn->digits+14 characters long                 */
1972 /* No error is possible.                                              */
1973 /*                                                                    */
1974 /* Note that this routine can generate a -0 or 0.000.  These are      */
1975 /* never generated in subset to-number or arithmetic, but can occur   */
1976 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234).              */
1977 /* ------------------------------------------------------------------ */
1978 /* If DECCHECK is enabled the string "?" is returned if a number is */
1979 /* invalid. */
1980
1981 /* TODIGIT -- macro to remove the leading digit from the unsigned */
1982 /* integer u at column cut (counting from the right, LSD=0) and place */
1983 /* it as an ASCII character into the character pointed to by c.  Note */
1984 /* that cut must be <= 9, and the maximum value for u is 2,000,000,000 */
1985 /* (as is needed for negative exponents of subnormals).  The unsigned */
1986 /* integer pow is used as a temporary variable. */
1987 #define TODIGIT(u, cut, c) {            \
1988   *(c)='0';                             \
1989   pow=powers[cut]*2;                    \
1990   if ((u)>pow) {                        \
1991     pow*=4;                             \
1992     if ((u)>=pow) {(u)-=pow; *(c)+=8;}  \
1993     pow/=2;                             \
1994     if ((u)>=pow) {(u)-=pow; *(c)+=4;}  \
1995     pow/=2;                             \
1996     }                                   \
1997   if ((u)>=pow) {(u)-=pow; *(c)+=2;}    \
1998   pow/=2;                               \
1999   if ((u)>=pow) {(u)-=pow; *(c)+=1;}    \
2000   }
2001
2002 static void
2003 decToString (const decNumber * dn, char *string, Flag eng)
2004 {
2005   Int exp = dn->exponent;       /* local copy */
2006   Int e;                        /* E-part value */
2007   Int pre;                      /* digits before the '.' */
2008   Int cut;                      /* for counting digits in a Unit */
2009   char *c = string;             /* work [output pointer] */
2010   const Unit *up = dn->lsu + D2U (dn->digits) - 1;      /* -> msu [input pointer] */
2011   uInt u, pow;                  /* work */
2012
2013 #if DECCHECK
2014   if (decCheckOperands (DECUNUSED, dn, DECUNUSED, DECUNUSED))
2015     {
2016       strcpy (string, "?");
2017       return;
2018     }
2019 #endif
2020
2021   if (decNumberIsNegative (dn))
2022     {                           /* Negatives get a minus (except */
2023       *c = '-';                 /* NaNs, which remove the '-' below) */
2024       c++;
2025     }
2026   if (dn->bits & DECSPECIAL)
2027     {                           /* Is a special value */
2028       if (decNumberIsInfinite (dn))
2029         {
2030           strcpy (c, "Infinity");
2031           return;
2032         }
2033       /* a NaN */
2034       if (dn->bits & DECSNAN)
2035         {                       /* signalling NaN */
2036           *c = 's';
2037           c++;
2038         }
2039       strcpy (c, "NaN");
2040       c += 3;                   /* step past */
2041       /* if not a clean non-zero coefficient, that's all we have in a */
2042       /* NaN string */
2043       if (exp != 0 || (*dn->lsu == 0 && dn->digits == 1))
2044         return;
2045       /* [drop through to add integer] */
2046     }
2047
2048   /* calculate how many digits in msu, and hence first cut */
2049   cut = dn->digits % DECDPUN;
2050   if (cut == 0)
2051     cut = DECDPUN;              /* msu is full */
2052   cut--;                        /* power of ten for digit */
2053
2054   if (exp == 0)
2055     {                           /* simple integer [common fastpath, */
2056       /*   used for NaNs, too] */
2057       for (; up >= dn->lsu; up--)
2058         {                       /* each Unit from msu */
2059           u = *up;              /* contains DECDPUN digits to lay out */
2060           for (; cut >= 0; c++, cut--)
2061             TODIGIT (u, cut, c);
2062           cut = DECDPUN - 1;    /* next Unit has all digits */
2063         }
2064       *c = '\0';                /* terminate the string */
2065       return;
2066     }
2067
2068   /* non-0 exponent -- assume plain form */
2069   pre = dn->digits + exp;       /* digits before '.' */
2070   e = 0;                        /* no E */
2071   if ((exp > 0) || (pre < -5))
2072     {                           /* need exponential form */
2073       e = exp + dn->digits - 1; /* calculate E value */
2074       pre = 1;                  /* assume one digit before '.' */
2075       if (eng && (e != 0))
2076         {                       /* may need to adjust */
2077           Int adj;              /* adjustment */
2078           /* The C remainder operator is undefined for negative numbers, so */
2079           /* we must use positive remainder calculation here */
2080           if (e < 0)
2081             {
2082               adj = (-e) % 3;
2083               if (adj != 0)
2084                 adj = 3 - adj;
2085             }
2086           else
2087             {                   /* e>0 */
2088               adj = e % 3;
2089             }
2090           e = e - adj;
2091           /* if we are dealing with zero we will use exponent which is a */
2092           /* multiple of three, as expected, but there will only be the */
2093           /* one zero before the E, still.  Otherwise note the padding. */
2094           if (!ISZERO (dn))
2095             pre += adj;
2096           else
2097             {                   /* is zero */
2098               if (adj != 0)
2099                 {               /* 0.00Esnn needed */
2100                   e = e + 3;
2101                   pre = -(2 - adj);
2102                 }
2103             }                   /* zero */
2104         }                       /* eng */
2105     }
2106
2107   /* lay out the digits of the coefficient, adding 0s and . as needed */
2108   u = *up;
2109   if (pre > 0)
2110     {                           /* xxx.xxx or xx00 (engineering) form */
2111       for (; pre > 0; pre--, c++, cut--)
2112         {
2113           if (cut < 0)
2114             {                   /* need new Unit */
2115               if (up == dn->lsu)
2116                 break;          /* out of input digits (pre>digits) */
2117               up--;
2118               cut = DECDPUN - 1;
2119               u = *up;
2120             }
2121           TODIGIT (u, cut, c);
2122         }
2123       if (up > dn->lsu || (up == dn->lsu && cut >= 0))
2124         {                       /* more to come, after '.' */
2125           *c = '.';
2126           c++;
2127           for (;; c++, cut--)
2128             {
2129               if (cut < 0)
2130                 {               /* need new Unit */
2131                   if (up == dn->lsu)
2132                     break;      /* out of input digits */
2133                   up--;
2134                   cut = DECDPUN - 1;
2135                   u = *up;
2136                 }
2137               TODIGIT (u, cut, c);
2138             }
2139         }
2140       else
2141         for (; pre > 0; pre--, c++)
2142           *c = '0';             /* 0 padding (for engineering) needed */
2143     }
2144   else
2145     {                           /* 0.xxx or 0.000xxx form */
2146       *c = '0';
2147       c++;
2148       *c = '.';
2149       c++;
2150       for (; pre < 0; pre++, c++)
2151         *c = '0';               /* add any 0's after '.' */
2152       for (;; c++, cut--)
2153         {
2154           if (cut < 0)
2155             {                   /* need new Unit */
2156               if (up == dn->lsu)
2157                 break;          /* out of input digits */
2158               up--;
2159               cut = DECDPUN - 1;
2160               u = *up;
2161             }
2162           TODIGIT (u, cut, c);
2163         }
2164     }
2165
2166   /* Finally add the E-part, if needed.  It will never be 0, has a
2167      base maximum and minimum of +999999999 through -999999999, but
2168      could range down to -1999999998 for subnormal numbers */
2169   if (e != 0)
2170     {
2171       Flag had = 0;             /* 1=had non-zero */
2172       *c = 'E';
2173       c++;
2174       *c = '+';
2175       c++;                      /* assume positive */
2176       u = e;                    /* .. */
2177       if (e < 0)
2178         {
2179           *(c - 1) = '-';       /* oops, need - */
2180           u = -e;               /* uInt, please */
2181         }
2182       /* layout the exponent (_itoa is not ANSI C) */
2183       for (cut = 9; cut >= 0; cut--)
2184         {
2185           TODIGIT (u, cut, c);
2186           if (*c == '0' && !had)
2187             continue;           /* skip leading zeros */
2188           had = 1;              /* had non-0 */
2189           c++;                  /* step for next */
2190         }                       /* cut */
2191     }
2192   *c = '\0';                    /* terminate the string (all paths) */
2193   return;
2194 }
2195
2196 /* ------------------------------------------------------------------ */
2197 /* decAddOp -- add/subtract operation                                 */
2198 /*                                                                    */
2199 /*   This computes C = A + B                                          */
2200 /*                                                                    */
2201 /*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
2202 /*   lhs is A                                                         */
2203 /*   rhs is B                                                         */
2204 /*   set is the context                                               */
2205 /*   negate is DECNEG if rhs should be negated, or 0 otherwise        */
2206 /*   status accumulates status for the caller                         */
2207 /*                                                                    */
2208 /* C must have space for set->digits digits.                          */
2209 /* ------------------------------------------------------------------ */
2210 /* If possible, we calculate the coefficient directly into C.         */
2211 /* However, if:                                                       */
2212 /*   -- we need a digits+1 calculation because numbers are unaligned  */
2213 /*      and span more than set->digits digits                         */
2214 /*   -- a carry to digits+1 digits looks possible                     */
2215 /*   -- C is the same as A or B, and the result would destructively   */
2216 /*      overlap the A or B coefficient                                */
2217 /* then we must calculate into a temporary buffer.  In this latter    */
2218 /* case we use the local (stack) buffer if possible, and only if too  */
2219 /* long for that do we resort to malloc.                              */
2220 /*                                                                    */
2221 /* Misalignment is handled as follows:                                */
2222 /*   Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp.    */
2223 /*   BPad: Apply the padding by a combination of shifting (whole      */
2224 /*         units) and multiplication (part units).                    */
2225 /*                                                                    */
2226 /* Addition, especially x=x+1, is speed-critical, so we take pains    */
2227 /* to make returning as fast as possible, by flagging any allocation. */
2228 /* ------------------------------------------------------------------ */
2229 static decNumber *
2230 decAddOp (decNumber * res, const decNumber * lhs,
2231           const decNumber * rhs, decContext * set, uByte negate, uInt * status)
2232 {
2233   decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
2234   decNumber *allocrhs = NULL;   /* .., rhs */
2235   Int rhsshift;                 /* working shift (in Units) */
2236   Int maxdigits;                /* longest logical length */
2237   Int mult;                     /* multiplier */
2238   Int residue;                  /* rounding accumulator */
2239   uByte bits;                   /* result bits */
2240   Flag diffsign;                /* non-0 if arguments have different sign */
2241   Unit *acc;                    /* accumulator for result */
2242   Unit accbuff[D2U (DECBUFFER + 1)];    /* local buffer [+1 is for possible */
2243   /* final carry digit or DECBUFFER=0] */
2244   Unit *allocacc = NULL;        /* -> allocated acc buffer, iff allocated */
2245   Flag alloced = 0;             /* set non-0 if any allocations */
2246   Int reqdigits = set->digits;  /* local copy; requested DIGITS */
2247   uByte merged;                 /* merged flags */
2248   Int padding;                  /* work */
2249
2250 #if DECCHECK
2251   if (decCheckOperands (res, lhs, rhs, set))
2252     return res;
2253 #endif
2254
2255   do
2256     {                           /* protect allocated storage */
2257 #if DECSUBSET
2258       if (!set->extended)
2259         {
2260           /* reduce operands and set lostDigits status, as needed */
2261           if (lhs->digits > reqdigits)
2262             {
2263               alloclhs = decRoundOperand (lhs, set, status);
2264               if (alloclhs == NULL)
2265                 break;
2266               lhs = alloclhs;
2267               alloced = 1;
2268             }
2269           if (rhs->digits > reqdigits)
2270             {
2271               allocrhs = decRoundOperand (rhs, set, status);
2272               if (allocrhs == NULL)
2273                 break;
2274               rhs = allocrhs;
2275               alloced = 1;
2276             }
2277         }
2278 #endif
2279       /* [following code does not require input rounding] */
2280
2281       /* note whether signs differ */
2282       diffsign = (Flag) ((lhs->bits ^ rhs->bits ^ negate) & DECNEG);
2283
2284       /* handle infinities and NaNs */
2285       merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2286       if (merged)
2287         {                       /* a special bit set */
2288           if (merged & (DECSNAN | DECNAN))      /* a NaN */
2289             decNaNs (res, lhs, rhs, status);
2290           else
2291             {                   /* one or two infinities */
2292               if (decNumberIsInfinite (lhs))
2293                 {               /* LHS is infinity */
2294                   /* two infinities with different signs is invalid */
2295                   if (decNumberIsInfinite (rhs) && diffsign)
2296                     {
2297                       *status |= DEC_Invalid_operation;
2298                       break;
2299                     }
2300                   bits = lhs->bits & DECNEG;    /* get sign from LHS */
2301                 }
2302               else
2303                 bits = (rhs->bits ^ negate) & DECNEG;   /* RHS must be Infinity */
2304               bits |= DECINF;
2305               decNumberZero (res);
2306               res->bits = bits; /* set +/- infinity */
2307             }                   /* an infinity */
2308           break;
2309         }
2310
2311       /* Quick exit for add 0s; return the non-0, modified as need be */
2312       if (ISZERO (lhs))
2313         {
2314           Int adjust;           /* work */
2315           Int lexp = lhs->exponent;     /* save in case LHS==RES */
2316           bits = lhs->bits;     /* .. */
2317           residue = 0;          /* clear accumulator */
2318           decCopyFit (res, rhs, set, &residue, status); /* copy (as needed) */
2319           res->bits ^= negate;  /* flip if rhs was negated */
2320 #if DECSUBSET
2321           if (set->extended)
2322             {                   /* exponents on zeros count */
2323 #endif
2324               /* exponent will be the lower of the two */
2325               adjust = lexp - res->exponent;    /* adjustment needed [if -ve] */
2326               if (ISZERO (res))
2327                 {               /* both 0: special IEEE 854 rules */
2328                   if (adjust < 0)
2329                     res->exponent = lexp;       /* set exponent */
2330                   /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
2331                   if (diffsign)
2332                     {
2333                       if (set->round != DEC_ROUND_FLOOR)
2334                         res->bits = 0;
2335                       else
2336                         res->bits = DECNEG;     /* preserve 0 sign */
2337                     }
2338                 }
2339               else
2340                 {               /* non-0 res */
2341                   if (adjust < 0)
2342                     {           /* 0-padding needed */
2343                       if ((res->digits - adjust) > set->digits)
2344                         {
2345                           adjust = res->digits - set->digits;   /* to fit exactly */
2346                           *status |= DEC_Rounded;       /* [but exact] */
2347                         }
2348                       res->digits =
2349                         decShiftToMost (res->lsu, res->digits, -adjust);
2350                       res->exponent += adjust;  /* set the exponent. */
2351                     }
2352                 }               /* non-0 res */
2353 #if DECSUBSET
2354             }                   /* extended */
2355 #endif
2356           decFinish (res, set, &residue, status);       /* clean and finalize */
2357           break;
2358         }
2359
2360       if (ISZERO (rhs))
2361         {                       /* [lhs is non-zero] */
2362           Int adjust;           /* work */
2363           Int rexp = rhs->exponent;     /* save in case RHS==RES */
2364           bits = rhs->bits;     /* be clean */
2365           residue = 0;          /* clear accumulator */
2366           decCopyFit (res, lhs, set, &residue, status); /* copy (as needed) */
2367 #if DECSUBSET
2368           if (set->extended)
2369             {                   /* exponents on zeros count */
2370 #endif
2371               /* exponent will be the lower of the two */
2372               /* [0-0 case handled above] */
2373               adjust = rexp - res->exponent;    /* adjustment needed [if -ve] */
2374               if (adjust < 0)
2375                 {               /* 0-padding needed */
2376                   if ((res->digits - adjust) > set->digits)
2377                     {
2378                       adjust = res->digits - set->digits;       /* to fit exactly */
2379                       *status |= DEC_Rounded;   /* [but exact] */
2380                     }
2381                   res->digits =
2382                     decShiftToMost (res->lsu, res->digits, -adjust);
2383                   res->exponent += adjust;      /* set the exponent. */
2384                 }
2385 #if DECSUBSET
2386             }                   /* extended */
2387 #endif
2388           decFinish (res, set, &residue, status);       /* clean and finalize */
2389           break;
2390         }
2391       /* [both fastpath and mainpath code below assume these cases */
2392       /* (notably 0-0) have already been handled] */
2393
2394       /* calculate the padding needed to align the operands */
2395       padding = rhs->exponent - lhs->exponent;
2396
2397       /* Fastpath cases where the numbers are aligned and normal, the RHS */
2398       /* is all in one unit, no operand rounding is needed, and no carry, */
2399       /* lengthening, or borrow is needed */
2400       if (rhs->digits <= DECDPUN && padding == 0 && rhs->exponent >= set->emin  /* [some normals drop through] */
2401           && rhs->digits <= reqdigits && lhs->digits <= reqdigits)
2402         {
2403           Int partial = *lhs->lsu;
2404           if (!diffsign)
2405             {                   /* adding */
2406               Int maxv = DECDPUNMAX;    /* highest no-overflow */
2407               if (lhs->digits < DECDPUN)
2408                 maxv = powers[lhs->digits] - 1;
2409               partial += *rhs->lsu;
2410               if (partial <= maxv)
2411                 {               /* no carry */
2412                   if (res != lhs)
2413                     decNumberCopy (res, lhs);   /* not in place */
2414                   *res->lsu = (Unit) partial;   /* [copy could have overwritten RHS] */
2415                   break;
2416                 }
2417               /* else drop out for careful add */
2418             }
2419           else
2420             {                   /* signs differ */
2421               partial -= *rhs->lsu;
2422               if (partial > 0)
2423                 {               /* no borrow needed, and non-0 result */
2424                   if (res != lhs)
2425                     decNumberCopy (res, lhs);   /* not in place */
2426                   *res->lsu = (Unit) partial;
2427                   /* this could have reduced digits [but result>0] */
2428                   res->digits = decGetDigits (res->lsu, D2U (res->digits));
2429                   break;
2430                 }
2431               /* else drop out for careful subtract */
2432             }
2433         }
2434
2435       /* Now align (pad) the lhs or rhs so we can add or subtract them, as
2436          necessary.  If one number is much larger than the other (that is,
2437          if in plain form there is a least one digit between the lowest
2438          digit or one and the highest of the other) we need to pad with up
2439          to DIGITS-1 trailing zeros, and then apply rounding (as exotic
2440          rounding modes may be affected by the residue).
2441        */
2442       rhsshift = 0;             /* rhs shift to left (padding) in Units */
2443       bits = lhs->bits;         /* assume sign is that of LHS */
2444       mult = 1;                 /* likely multiplier */
2445
2446       /* if padding==0 the operands are aligned; no padding needed */
2447       if (padding != 0)
2448         {
2449           /* some padding needed */
2450           /* We always pad the RHS, as we can then effect any required */
2451           /* padding by a combination of shifts and a multiply */
2452           Flag swapped = 0;
2453           if (padding < 0)
2454             {                   /* LHS needs the padding */
2455               const decNumber *t;
2456               padding = -padding;       /* will be +ve */
2457               bits = (uByte) (rhs->bits ^ negate);      /* assumed sign is now that of RHS */
2458               t = lhs;
2459               lhs = rhs;
2460               rhs = t;
2461               swapped = 1;
2462             }
2463
2464           /* If, after pad, rhs would be longer than lhs by digits+1 or */
2465           /* more then lhs cannot affect the answer, except as a residue, */
2466           /* so we only need to pad up to a length of DIGITS+1. */
2467           if (rhs->digits + padding > lhs->digits + reqdigits + 1)
2468             {
2469               /* The RHS is sufficient */
2470               /* for residue we use the relative sign indication... */
2471               Int shift = reqdigits - rhs->digits;      /* left shift needed */
2472               residue = 1;      /* residue for rounding */
2473               if (diffsign)
2474                 residue = -residue;     /* signs differ */
2475               /* copy, shortening if necessary */
2476               decCopyFit (res, rhs, set, &residue, status);
2477               /* if it was already shorter, then need to pad with zeros */
2478               if (shift > 0)
2479                 {
2480                   res->digits = decShiftToMost (res->lsu, res->digits, shift);
2481                   res->exponent -= shift;       /* adjust the exponent. */
2482                 }
2483               /* flip the result sign if unswapped and rhs was negated */
2484               if (!swapped)
2485                 res->bits ^= negate;
2486               decFinish (res, set, &residue, status);   /* done */
2487               break;
2488             }
2489
2490           /* LHS digits may affect result */
2491           rhsshift = D2U (padding + 1) - 1;     /* this much by Unit shift .. */
2492           mult = powers[padding - (rhsshift * DECDPUN)];        /* .. this by multiplication */
2493         }                       /* padding needed */
2494
2495       if (diffsign)
2496         mult = -mult;           /* signs differ */
2497
2498       /* determine the longer operand */
2499       maxdigits = rhs->digits + padding;        /* virtual length of RHS */
2500       if (lhs->digits > maxdigits)
2501         maxdigits = lhs->digits;
2502
2503       /* Decide on the result buffer to use; if possible place directly */
2504       /* into result. */
2505       acc = res->lsu;           /* assume build direct */
2506       /* If destructive overlap, or the number is too long, or a carry or */
2507       /* borrow to DIGITS+1 might be possible we must use a buffer. */
2508       /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
2509       if ((maxdigits >= reqdigits)      /* is, or could be, too large */
2510           || (res == rhs && rhsshift > 0))
2511         {                       /* destructive overlap */
2512           /* buffer needed; choose it */
2513           /* we'll need units for maxdigits digits, +1 Unit for carry or borrow */
2514           Int need = D2U (maxdigits) + 1;
2515           acc = accbuff;        /* assume use local buffer */
2516           if (need * sizeof (Unit) > sizeof (accbuff))
2517             {
2518               allocacc = (Unit *) malloc (need * sizeof (Unit));
2519               if (allocacc == NULL)
2520                 {               /* hopeless -- abandon */
2521                   *status |= DEC_Insufficient_storage;
2522                   break;
2523                 }
2524               acc = allocacc;
2525               alloced = 1;
2526             }
2527         }
2528
2529       res->bits = (uByte) (bits & DECNEG);      /* it's now safe to overwrite.. */
2530       res->exponent = lhs->exponent;    /* .. operands (even if aliased) */
2531
2532 #if DECTRACE
2533       decDumpAr ('A', lhs->lsu, D2U (lhs->digits));
2534       decDumpAr ('B', rhs->lsu, D2U (rhs->digits));
2535       printf ("  :h: %d %d\n", rhsshift, mult);
2536 #endif
2537
2538       /* add [A+B*m] or subtract [A+B*(-m)] */
2539       res->digits = decUnitAddSub (lhs->lsu, D2U (lhs->digits), rhs->lsu, D2U (rhs->digits), rhsshift, acc, mult) * DECDPUN;    /* [units -> digits] */
2540       if (res->digits < 0)
2541         {                       /* we borrowed */
2542           res->digits = -res->digits;
2543           res->bits ^= DECNEG;  /* flip the sign */
2544         }
2545 #if DECTRACE
2546       decDumpAr ('+', acc, D2U (res->digits));
2547 #endif
2548
2549       /* If we used a buffer we need to copy back, possibly shortening */
2550       /* (If we didn't use buffer it must have fit, so can't need rounding */
2551       /* and residue must be 0.) */
2552       residue = 0;              /* clear accumulator */
2553       if (acc != res->lsu)
2554         {
2555 #if DECSUBSET
2556           if (set->extended)
2557             {                   /* round from first significant digit */
2558 #endif
2559               /* remove leading zeros that we added due to rounding up to */
2560               /* integral Units -- before the test for rounding. */
2561               if (res->digits > reqdigits)
2562                 res->digits = decGetDigits (acc, D2U (res->digits));
2563               decSetCoeff (res, set, acc, res->digits, &residue, status);
2564 #if DECSUBSET
2565             }
2566           else
2567             {                   /* subset arithmetic rounds from original significant digit */
2568               /* We may have an underestimate.  This only occurs when both */
2569               /* numbers fit in DECDPUN digits and we are padding with a */
2570               /* negative multiple (-10, -100...) and the top digit(s) become */
2571               /* 0.  (This only matters if we are using X3.274 rules where the */
2572               /* leading zero could be included in the rounding.) */
2573               if (res->digits < maxdigits)
2574                 {
2575                   *(acc + D2U (res->digits)) = 0;       /* ensure leading 0 is there */
2576                   res->digits = maxdigits;
2577                 }
2578               else
2579                 {
2580                   /* remove leading zeros that we added due to rounding up to */
2581                   /* integral Units (but only those in excess of the original */
2582                   /* maxdigits length, unless extended) before test for rounding. */
2583                   if (res->digits > reqdigits)
2584                     {
2585                       res->digits = decGetDigits (acc, D2U (res->digits));
2586                       if (res->digits < maxdigits)
2587                         res->digits = maxdigits;
2588                     }
2589                 }
2590               decSetCoeff (res, set, acc, res->digits, &residue, status);
2591               /* Now apply rounding if needed before removing leading zeros. */
2592               /* This is safe because subnormals are not a possibility */
2593               if (residue != 0)
2594                 {
2595                   decApplyRound (res, set, residue, status);
2596                   residue = 0;  /* we did what we had to do */
2597                 }
2598             }                   /* subset */
2599 #endif
2600         }                       /* used buffer */
2601
2602       /* strip leading zeros [these were left on in case of subset subtract] */
2603       res->digits = decGetDigits (res->lsu, D2U (res->digits));
2604
2605       /* apply checks and rounding */
2606       decFinish (res, set, &residue, status);
2607
2608       /* "When the sum of two operands with opposite signs is exactly */
2609       /* zero, the sign of that sum shall be '+' in all rounding modes */
2610       /* except round toward -Infinity, in which mode that sign shall be */
2611       /* '-'."  [Subset zeros also never have '-', set by decFinish.] */
2612       if (ISZERO (res) && diffsign
2613 #if DECSUBSET
2614           && set->extended
2615 #endif
2616           && (*status & DEC_Inexact) == 0)
2617         {
2618           if (set->round == DEC_ROUND_FLOOR)
2619             res->bits |= DECNEG;        /* sign - */
2620           else
2621             res->bits &= ~DECNEG;       /* sign + */
2622         }
2623     }
2624   while (0);                    /* end protected */
2625
2626   if (alloced)
2627     {
2628       if (allocacc != NULL)
2629         free (allocacc);        /* drop any storage we used */
2630       if (allocrhs != NULL)
2631         free (allocrhs);        /* .. */
2632       if (alloclhs != NULL)
2633         free (alloclhs);        /* .. */
2634     }
2635   return res;
2636 }
2637
2638 /* ------------------------------------------------------------------ */
2639 /* decDivideOp -- division operation                                  */
2640 /*                                                                    */
2641 /*  This routine performs the calculations for all four division      */
2642 /*  operators (divide, divideInteger, remainder, remainderNear).      */
2643 /*                                                                    */
2644 /*  C=A op B                                                          */
2645 /*                                                                    */
2646 /*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
2647 /*   lhs is A                                                         */
2648 /*   rhs is B                                                         */
2649 /*   set is the context                                               */
2650 /*   op  is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively.    */
2651 /*   status is the usual accumulator                                  */
2652 /*                                                                    */
2653 /* C must have space for set->digits digits.                          */
2654 /*                                                                    */
2655 /* ------------------------------------------------------------------ */
2656 /*   The underlying algorithm of this routine is the same as in the   */
2657 /*   1981 S/370 implementation, that is, non-restoring long division  */
2658 /*   with bi-unit (rather than bi-digit) estimation for each unit     */
2659 /*   multiplier.  In this pseudocode overview, complications for the  */
2660 /*   Remainder operators and division residues for exact rounding are */
2661 /*   omitted for clarity.                                             */
2662 /*                                                                    */
2663 /*     Prepare operands and handle special values                     */
2664 /*     Test for x/0 and then 0/x                                      */
2665 /*     Exp =Exp1 - Exp2                                               */
2666 /*     Exp =Exp +len(var1) -len(var2)                                 */
2667 /*     Sign=Sign1 * Sign2                                             */
2668 /*     Pad accumulator (Var1) to double-length with 0's (pad1)        */
2669 /*     Pad Var2 to same length as Var1                                */
2670 /*     msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round  */
2671 /*     have=0                                                         */
2672 /*     Do until (have=digits+1 OR residue=0)                          */
2673 /*       if exp<0 then if integer divide/residue then leave           */
2674 /*       this_unit=0                                                  */
2675 /*       Do forever                                                   */
2676 /*          compare numbers                                           */
2677 /*          if <0 then leave inner_loop                               */
2678 /*          if =0 then (* quick exit without subtract *) do           */
2679 /*             this_unit=this_unit+1; output this_unit                */
2680 /*             leave outer_loop; end                                  */
2681 /*          Compare lengths of numbers (mantissae):                   */
2682 /*          If same then tops2=msu2pair -- {units 1&2 of var2}        */
2683 /*                  else tops2=msu2plus -- {0, unit 1 of var2}        */
2684 /*          tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
2685 /*          mult=tops1/tops2  -- Good and safe guess at divisor       */
2686 /*          if mult=0 then mult=1                                     */
2687 /*          this_unit=this_unit+mult                                  */
2688 /*          subtract                                                  */
2689 /*          end inner_loop                                            */
2690 /*        if have\=0 | this_unit\=0 then do                           */
2691 /*          output this_unit                                          */
2692 /*          have=have+1; end                                          */
2693 /*        var2=var2/10                                                */
2694 /*        exp=exp-1                                                   */
2695 /*        end outer_loop                                              */
2696 /*     exp=exp+1   -- set the proper exponent                         */
2697 /*     if have=0 then generate answer=0                               */
2698 /*     Return (Result is defined by Var1)                             */
2699 /*                                                                    */
2700 /* ------------------------------------------------------------------ */
2701 /* We need two working buffers during the long division; one (digits+ */
2702 /* 1) to accumulate the result, and the other (up to 2*digits+1) for  */
2703 /* long subtractions.  These are acc and var1 respectively.           */
2704 /* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
2705 /* ------------------------------------------------------------------ */
2706 static decNumber *
2707 decDivideOp (decNumber * res,
2708              const decNumber * lhs, const decNumber * rhs,
2709              decContext * set, Flag op, uInt * status)
2710 {
2711   decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
2712   decNumber *allocrhs = NULL;   /* .., rhs */
2713   Unit accbuff[D2U (DECBUFFER + DECDPUN)];      /* local buffer */
2714   Unit *acc = accbuff;          /* -> accumulator array for result */
2715   Unit *allocacc = NULL;        /* -> allocated buffer, iff allocated */
2716   Unit *accnext;                /* -> where next digit will go */
2717   Int acclength;                /* length of acc needed [Units] */
2718   Int accunits;                 /* count of units accumulated */
2719   Int accdigits;                /* count of digits accumulated */
2720
2721   Unit varbuff[D2U (DECBUFFER * 2 + DECDPUN) * sizeof (Unit)];  /* buffer for var1 */
2722   Unit *var1 = varbuff;         /* -> var1 array for long subtraction */
2723   Unit *varalloc = NULL;        /* -> allocated buffer, iff used */
2724
2725   const Unit *var2;             /* -> var2 array */
2726
2727   Int var1units, var2units;     /* actual lengths */
2728   Int var2ulen;                 /* logical length (units) */
2729   Int var1initpad = 0;          /* var1 initial padding (digits) */
2730   Unit *msu1;                   /* -> msu of each var */
2731   const Unit *msu2;             /* -> msu of each var */
2732   Int msu2plus;                 /* msu2 plus one [does not vary] */
2733   eInt msu2pair;                /* msu2 pair plus one [does not vary] */
2734   Int maxdigits;                /* longest LHS or required acc length */
2735   Int mult;                     /* multiplier for subtraction */
2736   Unit thisunit;                /* current unit being accumulated */
2737   Int residue;                  /* for rounding */
2738   Int reqdigits = set->digits;  /* requested DIGITS */
2739   Int exponent;                 /* working exponent */
2740   Int maxexponent = 0;          /* DIVIDE maximum exponent if unrounded */
2741   uByte bits;                   /* working sign */
2742   uByte merged;                 /* merged flags */
2743   Unit *target;                 /* work */
2744   const Unit *source;           /* work */
2745   uInt const *pow;              /* .. */
2746   Int shift, cut;               /* .. */
2747 #if DECSUBSET
2748   Int dropped;                  /* work */
2749 #endif
2750
2751 #if DECCHECK
2752   if (decCheckOperands (res, lhs, rhs, set))
2753     return res;
2754 #endif
2755
2756   do
2757     {                           /* protect allocated storage */
2758 #if DECSUBSET
2759       if (!set->extended)
2760         {
2761           /* reduce operands and set lostDigits status, as needed */
2762           if (lhs->digits > reqdigits)
2763             {
2764               alloclhs = decRoundOperand (lhs, set, status);
2765               if (alloclhs == NULL)
2766                 break;
2767               lhs = alloclhs;
2768             }
2769           if (rhs->digits > reqdigits)
2770             {
2771               allocrhs = decRoundOperand (rhs, set, status);
2772               if (allocrhs == NULL)
2773                 break;
2774               rhs = allocrhs;
2775             }
2776         }
2777 #endif
2778       /* [following code does not require input rounding] */
2779
2780       bits = (lhs->bits ^ rhs->bits) & DECNEG;  /* assumed sign for divisions */
2781
2782       /* handle infinities and NaNs */
2783       merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2784       if (merged)
2785         {                       /* a special bit set */
2786           if (merged & (DECSNAN | DECNAN))
2787             {                   /* one or two NaNs */
2788               decNaNs (res, lhs, rhs, status);
2789               break;
2790             }
2791           /* one or two infinities */
2792           if (decNumberIsInfinite (lhs))
2793             {                   /* LHS (dividend) is infinite */
2794               if (decNumberIsInfinite (rhs) ||  /* two infinities are invalid .. */
2795                   op & (REMAINDER | REMNEAR))
2796                 {               /* as is remainder of infinity */
2797                   *status |= DEC_Invalid_operation;
2798                   break;
2799                 }
2800               /* [Note that infinity/0 raises no exceptions] */
2801               decNumberZero (res);
2802               res->bits = bits | DECINF;        /* set +/- infinity */
2803               break;
2804             }
2805           else
2806             {                   /* RHS (divisor) is infinite */
2807               residue = 0;
2808               if (op & (REMAINDER | REMNEAR))
2809                 {
2810                   /* result is [finished clone of] lhs */
2811                   decCopyFit (res, lhs, set, &residue, status);
2812                 }
2813               else
2814                 {               /* a division */
2815                   decNumberZero (res);
2816                   res->bits = bits;     /* set +/- zero */
2817                   /* for DIVIDEINT the exponent is always 0.  For DIVIDE, result */
2818                   /* is a 0 with infinitely negative exponent, clamped to minimum */
2819                   if (op & DIVIDE)
2820                     {