OSDN Git Service

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