libdecnumber/decNumber.c

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