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