OSDN Git Service

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