libdecnumber/decNumber.c

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