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