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