libdecnumber/decNumber.c

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