qemu/libdecnumber/decNumber.c
<<
>>
Prefs
   1/* Decimal number arithmetic module for the decNumber C Library.
   2   Copyright (C) 2005, 2007 Free Software Foundation, Inc.
   3   Contributed by IBM Corporation.  Author Mike Cowlishaw.
   4
   5   This file is part of GCC.
   6
   7   GCC is free software; you can redistribute it and/or modify it under
   8   the terms of the GNU General Public License as published by the Free
   9   Software Foundation; either version 2, or (at your option) any later
  10   version.
  11
  12   In addition to the permissions in the GNU General Public License,
  13   the Free Software Foundation gives you unlimited permission to link
  14   the compiled version of this file into combinations with other
  15   programs, and to distribute those combinations without any
  16   restriction coming from the use of this file.  (The General Public
  17   License restrictions do apply in other respects; for example, they
  18   cover modification of the file, and distribution when not linked
  19   into a combine executable.)
  20
  21   GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  22   WARRANTY; without even the implied warranty of MERCHANTABILITY or
  23   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  24   for more details.
  25
  26   You should have received a copy of the GNU General Public License
  27   along with GCC; see the file COPYING.  If not, write to the Free
  28   Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  29   02110-1301, USA.  */
  30
  31/* ------------------------------------------------------------------ */
  32/* Decimal Number arithmetic module                                   */
  33/* ------------------------------------------------------------------ */
  34/* This module comprises the routines for General Decimal Arithmetic  */
  35/* as defined in the specification which may be found on the          */
  36/* http://www2.hursley.ibm.com/decimal web pages.  It implements both */
  37/* the full ('extended') arithmetic and the simpler ('subset')        */
  38/* arithmetic.                                                        */
  39/*                                                                    */
  40/* Usage notes:                                                       */
  41/*                                                                    */
  42/* 1. This code is ANSI C89 except:                                   */
  43/*                                                                    */
  44/*       If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and       */
  45/*       uint64_t types may be used.  To avoid these, set DECUSE64=0  */
  46/*       and DECDPUN<=4 (see documentation).                          */
  47/*                                                                    */
  48/* 2. The decNumber format which this library uses is optimized for   */
  49/*    efficient processing of relatively short numbers; in particular */
  50/*    it allows the use of fixed sized structures and minimizes copy  */
  51/*    and move operations.  It does, however, support arbitrary       */
  52/*    precision (up to 999,999,999 digits) and arbitrary exponent     */
  53/*    range (Emax in the range 0 through 999,999,999 and Emin in the  */
  54/*    range -999,999,999 through 0).  Mathematical functions (for     */
  55/*    example decNumberExp) as identified below are restricted more   */
  56/*    tightly: digits, emax, and -emin in the context must be <=      */
  57/*    DEC_MAX_MATH (999999), and their operand(s) must be within      */
  58/*    these bounds.                                                   */
  59/*                                                                    */
  60/* 3. Logical functions are further restricted; their operands must   */
  61/*    be finite, positive, have an exponent of zero, and all digits   */
  62/*    must be either 0 or 1.  The result will only contain digits     */
  63/*    which are 0 or 1 (and will have exponent=0 and a sign of 0).    */
  64/*                                                                    */
  65/* 4. Operands to operator functions are never modified unless they   */
  66/*    are also specified to be the result number (which is always     */
  67/*    permitted).  Other than that case, operands must not overlap.   */
  68/*                                                                    */
  69/* 5. Error handling: the type of the error is ORed into the status   */
  70/*    flags in the current context (decContext structure).  The       */
  71/*    SIGFPE signal is then raised if the corresponding trap-enabler  */
  72/*    flag in the decContext is set (is 1).                           */
  73/*                                                                    */
  74/*    It is the responsibility of the caller to clear the status      */
  75/*    flags as required.                                              */
  76/*                                                                    */
  77/*    The result of any routine which returns a number will always    */
  78/*    be a valid number (which may be a special value, such as an     */
  79/*    Infinity or NaN).                                               */
  80/*                                                                    */
  81/* 6. The decNumber format is not an exchangeable concrete            */
  82/*    representation as it comprises fields which may be machine-     */
  83/*    dependent (packed or unpacked, or special length, for example). */
  84/*    Canonical conversions to and from strings are provided; other   */
  85/*    conversions are available in separate modules.                  */
  86/*                                                                    */
  87/* 7. Normally, input operands are assumed to be valid.  Set DECCHECK */
  88/*    to 1 for extended operand checking (including NULL operands).   */
  89/*    Results are undefined if a badly-formed structure (or a NULL    */
  90/*    pointer to a structure) is provided, though with DECCHECK       */
  91/*    enabled the operator routines are protected against exceptions. */
  92/*    (Except if the result pointer is NULL, which is unrecoverable.) */
  93/*                                                                    */
  94/*    However, the routines will never cause exceptions if they are   */
  95/*    given well-formed operands, even if the value of the operands   */
  96/*    is inappropriate for the operation and DECCHECK is not set.     */
  97/*    (Except for SIGFPE, as and where documented.)                   */
  98/*                                                                    */
  99/* 8. Subset arithmetic is available only if DECSUBSET is set to 1.   */
 100/* ------------------------------------------------------------------ */
 101/* Implementation notes for maintenance of this module:               */
 102/*                                                                    */
 103/* 1. Storage leak protection:  Routines which use malloc are not     */
 104/*    permitted to use return for fastpath or error exits (i.e.,      */
 105/*    they follow strict structured programming conventions).         */
 106/*    Instead they have a do{}while(0); construct surrounding the     */
 107/*    code which is protected -- break may be used to exit this.      */
 108/*    Other routines can safely use the return statement inline.      */
 109/*                                                                    */
 110/*    Storage leak accounting can be enabled using DECALLOC.          */
 111/*                                                                    */
 112/* 2. All loops use the for(;;) construct.  Any do construct does     */
 113/*    not loop; it is for allocation protection as just described.    */
 114/*                                                                    */
 115/* 3. Setting status in the context must always be the very last      */
 116/*    action in a routine, as non-0 status may raise a trap and hence */
 117/*    the call to set status may not return (if the handler uses long */
 118/*    jump).  Therefore all cleanup must be done first.  In general,  */
 119/*    to achieve this status is accumulated and is only applied just  */
 120/*    before return by calling decContextSetStatus (via decStatus).   */
 121/*                                                                    */
 122/*    Routines which allocate storage cannot, in general, use the     */
 123/*    'top level' routines which could cause a non-returning          */
 124/*    transfer of control.  The decXxxxOp routines are safe (do not   */
 125/*    call decStatus even if traps are set in the context) and should */
 126/*    be used instead (they are also a little faster).                */
 127/*                                                                    */
 128/* 4. Exponent checking is minimized by allowing the exponent to      */
 129/*    grow outside its limits during calculations, provided that      */
 130/*    the decFinalize function is called later.  Multiplication and   */
 131/*    division, and intermediate calculations in exponentiation,      */
 132/*    require more careful checks because of the risk of 31-bit       */
 133/*    overflow (the most negative valid exponent is -1999999997, for  */
 134/*    a 999999999-digit number with adjusted exponent of -999999999). */
 135/*                                                                    */
 136/* 5. Rounding is deferred until finalization of results, with any    */
 137/*    'off to the right' data being represented as a single digit     */
 138/*    residue (in the range -1 through 9).  This avoids any double-   */
 139/*    rounding when more than one shortening takes place (for         */
 140/*    example, when a result is subnormal).                           */
 141/*                                                                    */
 142/* 6. The digits count is allowed to rise to a multiple of DECDPUN    */
 143/*    during many operations, so whole Units are handled and exact    */
 144/*    accounting of digits is not needed.  The correct digits value   */
 145/*    is found by decGetDigits, which accounts for leading zeros.     */
 146/*    This must be called before any rounding if the number of digits */
 147/*    is not known exactly.                                           */
 148/*                                                                    */
 149/* 7. The multiply-by-reciprocal 'trick' is used for partitioning     */
 150/*    numbers up to four digits, using appropriate constants.  This   */
 151/*    is not useful for longer numbers because overflow of 32 bits    */
 152/*    would lead to 4 multiplies, which is almost as expensive as     */
 153/*    a divide (unless a floating-point or 64-bit multiply is         */
 154/*    assumed to be available).                                       */
 155/*                                                                    */
 156/* 8. Unusual abbreviations that may be used in the commentary:       */
 157/*      lhs -- left hand side (operand, of an operation)              */
 158/*      lsd -- least significant digit (of coefficient)               */
 159/*      lsu -- least significant Unit (of coefficient)                */
 160/*      msd -- most significant digit (of coefficient)                */
 161/*      msi -- most significant item (in an array)                    */
 162/*      msu -- most significant Unit (of coefficient)                 */
 163/*      rhs -- right hand side (operand, of an operation)             */
 164/*      +ve -- positive                                               */
 165/*      -ve -- negative                                               */
 166/*      **  -- raise to the power                                     */
 167/* ------------------------------------------------------------------ */
 168
 169#include "qemu/osdep.h"
 170#include "libdecnumber/dconfig.h"
 171#include "libdecnumber/decNumber.h"
 172#include "libdecnumber/decNumberLocal.h"
 173
 174/* Constants */
 175/* Public lookup table used by the D2U macro */
 176const uByte d2utable[DECMAXD2U+1]=D2UTABLE;
 177
 178#define DECVERB     1              /* set to 1 for verbose DECCHECK */
 179#define powers      DECPOWERS      /* old internal name */
 180
 181/* Local constants */
 182#define DIVIDE      0x80           /* Divide operators */
 183#define REMAINDER   0x40           /* .. */
 184#define DIVIDEINT   0x20           /* .. */
 185#define REMNEAR     0x10           /* .. */
 186#define COMPARE     0x01           /* Compare operators */
 187#define COMPMAX     0x02           /* .. */
 188#define COMPMIN     0x03           /* .. */
 189#define COMPTOTAL   0x04           /* .. */
 190#define COMPNAN     0x05           /* .. [NaN processing] */
 191#define COMPSIG     0x06           /* .. [signaling COMPARE] */
 192#define COMPMAXMAG  0x07           /* .. */
 193#define COMPMINMAG  0x08           /* .. */
 194
 195#define DEC_sNaN     0x40000000    /* local status: sNaN signal */
 196#define BADINT  (Int)0x80000000    /* most-negative Int; error indicator */
 197/* Next two indicate an integer >= 10**6, and its parity (bottom bit) */
 198#define BIGEVEN (Int)0x80000002
 199#define BIGODD  (Int)0x80000003
 200
 201static Unit uarrone[1]={1};   /* Unit array of 1, used for incrementing */
 202
 203/* Granularity-dependent code */
 204#if DECDPUN<=4
 205  #define eInt  Int           /* extended integer */
 206  #define ueInt uInt          /* unsigned extended integer */
 207  /* Constant multipliers for divide-by-power-of five using reciprocal */
 208  /* multiply, after removing powers of 2 by shifting, and final shift */
 209  /* of 17 [we only need up to **4] */
 210  static const uInt multies[]={131073, 26215, 5243, 1049, 210};
 211  /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
 212  #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
 213#else
 214  /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
 215  #if !DECUSE64
 216    #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4
 217  #endif
 218  #define eInt  Long          /* extended integer */
 219  #define ueInt uLong         /* unsigned extended integer */
 220#endif
 221
 222/* Local routines */
 223static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *,
 224                              decContext *, uByte, uInt *);
 225static Flag        decBiStr(const char *, const char *, const char *);
 226static uInt        decCheckMath(const decNumber *, decContext *, uInt *);
 227static void        decApplyRound(decNumber *, decContext *, Int, uInt *);
 228static Int         decCompare(const decNumber *lhs, const decNumber *rhs, Flag);
 229static decNumber * decCompareOp(decNumber *, const decNumber *,
 230                              const decNumber *, decContext *,
 231                              Flag, uInt *);
 232static void        decCopyFit(decNumber *, const decNumber *, decContext *,
 233                              Int *, uInt *);
 234static decNumber * decDecap(decNumber *, Int);
 235static decNumber * decDivideOp(decNumber *, const decNumber *,
 236                              const decNumber *, decContext *, Flag, uInt *);
 237static decNumber * decExpOp(decNumber *, const decNumber *,
 238                              decContext *, uInt *);
 239static void        decFinalize(decNumber *, decContext *, Int *, uInt *);
 240static Int         decGetDigits(Unit *, Int);
 241static Int         decGetInt(const decNumber *);
 242static decNumber * decLnOp(decNumber *, const decNumber *,
 243                              decContext *, uInt *);
 244static decNumber * decMultiplyOp(decNumber *, const decNumber *,
 245                              const decNumber *, decContext *,
 246                              uInt *);
 247static decNumber * decNaNs(decNumber *, const decNumber *,
 248                              const decNumber *, decContext *, uInt *);
 249static decNumber * decQuantizeOp(decNumber *, const decNumber *,
 250                              const decNumber *, decContext *, Flag,
 251                              uInt *);
 252static void        decReverse(Unit *, Unit *);
 253static void        decSetCoeff(decNumber *, decContext *, const Unit *,
 254                              Int, Int *, uInt *);
 255static void        decSetMaxValue(decNumber *, decContext *);
 256static void        decSetOverflow(decNumber *, decContext *, uInt *);
 257static void        decSetSubnormal(decNumber *, decContext *, Int *, uInt *);
 258static Int         decShiftToLeast(Unit *, Int, Int);
 259static Int         decShiftToMost(Unit *, Int, Int);
 260static void        decStatus(decNumber *, uInt, decContext *);
 261static void        decToString(const decNumber *, char[], Flag);
 262static decNumber * decTrim(decNumber *, decContext *, Flag, Int *);
 263static Int         decUnitAddSub(const Unit *, Int, const Unit *, Int, Int,
 264                              Unit *, Int);
 265static Int         decUnitCompare(const Unit *, Int, const Unit *, Int, Int);
 266
 267#if !DECSUBSET
 268/* decFinish == decFinalize when no subset arithmetic needed */
 269#define decFinish(a,b,c,d) decFinalize(a,b,c,d)
 270#else
 271static void        decFinish(decNumber *, decContext *, Int *, uInt *);
 272static decNumber * decRoundOperand(const decNumber *, decContext *, uInt *);
 273#endif
 274
 275/* Local macros */
 276/* masked special-values bits */
 277#define SPECIALARG  (rhs->bits & DECSPECIAL)
 278#define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)
 279
 280/* Diagnostic macros, etc. */
 281#if DECALLOC
 282/* Handle malloc/free accounting.  If enabled, our accountable routines */
 283/* are used; otherwise the code just goes straight to the system malloc */
 284/* and free routines. */
 285#define malloc(a) decMalloc(a)
 286#define free(a) decFree(a)
 287#define DECFENCE 0x5a              /* corruption detector */
 288/* 'Our' malloc and free: */
 289static void *decMalloc(size_t);
 290static void  decFree(void *);
 291uInt decAllocBytes=0;              /* count of bytes allocated */
 292/* Note that DECALLOC code only checks for storage buffer overflow. */
 293/* To check for memory leaks, the decAllocBytes variable must be */
 294/* checked to be 0 at appropriate times (e.g., after the test */
 295/* harness completes a set of tests).  This checking may be unreliable */
 296/* if the testing is done in a multi-thread environment. */
 297#endif
 298
 299#if DECCHECK
 300/* Optional checking routines.  Enabling these means that decNumber */
 301/* and decContext operands to operator routines are checked for */
 302/* correctness.  This roughly doubles the execution time of the */
 303/* fastest routines (and adds 600+ bytes), so should not normally be */
 304/* used in 'production'. */
 305/* decCheckInexact is used to check that inexact results have a full */
 306/* complement of digits (where appropriate -- this is not the case */
 307/* for Quantize, for example) */
 308#define DECUNRESU ((decNumber *)(void *)0xffffffff)
 309#define DECUNUSED ((const decNumber *)(void *)0xffffffff)
 310#define DECUNCONT ((decContext *)(void *)(0xffffffff))
 311static Flag decCheckOperands(decNumber *, const decNumber *,
 312                             const decNumber *, decContext *);
 313static Flag decCheckNumber(const decNumber *);
 314static void decCheckInexact(const decNumber *, decContext *);
 315#endif
 316
 317#if DECTRACE || DECCHECK
 318/* Optional trace/debugging routines (may or may not be used) */
 319void decNumberShow(const decNumber *);  /* displays the components of a number */
 320static void decDumpAr(char, const Unit *, Int);
 321#endif
 322
 323/* ================================================================== */
 324/* Conversions                                                        */
 325/* ================================================================== */
 326
 327/* ------------------------------------------------------------------ */
 328/* from-int32 -- conversion from Int or uInt                          */
 329/*                                                                    */
 330/*  dn is the decNumber to receive the integer                        */
 331/*  in or uin is the integer to be converted                          */
 332/*  returns dn                                                        */
 333/*                                                                    */
 334/* No error is possible.                                              */
 335/* ------------------------------------------------------------------ */
 336decNumber * decNumberFromInt32(decNumber *dn, Int in) {
 337  uInt unsig;
 338  if (in>=0) unsig=in;
 339   else {                               /* negative (possibly BADINT) */
 340    if (in==BADINT) unsig=(uInt)1073741824*2; /* special case */
 341     else unsig=-in;                    /* invert */
 342    }
 343  /* in is now positive */
 344  decNumberFromUInt32(dn, unsig);
 345  if (in<0) dn->bits=DECNEG;            /* sign needed */
 346  return dn;
 347  } /* decNumberFromInt32 */
 348
 349decNumber * decNumberFromUInt32(decNumber *dn, uInt uin) {
 350  Unit *up;                             /* work pointer */
 351  decNumberZero(dn);                    /* clean */
 352  if (uin==0) return dn;                /* [or decGetDigits bad call] */
 353  for (up=dn->lsu; uin>0; up++) {
 354    *up=(Unit)(uin%(DECDPUNMAX+1));
 355    uin=uin/(DECDPUNMAX+1);
 356    }
 357  dn->digits=decGetDigits(dn->lsu, up-dn->lsu);
 358  return dn;
 359  } /* decNumberFromUInt32 */
 360
 361/* ------------------------------------------------------------------ */
 362/* to-int32 -- conversion to Int or uInt                              */
 363/*                                                                    */
 364/*  dn is the decNumber to convert                                    */
 365/*  set is the context for reporting errors                           */
 366/*  returns the converted decNumber, or 0 if Invalid is set           */
 367/*                                                                    */
 368/* Invalid is set if the decNumber does not have exponent==0 or if    */
 369/* it is a NaN, Infinite, or out-of-range.                            */
 370/* ------------------------------------------------------------------ */
 371Int decNumberToInt32(const decNumber *dn, decContext *set) {
 372  #if DECCHECK
 373  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
 374  #endif
 375
 376  /* special or too many digits, or bad exponent */
 377  if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0) ; /* bad */
 378   else { /* is a finite integer with 10 or fewer digits */
 379    Int d;                         /* work */
 380    const Unit *up;                /* .. */
 381    uInt hi=0, lo;                 /* .. */
 382    up=dn->lsu;                    /* -> lsu */
 383    lo=*up;                        /* get 1 to 9 digits */
 384    #if DECDPUN>1                  /* split to higher */
 385      hi=lo/10;
 386      lo=lo%10;
 387    #endif
 388    up++;
 389    /* collect remaining Units, if any, into hi */
 390    for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];
 391    /* now low has the lsd, hi the remainder */
 392    if (hi>214748364 || (hi==214748364 && lo>7)) { /* out of range? */
 393      /* most-negative is a reprieve */
 394      if (dn->bits&DECNEG && hi==214748364 && lo==8) return 0x80000000;
 395      /* bad -- drop through */
 396      }
 397     else { /* in-range always */
 398      Int i=X10(hi)+lo;
 399      if (dn->bits&DECNEG) return -i;
 400      return i;
 401      }
 402    } /* integer */
 403  decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
 404  return 0;
 405  } /* decNumberToInt32 */
 406
 407uInt decNumberToUInt32(const decNumber *dn, decContext *set) {
 408  #if DECCHECK
 409  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
 410  #endif
 411  /* special or too many digits, or bad exponent, or negative (<0) */
 412  if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0
 413    || (dn->bits&DECNEG && !ISZERO(dn)));                   /* bad */
 414   else { /* is a finite integer with 10 or fewer digits */
 415    Int d;                         /* work */
 416    const Unit *up;                /* .. */
 417    uInt hi=0, lo;                 /* .. */
 418    up=dn->lsu;                    /* -> lsu */
 419    lo=*up;                        /* get 1 to 9 digits */
 420    #if DECDPUN>1                  /* split to higher */
 421      hi=lo/10;
 422      lo=lo%10;
 423    #endif
 424    up++;
 425    /* collect remaining Units, if any, into hi */
 426    for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];
 427
 428    /* now low has the lsd, hi the remainder */
 429    if (hi>429496729 || (hi==429496729 && lo>5)) ; /* no reprieve possible */
 430     else return X10(hi)+lo;
 431    } /* integer */
 432  decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
 433  return 0;
 434  } /* decNumberToUInt32 */
 435
 436decNumber *decNumberFromInt64(decNumber *dn, int64_t in)
 437{
 438    uint64_t unsig = in;
 439    if (in < 0) {
 440        unsig = -unsig;
 441    }
 442
 443    decNumberFromUInt64(dn, unsig);
 444    if (in < 0) {
 445        dn->bits = DECNEG;        /* sign needed */
 446    }
 447    return dn;
 448} /* decNumberFromInt64 */
 449
 450decNumber *decNumberFromUInt64(decNumber *dn, uint64_t uin)
 451{
 452    Unit *up;                             /* work pointer */
 453    decNumberZero(dn);                    /* clean */
 454    if (uin == 0) {
 455        return dn;                /* [or decGetDigits bad call] */
 456    }
 457    for (up = dn->lsu; uin > 0; up++) {
 458        *up = (Unit)(uin % (DECDPUNMAX + 1));
 459        uin = uin / (DECDPUNMAX + 1);
 460    }
 461    dn->digits = decGetDigits(dn->lsu, up-dn->lsu);
 462    return dn;
 463} /* decNumberFromUInt64 */
 464
 465/* ------------------------------------------------------------------ */
 466/* to-int64 -- conversion to int64                                    */
 467/*                                                                    */
 468/*  dn is the decNumber to convert.  dn is assumed to have been       */
 469/*    rounded to a floating point integer value.                      */
 470/*  set is the context for reporting errors                           */
 471/*  returns the converted decNumber, or 0 if Invalid is set           */
 472/*                                                                    */
 473/* Invalid is set if the decNumber is a NaN, Infinite or is out of    */
 474/* range for a signed 64 bit integer.                                 */
 475/* ------------------------------------------------------------------ */
 476
 477int64_t decNumberIntegralToInt64(const decNumber *dn, decContext *set)
 478{
 479    if (decNumberIsSpecial(dn) || (dn->exponent < 0) ||
 480       (dn->digits + dn->exponent > 19)) {
 481        goto Invalid;
 482    } else {
 483        int64_t d;        /* work */
 484        const Unit *up;   /* .. */
 485        uint64_t hi = 0;
 486        up = dn->lsu;     /* -> lsu */
 487
 488        for (d = 1; d <= dn->digits; up++, d += DECDPUN) {
 489            uint64_t prev = hi;
 490            hi += *up * powers[d-1];
 491            if ((hi < prev) || (hi > INT64_MAX)) {
 492                goto Invalid;
 493            }
 494        }
 495
 496        uint64_t prev = hi;
 497        hi *= (uint64_t)powers[dn->exponent];
 498        if ((hi < prev) || (hi > INT64_MAX)) {
 499            goto Invalid;
 500        }
 501        return (decNumberIsNegative(dn)) ? -((int64_t)hi) : (int64_t)hi;
 502    }
 503
 504Invalid:
 505    decContextSetStatus(set, DEC_Invalid_operation);
 506    return 0;
 507} /* decNumberIntegralToInt64 */
 508
 509
 510/* ------------------------------------------------------------------ */
 511/* to-scientific-string -- conversion to numeric string               */
 512/* to-engineering-string -- conversion to numeric string              */
 513/*                                                                    */
 514/*   decNumberToString(dn, string);                                   */
 515/*   decNumberToEngString(dn, string);                                */
 516/*                                                                    */
 517/*  dn is the decNumber to convert                                    */
 518/*  string is the string where the result will be laid out            */
 519/*                                                                    */
 520/*  string must be at least dn->digits+14 characters long             */
 521/*                                                                    */
 522/*  No error is possible, and no status can be set.                   */
 523/* ------------------------------------------------------------------ */
 524char * decNumberToString(const decNumber *dn, char *string){
 525  decToString(dn, string, 0);
 526  return string;
 527  } /* DecNumberToString */
 528
 529char * decNumberToEngString(const decNumber *dn, char *string){
 530  decToString(dn, string, 1);
 531  return string;
 532  } /* DecNumberToEngString */
 533
 534/* ------------------------------------------------------------------ */
 535/* to-number -- conversion from numeric string                        */
 536/*                                                                    */
 537/* decNumberFromString -- convert string to decNumber                 */
 538/*   dn        -- the number structure to fill                        */
 539/*   chars[]   -- the string to convert ('\0' terminated)             */
 540/*   set       -- the context used for processing any error,          */
 541/*                determining the maximum precision available         */
 542/*                (set.digits), determining the maximum and minimum   */
 543/*                exponent (set.emax and set.emin), determining if    */
 544/*                extended values are allowed, and checking the       */
 545/*                rounding mode if overflow occurs or rounding is     */
 546/*                needed.                                             */
 547/*                                                                    */
 548/* The length of the coefficient and the size of the exponent are     */
 549/* checked by this routine, so the correct error (Underflow or        */
 550/* Overflow) can be reported or rounding applied, as necessary.       */
 551/*                                                                    */
 552/* If bad syntax is detected, the result will be a quiet NaN.         */
 553/* ------------------------------------------------------------------ */
 554decNumber * decNumberFromString(decNumber *dn, const char chars[],
 555                                decContext *set) {
 556  Int   exponent=0;                /* working exponent [assume 0] */
 557  uByte bits=0;                    /* working flags [assume +ve] */
 558  Unit  *res;                      /* where result will be built */
 559  Unit  resbuff[SD2U(DECBUFFER+9)];/* local buffer in case need temporary */
 560                                   /* [+9 allows for ln() constants] */
 561  Unit  *allocres=NULL;            /* -> allocated result, iff allocated */
 562  Int   d=0;                       /* count of digits found in decimal part */
 563  const char *dotchar=NULL;        /* where dot was found */
 564  const char *cfirst=chars;        /* -> first character of decimal part */
 565  const char *last=NULL;           /* -> last digit of decimal part */
 566  const char *c;                   /* work */
 567  Unit  *up;                       /* .. */
 568  #if DECDPUN>1
 569  Int   cut, out;                  /* .. */
 570  #endif
 571  Int   residue;                   /* rounding residue */
 572  uInt  status=0;                  /* error code */
 573
 574  #if DECCHECK
 575  if (decCheckOperands(DECUNRESU, DECUNUSED, DECUNUSED, set))
 576    return decNumberZero(dn);
 577  #endif
 578
 579  do {                             /* status & malloc protection */
 580    for (c=chars;; c++) {          /* -> input character */
 581      if (*c>='0' && *c<='9') {    /* test for Arabic digit */
 582        last=c;
 583        d++;                       /* count of real digits */
 584        continue;                  /* still in decimal part */
 585        }
 586      if (*c=='.' && dotchar==NULL) { /* first '.' */
 587        dotchar=c;                 /* record offset into decimal part */
 588        if (c==cfirst) cfirst++;   /* first digit must follow */
 589        continue;}
 590      if (c==chars) {              /* first in string... */
 591        if (*c=='-') {             /* valid - sign */
 592          cfirst++;
 593          bits=DECNEG;
 594          continue;}
 595        if (*c=='+') {             /* valid + sign */
 596          cfirst++;
 597          continue;}
 598        }
 599      /* *c is not a digit, or a valid +, -, or '.' */
 600      break;
 601      } /* c */
 602
 603    if (last==NULL) {              /* no digits yet */
 604      status=DEC_Conversion_syntax;/* assume the worst */
 605      if (*c=='\0') break;         /* and no more to come... */
 606      #if DECSUBSET
 607      /* if subset then infinities and NaNs are not allowed */
 608      if (!set->extended) break;   /* hopeless */
 609      #endif
 610      /* Infinities and NaNs are possible, here */
 611      if (dotchar!=NULL) break;    /* .. unless had a dot */
 612      decNumberZero(dn);           /* be optimistic */
 613      if (decBiStr(c, "infinity", "INFINITY")
 614       || decBiStr(c, "inf", "INF")) {
 615        dn->bits=bits | DECINF;
 616        status=0;                  /* is OK */
 617        break; /* all done */
 618        }
 619      /* a NaN expected */
 620      /* 2003.09.10 NaNs are now permitted to have a sign */
 621      dn->bits=bits | DECNAN;      /* assume simple NaN */
 622      if (*c=='s' || *c=='S') {    /* looks like an sNaN */
 623        c++;
 624        dn->bits=bits | DECSNAN;
 625        }
 626      if (*c!='n' && *c!='N') break;    /* check caseless "NaN" */
 627      c++;
 628      if (*c!='a' && *c!='A') break;    /* .. */
 629      c++;
 630      if (*c!='n' && *c!='N') break;    /* .. */
 631      c++;
 632      /* now either nothing, or nnnn payload, expected */
 633      /* -> start of integer and skip leading 0s [including plain 0] */
 634      for (cfirst=c; *cfirst=='0';) cfirst++;
 635      if (*cfirst=='\0') {         /* "NaN" or "sNaN", maybe with all 0s */
 636        status=0;                  /* it's good */
 637        break;                     /* .. */
 638        }
 639      /* something other than 0s; setup last and d as usual [no dots] */
 640      for (c=cfirst;; c++, d++) {
 641        if (*c<'0' || *c>'9') break; /* test for Arabic digit */
 642        last=c;
 643        }
 644      if (*c!='\0') break;         /* not all digits */
 645      if (d>set->digits-1) {
 646        /* [NB: payload in a decNumber can be full length unless */
 647        /* clamped, in which case can only be digits-1] */
 648        if (set->clamp) break;
 649        if (d>set->digits) break;
 650        } /* too many digits? */
 651      /* good; drop through to convert the integer to coefficient */
 652      status=0;                    /* syntax is OK */
 653      bits=dn->bits;               /* for copy-back */
 654      } /* last==NULL */
 655
 656     else if (*c!='\0') {          /* more to process... */
 657      /* had some digits; exponent is only valid sequence now */
 658      Flag nege;                   /* 1=negative exponent */
 659      const char *firstexp;        /* -> first significant exponent digit */
 660      status=DEC_Conversion_syntax;/* assume the worst */
 661      if (*c!='e' && *c!='E') break;
 662      /* Found 'e' or 'E' -- now process explicit exponent */
 663      /* 1998.07.11: sign no longer required */
 664      nege=0;
 665      c++;                         /* to (possible) sign */
 666      if (*c=='-') {nege=1; c++;}
 667       else if (*c=='+') c++;
 668      if (*c=='\0') break;
 669
 670      for (; *c=='0' && *(c+1)!='\0';) c++;  /* strip insignificant zeros */
 671      firstexp=c;                            /* save exponent digit place */
 672      for (; ;c++) {
 673        if (*c<'0' || *c>'9') break;         /* not a digit */
 674        exponent=X10(exponent)+(Int)*c-(Int)'0';
 675        } /* c */
 676      /* if not now on a '\0', *c must not be a digit */
 677      if (*c!='\0') break;
 678
 679      /* (this next test must be after the syntax checks) */
 680      /* if it was too long the exponent may have wrapped, so check */
 681      /* carefully and set it to a certain overflow if wrap possible */
 682      if (c>=firstexp+9+1) {
 683        if (c>firstexp+9+1 || *firstexp>'1') exponent=DECNUMMAXE*2;
 684        /* [up to 1999999999 is OK, for example 1E-1000000998] */
 685        }
 686      if (nege) exponent=-exponent;     /* was negative */
 687      status=0;                         /* is OK */
 688      } /* stuff after digits */
 689
 690    /* Here when whole string has been inspected; syntax is good */
 691    /* cfirst->first digit (never dot), last->last digit (ditto) */
 692
 693    /* strip leading zeros/dot [leave final 0 if all 0's] */
 694    if (*cfirst=='0') {                 /* [cfirst has stepped over .] */
 695      for (c=cfirst; c<last; c++, cfirst++) {
 696        if (*c=='.') continue;          /* ignore dots */
 697        if (*c!='0') break;             /* non-zero found */
 698        d--;                            /* 0 stripped */
 699        } /* c */
 700      #if DECSUBSET
 701      /* make a rapid exit for easy zeros if !extended */
 702      if (*cfirst=='0' && !set->extended) {
 703        decNumberZero(dn);              /* clean result */
 704        break;                          /* [could be return] */
 705        }
 706      #endif
 707      } /* at least one leading 0 */
 708
 709    /* Handle decimal point... */
 710    if (dotchar!=NULL && dotchar<last)  /* non-trailing '.' found? */
 711      exponent-=(last-dotchar);         /* adjust exponent */
 712    /* [we can now ignore the .] */
 713
 714    /* OK, the digits string is good.  Assemble in the decNumber, or in */
 715    /* a temporary units array if rounding is needed */
 716    if (d<=set->digits) res=dn->lsu;    /* fits into supplied decNumber */
 717     else {                             /* rounding needed */
 718      Int needbytes=D2U(d)*sizeof(Unit);/* bytes needed */
 719      res=resbuff;                      /* assume use local buffer */
 720      if (needbytes>(Int)sizeof(resbuff)) { /* too big for local */
 721        allocres=(Unit *)malloc(needbytes);
 722        if (allocres==NULL) {status|=DEC_Insufficient_storage; break;}
 723        res=allocres;
 724        }
 725      }
 726    /* res now -> number lsu, buffer, or allocated storage for Unit array */
 727
 728    /* Place the coefficient into the selected Unit array */
 729    /* [this is often 70% of the cost of this function when DECDPUN>1] */
 730    #if DECDPUN>1
 731    out=0;                         /* accumulator */
 732    up=res+D2U(d)-1;               /* -> msu */
 733    cut=d-(up-res)*DECDPUN;        /* digits in top unit */
 734    for (c=cfirst;; c++) {         /* along the digits */
 735      if (*c=='.') continue;       /* ignore '.' [don't decrement cut] */
 736      out=X10(out)+(Int)*c-(Int)'0';
 737      if (c==last) break;          /* done [never get to trailing '.'] */
 738      cut--;
 739      if (cut>0) continue;         /* more for this unit */
 740      *up=(Unit)out;               /* write unit */
 741      up--;                        /* prepare for unit below.. */
 742      cut=DECDPUN;                 /* .. */
 743      out=0;                       /* .. */
 744      } /* c */
 745    *up=(Unit)out;                 /* write lsu */
 746
 747    #else
 748    /* DECDPUN==1 */
 749    up=res;                        /* -> lsu */
 750    for (c=last; c>=cfirst; c--) { /* over each character, from least */
 751      if (*c=='.') continue;       /* ignore . [don't step up] */
 752      *up=(Unit)((Int)*c-(Int)'0');
 753      up++;
 754      } /* c */
 755    #endif
 756
 757    dn->bits=bits;
 758    dn->exponent=exponent;
 759    dn->digits=d;
 760
 761    /* if not in number (too long) shorten into the number */
 762    if (d>set->digits) {
 763      residue=0;
 764      decSetCoeff(dn, set, res, d, &residue, &status);
 765      /* always check for overflow or subnormal and round as needed */
 766      decFinalize(dn, set, &residue, &status);
 767      }
 768     else { /* no rounding, but may still have overflow or subnormal */
 769      /* [these tests are just for performance; finalize repeats them] */
 770      if ((dn->exponent-1<set->emin-dn->digits)
 771       || (dn->exponent-1>set->emax-set->digits)) {
 772        residue=0;
 773        decFinalize(dn, set, &residue, &status);
 774        }
 775      }
 776    /* decNumberShow(dn); */
 777    } while(0);                         /* [for break] */
 778
 779  if (allocres!=NULL) free(allocres);   /* drop any storage used */
 780  if (status!=0) decStatus(dn, status, set);
 781  return dn;
 782  } /* decNumberFromString */
 783
 784/* ================================================================== */
 785/* Operators                                                          */
 786/* ================================================================== */
 787
 788/* ------------------------------------------------------------------ */
 789/* decNumberAbs -- absolute value operator                            */
 790/*                                                                    */
 791/*   This computes C = abs(A)                                         */
 792/*                                                                    */
 793/*   res is C, the result.  C may be A                                */
 794/*   rhs is A                                                         */
 795/*   set is the context                                               */
 796/*                                                                    */
 797/* See also decNumberCopyAbs for a quiet bitwise version of this.     */
 798/* C must have space for set->digits digits.                          */
 799/* ------------------------------------------------------------------ */
 800/* This has the same effect as decNumberPlus unless A is negative,    */
 801/* in which case it has the same effect as decNumberMinus.            */
 802/* ------------------------------------------------------------------ */
 803decNumber * decNumberAbs(decNumber *res, const decNumber *rhs,
 804                         decContext *set) {
 805  decNumber dzero;                      /* for 0 */
 806  uInt status=0;                        /* accumulator */
 807
 808  #if DECCHECK
 809  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
 810  #endif
 811
 812  decNumberZero(&dzero);                /* set 0 */
 813  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
 814  decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status);
 815  if (status!=0) decStatus(res, status, set);
 816  #if DECCHECK
 817  decCheckInexact(res, set);
 818  #endif
 819  return res;
 820  } /* decNumberAbs */
 821
 822/* ------------------------------------------------------------------ */
 823/* decNumberAdd -- add two Numbers                                    */
 824/*                                                                    */
 825/*   This computes C = A + B                                          */
 826/*                                                                    */
 827/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
 828/*   lhs is A                                                         */
 829/*   rhs is B                                                         */
 830/*   set is the context                                               */
 831/*                                                                    */
 832/* C must have space for set->digits digits.                          */
 833/* ------------------------------------------------------------------ */
 834/* This just calls the routine shared with Subtract                   */
 835decNumber * decNumberAdd(decNumber *res, const decNumber *lhs,
 836                         const decNumber *rhs, decContext *set) {
 837  uInt status=0;                        /* accumulator */
 838  decAddOp(res, lhs, rhs, set, 0, &status);
 839  if (status!=0) decStatus(res, status, set);
 840  #if DECCHECK
 841  decCheckInexact(res, set);
 842  #endif
 843  return res;
 844  } /* decNumberAdd */
 845
 846/* ------------------------------------------------------------------ */
 847/* decNumberAnd -- AND two Numbers, digitwise                         */
 848/*                                                                    */
 849/*   This computes C = A & B                                          */
 850/*                                                                    */
 851/*   res is C, the result.  C may be A and/or B (e.g., X=X&X)         */
 852/*   lhs is A                                                         */
 853/*   rhs is B                                                         */
 854/*   set is the context (used for result length and error report)     */
 855/*                                                                    */
 856/* C must have space for set->digits digits.                          */
 857/*                                                                    */
 858/* Logical function restrictions apply (see above); a NaN is          */
 859/* returned with Invalid_operation if a restriction is violated.      */
 860/* ------------------------------------------------------------------ */
 861decNumber * decNumberAnd(decNumber *res, const decNumber *lhs,
 862                         const decNumber *rhs, decContext *set) {
 863  const Unit *ua, *ub;                  /* -> operands */
 864  const Unit *msua, *msub;              /* -> operand msus */
 865  Unit *uc,  *msuc;                     /* -> result and its msu */
 866  Int   msudigs;                        /* digits in res msu */
 867  #if DECCHECK
 868  if (decCheckOperands(res, lhs, rhs, set)) return res;
 869  #endif
 870
 871  if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
 872   || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
 873    decStatus(res, DEC_Invalid_operation, set);
 874    return res;
 875    }
 876
 877  /* operands are valid */
 878  ua=lhs->lsu;                          /* bottom-up */
 879  ub=rhs->lsu;                          /* .. */
 880  uc=res->lsu;                          /* .. */
 881  msua=ua+D2U(lhs->digits)-1;           /* -> msu of lhs */
 882  msub=ub+D2U(rhs->digits)-1;           /* -> msu of rhs */
 883  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
 884  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
 885  for (; uc<=msuc; ua++, ub++, uc++) {  /* Unit loop */
 886    Unit a, b;                          /* extract units */
 887    if (ua>msua) a=0;
 888     else a=*ua;
 889    if (ub>msub) b=0;
 890     else b=*ub;
 891    *uc=0;                              /* can now write back */
 892    if (a|b) {                          /* maybe 1 bits to examine */
 893      Int i, j;
 894      *uc=0;                            /* can now write back */
 895      /* This loop could be unrolled and/or use BIN2BCD tables */
 896      for (i=0; i<DECDPUN; i++) {
 897        if (a&b&1) *uc=*uc+(Unit)powers[i];  /* effect AND */
 898        j=a%10;
 899        a=a/10;
 900        j|=b%10;
 901        b=b/10;
 902        if (j>1) {
 903          decStatus(res, DEC_Invalid_operation, set);
 904          return res;
 905          }
 906        if (uc==msuc && i==msudigs-1) break; /* just did final digit */
 907        } /* each digit */
 908      } /* both OK */
 909    } /* each unit */
 910  /* [here uc-1 is the msu of the result] */
 911  res->digits=decGetDigits(res->lsu, uc-res->lsu);
 912  res->exponent=0;                      /* integer */
 913  res->bits=0;                          /* sign=0 */
 914  return res;  /* [no status to set] */
 915  } /* decNumberAnd */
 916
 917/* ------------------------------------------------------------------ */
 918/* decNumberCompare -- compare two Numbers                            */
 919/*                                                                    */
 920/*   This computes C = A ? B                                          */
 921/*                                                                    */
 922/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
 923/*   lhs is A                                                         */
 924/*   rhs is B                                                         */
 925/*   set is the context                                               */
 926/*                                                                    */
 927/* C must have space for one digit (or NaN).                          */
 928/* ------------------------------------------------------------------ */
 929decNumber * decNumberCompare(decNumber *res, const decNumber *lhs,
 930                             const decNumber *rhs, decContext *set) {
 931  uInt status=0;                        /* accumulator */
 932  decCompareOp(res, lhs, rhs, set, COMPARE, &status);
 933  if (status!=0) decStatus(res, status, set);
 934  return res;
 935  } /* decNumberCompare */
 936
 937/* ------------------------------------------------------------------ */
 938/* decNumberCompareSignal -- compare, signalling on all NaNs          */
 939/*                                                                    */
 940/*   This computes C = A ? B                                          */
 941/*                                                                    */
 942/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
 943/*   lhs is A                                                         */
 944/*   rhs is B                                                         */
 945/*   set is the context                                               */
 946/*                                                                    */
 947/* C must have space for one digit (or NaN).                          */
 948/* ------------------------------------------------------------------ */
 949decNumber * decNumberCompareSignal(decNumber *res, const decNumber *lhs,
 950                                   const decNumber *rhs, decContext *set) {
 951  uInt status=0;                        /* accumulator */
 952  decCompareOp(res, lhs, rhs, set, COMPSIG, &status);
 953  if (status!=0) decStatus(res, status, set);
 954  return res;
 955  } /* decNumberCompareSignal */
 956
 957/* ------------------------------------------------------------------ */
 958/* decNumberCompareTotal -- compare two Numbers, using total ordering */
 959/*                                                                    */
 960/*   This computes C = A ? B, under total ordering                    */
 961/*                                                                    */
 962/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
 963/*   lhs is A                                                         */
 964/*   rhs is B                                                         */
 965/*   set is the context                                               */
 966/*                                                                    */
 967/* C must have space for one digit; the result will always be one of  */
 968/* -1, 0, or 1.                                                       */
 969/* ------------------------------------------------------------------ */
 970decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs,
 971                                  const decNumber *rhs, decContext *set) {
 972  uInt status=0;                        /* accumulator */
 973  decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
 974  if (status!=0) decStatus(res, status, set);
 975  return res;
 976  } /* decNumberCompareTotal */
 977
 978/* ------------------------------------------------------------------ */
 979/* decNumberCompareTotalMag -- compare, total ordering of magnitudes  */
 980/*                                                                    */
 981/*   This computes C = |A| ? |B|, under total ordering                */
 982/*                                                                    */
 983/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
 984/*   lhs is A                                                         */
 985/*   rhs is B                                                         */
 986/*   set is the context                                               */
 987/*                                                                    */
 988/* C must have space for one digit; the result will always be one of  */
 989/* -1, 0, or 1.                                                       */
 990/* ------------------------------------------------------------------ */
 991decNumber * decNumberCompareTotalMag(decNumber *res, const decNumber *lhs,
 992                                     const decNumber *rhs, decContext *set) {
 993  uInt status=0;                   /* accumulator */
 994  uInt needbytes;                  /* for space calculations */
 995  decNumber bufa[D2N(DECBUFFER+1)];/* +1 in case DECBUFFER=0 */
 996  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
 997  decNumber bufb[D2N(DECBUFFER+1)];
 998  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
 999  decNumber *a, *b;                /* temporary pointers */
1000
1001  #if DECCHECK
1002  if (decCheckOperands(res, lhs, rhs, set)) return res;
1003  #endif
1004
1005  do {                                  /* protect allocated storage */
1006    /* if either is negative, take a copy and absolute */
1007    if (decNumberIsNegative(lhs)) {     /* lhs<0 */
1008      a=bufa;
1009      needbytes=sizeof(decNumber)+(D2U(lhs->digits)-1)*sizeof(Unit);
1010      if (needbytes>sizeof(bufa)) {     /* need malloc space */
1011        allocbufa=(decNumber *)malloc(needbytes);
1012        if (allocbufa==NULL) {          /* hopeless -- abandon */
1013          status|=DEC_Insufficient_storage;
1014          break;}
1015        a=allocbufa;                    /* use the allocated space */
1016        }
1017      decNumberCopy(a, lhs);            /* copy content */
1018      a->bits&=~DECNEG;                 /* .. and clear the sign */
1019      lhs=a;                            /* use copy from here on */
1020      }
1021    if (decNumberIsNegative(rhs)) {     /* rhs<0 */
1022      b=bufb;
1023      needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
1024      if (needbytes>sizeof(bufb)) {     /* need malloc space */
1025        allocbufb=(decNumber *)malloc(needbytes);
1026        if (allocbufb==NULL) {          /* hopeless -- abandon */
1027          status|=DEC_Insufficient_storage;
1028          break;}
1029        b=allocbufb;                    /* use the allocated space */
1030        }
1031      decNumberCopy(b, rhs);            /* copy content */
1032      b->bits&=~DECNEG;                 /* .. and clear the sign */
1033      rhs=b;                            /* use copy from here on */
1034      }
1035    decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
1036    } while(0);                         /* end protected */
1037
1038  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1039  if (allocbufb!=NULL) free(allocbufb); /* .. */
1040  if (status!=0) decStatus(res, status, set);
1041  return res;
1042  } /* decNumberCompareTotalMag */
1043
1044/* ------------------------------------------------------------------ */
1045/* decNumberDivide -- divide one number by another                    */
1046/*                                                                    */
1047/*   This computes C = A / B                                          */
1048/*                                                                    */
1049/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
1050/*   lhs is A                                                         */
1051/*   rhs is B                                                         */
1052/*   set is the context                                               */
1053/*                                                                    */
1054/* C must have space for set->digits digits.                          */
1055/* ------------------------------------------------------------------ */
1056decNumber * decNumberDivide(decNumber *res, const decNumber *lhs,
1057                            const decNumber *rhs, decContext *set) {
1058  uInt status=0;                        /* accumulator */
1059  decDivideOp(res, lhs, rhs, set, DIVIDE, &status);
1060  if (status!=0) decStatus(res, status, set);
1061  #if DECCHECK
1062  decCheckInexact(res, set);
1063  #endif
1064  return res;
1065  } /* decNumberDivide */
1066
1067/* ------------------------------------------------------------------ */
1068/* decNumberDivideInteger -- divide and return integer quotient       */
1069/*                                                                    */
1070/*   This computes C = A # B, where # is the integer divide operator  */
1071/*                                                                    */
1072/*   res is C, the result.  C may be A and/or B (e.g., X=X#X)         */
1073/*   lhs is A                                                         */
1074/*   rhs is B                                                         */
1075/*   set is the context                                               */
1076/*                                                                    */
1077/* C must have space for set->digits digits.                          */
1078/* ------------------------------------------------------------------ */
1079decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs,
1080                                   const decNumber *rhs, decContext *set) {
1081  uInt status=0;                        /* accumulator */
1082  decDivideOp(res, lhs, rhs, set, DIVIDEINT, &status);
1083  if (status!=0) decStatus(res, status, set);
1084  return res;
1085  } /* decNumberDivideInteger */
1086
1087/* ------------------------------------------------------------------ */
1088/* decNumberExp -- exponentiation                                     */
1089/*                                                                    */
1090/*   This computes C = exp(A)                                         */
1091/*                                                                    */
1092/*   res is C, the result.  C may be A                                */
1093/*   rhs is A                                                         */
1094/*   set is the context; note that rounding mode has no effect        */
1095/*                                                                    */
1096/* C must have space for set->digits digits.                          */
1097/*                                                                    */
1098/* Mathematical function restrictions apply (see above); a NaN is     */
1099/* returned with Invalid_operation if a restriction is violated.      */
1100/*                                                                    */
1101/* Finite results will always be full precision and Inexact, except   */
1102/* when A is a zero or -Infinity (giving 1 or 0 respectively).        */
1103/*                                                                    */
1104/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1105/* almost always be correctly rounded, but may be up to 1 ulp in      */
1106/* error in rare cases.                                               */
1107/* ------------------------------------------------------------------ */
1108/* This is a wrapper for decExpOp which can handle the slightly wider */
1109/* (double) range needed by Ln (which has to be able to calculate     */
1110/* exp(-a) where a can be the tiniest number (Ntiny).                 */
1111/* ------------------------------------------------------------------ */
1112decNumber * decNumberExp(decNumber *res, const decNumber *rhs,
1113                         decContext *set) {
1114  uInt status=0;                        /* accumulator */
1115  #if DECSUBSET
1116  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
1117  #endif
1118
1119  #if DECCHECK
1120  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1121  #endif
1122
1123  /* Check restrictions; these restrictions ensure that if h=8 (see */
1124  /* decExpOp) then the result will either overflow or underflow to 0. */
1125  /* Other math functions restrict the input range, too, for inverses. */
1126  /* If not violated then carry out the operation. */
1127  if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
1128    #if DECSUBSET
1129    if (!set->extended) {
1130      /* reduce operand and set lostDigits status, as needed */
1131      if (rhs->digits>set->digits) {
1132        allocrhs=decRoundOperand(rhs, set, &status);
1133        if (allocrhs==NULL) break;
1134        rhs=allocrhs;
1135        }
1136      }
1137    #endif
1138    decExpOp(res, rhs, set, &status);
1139    } while(0);                         /* end protected */
1140
1141  #if DECSUBSET
1142  if (allocrhs !=NULL) free(allocrhs);  /* drop any storage used */
1143  #endif
1144  /* apply significant status */
1145  if (status!=0) decStatus(res, status, set);
1146  #if DECCHECK
1147  decCheckInexact(res, set);
1148  #endif
1149  return res;
1150  } /* decNumberExp */
1151
1152/* ------------------------------------------------------------------ */
1153/* decNumberFMA -- fused multiply add                                 */
1154/*                                                                    */
1155/*   This computes D = (A * B) + C with only one rounding             */
1156/*                                                                    */
1157/*   res is D, the result.  D may be A or B or C (e.g., X=FMA(X,X,X)) */
1158/*   lhs is A                                                         */
1159/*   rhs is B                                                         */
1160/*   fhs is C [far hand side]                                         */
1161/*   set is the context                                               */
1162/*                                                                    */
1163/* Mathematical function restrictions apply (see above); a NaN is     */
1164/* returned with Invalid_operation if a restriction is violated.      */
1165/*                                                                    */
1166/* C must have space for set->digits digits.                          */
1167/* ------------------------------------------------------------------ */
1168decNumber * decNumberFMA(decNumber *res, const decNumber *lhs,
1169                         const decNumber *rhs, const decNumber *fhs,
1170                         decContext *set) {
1171  uInt status=0;                   /* accumulator */
1172  decContext dcmul;                /* context for the multiplication */
1173  uInt needbytes;                  /* for space calculations */
1174  decNumber bufa[D2N(DECBUFFER*2+1)];
1175  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
1176  decNumber *acc;                  /* accumulator pointer */
1177  decNumber dzero;                 /* work */
1178
1179  #if DECCHECK
1180  if (decCheckOperands(res, lhs, rhs, set)) return res;
1181  if (decCheckOperands(res, fhs, DECUNUSED, set)) return res;
1182  #endif
1183
1184  do {                                  /* protect allocated storage */
1185    #if DECSUBSET
1186    if (!set->extended) {               /* [undefined if subset] */
1187      status|=DEC_Invalid_operation;
1188      break;}
1189    #endif
1190    /* Check math restrictions [these ensure no overflow or underflow] */
1191    if ((!decNumberIsSpecial(lhs) && decCheckMath(lhs, set, &status))
1192     || (!decNumberIsSpecial(rhs) && decCheckMath(rhs, set, &status))
1193     || (!decNumberIsSpecial(fhs) && decCheckMath(fhs, set, &status))) break;
1194    /* set up context for multiply */
1195    dcmul=*set;
1196    dcmul.digits=lhs->digits+rhs->digits; /* just enough */
1197    /* [The above may be an over-estimate for subset arithmetic, but that's OK] */
1198    dcmul.emax=DEC_MAX_EMAX;            /* effectively unbounded .. */
1199    dcmul.emin=DEC_MIN_EMIN;            /* [thanks to Math restrictions] */
1200    /* set up decNumber space to receive the result of the multiply */
1201    acc=bufa;                           /* may fit */
1202    needbytes=sizeof(decNumber)+(D2U(dcmul.digits)-1)*sizeof(Unit);
1203    if (needbytes>sizeof(bufa)) {       /* need malloc space */
1204      allocbufa=(decNumber *)malloc(needbytes);
1205      if (allocbufa==NULL) {            /* hopeless -- abandon */
1206        status|=DEC_Insufficient_storage;
1207        break;}
1208      acc=allocbufa;                    /* use the allocated space */
1209      }
1210    /* multiply with extended range and necessary precision */
1211    /*printf("emin=%ld\n", dcmul.emin); */
1212    decMultiplyOp(acc, lhs, rhs, &dcmul, &status);
1213    /* Only Invalid operation (from sNaN or Inf * 0) is possible in */
1214    /* status; if either is seen than ignore fhs (in case it is */
1215    /* another sNaN) and set acc to NaN unless we had an sNaN */
1216    /* [decMultiplyOp leaves that to caller] */
1217    /* Note sNaN has to go through addOp to shorten payload if */
1218    /* necessary */
1219    if ((status&DEC_Invalid_operation)!=0) {
1220      if (!(status&DEC_sNaN)) {         /* but be true invalid */
1221        decNumberZero(res);             /* acc not yet set */
1222        res->bits=DECNAN;
1223        break;
1224        }
1225      decNumberZero(&dzero);            /* make 0 (any non-NaN would do) */
1226      fhs=&dzero;                       /* use that */
1227      }
1228    #if DECCHECK
1229     else { /* multiply was OK */
1230      if (status!=0) printf("Status=%08lx after FMA multiply\n", status);
1231      }
1232    #endif
1233    /* add the third operand and result -> res, and all is done */
1234    decAddOp(res, acc, fhs, set, 0, &status);
1235    } while(0);                         /* end protected */
1236
1237  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1238  if (status!=0) decStatus(res, status, set);
1239  #if DECCHECK
1240  decCheckInexact(res, set);
1241  #endif
1242  return res;
1243  } /* decNumberFMA */
1244
1245/* ------------------------------------------------------------------ */
1246/* decNumberInvert -- invert a Number, digitwise                      */
1247/*                                                                    */
1248/*   This computes C = ~A                                             */
1249/*                                                                    */
1250/*   res is C, the result.  C may be A (e.g., X=~X)                   */
1251/*   rhs is A                                                         */
1252/*   set is the context (used for result length and error report)     */
1253/*                                                                    */
1254/* C must have space for set->digits digits.                          */
1255/*                                                                    */
1256/* Logical function restrictions apply (see above); a NaN is          */
1257/* returned with Invalid_operation if a restriction is violated.      */
1258/* ------------------------------------------------------------------ */
1259decNumber * decNumberInvert(decNumber *res, const decNumber *rhs,
1260                            decContext *set) {
1261  const Unit *ua, *msua;                /* -> operand and its msu */
1262  Unit  *uc, *msuc;                     /* -> result and its msu */
1263  Int   msudigs;                        /* digits in res msu */
1264  #if DECCHECK
1265  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1266  #endif
1267
1268  if (rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
1269    decStatus(res, DEC_Invalid_operation, set);
1270    return res;
1271    }
1272  /* operand is valid */
1273  ua=rhs->lsu;                          /* bottom-up */
1274  uc=res->lsu;                          /* .. */
1275  msua=ua+D2U(rhs->digits)-1;           /* -> msu of rhs */
1276  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
1277  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
1278  for (; uc<=msuc; ua++, uc++) {        /* Unit loop */
1279    Unit a;                             /* extract unit */
1280    Int  i, j;                          /* work */
1281    if (ua>msua) a=0;
1282     else a=*ua;
1283    *uc=0;                              /* can now write back */
1284    /* always need to examine all bits in rhs */
1285    /* This loop could be unrolled and/or use BIN2BCD tables */
1286    for (i=0; i<DECDPUN; i++) {
1287      if ((~a)&1) *uc=*uc+(Unit)powers[i];   /* effect INVERT */
1288      j=a%10;
1289      a=a/10;
1290      if (j>1) {
1291        decStatus(res, DEC_Invalid_operation, set);
1292        return res;
1293        }
1294      if (uc==msuc && i==msudigs-1) break;   /* just did final digit */
1295      } /* each digit */
1296    } /* each unit */
1297  /* [here uc-1 is the msu of the result] */
1298  res->digits=decGetDigits(res->lsu, uc-res->lsu);
1299  res->exponent=0;                      /* integer */
1300  res->bits=0;                          /* sign=0 */
1301  return res;  /* [no status to set] */
1302  } /* decNumberInvert */
1303
1304/* ------------------------------------------------------------------ */
1305/* decNumberLn -- natural logarithm                                   */
1306/*                                                                    */
1307/*   This computes C = ln(A)                                          */
1308/*                                                                    */
1309/*   res is C, the result.  C may be A                                */
1310/*   rhs is A                                                         */
1311/*   set is the context; note that rounding mode has no effect        */
1312/*                                                                    */
1313/* C must have space for set->digits digits.                          */
1314/*                                                                    */
1315/* Notable cases:                                                     */
1316/*   A<0 -> Invalid                                                   */
1317/*   A=0 -> -Infinity (Exact)                                         */
1318/*   A=+Infinity -> +Infinity (Exact)                                 */
1319/*   A=1 exactly -> 0 (Exact)                                         */
1320/*                                                                    */
1321/* Mathematical function restrictions apply (see above); a NaN is     */
1322/* returned with Invalid_operation if a restriction is violated.      */
1323/*                                                                    */
1324/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1325/* almost always be correctly rounded, but may be up to 1 ulp in      */
1326/* error in rare cases.                                               */
1327/* ------------------------------------------------------------------ */
1328/* This is a wrapper for decLnOp which can handle the slightly wider  */
1329/* (+11) range needed by Ln, Log10, etc. (which may have to be able   */
1330/* to calculate at p+e+2).                                            */
1331/* ------------------------------------------------------------------ */
1332decNumber * decNumberLn(decNumber *res, const decNumber *rhs,
1333                        decContext *set) {
1334  uInt status=0;                   /* accumulator */
1335  #if DECSUBSET
1336  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
1337  #endif
1338
1339  #if DECCHECK
1340  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1341  #endif
1342
1343  /* Check restrictions; this is a math function; if not violated */
1344  /* then carry out the operation. */
1345  if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
1346    #if DECSUBSET
1347    if (!set->extended) {
1348      /* reduce operand and set lostDigits status, as needed */
1349      if (rhs->digits>set->digits) {
1350        allocrhs=decRoundOperand(rhs, set, &status);
1351        if (allocrhs==NULL) break;
1352        rhs=allocrhs;
1353        }
1354      /* special check in subset for rhs=0 */
1355      if (ISZERO(rhs)) {                /* +/- zeros -> error */
1356        status|=DEC_Invalid_operation;
1357        break;}
1358      } /* extended=0 */
1359    #endif
1360    decLnOp(res, rhs, set, &status);
1361    } while(0);                         /* end protected */
1362
1363  #if DECSUBSET
1364  if (allocrhs !=NULL) free(allocrhs);  /* drop any storage used */
1365  #endif
1366  /* apply significant status */
1367  if (status!=0) decStatus(res, status, set);
1368  #if DECCHECK
1369  decCheckInexact(res, set);
1370  #endif
1371  return res;
1372  } /* decNumberLn */
1373
1374/* ------------------------------------------------------------------ */
1375/* decNumberLogB - get adjusted exponent, by 754r rules               */
1376/*                                                                    */
1377/*   This computes C = adjustedexponent(A)                            */
1378/*                                                                    */
1379/*   res is C, the result.  C may be A                                */
1380/*   rhs is A                                                         */
1381/*   set is the context, used only for digits and status              */
1382/*                                                                    */
1383/* C must have space for 10 digits (A might have 10**9 digits and     */
1384/* an exponent of +999999999, or one digit and an exponent of         */
1385/* -1999999999).                                                      */
1386/*                                                                    */
1387/* This returns the adjusted exponent of A after (in theory) padding  */
1388/* with zeros on the right to set->digits digits while keeping the    */
1389/* same value.  The exponent is not limited by emin/emax.             */
1390/*                                                                    */
1391/* Notable cases:                                                     */
1392/*   A<0 -> Use |A|                                                   */
1393/*   A=0 -> -Infinity (Division by zero)                              */
1394/*   A=Infinite -> +Infinity (Exact)                                  */
1395/*   A=1 exactly -> 0 (Exact)                                         */
1396/*   NaNs are propagated as usual                                     */
1397/* ------------------------------------------------------------------ */
1398decNumber * decNumberLogB(decNumber *res, const decNumber *rhs,
1399                          decContext *set) {
1400  uInt status=0;                   /* accumulator */
1401
1402  #if DECCHECK
1403  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1404  #endif
1405
1406  /* NaNs as usual; Infinities return +Infinity; 0->oops */
1407  if (decNumberIsNaN(rhs)) decNaNs(res, rhs, NULL, set, &status);
1408   else if (decNumberIsInfinite(rhs)) decNumberCopyAbs(res, rhs);
1409   else if (decNumberIsZero(rhs)) {
1410    decNumberZero(res);                 /* prepare for Infinity */
1411    res->bits=DECNEG|DECINF;            /* -Infinity */
1412    status|=DEC_Division_by_zero;       /* as per 754r */
1413    }
1414   else { /* finite non-zero */
1415    Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */
1416    decNumberFromInt32(res, ae);        /* lay it out */
1417    }
1418
1419  if (status!=0) decStatus(res, status, set);
1420  return res;
1421  } /* decNumberLogB */
1422
1423/* ------------------------------------------------------------------ */
1424/* decNumberLog10 -- logarithm in base 10                             */
1425/*                                                                    */
1426/*   This computes C = log10(A)                                       */
1427/*                                                                    */
1428/*   res is C, the result.  C may be A                                */
1429/*   rhs is A                                                         */
1430/*   set is the context; note that rounding mode has no effect        */
1431/*                                                                    */
1432/* C must have space for set->digits digits.                          */
1433/*                                                                    */
1434/* Notable cases:                                                     */
1435/*   A<0 -> Invalid                                                   */
1436/*   A=0 -> -Infinity (Exact)                                         */
1437/*   A=+Infinity -> +Infinity (Exact)                                 */
1438/*   A=10**n (if n is an integer) -> n (Exact)                        */
1439/*                                                                    */
1440/* Mathematical function restrictions apply (see above); a NaN is     */
1441/* returned with Invalid_operation if a restriction is violated.      */
1442/*                                                                    */
1443/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1444/* almost always be correctly rounded, but may be up to 1 ulp in      */
1445/* error in rare cases.                                               */
1446/* ------------------------------------------------------------------ */
1447/* This calculates ln(A)/ln(10) using appropriate precision.  For     */
1448/* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the      */
1449/* requested digits and t is the number of digits in the exponent     */
1450/* (maximum 6).  For ln(10) it is p + 3; this is often handled by the */
1451/* fastpath in decLnOp.  The final division is done to the requested  */
1452/* precision.                                                         */
1453/* ------------------------------------------------------------------ */
1454decNumber * decNumberLog10(decNumber *res, const decNumber *rhs,
1455                          decContext *set) {
1456  uInt status=0, ignore=0;         /* status accumulators */
1457  uInt needbytes;                  /* for space calculations */
1458  Int p;                           /* working precision */
1459  Int t;                           /* digits in exponent of A */
1460
1461  /* buffers for a and b working decimals */
1462  /* (adjustment calculator, same size) */
1463  decNumber bufa[D2N(DECBUFFER+2)];
1464  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
1465  decNumber *a=bufa;               /* temporary a */
1466  decNumber bufb[D2N(DECBUFFER+2)];
1467  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
1468  decNumber *b=bufb;               /* temporary b */
1469  decNumber bufw[D2N(10)];         /* working 2-10 digit number */
1470  decNumber *w=bufw;               /* .. */
1471  #if DECSUBSET
1472  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
1473  #endif
1474
1475  decContext aset;                 /* working context */
1476
1477  #if DECCHECK
1478  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1479  #endif
1480
1481  /* Check restrictions; this is a math function; if not violated */
1482  /* then carry out the operation. */
1483  if (!decCheckMath(rhs, set, &status)) do { /* protect malloc */
1484    #if DECSUBSET
1485    if (!set->extended) {
1486      /* reduce operand and set lostDigits status, as needed */
1487      if (rhs->digits>set->digits) {
1488        allocrhs=decRoundOperand(rhs, set, &status);
1489        if (allocrhs==NULL) break;
1490        rhs=allocrhs;
1491        }
1492      /* special check in subset for rhs=0 */
1493      if (ISZERO(rhs)) {                /* +/- zeros -> error */
1494        status|=DEC_Invalid_operation;
1495        break;}
1496      } /* extended=0 */
1497    #endif
1498
1499    decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */
1500
1501    /* handle exact powers of 10; only check if +ve finite */
1502    if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) {
1503      Int residue=0;               /* (no residue) */
1504      uInt copystat=0;             /* clean status */
1505
1506      /* round to a single digit... */
1507      aset.digits=1;
1508      decCopyFit(w, rhs, &aset, &residue, &copystat); /* copy & shorten */
1509      /* if exact and the digit is 1, rhs is a power of 10 */
1510      if (!(copystat&DEC_Inexact) && w->lsu[0]==1) {
1511        /* the exponent, conveniently, is the power of 10; making */
1512        /* this the result needs a little care as it might not fit, */
1513        /* so first convert it into the working number, and then move */
1514        /* to res */
1515        decNumberFromInt32(w, w->exponent);
1516        residue=0;
1517        decCopyFit(res, w, set, &residue, &status); /* copy & round */
1518        decFinish(res, set, &residue, &status);     /* cleanup/set flags */
1519        break;
1520        } /* not a power of 10 */
1521      } /* not a candidate for exact */
1522
1523    /* simplify the information-content calculation to use 'total */
1524    /* number of digits in a, including exponent' as compared to the */
1525    /* requested digits, as increasing this will only rarely cost an */
1526    /* iteration in ln(a) anyway */
1527    t=6;                                /* it can never be >6 */
1528
1529    /* allocate space when needed... */
1530    p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3;
1531    needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
1532    if (needbytes>sizeof(bufa)) {       /* need malloc space */
1533      allocbufa=(decNumber *)malloc(needbytes);
1534      if (allocbufa==NULL) {            /* hopeless -- abandon */
1535        status|=DEC_Insufficient_storage;
1536        break;}
1537      a=allocbufa;                      /* use the allocated space */
1538      }
1539    aset.digits=p;                      /* as calculated */
1540    aset.emax=DEC_MAX_MATH;             /* usual bounds */
1541    aset.emin=-DEC_MAX_MATH;            /* .. */
1542    aset.clamp=0;                       /* and no concrete format */
1543    decLnOp(a, rhs, &aset, &status);    /* a=ln(rhs) */
1544
1545    /* skip the division if the result so far is infinite, NaN, or */
1546    /* zero, or there was an error; note NaN from sNaN needs copy */
1547    if (status&DEC_NaNs && !(status&DEC_sNaN)) break;
1548    if (a->bits&DECSPECIAL || ISZERO(a)) {
1549      decNumberCopy(res, a);            /* [will fit] */
1550      break;}
1551
1552    /* for ln(10) an extra 3 digits of precision are needed */
1553    p=set->digits+3;
1554    needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
1555    if (needbytes>sizeof(bufb)) {       /* need malloc space */
1556      allocbufb=(decNumber *)malloc(needbytes);
1557      if (allocbufb==NULL) {            /* hopeless -- abandon */
1558        status|=DEC_Insufficient_storage;
1559        break;}
1560      b=allocbufb;                      /* use the allocated space */
1561      }
1562    decNumberZero(w);                   /* set up 10... */
1563    #if DECDPUN==1
1564    w->lsu[1]=1; w->lsu[0]=0;           /* .. */
1565    #else
1566    w->lsu[0]=10;                       /* .. */
1567    #endif
1568    w->digits=2;                        /* .. */
1569
1570    aset.digits=p;
1571    decLnOp(b, w, &aset, &ignore);      /* b=ln(10) */
1572
1573    aset.digits=set->digits;            /* for final divide */
1574    decDivideOp(res, a, b, &aset, DIVIDE, &status); /* into result */
1575    } while(0);                         /* [for break] */
1576
1577  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1578  if (allocbufb!=NULL) free(allocbufb); /* .. */
1579  #if DECSUBSET
1580  if (allocrhs !=NULL) free(allocrhs);  /* .. */
1581  #endif
1582  /* apply significant status */
1583  if (status!=0) decStatus(res, status, set);
1584  #if DECCHECK
1585  decCheckInexact(res, set);
1586  #endif
1587  return res;
1588  } /* decNumberLog10 */
1589
1590/* ------------------------------------------------------------------ */
1591/* decNumberMax -- compare two Numbers and return the maximum         */
1592/*                                                                    */
1593/*   This computes C = A ? B, returning the maximum by 754R rules     */
1594/*                                                                    */
1595/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
1596/*   lhs is A                                                         */
1597/*   rhs is B                                                         */
1598/*   set is the context                                               */
1599/*                                                                    */
1600/* C must have space for set->digits digits.                          */
1601/* ------------------------------------------------------------------ */
1602decNumber * decNumberMax(decNumber *res, const decNumber *lhs,
1603                         const decNumber *rhs, decContext *set) {
1604  uInt status=0;                        /* accumulator */
1605  decCompareOp(res, lhs, rhs, set, COMPMAX, &status);
1606  if (status!=0) decStatus(res, status, set);
1607  #if DECCHECK
1608  decCheckInexact(res, set);
1609  #endif
1610  return res;
1611  } /* decNumberMax */
1612
1613/* ------------------------------------------------------------------ */
1614/* decNumberMaxMag -- compare and return the maximum by magnitude     */
1615/*                                                                    */
1616/*   This computes C = A ? B, returning the maximum by 754R rules     */
1617/*                                                                    */
1618/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
1619/*   lhs is A                                                         */
1620/*   rhs is B                                                         */
1621/*   set is the context                                               */
1622/*                                                                    */
1623/* C must have space for set->digits digits.                          */
1624/* ------------------------------------------------------------------ */
1625decNumber * decNumberMaxMag(decNumber *res, const decNumber *lhs,
1626                         const decNumber *rhs, decContext *set) {
1627  uInt status=0;                        /* accumulator */
1628  decCompareOp(res, lhs, rhs, set, COMPMAXMAG, &status);
1629  if (status!=0) decStatus(res, status, set);
1630  #if DECCHECK
1631  decCheckInexact(res, set);
1632  #endif
1633  return res;
1634  } /* decNumberMaxMag */
1635
1636/* ------------------------------------------------------------------ */
1637/* decNumberMin -- compare two Numbers and return the minimum         */
1638/*                                                                    */
1639/*   This computes C = A ? B, returning the minimum by 754R rules     */
1640/*                                                                    */
1641/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
1642/*   lhs is A                                                         */
1643/*   rhs is B                                                         */
1644/*   set is the context                                               */
1645/*                                                                    */
1646/* C must have space for set->digits digits.                          */
1647/* ------------------------------------------------------------------ */
1648decNumber * decNumberMin(decNumber *res, const decNumber *lhs,
1649                         const decNumber *rhs, decContext *set) {
1650  uInt status=0;                        /* accumulator */
1651  decCompareOp(res, lhs, rhs, set, COMPMIN, &status);
1652  if (status!=0) decStatus(res, status, set);
1653  #if DECCHECK
1654  decCheckInexact(res, set);
1655  #endif
1656  return res;
1657  } /* decNumberMin */
1658
1659/* ------------------------------------------------------------------ */
1660/* decNumberMinMag -- compare and return the minimum by magnitude     */
1661/*                                                                    */
1662/*   This computes C = A ? B, returning the minimum by 754R rules     */
1663/*                                                                    */
1664/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
1665/*   lhs is A                                                         */
1666/*   rhs is B                                                         */
1667/*   set is the context                                               */
1668/*                                                                    */
1669/* C must have space for set->digits digits.                          */
1670/* ------------------------------------------------------------------ */
1671decNumber * decNumberMinMag(decNumber *res, const decNumber *lhs,
1672                         const decNumber *rhs, decContext *set) {
1673  uInt status=0;                        /* accumulator */
1674  decCompareOp(res, lhs, rhs, set, COMPMINMAG, &status);
1675  if (status!=0) decStatus(res, status, set);
1676  #if DECCHECK
1677  decCheckInexact(res, set);
1678  #endif
1679  return res;
1680  } /* decNumberMinMag */
1681
1682/* ------------------------------------------------------------------ */
1683/* decNumberMinus -- prefix minus operator                            */
1684/*                                                                    */
1685/*   This computes C = 0 - A                                          */
1686/*                                                                    */
1687/*   res is C, the result.  C may be A                                */
1688/*   rhs is A                                                         */
1689/*   set is the context                                               */
1690/*                                                                    */
1691/* See also decNumberCopyNegate for a quiet bitwise version of this.  */
1692/* C must have space for set->digits digits.                          */
1693/* ------------------------------------------------------------------ */
1694/* Simply use AddOp for the subtract, which will do the necessary.    */
1695/* ------------------------------------------------------------------ */
1696decNumber * decNumberMinus(decNumber *res, const decNumber *rhs,
1697                           decContext *set) {
1698  decNumber dzero;
1699  uInt status=0;                        /* accumulator */
1700
1701  #if DECCHECK
1702  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1703  #endif
1704
1705  decNumberZero(&dzero);                /* make 0 */
1706  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
1707  decAddOp(res, &dzero, rhs, set, DECNEG, &status);
1708  if (status!=0) decStatus(res, status, set);
1709  #if DECCHECK
1710  decCheckInexact(res, set);
1711  #endif
1712  return res;
1713  } /* decNumberMinus */
1714
1715/* ------------------------------------------------------------------ */
1716/* decNumberNextMinus -- next towards -Infinity                       */
1717/*                                                                    */
1718/*   This computes C = A - infinitesimal, rounded towards -Infinity   */
1719/*                                                                    */
1720/*   res is C, the result.  C may be A                                */
1721/*   rhs is A                                                         */
1722/*   set is the context                                               */
1723/*                                                                    */
1724/* This is a generalization of 754r NextDown.                         */
1725/* ------------------------------------------------------------------ */
1726decNumber * decNumberNextMinus(decNumber *res, const decNumber *rhs,
1727                               decContext *set) {
1728  decNumber dtiny;                           /* constant */
1729  decContext workset=*set;                   /* work */
1730  uInt status=0;                             /* accumulator */
1731  #if DECCHECK
1732  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1733  #endif
1734
1735  /* +Infinity is the special case */
1736  if ((rhs->bits&(DECINF|DECNEG))==DECINF) {
1737    decSetMaxValue(res, set);                /* is +ve */
1738    /* there is no status to set */
1739    return res;
1740    }
1741  decNumberZero(&dtiny);                     /* start with 0 */
1742  dtiny.lsu[0]=1;                            /* make number that is .. */
1743  dtiny.exponent=DEC_MIN_EMIN-1;             /* .. smaller than tiniest */
1744  workset.round=DEC_ROUND_FLOOR;
1745  decAddOp(res, rhs, &dtiny, &workset, DECNEG, &status);
1746  status&=DEC_Invalid_operation|DEC_sNaN;    /* only sNaN Invalid please */
1747  if (status!=0) decStatus(res, status, set);
1748  return res;
1749  } /* decNumberNextMinus */
1750
1751/* ------------------------------------------------------------------ */
1752/* decNumberNextPlus -- next towards +Infinity                        */
1753/*                                                                    */
1754/*   This computes C = A + infinitesimal, rounded towards +Infinity   */
1755/*                                                                    */
1756/*   res is C, the result.  C may be A                                */
1757/*   rhs is A                                                         */
1758/*   set is the context                                               */
1759/*                                                                    */
1760/* This is a generalization of 754r NextUp.                           */
1761/* ------------------------------------------------------------------ */
1762decNumber * decNumberNextPlus(decNumber *res, const decNumber *rhs,
1763                              decContext *set) {
1764  decNumber dtiny;                           /* constant */
1765  decContext workset=*set;                   /* work */
1766  uInt status=0;                             /* accumulator */
1767  #if DECCHECK
1768  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1769  #endif
1770
1771  /* -Infinity is the special case */
1772  if ((rhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
1773    decSetMaxValue(res, set);
1774    res->bits=DECNEG;                        /* negative */
1775    /* there is no status to set */
1776    return res;
1777    }
1778  decNumberZero(&dtiny);                     /* start with 0 */
1779  dtiny.lsu[0]=1;                            /* make number that is .. */
1780  dtiny.exponent=DEC_MIN_EMIN-1;             /* .. smaller than tiniest */
1781  workset.round=DEC_ROUND_CEILING;
1782  decAddOp(res, rhs, &dtiny, &workset, 0, &status);
1783  status&=DEC_Invalid_operation|DEC_sNaN;    /* only sNaN Invalid please */
1784  if (status!=0) decStatus(res, status, set);
1785  return res;
1786  } /* decNumberNextPlus */
1787
1788/* ------------------------------------------------------------------ */
1789/* decNumberNextToward -- next towards rhs                            */
1790/*                                                                    */
1791/*   This computes C = A +/- infinitesimal, rounded towards           */
1792/*   +/-Infinity in the direction of B, as per 754r nextafter rules   */
1793/*                                                                    */
1794/*   res is C, the result.  C may be A or B.                          */
1795/*   lhs is A                                                         */
1796/*   rhs is B                                                         */
1797/*   set is the context                                               */
1798/*                                                                    */
1799/* This is a generalization of 754r NextAfter.                        */
1800/* ------------------------------------------------------------------ */
1801decNumber * decNumberNextToward(decNumber *res, const decNumber *lhs,
1802                                const decNumber *rhs, decContext *set) {
1803  decNumber dtiny;                           /* constant */
1804  decContext workset=*set;                   /* work */
1805  Int result;                                /* .. */
1806  uInt status=0;                             /* accumulator */
1807  #if DECCHECK
1808  if (decCheckOperands(res, lhs, rhs, set)) return res;
1809  #endif
1810
1811  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) {
1812    decNaNs(res, lhs, rhs, set, &status);
1813    }
1814   else { /* Is numeric, so no chance of sNaN Invalid, etc. */
1815    result=decCompare(lhs, rhs, 0);     /* sign matters */
1816    if (result==BADINT) status|=DEC_Insufficient_storage; /* rare */
1817     else { /* valid compare */
1818      if (result==0) decNumberCopySign(res, lhs, rhs); /* easy */
1819       else { /* differ: need NextPlus or NextMinus */
1820        uByte sub;                      /* add or subtract */
1821        if (result<0) {                 /* lhs<rhs, do nextplus */
1822          /* -Infinity is the special case */
1823          if ((lhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
1824            decSetMaxValue(res, set);
1825            res->bits=DECNEG;           /* negative */
1826            return res;                 /* there is no status to set */
1827            }
1828          workset.round=DEC_ROUND_CEILING;
1829          sub=0;                        /* add, please */
1830          } /* plus */
1831         else {                         /* lhs>rhs, do nextminus */
1832          /* +Infinity is the special case */
1833          if ((lhs->bits&(DECINF|DECNEG))==DECINF) {
1834            decSetMaxValue(res, set);
1835            return res;                 /* there is no status to set */
1836            }
1837          workset.round=DEC_ROUND_FLOOR;
1838          sub=DECNEG;                   /* subtract, please */
1839          } /* minus */
1840        decNumberZero(&dtiny);          /* start with 0 */
1841        dtiny.lsu[0]=1;                 /* make number that is .. */
1842        dtiny.exponent=DEC_MIN_EMIN-1;  /* .. smaller than tiniest */
1843        decAddOp(res, lhs, &dtiny, &workset, sub, &status); /* + or - */
1844        /* turn off exceptions if the result is a normal number */
1845        /* (including Nmin), otherwise let all status through */
1846        if (decNumberIsNormal(res, set)) status=0;
1847        } /* unequal */
1848      } /* compare OK */
1849    } /* numeric */
1850  if (status!=0) decStatus(res, status, set);
1851  return res;
1852  } /* decNumberNextToward */
1853
1854/* ------------------------------------------------------------------ */
1855/* decNumberOr -- OR two Numbers, digitwise                           */
1856/*                                                                    */
1857/*   This computes C = A | B                                          */
1858/*                                                                    */
1859/*   res is C, the result.  C may be A and/or B (e.g., X=X|X)         */
1860/*   lhs is A                                                         */
1861/*   rhs is B                                                         */
1862/*   set is the context (used for result length and error report)     */
1863/*                                                                    */
1864/* C must have space for set->digits digits.                          */
1865/*                                                                    */
1866/* Logical function restrictions apply (see above); a NaN is          */
1867/* returned with Invalid_operation if a restriction is violated.      */
1868/* ------------------------------------------------------------------ */
1869decNumber * decNumberOr(decNumber *res, const decNumber *lhs,
1870                        const decNumber *rhs, decContext *set) {
1871  const Unit *ua, *ub;                  /* -> operands */
1872  const Unit *msua, *msub;              /* -> operand msus */
1873  Unit  *uc, *msuc;                     /* -> result and its msu */
1874  Int   msudigs;                        /* digits in res msu */
1875  #if DECCHECK
1876  if (decCheckOperands(res, lhs, rhs, set)) return res;
1877  #endif
1878
1879  if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
1880   || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
1881    decStatus(res, DEC_Invalid_operation, set);
1882    return res;
1883    }
1884  /* operands are valid */
1885  ua=lhs->lsu;                          /* bottom-up */
1886  ub=rhs->lsu;                          /* .. */
1887  uc=res->lsu;                          /* .. */
1888  msua=ua+D2U(lhs->digits)-1;           /* -> msu of lhs */
1889  msub=ub+D2U(rhs->digits)-1;           /* -> msu of rhs */
1890  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
1891  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
1892  for (; uc<=msuc; ua++, ub++, uc++) {  /* Unit loop */
1893    Unit a, b;                          /* extract units */
1894    if (ua>msua) a=0;
1895     else a=*ua;
1896    if (ub>msub) b=0;
1897     else b=*ub;
1898    *uc=0;                              /* can now write back */
1899    if (a|b) {                          /* maybe 1 bits to examine */
1900      Int i, j;
1901      /* This loop could be unrolled and/or use BIN2BCD tables */
1902      for (i=0; i<DECDPUN; i++) {
1903        if ((a|b)&1) *uc=*uc+(Unit)powers[i];     /* effect OR */
1904        j=a%10;
1905        a=a/10;
1906        j|=b%10;
1907        b=b/10;
1908        if (j>1) {
1909          decStatus(res, DEC_Invalid_operation, set);
1910          return res;
1911          }
1912        if (uc==msuc && i==msudigs-1) break;      /* just did final digit */
1913        } /* each digit */
1914      } /* non-zero */
1915    } /* each unit */
1916  /* [here uc-1 is the msu of the result] */
1917  res->digits=decGetDigits(res->lsu, uc-res->lsu);
1918  res->exponent=0;                      /* integer */
1919  res->bits=0;                          /* sign=0 */
1920  return res;  /* [no status to set] */
1921  } /* decNumberOr */
1922
1923/* ------------------------------------------------------------------ */
1924/* decNumberPlus -- prefix plus operator                              */
1925/*                                                                    */
1926/*   This computes C = 0 + A                                          */
1927/*                                                                    */
1928/*   res is C, the result.  C may be A                                */
1929/*   rhs is A                                                         */
1930/*   set is the context                                               */
1931/*                                                                    */
1932/* See also decNumberCopy for a quiet bitwise version of this.        */
1933/* C must have space for set->digits digits.                          */
1934/* ------------------------------------------------------------------ */
1935/* This simply uses AddOp; Add will take fast path after preparing A. */
1936/* Performance is a concern here, as this routine is often used to    */
1937/* check operands and apply rounding and overflow/underflow testing.  */
1938/* ------------------------------------------------------------------ */
1939decNumber * decNumberPlus(decNumber *res, const decNumber *rhs,
1940                          decContext *set) {
1941  decNumber dzero;
1942  uInt status=0;                        /* accumulator */
1943  #if DECCHECK
1944  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1945  #endif
1946
1947  decNumberZero(&dzero);                /* make 0 */
1948  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
1949  decAddOp(res, &dzero, rhs, set, 0, &status);
1950  if (status!=0) decStatus(res, status, set);
1951  #if DECCHECK
1952  decCheckInexact(res, set);
1953  #endif
1954  return res;
1955  } /* decNumberPlus */
1956
1957/* ------------------------------------------------------------------ */
1958/* decNumberMultiply -- multiply two Numbers                          */
1959/*                                                                    */
1960/*   This computes C = A x B                                          */
1961/*                                                                    */
1962/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
1963/*   lhs is A                                                         */
1964/*   rhs is B                                                         */
1965/*   set is the context                                               */
1966/*                                                                    */
1967/* C must have space for set->digits digits.                          */
1968/* ------------------------------------------------------------------ */
1969decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs,
1970                              const decNumber *rhs, decContext *set) {
1971  uInt status=0;                   /* accumulator */
1972  decMultiplyOp(res, lhs, rhs, set, &status);
1973  if (status!=0) decStatus(res, status, set);
1974  #if DECCHECK
1975  decCheckInexact(res, set);
1976  #endif
1977  return res;
1978  } /* decNumberMultiply */
1979
1980/* ------------------------------------------------------------------ */
1981/* decNumberPower -- raise a number to a power                        */
1982/*                                                                    */
1983/*   This computes C = A ** B                                         */
1984/*                                                                    */
1985/*   res is C, the result.  C may be A and/or B (e.g., X=X**X)        */
1986/*   lhs is A                                                         */
1987/*   rhs is B                                                         */
1988/*   set is the context                                               */
1989/*                                                                    */
1990/* C must have space for set->digits digits.                          */
1991/*                                                                    */
1992/* Mathematical function restrictions apply (see above); a NaN is     */
1993/* returned with Invalid_operation if a restriction is violated.      */
1994/*                                                                    */
1995/* However, if 1999999997<=B<=999999999 and B is an integer then the  */
1996/* restrictions on A and the context are relaxed to the usual bounds, */
1997/* for compatibility with the earlier (integer power only) version    */
1998/* of this function.                                                  */
1999/*                                                                    */
2000/* When B is an integer, the result may be exact, even if rounded.    */
2001/*                                                                    */
2002/* The final result is rounded according to the context; it will      */
2003/* almost always be correctly rounded, but may be up to 1 ulp in      */
2004/* error in rare cases.                                               */
2005/* ------------------------------------------------------------------ */
2006decNumber * decNumberPower(decNumber *res, const decNumber *lhs,
2007                           const decNumber *rhs, decContext *set) {
2008  #if DECSUBSET
2009  decNumber *alloclhs=NULL;        /* non-NULL if rounded lhs allocated */
2010  decNumber *allocrhs=NULL;        /* .., rhs */
2011  #endif
2012  decNumber *allocdac=NULL;        /* -> allocated acc buffer, iff used */
2013  decNumber *allocinv=NULL;        /* -> allocated 1/x buffer, iff used */
2014  Int   reqdigits=set->digits;     /* requested DIGITS */
2015  Int   n;                         /* rhs in binary */
2016  Flag  rhsint=0;                  /* 1 if rhs is an integer */
2017  Flag  useint=0;                  /* 1 if can use integer calculation */
2018  Flag  isoddint=0;                /* 1 if rhs is an integer and odd */
2019  Int   i;                         /* work */
2020  #if DECSUBSET
2021  Int   dropped;                   /* .. */
2022  #endif
2023  uInt  needbytes;                 /* buffer size needed */
2024  Flag  seenbit;                   /* seen a bit while powering */
2025  Int   residue=0;                 /* rounding residue */
2026  uInt  status=0;                  /* accumulators */
2027  uByte bits=0;                    /* result sign if errors */
2028  decContext aset;                 /* working context */
2029  decNumber dnOne;                 /* work value 1... */
2030  /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
2031  decNumber dacbuff[D2N(DECBUFFER+9)];
2032  decNumber *dac=dacbuff;          /* -> result accumulator */
2033  /* same again for possible 1/lhs calculation */
2034  decNumber invbuff[D2N(DECBUFFER+9)];
2035
2036  #if DECCHECK
2037  if (decCheckOperands(res, lhs, rhs, set)) return res;
2038  #endif
2039
2040  do {                             /* protect allocated storage */
2041    #if DECSUBSET
2042    if (!set->extended) { /* reduce operands and set status, as needed */
2043      if (lhs->digits>reqdigits) {
2044        alloclhs=decRoundOperand(lhs, set, &status);
2045        if (alloclhs==NULL) break;
2046        lhs=alloclhs;
2047        }
2048      if (rhs->digits>reqdigits) {
2049        allocrhs=decRoundOperand(rhs, set, &status);
2050        if (allocrhs==NULL) break;
2051        rhs=allocrhs;
2052        }
2053      }
2054    #endif
2055    /* [following code does not require input rounding] */
2056
2057    /* handle NaNs and rhs Infinity (lhs infinity is harder) */
2058    if (SPECIALARGS) {
2059      if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { /* NaNs */
2060        decNaNs(res, lhs, rhs, set, &status);
2061        break;}
2062      if (decNumberIsInfinite(rhs)) {   /* rhs Infinity */
2063        Flag rhsneg=rhs->bits&DECNEG;   /* save rhs sign */
2064        if (decNumberIsNegative(lhs)    /* lhs<0 */
2065         && !decNumberIsZero(lhs))      /* .. */
2066          status|=DEC_Invalid_operation;
2067         else {                         /* lhs >=0 */
2068          decNumberZero(&dnOne);        /* set up 1 */
2069          dnOne.lsu[0]=1;
2070          decNumberCompare(dac, lhs, &dnOne, set); /* lhs ? 1 */
2071          decNumberZero(res);           /* prepare for 0/1/Infinity */
2072          if (decNumberIsNegative(dac)) {    /* lhs<1 */
2073            if (rhsneg) res->bits|=DECINF;   /* +Infinity [else is +0] */
2074            }
2075           else if (dac->lsu[0]==0) {        /* lhs=1 */
2076            /* 1**Infinity is inexact, so return fully-padded 1.0000 */
2077            Int shift=set->digits-1;
2078            *res->lsu=1;                     /* was 0, make int 1 */
2079            res->digits=decShiftToMost(res->lsu, 1, shift);
2080            res->exponent=-shift;            /* make 1.0000... */
2081            status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */
2082            }
2083           else {                            /* lhs>1 */
2084            if (!rhsneg) res->bits|=DECINF;  /* +Infinity [else is +0] */
2085            }
2086          } /* lhs>=0 */
2087        break;}
2088      /* [lhs infinity drops through] */
2089      } /* specials */
2090
2091    /* Original rhs may be an integer that fits and is in range */
2092    n=decGetInt(rhs);
2093    if (n!=BADINT) {                    /* it is an integer */
2094      rhsint=1;                         /* record the fact for 1**n */
2095      isoddint=(Flag)n&1;               /* [works even if big] */
2096      if (n!=BIGEVEN && n!=BIGODD)      /* can use integer path? */
2097        useint=1;                       /* looks good */
2098      }
2099
2100    if (decNumberIsNegative(lhs)        /* -x .. */
2101      && isoddint) bits=DECNEG;         /* .. to an odd power */
2102
2103    /* handle LHS infinity */
2104    if (decNumberIsInfinite(lhs)) {     /* [NaNs already handled] */
2105      uByte rbits=rhs->bits;            /* save */
2106      decNumberZero(res);               /* prepare */
2107      if (n==0) *res->lsu=1;            /* [-]Inf**0 => 1 */
2108       else {
2109        /* -Inf**nonint -> error */
2110        if (!rhsint && decNumberIsNegative(lhs)) {
2111          status|=DEC_Invalid_operation;     /* -Inf**nonint is error */
2112          break;}
2113        if (!(rbits & DECNEG)) bits|=DECINF; /* was not a **-n */
2114        /* [otherwise will be 0 or -0] */
2115        res->bits=bits;
2116        }
2117      break;}
2118
2119    /* similarly handle LHS zero */
2120    if (decNumberIsZero(lhs)) {
2121      if (n==0) {                            /* 0**0 => Error */
2122        #if DECSUBSET
2123        if (!set->extended) {                /* [unless subset] */
2124          decNumberZero(res);
2125          *res->lsu=1;                       /* return 1 */
2126          break;}
2127        #endif
2128        status|=DEC_Invalid_operation;
2129        }
2130       else {                                /* 0**x */
2131        uByte rbits=rhs->bits;               /* save */
2132        if (rbits & DECNEG) {                /* was a 0**(-n) */
2133          #if DECSUBSET
2134          if (!set->extended) {              /* [bad if subset] */
2135            status|=DEC_Invalid_operation;
2136            break;}
2137          #endif
2138          bits|=DECINF;
2139          }
2140        decNumberZero(res);                  /* prepare */
2141        /* [otherwise will be 0 or -0] */
2142        res->bits=bits;
2143        }
2144      break;}
2145
2146    /* here both lhs and rhs are finite; rhs==0 is handled in the */
2147    /* integer path.  Next handle the non-integer cases */
2148    if (!useint) {                      /* non-integral rhs */
2149      /* any -ve lhs is bad, as is either operand or context out of */
2150      /* bounds */
2151      if (decNumberIsNegative(lhs)) {
2152        status|=DEC_Invalid_operation;
2153        break;}
2154      if (decCheckMath(lhs, set, &status)
2155       || decCheckMath(rhs, set, &status)) break; /* variable status */
2156
2157      decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */
2158      aset.emax=DEC_MAX_MATH;           /* usual bounds */
2159      aset.emin=-DEC_MAX_MATH;          /* .. */
2160      aset.clamp=0;                     /* and no concrete format */
2161
2162      /* calculate the result using exp(ln(lhs)*rhs), which can */
2163      /* all be done into the accumulator, dac.  The precision needed */
2164      /* is enough to contain the full information in the lhs (which */
2165      /* is the total digits, including exponent), or the requested */
2166      /* precision, if larger, + 4; 6 is used for the exponent */
2167      /* maximum length, and this is also used when it is shorter */
2168      /* than the requested digits as it greatly reduces the >0.5 ulp */
2169      /* cases at little cost (because Ln doubles digits each */
2170      /* iteration so a few extra digits rarely causes an extra */
2171      /* iteration) */
2172      aset.digits=MAXI(lhs->digits, set->digits)+6+4;
2173      } /* non-integer rhs */
2174
2175     else { /* rhs is in-range integer */
2176      if (n==0) {                       /* x**0 = 1 */
2177        /* (0**0 was handled above) */
2178        decNumberZero(res);             /* result=1 */
2179        *res->lsu=1;                    /* .. */
2180        break;}
2181      /* rhs is a non-zero integer */
2182      if (n<0) n=-n;                    /* use abs(n) */
2183
2184      aset=*set;                        /* clone the context */
2185      aset.round=DEC_ROUND_HALF_EVEN;   /* internally use balanced */
2186      /* calculate the working DIGITS */
2187      aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2;
2188      #if DECSUBSET
2189      if (!set->extended) aset.digits--;     /* use classic precision */
2190      #endif
2191      /* it's an error if this is more than can be handled */
2192      if (aset.digits>DECNUMMAXP) {status|=DEC_Invalid_operation; break;}
2193      } /* integer path */
2194
2195    /* aset.digits is the count of digits for the accumulator needed */
2196    /* if accumulator is too long for local storage, then allocate */
2197    needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit);
2198    /* [needbytes also used below if 1/lhs needed] */
2199    if (needbytes>sizeof(dacbuff)) {
2200      allocdac=(decNumber *)malloc(needbytes);
2201      if (allocdac==NULL) {   /* hopeless -- abandon */
2202        status|=DEC_Insufficient_storage;
2203        break;}
2204      dac=allocdac;           /* use the allocated space */
2205      }
2206    /* here, aset is set up and accumulator is ready for use */
2207
2208    if (!useint) {                           /* non-integral rhs */
2209      /* x ** y; special-case x=1 here as it will otherwise always */
2210      /* reduce to integer 1; decLnOp has a fastpath which detects */
2211      /* the case of x=1 */
2212      decLnOp(dac, lhs, &aset, &status);     /* dac=ln(lhs) */
2213      /* [no error possible, as lhs 0 already handled] */
2214      if (ISZERO(dac)) {                     /* x==1, 1.0, etc. */
2215        /* need to return fully-padded 1.0000 etc., but rhsint->1 */
2216        *dac->lsu=1;                         /* was 0, make int 1 */
2217        if (!rhsint) {                       /* add padding */
2218          Int shift=set->digits-1;
2219          dac->digits=decShiftToMost(dac->lsu, 1, shift);
2220          dac->exponent=-shift;              /* make 1.0000... */
2221          status|=DEC_Inexact|DEC_Rounded;   /* deemed inexact */
2222          }
2223        }
2224       else {
2225        decMultiplyOp(dac, dac, rhs, &aset, &status);  /* dac=dac*rhs */
2226        decExpOp(dac, dac, &aset, &status);            /* dac=exp(dac) */
2227        }
2228      /* and drop through for final rounding */
2229      } /* non-integer rhs */
2230
2231     else {                             /* carry on with integer */
2232      decNumberZero(dac);               /* acc=1 */
2233      *dac->lsu=1;                      /* .. */
2234
2235      /* if a negative power the constant 1 is needed, and if not subset */
2236      /* invert the lhs now rather than inverting the result later */
2237      if (decNumberIsNegative(rhs)) {   /* was a **-n [hence digits>0] */
2238        decNumber *inv=invbuff;         /* assume use fixed buffer */
2239        decNumberCopy(&dnOne, dac);     /* dnOne=1;  [needed now or later] */
2240        #if DECSUBSET
2241        if (set->extended) {            /* need to calculate 1/lhs */
2242        #endif
2243          /* divide lhs into 1, putting result in dac [dac=1/dac] */
2244          decDivideOp(dac, &dnOne, lhs, &aset, DIVIDE, &status);
2245          /* now locate or allocate space for the inverted lhs */
2246          if (needbytes>sizeof(invbuff)) {
2247            allocinv=(decNumber *)malloc(needbytes);
2248            if (allocinv==NULL) {       /* hopeless -- abandon */
2249              status|=DEC_Insufficient_storage;
2250              break;}
2251            inv=allocinv;               /* use the allocated space */
2252            }
2253          /* [inv now points to big-enough buffer or allocated storage] */
2254          decNumberCopy(inv, dac);      /* copy the 1/lhs */
2255          decNumberCopy(dac, &dnOne);   /* restore acc=1 */
2256          lhs=inv;                      /* .. and go forward with new lhs */
2257        #if DECSUBSET
2258          }
2259        #endif
2260        }
2261
2262      /* Raise-to-the-power loop... */
2263      seenbit=0;                   /* set once a 1-bit is encountered */
2264      for (i=1;;i++){              /* for each bit [top bit ignored] */
2265        /* abandon if had overflow or terminal underflow */
2266        if (status & (DEC_Overflow|DEC_Underflow)) { /* interesting? */
2267          if (status&DEC_Overflow || ISZERO(dac)) break;
2268          }
2269        /* [the following two lines revealed an optimizer bug in a C++ */
2270        /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
2271        n=n<<1;                    /* move next bit to testable position */
2272        if (n<0) {                 /* top bit is set */
2273          seenbit=1;               /* OK, significant bit seen */
2274          decMultiplyOp(dac, dac, lhs, &aset, &status); /* dac=dac*x */
2275          }
2276        if (i==31) break;          /* that was the last bit */
2277        if (!seenbit) continue;    /* no need to square 1 */
2278        decMultiplyOp(dac, dac, dac, &aset, &status); /* dac=dac*dac [square] */
2279        } /*i*/ /* 32 bits */
2280
2281      /* complete internal overflow or underflow processing */
2282      if (status & (DEC_Overflow|DEC_Underflow)) {
2283        #if DECSUBSET
2284        /* If subset, and power was negative, reverse the kind of -erflow */
2285        /* [1/x not yet done] */
2286        if (!set->extended && decNumberIsNegative(rhs)) {
2287          if (status & DEC_Overflow)
2288            status^=DEC_Overflow | DEC_Underflow | DEC_Subnormal;
2289           else { /* trickier -- Underflow may or may not be set */
2290            status&=~(DEC_Underflow | DEC_Subnormal); /* [one or both] */
2291            status|=DEC_Overflow;
2292            }
2293          }
2294        #endif
2295        dac->bits=(dac->bits & ~DECNEG) | bits; /* force correct sign */
2296        /* round subnormals [to set.digits rather than aset.digits] */
2297        /* or set overflow result similarly as required */
2298        decFinalize(dac, set, &residue, &status);
2299        decNumberCopy(res, dac);   /* copy to result (is now OK length) */
2300        break;
2301        }
2302
2303      #if DECSUBSET
2304      if (!set->extended &&                  /* subset math */
2305          decNumberIsNegative(rhs)) {        /* was a **-n [hence digits>0] */
2306        /* so divide result into 1 [dac=1/dac] */
2307        decDivideOp(dac, &dnOne, dac, &aset, DIVIDE, &status);
2308        }
2309      #endif
2310      } /* rhs integer path */
2311
2312    /* reduce result to the requested length and copy to result */
2313    decCopyFit(res, dac, set, &residue, &status);
2314    decFinish(res, set, &residue, &status);  /* final cleanup */
2315    #if DECSUBSET
2316    if (!set->extended) decTrim(res, set, 0, &dropped); /* trailing zeros */
2317    #endif
2318    } while(0);                         /* end protected */
2319
2320  if (allocdac!=NULL) free(allocdac);   /* drop any storage used */
2321  if (allocinv!=NULL) free(allocinv);   /* .. */
2322  #if DECSUBSET
2323  if (alloclhs!=NULL) free(alloclhs);   /* .. */
2324  if (allocrhs!=NULL) free(allocrhs);   /* .. */
2325  #endif
2326  if (status!=0) decStatus(res, status, set);
2327  #if DECCHECK
2328  decCheckInexact(res, set);
2329  #endif
2330  return res;
2331  } /* decNumberPower */
2332
2333/* ------------------------------------------------------------------ */
2334/* decNumberQuantize -- force exponent to requested value             */
2335/*                                                                    */
2336/*   This computes C = op(A, B), where op adjusts the coefficient     */
2337/*   of C (by rounding or shifting) such that the exponent (-scale)   */
2338/*   of C has exponent of B.  The numerical value of C will equal A,  */
2339/*   except for the effects of any rounding that occurred.            */
2340/*                                                                    */
2341/*   res is C, the result.  C may be A or B                           */
2342/*   lhs is A, the number to adjust                                   */
2343/*   rhs is B, the number with exponent to match                      */
2344/*   set is the context                                               */
2345/*                                                                    */
2346/* C must have space for set->digits digits.                          */
2347/*                                                                    */
2348/* Unless there is an error or the result is infinite, the exponent   */
2349/* after the operation is guaranteed to be equal to that of B.        */
2350/* ------------------------------------------------------------------ */
2351decNumber * decNumberQuantize(decNumber *res, const decNumber *lhs,
2352                              const decNumber *rhs, decContext *set) {
2353  uInt status=0;                        /* accumulator */
2354  decQuantizeOp(res, lhs, rhs, set, 1, &status);
2355  if (status!=0) decStatus(res, status, set);
2356  return res;
2357  } /* decNumberQuantize */
2358
2359/* ------------------------------------------------------------------ */
2360/* decNumberReduce -- remove trailing zeros                           */
2361/*                                                                    */
2362/*   This computes C = 0 + A, and normalizes the result               */
2363/*                                                                    */
2364/*   res is C, the result.  C may be A                                */
2365/*   rhs is A                                                         */
2366/*   set is the context                                               */
2367/*                                                                    */
2368/* C must have space for set->digits digits.                          */
2369/* ------------------------------------------------------------------ */
2370/* Previously known as Normalize */
2371decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs,
2372                               decContext *set) {
2373  return decNumberReduce(res, rhs, set);
2374  } /* decNumberNormalize */
2375
2376decNumber * decNumberReduce(decNumber *res, const decNumber *rhs,
2377                            decContext *set) {
2378  #if DECSUBSET
2379  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
2380  #endif
2381  uInt status=0;                   /* as usual */
2382  Int  residue=0;                  /* as usual */
2383  Int  dropped;                    /* work */
2384
2385  #if DECCHECK
2386  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
2387  #endif
2388
2389  do {                             /* protect allocated storage */
2390    #if DECSUBSET
2391    if (!set->extended) {
2392      /* reduce operand and set lostDigits status, as needed */
2393      if (rhs->digits>set->digits) {
2394        allocrhs=decRoundOperand(rhs, set, &status);
2395        if (allocrhs==NULL) break;
2396        rhs=allocrhs;
2397        }
2398      }
2399    #endif
2400    /* [following code does not require input rounding] */
2401
2402    /* Infinities copy through; NaNs need usual treatment */
2403    if (decNumberIsNaN(rhs)) {
2404      decNaNs(res, rhs, NULL, set, &status);
2405      break;
2406      }
2407
2408    /* reduce result to the requested length and copy to result */
2409    decCopyFit(res, rhs, set, &residue, &status); /* copy & round */
2410    decFinish(res, set, &residue, &status);       /* cleanup/set flags */
2411    decTrim(res, set, 1, &dropped);               /* normalize in place */
2412    } while(0);                              /* end protected */
2413
2414  #if DECSUBSET
2415  if (allocrhs !=NULL) free(allocrhs);       /* .. */
2416  #endif
2417  if (status!=0) decStatus(res, status, set);/* then report status */
2418  return res;
2419  } /* decNumberReduce */
2420
2421/* ------------------------------------------------------------------ */
2422/* decNumberRescale -- force exponent to requested value              */
2423/*                                                                    */
2424/*   This computes C = op(A, B), where op adjusts the coefficient     */
2425/*   of C (by rounding or shifting) such that the exponent (-scale)   */
2426/*   of C has the value B.  The numerical value of C will equal A,    */
2427/*   except for the effects of any rounding that occurred.            */
2428/*                                                                    */
2429/*   res is C, the result.  C may be A or B                           */
2430/*   lhs is A, the number to adjust                                   */
2431/*   rhs is B, the requested exponent                                 */
2432/*   set is the context                                               */
2433/*                                                                    */
2434/* C must have space for set->digits digits.                          */
2435/*                                                                    */
2436/* Unless there is an error or the result is infinite, the exponent   */
2437/* after the operation is guaranteed to be equal to B.                */
2438/* ------------------------------------------------------------------ */
2439decNumber * decNumberRescale(decNumber *res, const decNumber *lhs,
2440                             const decNumber *rhs, decContext *set) {
2441  uInt status=0;                        /* accumulator */
2442  decQuantizeOp(res, lhs, rhs, set, 0, &status);
2443  if (status!=0) decStatus(res, status, set);
2444  return res;
2445  } /* decNumberRescale */
2446
2447/* ------------------------------------------------------------------ */
2448/* decNumberRemainder -- divide and return remainder                  */
2449/*                                                                    */
2450/*   This computes C = A % B                                          */
2451/*                                                                    */
2452/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
2453/*   lhs is A                                                         */
2454/*   rhs is B                                                         */
2455/*   set is the context                                               */
2456/*                                                                    */
2457/* C must have space for set->digits digits.                          */
2458/* ------------------------------------------------------------------ */
2459decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs,
2460                               const decNumber *rhs, decContext *set) {
2461  uInt status=0;                        /* accumulator */
2462  decDivideOp(res, lhs, rhs, set, REMAINDER, &status);
2463  if (status!=0) decStatus(res, status, set);
2464  #if DECCHECK
2465  decCheckInexact(res, set);
2466  #endif
2467  return res;
2468  } /* decNumberRemainder */
2469
2470/* ------------------------------------------------------------------ */
2471/* decNumberRemainderNear -- divide and return remainder from nearest */
2472/*                                                                    */
2473/*   This computes C = A % B, where % is the IEEE remainder operator  */
2474/*                                                                    */
2475/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
2476/*   lhs is A                                                         */
2477/*   rhs is B                                                         */
2478/*   set is the context                                               */
2479/*                                                                    */
2480/* C must have space for set->digits digits.                          */
2481/* ------------------------------------------------------------------ */
2482decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs,
2483                                   const decNumber *rhs, decContext *set) {
2484  uInt status=0;                        /* accumulator */
2485  decDivideOp(res, lhs, rhs, set, REMNEAR, &status);
2486  if (status!=0) decStatus(res, status, set);
2487  #if DECCHECK
2488  decCheckInexact(res, set);
2489  #endif
2490  return res;
2491  } /* decNumberRemainderNear */
2492
2493/* ------------------------------------------------------------------ */
2494/* decNumberRotate -- rotate the coefficient of a Number left/right   */
2495/*                                                                    */
2496/*   This computes C = A rot B  (in base ten and rotating set->digits */
2497/*   digits).                                                         */
2498/*                                                                    */
2499/*   res is C, the result.  C may be A and/or B (e.g., X=XrotX)       */
2500/*   lhs is A                                                         */
2501/*   rhs is B, the number of digits to rotate (-ve to right)          */
2502/*   set is the context                                               */
2503/*                                                                    */
2504/* The digits of the coefficient of A are rotated to the left (if B   */
2505/* is positive) or to the right (if B is negative) without adjusting  */
2506/* the exponent or the sign of A.  If lhs->digits is less than        */
2507/* set->digits the coefficient is padded with zeros on the left       */
2508/* before the rotate.  Any leading zeros in the result are removed    */
2509/* as usual.                                                          */
2510/*                                                                    */
2511/* B must be an integer (q=0) and in the range -set->digits through   */
2512/* +set->digits.                                                      */
2513/* C must have space for set->digits digits.                          */
2514/* NaNs are propagated as usual.  Infinities are unaffected (but      */
2515/* B must be valid).  No status is set unless B is invalid or an      */
2516/* operand is an sNaN.                                                */
2517/* ------------------------------------------------------------------ */
2518decNumber * decNumberRotate(decNumber *res, const decNumber *lhs,
2519                           const decNumber *rhs, decContext *set) {
2520  uInt status=0;              /* accumulator */
2521  Int  rotate;                /* rhs as an Int */
2522
2523  #if DECCHECK
2524  if (decCheckOperands(res, lhs, rhs, set)) return res;
2525  #endif
2526
2527  /* NaNs propagate as normal */
2528  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2529    decNaNs(res, lhs, rhs, set, &status);
2530   /* rhs must be an integer */
2531   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2532    status=DEC_Invalid_operation;
2533   else { /* both numeric, rhs is an integer */
2534    rotate=decGetInt(rhs);                   /* [cannot fail] */
2535    if (rotate==BADINT                       /* something bad .. */
2536     || rotate==BIGODD || rotate==BIGEVEN    /* .. very big .. */
2537     || abs(rotate)>set->digits)             /* .. or out of range */
2538      status=DEC_Invalid_operation;
2539     else {                                  /* rhs is OK */
2540      decNumberCopy(res, lhs);
2541      /* convert -ve rotate to equivalent positive rotation */
2542      if (rotate<0) rotate=set->digits+rotate;
2543      if (rotate!=0 && rotate!=set->digits   /* zero or full rotation */
2544       && !decNumberIsInfinite(res)) {       /* lhs was infinite */
2545        /* left-rotate to do; 0 < rotate < set->digits */
2546        uInt units, shift;                   /* work */
2547        uInt msudigits;                      /* digits in result msu */
2548        Unit *msu=res->lsu+D2U(res->digits)-1;    /* current msu */
2549        Unit *msumax=res->lsu+D2U(set->digits)-1; /* rotation msu */
2550        for (msu++; msu<=msumax; msu++) *msu=0;   /* ensure high units=0 */
2551        res->digits=set->digits;                  /* now full-length */
2552        msudigits=MSUDIGITS(res->digits);         /* actual digits in msu */
2553
2554        /* rotation here is done in-place, in three steps */
2555        /* 1. shift all to least up to one unit to unit-align final */
2556        /*    lsd [any digits shifted out are rotated to the left, */
2557        /*    abutted to the original msd (which may require split)] */
2558        /* */
2559        /*    [if there are no whole units left to rotate, the */
2560        /*    rotation is now complete] */
2561        /* */
2562        /* 2. shift to least, from below the split point only, so that */
2563        /*    the final msd is in the right place in its Unit [any */
2564        /*    digits shifted out will fit exactly in the current msu, */
2565        /*    left aligned, no split required] */
2566        /* */
2567        /* 3. rotate all the units by reversing left part, right */
2568        /*    part, and then whole */
2569        /* */
2570        /* example: rotate right 8 digits (2 units + 2), DECDPUN=3. */
2571        /* */
2572        /*   start: 00a bcd efg hij klm npq */
2573        /* */
2574        /*      1a  000 0ab cde fgh|ijk lmn [pq saved] */
2575        /*      1b  00p qab cde fgh|ijk lmn */
2576        /* */
2577        /*      2a  00p qab cde fgh|00i jkl [mn saved] */
2578        /*      2b  mnp qab cde fgh|00i jkl */
2579        /* */
2580        /*      3a  fgh cde qab mnp|00i jkl */
2581        /*      3b  fgh cde qab mnp|jkl 00i */
2582        /*      3c  00i jkl mnp qab cde fgh */
2583
2584        /* Step 1: amount to shift is the partial right-rotate count */
2585        rotate=set->digits-rotate;      /* make it right-rotate */
2586        units=rotate/DECDPUN;           /* whole units to rotate */
2587        shift=rotate%DECDPUN;           /* left-over digits count */
2588        if (shift>0) {                  /* not an exact number of units */
2589          uInt save=res->lsu[0]%powers[shift];    /* save low digit(s) */
2590          decShiftToLeast(res->lsu, D2U(res->digits), shift);
2591          if (shift>msudigits) {        /* msumax-1 needs >0 digits */
2592            uInt rem=save%powers[shift-msudigits];/* split save */
2593            *msumax=(Unit)(save/powers[shift-msudigits]); /* and insert */
2594            *(msumax-1)=*(msumax-1)
2595                       +(Unit)(rem*powers[DECDPUN-(shift-msudigits)]); /* .. */
2596            }
2597           else { /* all fits in msumax */
2598            *msumax=*msumax+(Unit)(save*powers[msudigits-shift]); /* [maybe *1] */
2599            }
2600          } /* digits shift needed */
2601
2602        /* If whole units to rotate... */
2603        if (units>0) {                  /* some to do */
2604          /* Step 2: the units to touch are the whole ones in rotate, */
2605          /*   if any, and the shift is DECDPUN-msudigits (which may be */
2606          /*   0, again) */
2607          shift=DECDPUN-msudigits;
2608          if (shift>0) {                /* not an exact number of units */
2609            uInt save=res->lsu[0]%powers[shift];  /* save low digit(s) */
2610            decShiftToLeast(res->lsu, units, shift);
2611            *msumax=*msumax+(Unit)(save*powers[msudigits]);
2612            } /* partial shift needed */
2613
2614          /* Step 3: rotate the units array using triple reverse */
2615          /* (reversing is easy and fast) */
2616          decReverse(res->lsu+units, msumax);     /* left part */
2617          decReverse(res->lsu, res->lsu+units-1); /* right part */
2618          decReverse(res->lsu, msumax);           /* whole */
2619          } /* whole units to rotate */
2620        /* the rotation may have left an undetermined number of zeros */
2621        /* on the left, so true length needs to be calculated */
2622        res->digits=decGetDigits(res->lsu, msumax-res->lsu+1);
2623        } /* rotate needed */
2624      } /* rhs OK */
2625    } /* numerics */
2626  if (status!=0) decStatus(res, status, set);
2627  return res;
2628  } /* decNumberRotate */
2629
2630/* ------------------------------------------------------------------ */
2631/* decNumberSameQuantum -- test for equal exponents                   */
2632/*                                                                    */
2633/*   res is the result number, which will contain either 0 or 1       */
2634/*   lhs is a number to test                                          */
2635/*   rhs is the second (usually a pattern)                            */
2636/*                                                                    */
2637/* No errors are possible and no context is needed.                   */
2638/* ------------------------------------------------------------------ */
2639decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs,
2640                                 const decNumber *rhs) {
2641  Unit ret=0;                      /* return value */
2642
2643  #if DECCHECK
2644  if (decCheckOperands(res, lhs, rhs, DECUNCONT)) return res;
2645  #endif
2646
2647  if (SPECIALARGS) {
2648    if (decNumberIsNaN(lhs) && decNumberIsNaN(rhs)) ret=1;
2649     else if (decNumberIsInfinite(lhs) && decNumberIsInfinite(rhs)) ret=1;
2650     /* [anything else with a special gives 0] */
2651    }
2652   else if (lhs->exponent==rhs->exponent) ret=1;
2653
2654  decNumberZero(res);              /* OK to overwrite an operand now */
2655  *res->lsu=ret;
2656  return res;
2657  } /* decNumberSameQuantum */
2658
2659/* ------------------------------------------------------------------ */
2660/* decNumberScaleB -- multiply by a power of 10                       */
2661/*                                                                    */
2662/* This computes C = A x 10**B where B is an integer (q=0) with       */
2663/* maximum magnitude 2*(emax+digits)                                  */
2664/*                                                                    */
2665/*   res is C, the result.  C may be A or B                           */
2666/*   lhs is A, the number to adjust                                   */
2667/*   rhs is B, the requested power of ten to use                      */
2668/*   set is the context                                               */
2669/*                                                                    */
2670/* C must have space for set->digits digits.                          */
2671/*                                                                    */
2672/* The result may underflow or overflow.                              */
2673/* ------------------------------------------------------------------ */
2674decNumber * decNumberScaleB(decNumber *res, const decNumber *lhs,
2675                            const decNumber *rhs, decContext *set) {
2676  Int  reqexp;                /* requested exponent change [B] */
2677  uInt status=0;              /* accumulator */
2678  Int  residue;               /* work */
2679
2680  #if DECCHECK
2681  if (decCheckOperands(res, lhs, rhs, set)) return res;
2682  #endif
2683
2684  /* Handle special values except lhs infinite */
2685  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2686    decNaNs(res, lhs, rhs, set, &status);
2687    /* rhs must be an integer */
2688   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2689    status=DEC_Invalid_operation;
2690   else {
2691    /* lhs is a number; rhs is a finite with q==0 */
2692    reqexp=decGetInt(rhs);                   /* [cannot fail] */
2693    if (reqexp==BADINT                       /* something bad .. */
2694     || reqexp==BIGODD || reqexp==BIGEVEN    /* .. very big .. */
2695     || abs(reqexp)>(2*(set->digits+set->emax))) /* .. or out of range */
2696      status=DEC_Invalid_operation;
2697     else {                                  /* rhs is OK */
2698      decNumberCopy(res, lhs);               /* all done if infinite lhs */
2699      if (!decNumberIsInfinite(res)) {       /* prepare to scale */
2700        res->exponent+=reqexp;               /* adjust the exponent */
2701        residue=0;
2702        decFinalize(res, set, &residue, &status); /* .. and check */
2703        } /* finite LHS */
2704      } /* rhs OK */
2705    } /* rhs finite */
2706  if (status!=0) decStatus(res, status, set);
2707  return res;
2708  } /* decNumberScaleB */
2709
2710/* ------------------------------------------------------------------ */
2711/* decNumberShift -- shift the coefficient of a Number left or right  */
2712/*                                                                    */
2713/*   This computes C = A << B or C = A >> -B  (in base ten).          */
2714/*                                                                    */
2715/*   res is C, the result.  C may be A and/or B (e.g., X=X<<X)        */
2716/*   lhs is A                                                         */
2717/*   rhs is B, the number of digits to shift (-ve to right)           */
2718/*   set is the context                                               */
2719/*                                                                    */
2720/* The digits of the coefficient of A are shifted to the left (if B   */
2721/* is positive) or to the right (if B is negative) without adjusting  */
2722/* the exponent or the sign of A.                                     */
2723/*                                                                    */
2724/* B must be an integer (q=0) and in the range -set->digits through   */
2725/* +set->digits.                                                      */
2726/* C must have space for set->digits digits.                          */
2727/* NaNs are propagated as usual.  Infinities are unaffected (but      */
2728/* B must be valid).  No status is set unless B is invalid or an      */
2729/* operand is an sNaN.                                                */
2730/* ------------------------------------------------------------------ */
2731decNumber * decNumberShift(decNumber *res, const decNumber *lhs,
2732                           const decNumber *rhs, decContext *set) {
2733  uInt status=0;              /* accumulator */
2734  Int  shift;                 /* rhs as an Int */
2735
2736  #if DECCHECK
2737  if (decCheckOperands(res, lhs, rhs, set)) return res;
2738  #endif
2739
2740  /* NaNs propagate as normal */
2741  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2742    decNaNs(res, lhs, rhs, set, &status);
2743   /* rhs must be an integer */
2744   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2745    status=DEC_Invalid_operation;
2746   else { /* both numeric, rhs is an integer */
2747    shift=decGetInt(rhs);                    /* [cannot fail] */
2748    if (shift==BADINT                        /* something bad .. */
2749     || shift==BIGODD || shift==BIGEVEN      /* .. very big .. */
2750     || abs(shift)>set->digits)              /* .. or out of range */
2751      status=DEC_Invalid_operation;
2752     else {                                  /* rhs is OK */
2753      decNumberCopy(res, lhs);
2754      if (shift!=0 && !decNumberIsInfinite(res)) { /* something to do */
2755        if (shift>0) {                       /* to left */
2756          if (shift==set->digits) {          /* removing all */
2757            *res->lsu=0;                     /* so place 0 */
2758            res->digits=1;                   /* .. */
2759            }
2760           else {                            /* */
2761            /* first remove leading digits if necessary */
2762            if (res->digits+shift>set->digits) {
2763              decDecap(res, res->digits+shift-set->digits);
2764              /* that updated res->digits; may have gone to 1 (for a */
2765              /* single digit or for zero */
2766              }
2767            if (res->digits>1 || *res->lsu)  /* if non-zero.. */
2768              res->digits=decShiftToMost(res->lsu, res->digits, shift);
2769            } /* partial left */
2770          } /* left */
2771         else { /* to right */
2772          if (-shift>=res->digits) {         /* discarding all */
2773            *res->lsu=0;                     /* so place 0 */
2774            res->digits=1;                   /* .. */
2775            }
2776           else {
2777            decShiftToLeast(res->lsu, D2U(res->digits), -shift);
2778            res->digits-=(-shift);
2779            }
2780          } /* to right */
2781        } /* non-0 non-Inf shift */
2782      } /* rhs OK */
2783    } /* numerics */
2784  if (status!=0) decStatus(res, status, set);
2785  return res;
2786  } /* decNumberShift */
2787
2788/* ------------------------------------------------------------------ */
2789/* decNumberSquareRoot -- square root operator                        */
2790/*                                                                    */
2791/*   This computes C = squareroot(A)                                  */
2792/*                                                                    */
2793/*   res is C, the result.  C may be A                                */
2794/*   rhs is A                                                         */
2795/*   set is the context; note that rounding mode has no effect        */
2796/*                                                                    */
2797/* C must have space for set->digits digits.                          */
2798/* ------------------------------------------------------------------ */
2799/* This uses the following varying-precision algorithm in:            */
2800/*                                                                    */
2801/*   Properly Rounded Variable Precision Square Root, T. E. Hull and  */
2802/*   A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
2803/*   pp229-237, ACM, September 1985.                                  */
2804/*                                                                    */
2805/* The square-root is calculated using Newton's method, after which   */
2806/* a check is made to ensure the result is correctly rounded.         */
2807/*                                                                    */
2808/* % [Reformatted original Numerical Turing source code follows.]     */
2809/* function sqrt(x : real) : real                                     */
2810/* % sqrt(x) returns the properly rounded approximation to the square */
2811/* % root of x, in the precision of the calling environment, or it    */
2812/* % fails if x < 0.                                                  */
2813/* % t e hull and a abrham, august, 1984                              */
2814/* if x <= 0 then                                                     */
2815/*   if x < 0 then                                                    */
2816/*     assert false                                                   */
2817/*   else                                                             */
2818/*     result 0                                                       */
2819/*   end if                                                           */
2820/* end if                                                             */
2821/* var f := setexp(x, 0)  % fraction part of x   [0.1 <= x < 1]       */
2822/* var e := getexp(x)     % exponent part of x                        */
2823/* var approx : real                                                  */
2824/* if e mod 2 = 0  then                                               */
2825/*   approx := .259 + .819 * f   % approx to root of f                */
2826/* else                                                               */
2827/*   f := f/l0                   % adjustments                        */
2828/*   e := e + 1                  %   for odd                          */
2829/*   approx := .0819 + 2.59 * f  %   exponent                         */
2830/* end if                                                             */
2831/*                                                                    */
2832/* var p:= 3                                                          */
2833/* const maxp := currentprecision + 2                                 */
2834/* loop                                                               */
2835/*   p := min(2*p - 2, maxp)     % p = 4,6,10, . . . , maxp           */
2836/*   precision p                                                      */
2837/*   approx := .5 * (approx + f/approx)                               */
2838/*   exit when p = maxp                                               */
2839/* end loop                                                           */
2840/*                                                                    */
2841/* % approx is now within 1 ulp of the properly rounded square root   */
2842/* % of f; to ensure proper rounding, compare squares of (approx -    */
2843/* % l/2 ulp) and (approx + l/2 ulp) with f.                          */
2844/* p := currentprecision                                              */
2845/* begin                                                              */
2846/*   precision p + 2                                                  */
2847/*   const approxsubhalf := approx - setexp(.5, -p)                   */
2848/*   if mulru(approxsubhalf, approxsubhalf) > f then                  */
2849/*     approx := approx - setexp(.l, -p + 1)                          */
2850/*   else                                                             */
2851/*     const approxaddhalf := approx + setexp(.5, -p)                 */
2852/*     if mulrd(approxaddhalf, approxaddhalf) < f then                */
2853/*       approx := approx + setexp(.l, -p + 1)                        */
2854/*     end if                                                         */
2855/*   end if                                                           */
2856/* end                                                                */
2857/* result setexp(approx, e div 2)  % fix exponent                     */
2858/* end sqrt                                                           */
2859/* ------------------------------------------------------------------ */
2860decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs,
2861                                decContext *set) {
2862  decContext workset, approxset;   /* work contexts */
2863  decNumber dzero;                 /* used for constant zero */
2864  Int  maxp;                       /* largest working precision */
2865  Int  workp;                      /* working precision */
2866  Int  residue=0;                  /* rounding residue */
2867  uInt status=0, ignore=0;         /* status accumulators */
2868  uInt rstatus;                    /* .. */
2869  Int  exp;                        /* working exponent */
2870  Int  ideal;                      /* ideal (preferred) exponent */
2871  Int  needbytes;                  /* work */
2872  Int  dropped;                    /* .. */
2873
2874  #if DECSUBSET
2875  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
2876  #endif
2877  /* buffer for f [needs +1 in case DECBUFFER 0] */
2878  decNumber buff[D2N(DECBUFFER+1)];
2879  /* buffer for a [needs +2 to match likely maxp] */
2880  decNumber bufa[D2N(DECBUFFER+2)];
2881  /* buffer for temporary, b [must be same size as a] */
2882  decNumber bufb[D2N(DECBUFFER+2)];
2883  decNumber *allocbuff=NULL;       /* -> allocated buff, iff allocated */
2884  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
2885  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
2886  decNumber *f=buff;               /* reduced fraction */
2887  decNumber *a=bufa;               /* approximation to result */
2888  decNumber *b=bufb;               /* intermediate result */
2889  /* buffer for temporary variable, up to 3 digits */
2890  decNumber buft[D2N(3)];
2891  decNumber *t=buft;               /* up-to-3-digit constant or work */
2892
2893  #if DECCHECK
2894  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
2895  #endif
2896
2897  do {                             /* protect allocated storage */
2898    #if DECSUBSET
2899    if (!set->extended) {
2900      /* reduce operand and set lostDigits status, as needed */
2901      if (rhs->digits>set->digits) {
2902        allocrhs=decRoundOperand(rhs, set, &status);
2903        if (allocrhs==NULL) break;
2904        /* [Note: 'f' allocation below could reuse this buffer if */
2905        /* used, but as this is rare they are kept separate for clarity.] */
2906        rhs=allocrhs;
2907        }
2908      }
2909    #endif
2910    /* [following code does not require input rounding] */
2911
2912    /* handle infinities and NaNs */
2913    if (SPECIALARG) {
2914      if (decNumberIsInfinite(rhs)) {         /* an infinity */
2915        if (decNumberIsNegative(rhs)) status|=DEC_Invalid_operation;
2916         else decNumberCopy(res, rhs);        /* +Infinity */
2917        }
2918       else decNaNs(res, rhs, NULL, set, &status); /* a NaN */
2919      break;
2920      }
2921
2922    /* calculate the ideal (preferred) exponent [floor(exp/2)] */
2923    /* [We would like to write: ideal=rhs->exponent>>1, but this */
2924    /* generates a compiler warning.  Generated code is the same.] */
2925    ideal=(rhs->exponent&~1)/2;         /* target */
2926
2927    /* handle zeros */
2928    if (ISZERO(rhs)) {
2929      decNumberCopy(res, rhs);          /* could be 0 or -0 */
2930      res->exponent=ideal;              /* use the ideal [safe] */
2931      /* use decFinish to clamp any out-of-range exponent, etc. */
2932      decFinish(res, set, &residue, &status);
2933      break;
2934      }
2935
2936    /* any other -x is an oops */
2937    if (decNumberIsNegative(rhs)) {
2938      status|=DEC_Invalid_operation;
2939      break;
2940      }
2941
2942    /* space is needed for three working variables */
2943    /*   f -- the same precision as the RHS, reduced to 0.01->0.99... */
2944    /*   a -- Hull's approximation -- precision, when assigned, is */
2945    /*        currentprecision+1 or the input argument precision, */
2946    /*        whichever is larger (+2 for use as temporary) */
2947    /*   b -- intermediate temporary result (same size as a) */
2948    /* if any is too long for local storage, then allocate */
2949    workp=MAXI(set->digits+1, rhs->digits);  /* actual rounding precision */
2950    maxp=workp+2;                            /* largest working precision */
2951
2952    needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
2953    if (needbytes>(Int)sizeof(buff)) {
2954      allocbuff=(decNumber *)malloc(needbytes);
2955      if (allocbuff==NULL) {  /* hopeless -- abandon */
2956        status|=DEC_Insufficient_storage;
2957        break;}
2958      f=allocbuff;            /* use the allocated space */
2959      }
2960    /* a and b both need to be able to hold a maxp-length number */
2961    needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit);
2962    if (needbytes>(Int)sizeof(bufa)) {            /* [same applies to b] */
2963      allocbufa=(decNumber *)malloc(needbytes);
2964      allocbufb=(decNumber *)malloc(needbytes);
2965      if (allocbufa==NULL || allocbufb==NULL) {   /* hopeless */
2966        status|=DEC_Insufficient_storage;
2967        break;}
2968      a=allocbufa;            /* use the allocated spaces */
2969      b=allocbufb;            /* .. */
2970      }
2971
2972    /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
2973    decNumberCopy(f, rhs);
2974    exp=f->exponent+f->digits;               /* adjusted to Hull rules */
2975    f->exponent=-(f->digits);                /* to range */
2976
2977    /* set up working context */
2978    decContextDefault(&workset, DEC_INIT_DECIMAL64);
2979
2980    /* [Until further notice, no error is possible and status bits */
2981    /* (Rounded, etc.) should be ignored, not accumulated.] */
2982
2983    /* Calculate initial approximation, and allow for odd exponent */
2984    workset.digits=workp;                    /* p for initial calculation */
2985    t->bits=0; t->digits=3;
2986    a->bits=0; a->digits=3;
2987    if ((exp & 1)==0) {                      /* even exponent */
2988      /* Set t=0.259, a=0.819 */
2989      t->exponent=-3;
2990      a->exponent=-3;
2991      #if DECDPUN>=3
2992        t->lsu[0]=259;
2993        a->lsu[0]=819;
2994      #elif DECDPUN==2
2995        t->lsu[0]=59; t->lsu[1]=2;
2996        a->lsu[0]=19; a->lsu[1]=8;
2997      #else
2998        t->lsu[0]=9; t->lsu[1]=5; t->lsu[2]=2;
2999        a->lsu[0]=9; a->lsu[1]=1; a->lsu[2]=8;
3000      #endif
3001      }
3002     else {                                  /* odd exponent */
3003      /* Set t=0.0819, a=2.59 */
3004      f->exponent--;                         /* f=f/10 */
3005      exp++;                                 /* e=e+1 */
3006      t->exponent=-4;
3007      a->exponent=-2;
3008      #if DECDPUN>=3
3009        t->lsu[0]=819;
3010        a->lsu[0]=259;
3011      #elif DECDPUN==2
3012        t->lsu[0]=19; t->lsu[1]=8;
3013        a->lsu[0]=59; a->lsu[1]=2;
3014      #else
3015        t->lsu[0]=9; t->lsu[1]=1; t->lsu[2]=8;
3016        a->lsu[0]=9; a->lsu[1]=5; a->lsu[2]=2;
3017      #endif
3018      }
3019    decMultiplyOp(a, a, f, &workset, &ignore);    /* a=a*f */
3020    decAddOp(a, a, t, &workset, 0, &ignore);      /* ..+t */
3021    /* [a is now the initial approximation for sqrt(f), calculated with */
3022    /* currentprecision, which is also a's precision.] */
3023
3024    /* the main calculation loop */
3025    decNumberZero(&dzero);                   /* make 0 */
3026    decNumberZero(t);                        /* set t = 0.5 */
3027    t->lsu[0]=5;                             /* .. */
3028    t->exponent=-1;                          /* .. */
3029    workset.digits=3;                        /* initial p */
3030    for (;;) {
3031      /* set p to min(2*p - 2, maxp)  [hence 3; or: 4, 6, 10, ... , maxp] */
3032      workset.digits=workset.digits*2-2;
3033      if (workset.digits>maxp) workset.digits=maxp;
3034      /* a = 0.5 * (a + f/a) */
3035      /* [calculated at p then rounded to currentprecision] */
3036      decDivideOp(b, f, a, &workset, DIVIDE, &ignore); /* b=f/a */
3037      decAddOp(b, b, a, &workset, 0, &ignore);    /* b=b+a */
3038      decMultiplyOp(a, b, t, &workset, &ignore);  /* a=b*0.5 */
3039      if (a->digits==maxp) break;            /* have required digits */
3040      } /* loop */
3041
3042    /* Here, 0.1 <= a < 1 [Hull], and a has maxp digits */
3043    /* now reduce to length, etc.; this needs to be done with a */
3044    /* having the correct exponent so as to handle subnormals */
3045    /* correctly */
3046    approxset=*set;                          /* get emin, emax, etc. */
3047    approxset.round=DEC_ROUND_HALF_EVEN;
3048    a->exponent+=exp/2;                      /* set correct exponent */
3049
3050    rstatus=0;                               /* clear status */
3051    residue=0;                               /* .. and accumulator */
3052    decCopyFit(a, a, &approxset, &residue, &rstatus);  /* reduce (if needed) */
3053    decFinish(a, &approxset, &residue, &rstatus);      /* clean and finalize */
3054
3055    /* Overflow was possible if the input exponent was out-of-range, */
3056    /* in which case quit */
3057    if (rstatus&DEC_Overflow) {
3058      status=rstatus;                        /* use the status as-is */
3059      decNumberCopy(res, a);                 /* copy to result */
3060      break;
3061      }
3062
3063    /* Preserve status except Inexact/Rounded */
3064    status|=(rstatus & ~(DEC_Rounded|DEC_Inexact));
3065
3066    /* Carry out the Hull correction */
3067    a->exponent-=exp/2;                      /* back to 0.1->1 */
3068
3069    /* a is now at final precision and within 1 ulp of the properly */
3070    /* rounded square root of f; to ensure proper rounding, compare */
3071    /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
3072    /* Here workset.digits=maxp and t=0.5, and a->digits determines */
3073    /* the ulp */
3074    workset.digits--;                             /* maxp-1 is OK now */
3075    t->exponent=-a->digits-1;                     /* make 0.5 ulp */
3076    decAddOp(b, a, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */
3077    workset.round=DEC_ROUND_UP;
3078    decMultiplyOp(b, b, b, &workset, &ignore);    /* b = mulru(b, b) */
3079    decCompareOp(b, f, b, &workset, COMPARE, &ignore); /* b ? f, reversed */
3080    if (decNumberIsNegative(b)) {                 /* f < b [i.e., b > f] */
3081      /* this is the more common adjustment, though both are rare */
3082      t->exponent++;                              /* make 1.0 ulp */
3083      t->lsu[0]=1;                                /* .. */
3084      decAddOp(a, a, t, &workset, DECNEG, &ignore); /* a = a - 1 ulp */
3085      /* assign to approx [round to length] */
3086      approxset.emin-=exp/2;                      /* adjust to match a */
3087      approxset.emax-=exp/2;
3088      decAddOp(a, &dzero, a, &approxset, 0, &ignore);
3089      }
3090     else {
3091      decAddOp(b, a, t, &workset, 0, &ignore);    /* b = a + 0.5 ulp */
3092      workset.round=DEC_ROUND_DOWN;
3093      decMultiplyOp(b, b, b, &workset, &ignore);  /* b = mulrd(b, b) */
3094      decCompareOp(b, b, f, &workset, COMPARE, &ignore);   /* b ? f */
3095      if (decNumberIsNegative(b)) {               /* b < f */
3096        t->exponent++;                            /* make 1.0 ulp */
3097        t->lsu[0]=1;                              /* .. */
3098        decAddOp(a, a, t, &workset, 0, &ignore);  /* a = a + 1 ulp */
3099        /* assign to approx [round to length] */
3100        approxset.emin-=exp/2;                    /* adjust to match a */
3101        approxset.emax-=exp/2;
3102        decAddOp(a, &dzero, a, &approxset, 0, &ignore);
3103        }
3104      }
3105    /* [no errors are possible in the above, and rounding/inexact during */
3106    /* estimation are irrelevant, so status was not accumulated] */
3107
3108    /* Here, 0.1 <= a < 1  (still), so adjust back */
3109    a->exponent+=exp/2;                      /* set correct exponent */
3110
3111    /* count droppable zeros [after any subnormal rounding] by */
3112    /* trimming a copy */
3113    decNumberCopy(b, a);
3114    decTrim(b, set, 1, &dropped);            /* [drops trailing zeros] */
3115
3116    /* Set Inexact and Rounded.  The answer can only be exact if */
3117    /* it is short enough so that squaring it could fit in workp digits, */
3118    /* and it cannot have trailing zeros due to clamping, so these are */
3119    /* the only (relatively rare) conditions a careful check is needed */
3120    if (b->digits*2-1 > workp && !set->clamp) { /* cannot fit */
3121      status|=DEC_Inexact|DEC_Rounded;
3122      }
3123     else {                                  /* could be exact/unrounded */
3124      uInt mstatus=0;                        /* local status */
3125      decMultiplyOp(b, b, b, &workset, &mstatus); /* try the multiply */
3126      if (mstatus&DEC_Overflow) {            /* result just won't fit */
3127        status|=DEC_Inexact|DEC_Rounded;
3128        }
3129       else {                                /* plausible */
3130        decCompareOp(t, b, rhs, &workset, COMPARE, &mstatus); /* b ? rhs */
3131        if (!ISZERO(t)) status|=DEC_Inexact|DEC_Rounded; /* not equal */
3132         else {                              /* is Exact */
3133          /* here, dropped is the count of trailing zeros in 'a' */
3134          /* use closest exponent to ideal... */
3135          Int todrop=ideal-a->exponent;      /* most that can be dropped */
3136          if (todrop<0) status|=DEC_Rounded; /* ideally would add 0s */
3137           else {                            /* unrounded */
3138            if (dropped<todrop) {            /* clamp to those available */
3139              todrop=dropped;
3140              status|=DEC_Clamped;
3141              }
3142            if (todrop>0) {                  /* have some to drop */
3143              decShiftToLeast(a->lsu, D2U(a->digits), todrop);
3144              a->exponent+=todrop;           /* maintain numerical value */
3145              a->digits-=todrop;             /* new length */
3146              }
3147            }
3148          }
3149        }
3150      }
3151
3152    /* double-check Underflow, as perhaps the result could not have */
3153    /* been subnormal (initial argument too big), or it is now Exact */
3154    if (status&DEC_Underflow) {
3155      Int ae=rhs->exponent+rhs->digits-1;    /* adjusted exponent */
3156      /* check if truly subnormal */
3157      #if DECEXTFLAG                         /* DEC_Subnormal too */
3158        if (ae>=set->emin*2) status&=~(DEC_Subnormal|DEC_Underflow);
3159      #else
3160        if (ae>=set->emin*2) status&=~DEC_Underflow;
3161      #endif
3162      /* check if truly inexact */
3163      if (!(status&DEC_Inexact)) status&=~DEC_Underflow;
3164      }
3165
3166    decNumberCopy(res, a);                   /* a is now the result */
3167    } while(0);                              /* end protected */
3168
3169  if (allocbuff!=NULL) free(allocbuff);      /* drop any storage used */
3170  if (allocbufa!=NULL) free(allocbufa);      /* .. */
3171  if (allocbufb!=NULL) free(allocbufb);      /* .. */
3172  #if DECSUBSET
3173  if (allocrhs !=NULL) free(allocrhs);       /* .. */
3174  #endif
3175  if (status!=0) decStatus(res, status, set);/* then report status */
3176  #if DECCHECK
3177  decCheckInexact(res, set);
3178  #endif
3179  return res;
3180  } /* decNumberSquareRoot */
3181
3182/* ------------------------------------------------------------------ */
3183/* decNumberSubtract -- subtract two Numbers                          */
3184/*                                                                    */
3185/*   This computes C = A - B                                          */
3186/*                                                                    */
3187/*   res is C, the result.  C may be A and/or B (e.g., X=X-X)         */
3188/*   lhs is A                                                         */
3189/*   rhs is B                                                         */
3190/*   set is the context                                               */
3191/*                                                                    */
3192/* C must have space for set->digits digits.                          */
3193/* ------------------------------------------------------------------ */
3194decNumber * decNumberSubtract(decNumber *res, const decNumber *lhs,
3195                              const decNumber *rhs, decContext *set) {
3196  uInt status=0;                        /* accumulator */
3197
3198  decAddOp(res, lhs, rhs, set, DECNEG, &status);
3199  if (status!=0) decStatus(res, status, set);
3200  #if DECCHECK
3201  decCheckInexact(res, set);
3202  #endif
3203  return res;
3204  } /* decNumberSubtract */
3205
3206/* ------------------------------------------------------------------ */
3207/* decNumberToIntegralExact -- round-to-integral-value with InExact   */
3208/* decNumberToIntegralValue -- round-to-integral-value                */
3209/*                                                                    */
3210/*   res is the result                                                */
3211/*   rhs is input number                                              */
3212/*   set is the context                                               */
3213/*                                                                    */
3214/* res must have space for any value of rhs.                          */
3215/*                                                                    */
3216/* This implements the IEEE special operators and therefore treats    */
3217/* special values as valid.  For finite numbers it returns            */
3218/* rescale(rhs, 0) if rhs->exponent is <0.                            */
3219/* Otherwise the result is rhs (so no error is possible, except for   */
3220/* sNaN).                                                             */
3221/*                                                                    */
3222/* The context is used for rounding mode and status after sNaN, but   */
3223/* the digits setting is ignored.  The Exact version will signal      */
3224/* Inexact if the result differs numerically from rhs; the other      */
3225/* never signals Inexact.                                             */
3226/* ------------------------------------------------------------------ */
3227decNumber * decNumberToIntegralExact(decNumber *res, const decNumber *rhs,
3228                                     decContext *set) {
3229  decNumber dn;
3230  decContext workset;              /* working context */
3231  uInt status=0;                   /* accumulator */
3232
3233  #if DECCHECK
3234  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
3235  #endif
3236
3237  /* handle infinities and NaNs */
3238  if (SPECIALARG) {
3239    if (decNumberIsInfinite(rhs)) decNumberCopy(res, rhs); /* an Infinity */
3240     else decNaNs(res, rhs, NULL, set, &status); /* a NaN */
3241    }
3242   else { /* finite */
3243    /* have a finite number; no error possible (res must be big enough) */
3244    if (rhs->exponent>=0) return decNumberCopy(res, rhs);
3245    /* that was easy, but if negative exponent there is work to do... */
3246    workset=*set;                  /* clone rounding, etc. */
3247    workset.digits=rhs->digits;    /* no length rounding */
3248    workset.traps=0;               /* no traps */
3249    decNumberZero(&dn);            /* make a number with exponent 0 */
3250    decNumberQuantize(res, rhs, &dn, &workset);
3251    status|=workset.status;
3252    }
3253  if (status!=0) decStatus(res, status, set);
3254  return res;
3255  } /* decNumberToIntegralExact */
3256
3257decNumber * decNumberToIntegralValue(decNumber *res, const decNumber *rhs,
3258                                     decContext *set) {
3259  decContext workset=*set;         /* working context */
3260  workset.traps=0;                 /* no traps */
3261  decNumberToIntegralExact(res, rhs, &workset);
3262  /* this never affects set, except for sNaNs; NaN will have been set */
3263  /* or propagated already, so no need to call decStatus */
3264  set->status|=workset.status&DEC_Invalid_operation;
3265  return res;
3266  } /* decNumberToIntegralValue */
3267
3268/* ------------------------------------------------------------------ */
3269/* decNumberXor -- XOR two Numbers, digitwise                         */
3270/*                                                                    */
3271/*   This computes C = A ^ B                                          */
3272/*                                                                    */
3273/*   res is C, the result.  C may be A and/or B (e.g., X=X^X)         */
3274/*   lhs is A                                                         */
3275/*   rhs is B                                                         */
3276/*   set is the context (used for result length and error report)     */
3277/*                                                                    */
3278/* C must have space for set->digits digits.                          */
3279/*                                                                    */
3280/* Logical function restrictions apply (see above); a NaN is          */
3281/* returned with Invalid_operation if a restriction is violated.      */
3282/* ------------------------------------------------------------------ */
3283decNumber * decNumberXor(decNumber *res, const decNumber *lhs,
3284                         const decNumber *rhs, decContext *set) {
3285  const Unit *ua, *ub;                  /* -> operands */
3286  const Unit *msua, *msub;              /* -> operand msus */
3287  Unit  *uc, *msuc;                     /* -> result and its msu */
3288  Int   msudigs;                        /* digits in res msu */
3289  #if DECCHECK
3290  if (decCheckOperands(res, lhs, rhs, set)) return res;
3291  #endif
3292
3293  if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
3294   || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
3295    decStatus(res, DEC_Invalid_operation, set);
3296    return res;
3297    }
3298  /* operands are valid */
3299  ua=lhs->lsu;                          /* bottom-up */
3300  ub=rhs->lsu;                          /* .. */
3301  uc=res->lsu;                          /* .. */
3302  msua=ua+D2U(lhs->digits)-1;           /* -> msu of lhs */
3303  msub=ub+D2U(rhs->digits)-1;           /* -> msu of rhs */
3304  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
3305  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
3306  for (; uc<=msuc; ua++, ub++, uc++) {  /* Unit loop */
3307    Unit a, b;                          /* extract units */
3308    if (ua>msua) a=0;
3309     else a=*ua;
3310    if (ub>msub) b=0;
3311     else b=*ub;
3312    *uc=0;                              /* can now write back */
3313    if (a|b) {                          /* maybe 1 bits to examine */
3314      Int i, j;
3315      /* This loop could be unrolled and/or use BIN2BCD tables */
3316      for (i=0; i<DECDPUN; i++) {
3317        if ((a^b)&1) *uc=*uc+(Unit)powers[i];     /* effect XOR */
3318        j=a%10;
3319        a=a/10;
3320        j|=b%10;
3321        b=b/10;
3322        if (j>1) {
3323          decStatus(res, DEC_Invalid_operation, set);
3324          return res;
3325          }
3326        if (uc==msuc && i==msudigs-1) break;      /* just did final digit */
3327        } /* each digit */
3328      } /* non-zero */
3329    } /* each unit */
3330  /* [here uc-1 is the msu of the result] */
3331  res->digits=decGetDigits(res->lsu, uc-res->lsu);
3332  res->exponent=0;                      /* integer */
3333  res->bits=0;                          /* sign=0 */
3334  return res;  /* [no status to set] */
3335  } /* decNumberXor */
3336
3337
3338/* ================================================================== */
3339/* Utility routines                                                   */
3340/* ================================================================== */
3341
3342/* ------------------------------------------------------------------ */
3343/* decNumberClass -- return the decClass of a decNumber               */
3344/*   dn -- the decNumber to test                                      */
3345/*   set -- the context to use for Emin                               */
3346/*   returns the decClass enum                                        */
3347/* ------------------------------------------------------------------ */
3348enum decClass decNumberClass(const decNumber *dn, decContext *set) {
3349  if (decNumberIsSpecial(dn)) {
3350    if (decNumberIsQNaN(dn)) return DEC_CLASS_QNAN;
3351    if (decNumberIsSNaN(dn)) return DEC_CLASS_SNAN;
3352    /* must be an infinity */
3353    if (decNumberIsNegative(dn)) return DEC_CLASS_NEG_INF;
3354    return DEC_CLASS_POS_INF;
3355    }
3356  /* is finite */
3357  if (decNumberIsNormal(dn, set)) { /* most common */
3358    if (decNumberIsNegative(dn)) return DEC_CLASS_NEG_NORMAL;
3359    return DEC_CLASS_POS_NORMAL;
3360    }
3361  /* is subnormal or zero */
3362  if (decNumberIsZero(dn)) {    /* most common */
3363    if (decNumberIsNegative(dn)) return DEC_CLASS_NEG_ZERO;
3364    return DEC_CLASS_POS_ZERO;
3365    }
3366  if (decNumberIsNegative(dn)) return DEC_CLASS_NEG_SUBNORMAL;
3367  return DEC_CLASS_POS_SUBNORMAL;
3368  } /* decNumberClass */
3369
3370/* ------------------------------------------------------------------ */
3371/* decNumberClassToString -- convert decClass to a string             */
3372/*                                                                    */
3373/*  eclass is a valid decClass                                        */
3374/*  returns a constant string describing the class (max 13+1 chars)   */
3375/* ------------------------------------------------------------------ */
3376const char *decNumberClassToString(enum decClass eclass) {
3377  if (eclass==DEC_CLASS_POS_NORMAL)    return DEC_ClassString_PN;
3378  if (eclass==DEC_CLASS_NEG_NORMAL)    return DEC_ClassString_NN;
3379  if (eclass==DEC_CLASS_POS_ZERO)      return DEC_ClassString_PZ;
3380  if (eclass==DEC_CLASS_NEG_ZERO)      return DEC_ClassString_NZ;
3381  if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
3382  if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
3383  if (eclass==DEC_CLASS_POS_INF)       return DEC_ClassString_PI;
3384  if (eclass==DEC_CLASS_NEG_INF)       return DEC_ClassString_NI;
3385  if (eclass==DEC_CLASS_QNAN)          return DEC_ClassString_QN;
3386  if (eclass==DEC_CLASS_SNAN)          return DEC_ClassString_SN;
3387  return DEC_ClassString_UN;           /* Unknown */
3388  } /* decNumberClassToString */
3389
3390/* ------------------------------------------------------------------ */
3391/* decNumberCopy -- copy a number                                     */
3392/*                                                                    */
3393/*   dest is the target decNumber                                     */
3394/*   src  is the source decNumber                                     */
3395/*   returns dest                                                     */
3396/*                                                                    */
3397/* (dest==src is allowed and is a no-op)                              */
3398/* All fields are updated as required.  This is a utility operation,  */
3399/* so special values are unchanged and no error is possible.          */
3400/* ------------------------------------------------------------------ */
3401decNumber * decNumberCopy(decNumber *dest, const decNumber *src) {
3402
3403  #if DECCHECK
3404  if (src==NULL) return decNumberZero(dest);
3405  #endif
3406
3407  if (dest==src) return dest;                /* no copy required */
3408
3409  /* Use explicit assignments here as structure assignment could copy */
3410  /* more than just the lsu (for small DECDPUN).  This would not affect */
3411  /* the value of the results, but could disturb test harness spill */
3412  /* checking. */
3413  dest->bits=src->bits;
3414  dest->exponent=src->exponent;
3415  dest->digits=src->digits;
3416  dest->lsu[0]=src->lsu[0];
3417  if (src->digits>DECDPUN) {                 /* more Units to come */
3418    const Unit *smsup, *s;                   /* work */
3419    Unit  *d;                                /* .. */
3420    /* memcpy for the remaining Units would be safe as they cannot */
3421    /* overlap.  However, this explicit loop is faster in short cases. */
3422    d=dest->lsu+1;                           /* -> first destination */
3423    smsup=src->lsu+D2U(src->digits);         /* -> source msu+1 */
3424    for (s=src->lsu+1; s<smsup; s++, d++) *d=*s;
3425    }
3426  return dest;
3427  } /* decNumberCopy */
3428
3429/* ------------------------------------------------------------------ */
3430/* decNumberCopyAbs -- quiet absolute value operator                  */
3431/*                                                                    */
3432/*   This sets C = abs(A)                                             */
3433/*                                                                    */
3434/*   res is C, the result.  C may be A                                */
3435/*   rhs is A                                                         */
3436/*                                                                    */
3437/* C must have space for set->digits digits.                          */
3438/* No exception or error can occur; this is a quiet bitwise operation.*/
3439/* See also decNumberAbs for a checking version of this.              */
3440/* ------------------------------------------------------------------ */
3441decNumber * decNumberCopyAbs(decNumber *res, const decNumber *rhs) {
3442  #if DECCHECK
3443  if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
3444  #endif
3445  decNumberCopy(res, rhs);
3446  res->bits&=~DECNEG;                   /* turn off sign */
3447  return res;
3448  } /* decNumberCopyAbs */
3449
3450/* ------------------------------------------------------------------ */
3451/* decNumberCopyNegate -- quiet negate value operator                 */
3452/*                                                                    */
3453/*   This sets C = negate(A)                                          */
3454/*                                                                    */
3455/*   res is C, the result.  C may be A                                */
3456/*   rhs is A                                                         */
3457/*                                                                    */
3458/* C must have space for set->digits digits.                          */
3459/* No exception or error can occur; this is a quiet bitwise operation.*/
3460/* See also decNumberMinus for a checking version of this.            */
3461/* ------------------------------------------------------------------ */
3462decNumber * decNumberCopyNegate(decNumber *res, const decNumber *rhs) {
3463  #if DECCHECK
3464  if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
3465  #endif
3466  decNumberCopy(res, rhs);
3467  res->bits^=DECNEG;                    /* invert the sign */
3468  return res;
3469  } /* decNumberCopyNegate */
3470
3471/* ------------------------------------------------------------------ */
3472/* decNumberCopySign -- quiet copy and set sign operator              */
3473/*                                                                    */
3474/*   This sets C = A with the sign of B                               */
3475/*                                                                    */
3476/*   res is C, the result.  C may be A                                */
3477/*   lhs is A                                                         */
3478/*   rhs is B                                                         */
3479/*                                                                    */
3480/* C must have space for set->digits digits.                          */
3481/* No exception or error can occur; this is a quiet bitwise operation.*/
3482/* ------------------------------------------------------------------ */
3483decNumber * decNumberCopySign(decNumber *res, const decNumber *lhs,
3484                              const decNumber *rhs) {
3485  uByte sign;                           /* rhs sign */
3486  #if DECCHECK
3487  if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
3488  #endif
3489  sign=rhs->bits & DECNEG;              /* save sign bit */
3490  decNumberCopy(res, lhs);
3491  res->bits&=~DECNEG;                   /* clear the sign */
3492  res->bits|=sign;                      /* set from rhs */
3493  return res;
3494  } /* decNumberCopySign */
3495
3496/* ------------------------------------------------------------------ */
3497/* decNumberGetBCD -- get the coefficient in BCD8                     */
3498/*   dn is the source decNumber                                       */
3499/*   bcd is the uInt array that will receive dn->digits BCD bytes,    */
3500/*     most-significant at offset 0                                   */
3501/*   returns bcd                                                      */
3502/*                                                                    */
3503/* bcd must have at least dn->digits bytes.  No error is possible; if */
3504/* dn is a NaN or Infinite, digits must be 1 and the coefficient 0.   */
3505/* ------------------------------------------------------------------ */
3506uByte * decNumberGetBCD(const decNumber *dn, uint8_t *bcd) {
3507  uByte *ub=bcd+dn->digits-1;      /* -> lsd */
3508  const Unit *up=dn->lsu;          /* Unit pointer, -> lsu */
3509
3510  #if DECDPUN==1                   /* trivial simple copy */
3511    for (; ub>=bcd; ub--, up++) *ub=*up;
3512  #else                            /* chopping needed */
3513    uInt u=*up;                    /* work */
3514    uInt cut=DECDPUN;              /* downcounter through unit */
3515    for (; ub>=bcd; ub--) {
3516      *ub=(uByte)(u%10);           /* [*6554 trick inhibits, here] */
3517      u=u/10;
3518      cut--;
3519      if (cut>0) continue;         /* more in this unit */
3520      up++;
3521      u=*up;
3522      cut=DECDPUN;
3523      }
3524  #endif
3525  return bcd;
3526  } /* decNumberGetBCD */
3527
3528/* ------------------------------------------------------------------ */
3529/* decNumberSetBCD -- set (replace) the coefficient from BCD8         */
3530/*   dn is the target decNumber                                       */
3531/*   bcd is the uInt array that will source n BCD bytes, most-        */
3532/*     significant at offset 0                                        */
3533/*   n is the number of digits in the source BCD array (bcd)          */
3534/*   returns dn                                                       */
3535/*                                                                    */
3536/* dn must have space for at least n digits.  No error is possible;   */
3537/* if dn is a NaN, or Infinite, or is to become a zero, n must be 1   */
3538/* and bcd[0] zero.                                                   */
3539/* ------------------------------------------------------------------ */
3540decNumber * decNumberSetBCD(decNumber *dn, const uByte *bcd, uInt n) {
3541  Unit *up = dn->lsu + D2U(n) - 1;      /* -> msu [target pointer] */
3542  const uByte *ub=bcd;                  /* -> source msd */
3543
3544  #if DECDPUN==1                        /* trivial simple copy */
3545    for (; ub<bcd+n; ub++, up--) *up=*ub;
3546  #else                                 /* some assembly needed */
3547    /* calculate how many digits in msu, and hence first cut */
3548    Int cut=MSUDIGITS(n);               /* [faster than remainder] */
3549    for (;up>=dn->lsu; up--) {          /* each Unit from msu */
3550      *up=0;                            /* will take <=DECDPUN digits */
3551      for (; cut>0; ub++, cut--) *up=X10(*up)+*ub;
3552      cut=DECDPUN;                      /* next Unit has all digits */
3553      }
3554  #endif
3555  dn->digits=n;                         /* set digit count */
3556  return dn;
3557  } /* decNumberSetBCD */
3558
3559/* ------------------------------------------------------------------ */
3560/* decNumberIsNormal -- test normality of a decNumber                 */
3561/*   dn is the decNumber to test                                      */
3562/*   set is the context to use for Emin                               */
3563/*   returns 1 if |dn| is finite and >=Nmin, 0 otherwise              */
3564/* ------------------------------------------------------------------ */
3565Int decNumberIsNormal(const decNumber *dn, decContext *set) {
3566  Int ae;                               /* adjusted exponent */
3567  #if DECCHECK
3568  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
3569  #endif
3570
3571  if (decNumberIsSpecial(dn)) return 0; /* not finite */
3572  if (decNumberIsZero(dn)) return 0;    /* not non-zero */
3573
3574  ae=dn->exponent+dn->digits-1;         /* adjusted exponent */
3575  if (ae<set->emin) return 0;           /* is subnormal */
3576  return 1;
3577  } /* decNumberIsNormal */
3578
3579/* ------------------------------------------------------------------ */
3580/* decNumberIsSubnormal -- test subnormality of a decNumber           */
3581/*   dn is the decNumber to test                                      */
3582/*   set is the context to use for Emin                               */
3583/*   returns 1 if |dn| is finite, non-zero, and <Nmin, 0 otherwise    */
3584/* ------------------------------------------------------------------ */
3585Int decNumberIsSubnormal(const decNumber *dn, decContext *set) {
3586  Int ae;                               /* adjusted exponent */
3587  #if DECCHECK
3588  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
3589  #endif
3590
3591  if (decNumberIsSpecial(dn)) return 0; /* not finite */
3592  if (decNumberIsZero(dn)) return 0;    /* not non-zero */
3593
3594  ae=dn->exponent+dn->digits-1;         /* adjusted exponent */
3595  if (ae<set->emin) return 1;           /* is subnormal */
3596  return 0;
3597  } /* decNumberIsSubnormal */
3598
3599/* ------------------------------------------------------------------ */
3600/* decNumberTrim -- remove insignificant zeros                        */
3601/*                                                                    */
3602/*   dn is the number to trim                                         */
3603/*   returns dn                                                       */
3604/*                                                                    */
3605/* All fields are updated as required.  This is a utility operation,  */
3606/* so special values are unchanged and no error is possible.          */
3607/* ------------------------------------------------------------------ */
3608decNumber * decNumberTrim(decNumber *dn) {
3609  Int  dropped;                    /* work */
3610  decContext set;                  /* .. */
3611  #if DECCHECK
3612  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, DECUNCONT)) return dn;
3613  #endif
3614  decContextDefault(&set, DEC_INIT_BASE);    /* clamp=0 */
3615  return decTrim(dn, &set, 0, &dropped);
3616  } /* decNumberTrim */
3617
3618/* ------------------------------------------------------------------ */
3619/* decNumberVersion -- return the name and version of this module     */
3620/*                                                                    */
3621/* No error is possible.                                              */
3622/* ------------------------------------------------------------------ */
3623const char * decNumberVersion(void) {
3624  return DECVERSION;
3625  } /* decNumberVersion */
3626
3627/* ------------------------------------------------------------------ */
3628/* decNumberZero -- set a number to 0                                 */
3629/*                                                                    */
3630/*   dn is the number to set, with space for one digit                */
3631/*   returns dn                                                       */
3632/*                                                                    */
3633/* No error is possible.                                              */
3634/* ------------------------------------------------------------------ */
3635/* Memset is not used as it is much slower in some environments. */
3636decNumber * decNumberZero(decNumber *dn) {
3637
3638  #if DECCHECK
3639  if (decCheckOperands(dn, DECUNUSED, DECUNUSED, DECUNCONT)) return dn;
3640  #endif
3641
3642  dn->bits=0;
3643  dn->exponent=0;
3644  dn->digits=1;
3645  dn->lsu[0]=0;
3646  return dn;
3647  } /* decNumberZero */
3648
3649/* ================================================================== */
3650/* Local routines                                                     */
3651/* ================================================================== */
3652
3653/* ------------------------------------------------------------------ */
3654/* decToString -- lay out a number into a string                      */
3655/*                                                                    */
3656/*   dn     is the number to lay out                                  */
3657/*   string is where to lay out the number                            */
3658/*   eng    is 1 if Engineering, 0 if Scientific                      */
3659/*                                                                    */
3660/* string must be at least dn->digits+14 characters long              */
3661/* No error is possible.                                              */
3662/*                                                                    */
3663/* Note that this routine can generate a -0 or 0.000.  These are      */
3664/* never generated in subset to-number or arithmetic, but can occur   */
3665/* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234).              */
3666/* ------------------------------------------------------------------ */
3667/* If DECCHECK is enabled the string "?" is returned if a number is */
3668/* invalid. */
3669static void decToString(const decNumber *dn, char *string, Flag eng) {
3670  Int exp=dn->exponent;       /* local copy */
3671  Int e;                      /* E-part value */
3672  Int pre;                    /* digits before the '.' */
3673  Int cut;                    /* for counting digits in a Unit */
3674  char *c=string;             /* work [output pointer] */
3675  const Unit *up=dn->lsu+D2U(dn->digits)-1; /* -> msu [input pointer] */
3676  uInt u, pow;                /* work */
3677
3678  #if DECCHECK
3679  if (decCheckOperands(DECUNRESU, dn, DECUNUSED, DECUNCONT)) {
3680    strcpy(string, "?");
3681    return;}
3682  #endif
3683
3684  if (decNumberIsNegative(dn)) {   /* Negatives get a minus */
3685    *c='-';
3686    c++;
3687    }
3688  if (dn->bits&DECSPECIAL) {       /* Is a special value */
3689    if (decNumberIsInfinite(dn)) {
3690      strcpy(c,   "Inf");
3691      strcpy(c+3, "inity");
3692      return;}
3693    /* a NaN */
3694    if (dn->bits&DECSNAN) {        /* signalling NaN */
3695      *c='s';
3696      c++;
3697      }
3698    strcpy(c, "NaN");
3699    c+=3;                          /* step past */
3700    /* if not a clean non-zero coefficient, that's all there is in a */
3701    /* NaN string */
3702    if (exp!=0 || (*dn->lsu==0 && dn->digits==1)) return;
3703    /* [drop through to add integer] */
3704    }
3705
3706  /* calculate how many digits in msu, and hence first cut */
3707  cut=MSUDIGITS(dn->digits);       /* [faster than remainder] */
3708  cut--;                           /* power of ten for digit */
3709
3710  if (exp==0) {                    /* simple integer [common fastpath] */
3711    for (;up>=dn->lsu; up--) {     /* each Unit from msu */
3712      u=*up;                       /* contains DECDPUN digits to lay out */
3713      for (; cut>=0; c++, cut--) TODIGIT(u, cut, c, pow);
3714      cut=DECDPUN-1;               /* next Unit has all digits */
3715      }
3716    *c='\0';                       /* terminate the string */
3717    return;}
3718
3719  /* non-0 exponent -- assume plain form */
3720  pre=dn->digits+exp;              /* digits before '.' */
3721  e=0;                             /* no E */
3722  if ((exp>0) || (pre<-5)) {       /* need exponential form */
3723    e=exp+dn->digits-1;            /* calculate E value */
3724    pre=1;                         /* assume one digit before '.' */
3725    if (eng && (e!=0)) {           /* engineering: may need to adjust */
3726      Int adj;                     /* adjustment */
3727      /* The C remainder operator is undefined for negative numbers, so */
3728      /* a positive remainder calculation must be used here */
3729      if (e<0) {
3730        adj=(-e)%3;
3731        if (adj!=0) adj=3-adj;
3732        }
3733       else { /* e>0 */
3734        adj=e%3;
3735        }
3736      e=e-adj;
3737      /* if dealing with zero still produce an exponent which is a */
3738      /* multiple of three, as expected, but there will only be the */
3739      /* one zero before the E, still.  Otherwise note the padding. */
3740      if (!ISZERO(dn)) pre+=adj;
3741       else {  /* is zero */
3742        if (adj!=0) {              /* 0.00Esnn needed */
3743          e=e+3;
3744          pre=-(2-adj);
3745          }
3746        } /* zero */
3747      } /* eng */
3748    } /* need exponent */
3749
3750  /* lay out the digits of the coefficient, adding 0s and . as needed */
3751  u=*up;
3752  if (pre>0) {                     /* xxx.xxx or xx00 (engineering) form */
3753    Int n=pre;
3754    for (; pre>0; pre--, c++, cut--) {
3755      if (cut<0) {                 /* need new Unit */
3756        if (up==dn->lsu) break;    /* out of input digits (pre>digits) */
3757        up--;
3758        cut=DECDPUN-1;
3759        u=*up;
3760        }
3761      TODIGIT(u, cut, c, pow);
3762      }
3763    if (n<dn->digits) {            /* more to come, after '.' */
3764      *c='.'; c++;
3765      for (;; c++, cut--) {
3766        if (cut<0) {               /* need new Unit */
3767          if (up==dn->lsu) break;  /* out of input digits */
3768          up--;
3769          cut=DECDPUN-1;
3770          u=*up;
3771          }
3772        TODIGIT(u, cut, c, pow);
3773        }
3774      }
3775     else for (; pre>0; pre--, c++) *c='0'; /* 0 padding (for engineering) needed */
3776    }
3777   else {                          /* 0.xxx or 0.000xxx form */
3778    *c='0'; c++;
3779    *c='.'; c++;
3780    for (; pre<0; pre++, c++) *c='0';   /* add any 0's after '.' */
3781    for (; ; c++, cut--) {
3782      if (cut<0) {                 /* need new Unit */
3783        if (up==dn->lsu) break;    /* out of input digits */
3784        up--;
3785        cut=DECDPUN-1;
3786        u=*up;
3787        }
3788      TODIGIT(u, cut, c, pow);
3789      }
3790    }
3791
3792  /* Finally add the E-part, if needed.  It will never be 0, has a
3793     base maximum and minimum of +999999999 through -999999999, but
3794     could range down to -1999999998 for anormal numbers */
3795  if (e!=0) {
3796    Flag had=0;               /* 1=had non-zero */
3797    *c='E'; c++;
3798    *c='+'; c++;              /* assume positive */
3799    u=e;                      /* .. */
3800    if (e<0) {
3801      *(c-1)='-';             /* oops, need - */
3802      u=-e;                   /* uInt, please */
3803      }
3804    /* lay out the exponent [_itoa or equivalent is not ANSI C] */
3805    for (cut=9; cut>=0; cut--) {
3806      TODIGIT(u, cut, c, pow);
3807      if (*c=='0' && !had) continue;    /* skip leading zeros */
3808      had=1;                            /* had non-0 */
3809      c++;                              /* step for next */
3810      } /* cut */
3811    }
3812  *c='\0';          /* terminate the string (all paths) */
3813  return;
3814  } /* decToString */
3815
3816/* ------------------------------------------------------------------ */
3817/* decAddOp -- add/subtract operation                                 */
3818/*                                                                    */
3819/*   This computes C = A + B                                          */
3820/*                                                                    */
3821/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
3822/*   lhs is A                                                         */
3823/*   rhs is B                                                         */
3824/*   set is the context                                               */
3825/*   negate is DECNEG if rhs should be negated, or 0 otherwise        */
3826/*   status accumulates status for the caller                         */
3827/*                                                                    */
3828/* C must have space for set->digits digits.                          */
3829/* Inexact in status must be 0 for correct Exact zero sign in result  */
3830/* ------------------------------------------------------------------ */
3831/* If possible, the coefficient is calculated directly into C.        */
3832/* However, if:                                                       */
3833/*   -- a digits+1 calculation is needed because the numbers are      */
3834/*      unaligned and span more than set->digits digits               */
3835/*   -- a carry to digits+1 digits looks possible                     */
3836/*   -- C is the same as A or B, and the result would destructively   */
3837/*      overlap the A or B coefficient                                */
3838/* then the result must be calculated into a temporary buffer.  In    */
3839/* this case a local (stack) buffer is used if possible, and only if  */
3840/* too long for that does malloc become the final resort.             */
3841/*                                                                    */
3842/* Misalignment is handled as follows:                                */
3843/*   Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp.    */
3844/*   BPad: Apply the padding by a combination of shifting (whole      */
3845/*         units) and multiplication (part units).                    */
3846/*                                                                    */
3847/* Addition, especially x=x+1, is speed-critical.                     */
3848/* The static buffer is larger than might be expected to allow for    */
3849/* calls from higher-level functions (notably exp).                   */
3850/* ------------------------------------------------------------------ */
3851static decNumber * decAddOp(decNumber *res, const decNumber *lhs,
3852                            const decNumber *rhs, decContext *set,
3853                            uByte negate, uInt *status) {
3854  #if DECSUBSET
3855  decNumber *alloclhs=NULL;        /* non-NULL if rounded lhs allocated */
3856  decNumber *allocrhs=NULL;        /* .., rhs */
3857  #endif
3858  Int   rhsshift;                  /* working shift (in Units) */
3859  Int   maxdigits;                 /* longest logical length */
3860  Int   mult;                      /* multiplier */
3861  Int   residue;                   /* rounding accumulator */
3862  uByte bits;                      /* result bits */
3863  Flag  diffsign;                  /* non-0 if arguments have different sign */
3864  Unit  *acc;                      /* accumulator for result */
3865  Unit  accbuff[SD2U(DECBUFFER*2+20)]; /* local buffer [*2+20 reduces many */
3866                                   /* allocations when called from */
3867                                   /* other operations, notable exp] */
3868  Unit  *allocacc=NULL;            /* -> allocated acc buffer, iff allocated */
3869  Int   reqdigits=set->digits;     /* local copy; requested DIGITS */
3870  Int   padding;                   /* work */
3871
3872  #if DECCHECK
3873  if (decCheckOperands(res, lhs, rhs, set)) return res;
3874  #endif
3875
3876  do {                             /* protect allocated storage */
3877    #if DECSUBSET
3878    if (!set->extended) {
3879      /* reduce operands and set lostDigits status, as needed */
3880      if (lhs->digits>reqdigits) {
3881        alloclhs=decRoundOperand(lhs, set, status);
3882        if (alloclhs==NULL) break;
3883        lhs=alloclhs;
3884        }
3885      if (rhs->digits>reqdigits) {
3886        allocrhs=decRoundOperand(rhs, set, status);
3887        if (allocrhs==NULL) break;
3888        rhs=allocrhs;
3889        }
3890      }
3891    #endif
3892    /* [following code does not require input rounding] */
3893
3894    /* note whether signs differ [used all paths] */
3895    diffsign=(Flag)((lhs->bits^rhs->bits^negate)&DECNEG);
3896
3897    /* handle infinities and NaNs */
3898    if (SPECIALARGS) {                  /* a special bit set */
3899      if (SPECIALARGS & (DECSNAN | DECNAN))  /* a NaN */
3900        decNaNs(res, lhs, rhs, set, status);
3901       else { /* one or two infinities */
3902        if (decNumberIsInfinite(lhs)) { /* LHS is infinity */
3903          /* two infinities with different signs is invalid */
3904          if (decNumberIsInfinite(rhs) && diffsign) {
3905            *status|=DEC_Invalid_operation;
3906            break;
3907            }
3908          bits=lhs->bits & DECNEG;      /* get sign from LHS */
3909          }
3910         else bits=(rhs->bits^negate) & DECNEG;/* RHS must be Infinity */
3911        bits|=DECINF;
3912        decNumberZero(res);
3913        res->bits=bits;                 /* set +/- infinity */
3914        } /* an infinity */
3915      break;
3916      }
3917
3918    /* Quick exit for add 0s; return the non-0, modified as need be */
3919    if (ISZERO(lhs)) {
3920      Int adjust;                       /* work */
3921      Int lexp=lhs->exponent;           /* save in case LHS==RES */
3922      bits=lhs->bits;                   /* .. */
3923      residue=0;                        /* clear accumulator */
3924      decCopyFit(res, rhs, set, &residue, status); /* copy (as needed) */
3925      res->bits^=negate;                /* flip if rhs was negated */
3926      #if DECSUBSET
3927      if (set->extended) {              /* exponents on zeros count */
3928      #endif
3929        /* exponent will be the lower of the two */
3930        adjust=lexp-res->exponent;      /* adjustment needed [if -ve] */
3931        if (ISZERO(res)) {              /* both 0: special IEEE 854 rules */
3932          if (adjust<0) res->exponent=lexp;  /* set exponent */
3933          /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
3934          if (diffsign) {
3935            if (set->round!=DEC_ROUND_FLOOR) res->bits=0;
3936             else res->bits=DECNEG;     /* preserve 0 sign */
3937            }
3938          }
3939         else { /* non-0 res */
3940          if (adjust<0) {     /* 0-padding needed */
3941            if ((res->digits-adjust)>set->digits) {
3942              adjust=res->digits-set->digits;     /* to fit exactly */
3943              *status|=DEC_Rounded;               /* [but exact] */
3944              }
3945            res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
3946            res->exponent+=adjust;                /* set the exponent. */
3947            }
3948          } /* non-0 res */
3949      #if DECSUBSET
3950        } /* extended */
3951      #endif
3952      decFinish(res, set, &residue, status);      /* clean and finalize */
3953      break;}
3954
3955    if (ISZERO(rhs)) {                  /* [lhs is non-zero] */
3956      Int adjust;                       /* work */
3957      Int rexp=rhs->exponent;           /* save in case RHS==RES */
3958      bits=rhs->bits;                   /* be clean */
3959      residue=0;                        /* clear accumulator */
3960      decCopyFit(res, lhs, set, &residue, status); /* copy (as needed) */
3961      #if DECSUBSET
3962      if (set->extended) {              /* exponents on zeros count */
3963      #endif
3964        /* exponent will be the lower of the two */
3965        /* [0-0 case handled above] */
3966        adjust=rexp-res->exponent;      /* adjustment needed [if -ve] */
3967        if (adjust<0) {     /* 0-padding needed */
3968          if ((res->digits-adjust)>set->digits) {
3969            adjust=res->digits-set->digits;     /* to fit exactly */
3970            *status|=DEC_Rounded;               /* [but exact] */
3971            }
3972          res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
3973          res->exponent+=adjust;                /* set the exponent. */
3974          }
3975      #if DECSUBSET
3976        } /* extended */
3977      #endif
3978      decFinish(res, set, &residue, status);      /* clean and finalize */
3979      break;}
3980
3981    /* [NB: both fastpath and mainpath code below assume these cases */
3982    /* (notably 0-0) have already been handled] */
3983
3984    /* calculate the padding needed to align the operands */
3985    padding=rhs->exponent-lhs->exponent;
3986
3987    /* Fastpath cases where the numbers are aligned and normal, the RHS */
3988    /* is all in one unit, no operand rounding is needed, and no carry, */
3989    /* lengthening, or borrow is needed */
3990    if (padding==0
3991        && rhs->digits<=DECDPUN
3992        && rhs->exponent>=set->emin     /* [some normals drop through] */
3993        && rhs->exponent<=set->emax-set->digits+1 /* [could clamp] */
3994        && rhs->digits<=reqdigits
3995        && lhs->digits<=reqdigits) {
3996      Int partial=*lhs->lsu;
3997      if (!diffsign) {                  /* adding */
3998        partial+=*rhs->lsu;
3999        if ((partial<=DECDPUNMAX)       /* result fits in unit */
4000         && (lhs->digits>=DECDPUN ||    /* .. and no digits-count change */
4001             partial<(Int)powers[lhs->digits])) { /* .. */
4002          if (res!=lhs) decNumberCopy(res, lhs);  /* not in place */
4003          *res->lsu=(Unit)partial;      /* [copy could have overwritten RHS] */
4004          break;
4005          }
4006        /* else drop out for careful add */
4007        }
4008       else {                           /* signs differ */
4009        partial-=*rhs->lsu;
4010        if (partial>0) { /* no borrow needed, and non-0 result */
4011          if (res!=lhs) decNumberCopy(res, lhs);  /* not in place */
4012          *res->lsu=(Unit)partial;
4013          /* this could have reduced digits [but result>0] */
4014          res->digits=decGetDigits(res->lsu, D2U(res->digits));
4015          break;
4016          }
4017        /* else drop out for careful subtract */
4018        }
4019      }
4020
4021    /* Now align (pad) the lhs or rhs so they can be added or */
4022    /* subtracted, as necessary.  If one number is much larger than */
4023    /* the other (that is, if in plain form there is a least one */
4024    /* digit between the lowest digit of one and the highest of the */
4025    /* other) padding with up to DIGITS-1 trailing zeros may be */
4026    /* needed; then apply rounding (as exotic rounding modes may be */
4027    /* affected by the residue). */
4028    rhsshift=0;               /* rhs shift to left (padding) in Units */
4029    bits=lhs->bits;           /* assume sign is that of LHS */
4030    mult=1;                   /* likely multiplier */
4031
4032    /* [if padding==0 the operands are aligned; no padding is needed] */
4033    if (padding!=0) {
4034      /* some padding needed; always pad the RHS, as any required */
4035      /* padding can then be effected by a simple combination of */
4036      /* shifts and a multiply */
4037      Flag swapped=0;
4038      if (padding<0) {                  /* LHS needs the padding */
4039        const decNumber *t;
4040        padding=-padding;               /* will be +ve */
4041        bits=(uByte)(rhs->bits^negate); /* assumed sign is now that of RHS */
4042        t=lhs; lhs=rhs; rhs=t;
4043        swapped=1;
4044        }
4045
4046      /* If, after pad, rhs would be longer than lhs by digits+1 or */
4047      /* more then lhs cannot affect the answer, except as a residue, */
4048      /* so only need to pad up to a length of DIGITS+1. */
4049      if (rhs->digits+padding > lhs->digits+reqdigits+1) {
4050        /* The RHS is sufficient */
4051        /* for residue use the relative sign indication... */
4052        Int shift=reqdigits-rhs->digits;     /* left shift needed */
4053        residue=1;                           /* residue for rounding */
4054        if (diffsign) residue=-residue;      /* signs differ */
4055        /* copy, shortening if necessary */
4056        decCopyFit(res, rhs, set, &residue, status);
4057        /* if it was already shorter, then need to pad with zeros */
4058        if (shift>0) {
4059          res->digits=decShiftToMost(res->lsu, res->digits, shift);
4060          res->exponent-=shift;              /* adjust the exponent. */
4061          }
4062        /* flip the result sign if unswapped and rhs was negated */
4063        if (!swapped) res->bits^=negate;
4064        decFinish(res, set, &residue, status);    /* done */
4065        break;}
4066
4067      /* LHS digits may affect result */
4068      rhsshift=D2U(padding+1)-1;        /* this much by Unit shift .. */
4069      mult=powers[padding-(rhsshift*DECDPUN)]; /* .. this by multiplication */
4070      } /* padding needed */
4071
4072    if (diffsign) mult=-mult;           /* signs differ */
4073
4074    /* determine the longer operand */
4075    maxdigits=rhs->digits+padding;      /* virtual length of RHS */
4076    if (lhs->digits>maxdigits) maxdigits=lhs->digits;
4077
4078    /* Decide on the result buffer to use; if possible place directly */
4079    /* into result. */
4080    acc=res->lsu;                       /* assume add direct to result */
4081    /* If destructive overlap, or the number is too long, or a carry or */
4082    /* borrow to DIGITS+1 might be possible, a buffer must be used. */
4083    /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
4084    if ((maxdigits>=reqdigits)          /* is, or could be, too large */
4085     || (res==rhs && rhsshift>0)) {     /* destructive overlap */
4086      /* buffer needed, choose it; units for maxdigits digits will be */
4087      /* needed, +1 Unit for carry or borrow */
4088      Int need=D2U(maxdigits)+1;
4089      acc=accbuff;                      /* assume use local buffer */
4090      if (need*sizeof(Unit)>sizeof(accbuff)) {
4091        /* printf("malloc add %ld %ld\n", need, sizeof(accbuff)); */
4092        allocacc=(Unit *)malloc(need*sizeof(Unit));
4093        if (allocacc==NULL) {           /* hopeless -- abandon */
4094          *status|=DEC_Insufficient_storage;
4095          break;}
4096        acc=allocacc;
4097        }
4098      }
4099
4100    res->bits=(uByte)(bits&DECNEG);     /* it's now safe to overwrite.. */
4101    res->exponent=lhs->exponent;        /* .. operands (even if aliased) */
4102
4103    #if DECTRACE
4104      decDumpAr('A', lhs->lsu, D2U(lhs->digits));
4105      decDumpAr('B', rhs->lsu, D2U(rhs->digits));
4106      printf("  :h: %ld %ld\n", rhsshift, mult);
4107    #endif
4108
4109    /* add [A+B*m] or subtract [A+B*(-m)] */
4110    res->digits=decUnitAddSub(lhs->lsu, D2U(lhs->digits),
4111                              rhs->lsu, D2U(rhs->digits),
4112                              rhsshift, acc, mult)
4113               *DECDPUN;           /* [units -> digits] */
4114    if (res->digits<0) {           /* borrowed... */
4115      res->digits=-res->digits;
4116      res->bits^=DECNEG;           /* flip the sign */
4117      }
4118    #if DECTRACE
4119      decDumpAr('+', acc, D2U(res->digits));
4120    #endif
4121
4122    /* If a buffer was used the result must be copied back, possibly */
4123    /* shortening.  (If no buffer was used then the result must have */
4124    /* fit, so can't need rounding and residue must be 0.) */
4125    residue=0;                     /* clear accumulator */
4126    if (acc!=res->lsu) {
4127      #if DECSUBSET
4128      if (set->extended) {         /* round from first significant digit */
4129      #endif
4130        /* remove leading zeros that were added due to rounding up to */
4131        /* integral Units -- before the test for rounding. */
4132        if (res->digits>reqdigits)
4133          res->digits=decGetDigits(acc, D2U(res->digits));
4134        decSetCoeff(res, set, acc, res->digits, &residue, status);
4135      #if DECSUBSET
4136        }
4137       else { /* subset arithmetic rounds from original significant digit */
4138        /* May have an underestimate.  This only occurs when both */
4139        /* numbers fit in DECDPUN digits and are padding with a */
4140        /* negative multiple (-10, -100...) and the top digit(s) become */
4141        /* 0.  (This only matters when using X3.274 rules where the */
4142        /* leading zero could be included in the rounding.) */
4143        if (res->digits<maxdigits) {
4144          *(acc+D2U(res->digits))=0; /* ensure leading 0 is there */
4145          res->digits=maxdigits;
4146          }
4147         else {
4148          /* remove leading zeros that added due to rounding up to */
4149          /* integral Units (but only those in excess of the original */
4150          /* maxdigits length, unless extended) before test for rounding. */
4151          if (res->digits>reqdigits) {
4152            res->digits=decGetDigits(acc, D2U(res->digits));
4153            if (res->digits<maxdigits) res->digits=maxdigits;
4154            }
4155          }
4156        decSetCoeff(res, set, acc, res->digits, &residue, status);
4157        /* Now apply rounding if needed before removing leading zeros. */
4158        /* This is safe because subnormals are not a possibility */
4159        if (residue!=0) {
4160          decApplyRound(res, set, residue, status);
4161          residue=0;                 /* did what needed to be done */
4162          }
4163        } /* subset */
4164      #endif
4165      } /* used buffer */
4166
4167    /* strip leading zeros [these were left on in case of subset subtract] */
4168    res->digits=decGetDigits(res->lsu, D2U(res->digits));
4169
4170    /* apply checks and rounding */
4171    decFinish(res, set, &residue, status);
4172
4173    /* "When the sum of two operands with opposite signs is exactly */
4174    /* zero, the sign of that sum shall be '+' in all rounding modes */
4175    /* except round toward -Infinity, in which mode that sign shall be */
4176    /* '-'."  [Subset zeros also never have '-', set by decFinish.] */
4177    if (ISZERO(res) && diffsign
4178     #if DECSUBSET
4179     && set->extended
4180     #endif
4181     && (*status&DEC_Inexact)==0) {
4182      if (set->round==DEC_ROUND_FLOOR) res->bits|=DECNEG;   /* sign - */
4183                                  else res->bits&=~DECNEG;  /* sign + */
4184      }
4185    } while(0);                              /* end protected */
4186
4187  if (allocacc!=NULL) free(allocacc);        /* drop any storage used */
4188  #if DECSUBSET
4189  if (allocrhs!=NULL) free(allocrhs);        /* .. */
4190  if (alloclhs!=NULL) free(alloclhs);        /* .. */
4191  #endif
4192  return res;
4193  } /* decAddOp */
4194
4195/* ------------------------------------------------------------------ */
4196/* decDivideOp -- division operation                                  */
4197/*                                                                    */
4198/*  This routine performs the calculations for all four division      */
4199/*  operators (divide, divideInteger, remainder, remainderNear).      */
4200/*                                                                    */
4201/*  C=A op B                                                          */
4202/*                                                                    */
4203/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
4204/*   lhs is A                                                         */
4205/*   rhs is B                                                         */
4206/*   set is the context                                               */
4207/*   op  is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively.    */
4208/*   status is the usual accumulator                                  */
4209/*                                                                    */
4210/* C must have space for set->digits digits.                          */
4211/*                                                                    */
4212/* ------------------------------------------------------------------ */
4213/*   The underlying algorithm of this routine is the same as in the   */
4214/*   1981 S/370 implementation, that is, non-restoring long division  */
4215/*   with bi-unit (rather than bi-digit) estimation for each unit     */
4216/*   multiplier.  In this pseudocode overview, complications for the  */
4217/*   Remainder operators and division residues for exact rounding are */
4218/*   omitted for clarity.                                             */
4219/*                                                                    */
4220/*     Prepare operands and handle special values                     */
4221/*     Test for x/0 and then 0/x                                      */
4222/*     Exp =Exp1 - Exp2                                               */
4223/*     Exp =Exp +len(var1) -len(var2)                                 */
4224/*     Sign=Sign1 * Sign2                                             */
4225/*     Pad accumulator (Var1) to double-length with 0's (pad1)        */
4226/*     Pad Var2 to same length as Var1                                */
4227/*     msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round  */
4228/*     have=0                                                         */
4229/*     Do until (have=digits+1 OR residue=0)                          */
4230/*       if exp<0 then if integer divide/residue then leave           */
4231/*       this_unit=0                                                  */
4232/*       Do forever                                                   */
4233/*          compare numbers                                           */
4234/*          if <0 then leave inner_loop                               */
4235/*          if =0 then (* quick exit without subtract *) do           */
4236/*             this_unit=this_unit+1; output this_unit                */
4237/*             leave outer_loop; end                                  */
4238/*          Compare lengths of numbers (mantissae):                   */
4239/*          If same then tops2=msu2pair -- {units 1&2 of var2}        */
4240/*                  else tops2=msu2plus -- {0, unit 1 of var2}        */
4241/*          tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
4242/*          mult=tops1/tops2  -- Good and safe guess at divisor       */
4243/*          if mult=0 then mult=1                                     */
4244/*          this_unit=this_unit+mult                                  */
4245/*          subtract                                                  */
4246/*          end inner_loop                                            */
4247/*        if have\=0 | this_unit\=0 then do                           */
4248/*          output this_unit                                          */
4249/*          have=have+1; end                                          */
4250/*        var2=var2/10                                                */
4251/*        exp=exp-1                                                   */
4252/*        end outer_loop                                              */
4253/*     exp=exp+1   -- set the proper exponent                         */
4254/*     if have=0 then generate answer=0                               */
4255/*     Return (Result is defined by Var1)                             */
4256/*                                                                    */
4257/* ------------------------------------------------------------------ */
4258/* Two working buffers are needed during the division; one (digits+   */
4259/* 1) to accumulate the result, and the other (up to 2*digits+1) for  */
4260/* long subtractions.  These are acc and var1 respectively.           */
4261/* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
4262/* The static buffers may be larger than might be expected to allow   */
4263/* for calls from higher-level functions (notably exp).               */
4264/* ------------------------------------------------------------------ */
4265static decNumber * decDivideOp(decNumber *res,
4266                               const decNumber *lhs, const decNumber *rhs,
4267                               decContext *set, Flag op, uInt *status) {
4268  #if DECSUBSET
4269  decNumber *alloclhs=NULL;        /* non-NULL if rounded lhs allocated */
4270  decNumber *allocrhs=NULL;        /* .., rhs */
4271  #endif
4272  Unit  accbuff[SD2U(DECBUFFER+DECDPUN+10)]; /* local buffer */
4273  Unit  *acc=accbuff;              /* -> accumulator array for result */
4274  Unit  *allocacc=NULL;            /* -> allocated buffer, iff allocated */
4275  Unit  *accnext;                  /* -> where next digit will go */
4276  Int   acclength;                 /* length of acc needed [Units] */
4277  Int   accunits;                  /* count of units accumulated */
4278  Int   accdigits;                 /* count of digits accumulated */
4279
4280  Unit  varbuff[SD2U(DECBUFFER*2+DECDPUN)*sizeof(Unit)]; /* buffer for var1 */
4281  Unit  *var1=varbuff;             /* -> var1 array for long subtraction */
4282  Unit  *varalloc=NULL;            /* -> allocated buffer, iff used */
4283  Unit  *msu1;                     /* -> msu of var1 */
4284
4285  const Unit *var2;                /* -> var2 array */
4286  const Unit *msu2;                /* -> msu of var2 */
4287  Int   msu2plus;                  /* msu2 plus one [does not vary] */
4288  eInt  msu2pair;                  /* msu2 pair plus one [does not vary] */
4289
4290  Int   var1units, var2units;      /* actual lengths */
4291  Int   var2ulen;                  /* logical length (units) */
4292  Int   var1initpad=0;             /* var1 initial padding (digits) */
4293  Int   maxdigits;                 /* longest LHS or required acc length */
4294  Int   mult;                      /* multiplier for subtraction */
4295  Unit  thisunit;                  /* current unit being accumulated */
4296  Int   residue;                   /* for rounding */
4297  Int   reqdigits=set->digits;     /* requested DIGITS */
4298  Int   exponent;                  /* working exponent */
4299  Int   maxexponent=0;             /* DIVIDE maximum exponent if unrounded */
4300  uByte bits;                      /* working sign */
4301  Unit  *target;                   /* work */
4302  const Unit *source;              /* .. */
4303  uLong const *pow;                /* .. */
4304  Int   shift, cut;                /* .. */
4305  #if DECSUBSET
4306  Int   dropped;                   /* work */
4307  #endif
4308
4309  #if DECCHECK
4310  if (decCheckOperands(res, lhs, rhs, set)) return res;
4311  #endif
4312
4313  do {                             /* protect allocated storage */
4314    #if DECSUBSET
4315    if (!set->extended) {
4316      /* reduce operands and set lostDigits status, as needed */
4317      if (lhs->digits>reqdigits) {
4318        alloclhs=decRoundOperand(lhs, set, status);
4319        if (alloclhs==NULL) break;
4320        lhs=alloclhs;
4321        }
4322      if (rhs->digits>reqdigits) {
4323        allocrhs=decRoundOperand(rhs, set, status);
4324        if (allocrhs==NULL) break;
4325        rhs=allocrhs;
4326        }
4327      }
4328    #endif
4329    /* [following code does not require input rounding] */
4330
4331    bits=(lhs->bits^rhs->bits)&DECNEG;  /* assumed sign for divisions */
4332
4333    /* handle infinities and NaNs */
4334    if (SPECIALARGS) {                  /* a special bit set */
4335      if (SPECIALARGS & (DECSNAN | DECNAN)) { /* one or two NaNs */
4336        decNaNs(res, lhs, rhs, set, status);
4337        break;
4338        }
4339      /* one or two infinities */
4340      if (decNumberIsInfinite(lhs)) {   /* LHS (dividend) is infinite */
4341        if (decNumberIsInfinite(rhs) || /* two infinities are invalid .. */
4342            op & (REMAINDER | REMNEAR)) { /* as is remainder of infinity */
4343          *status|=DEC_Invalid_operation;
4344          break;
4345          }
4346        /* [Note that infinity/0 raises no exceptions] */
4347        decNumberZero(res);
4348        res->bits=bits|DECINF;          /* set +/- infinity */
4349        break;
4350        }
4351       else {                           /* RHS (divisor) is infinite */
4352        residue=0;
4353        if (op&(REMAINDER|REMNEAR)) {
4354          /* result is [finished clone of] lhs */
4355          decCopyFit(res, lhs, set, &residue, status);
4356          }
4357         else {  /* a division */
4358          decNumberZero(res);
4359          res->bits=bits;               /* set +/- zero */
4360          /* for DIVIDEINT the exponent is always 0.  For DIVIDE, result */
4361          /* is a 0 with infinitely negative exponent, clamped to minimum */
4362          if (op&DIVIDE) {
4363            res->exponent=set->emin-set->digits+1;
4364            *status|=DEC_Clamped;
4365            }
4366          }
4367        decFinish(res, set, &residue, status);
4368        break;
4369        }
4370      }
4371
4372    /* handle 0 rhs (x/0) */
4373    if (ISZERO(rhs)) {                  /* x/0 is always exceptional */
4374      if (ISZERO(lhs)) {
4375        decNumberZero(res);             /* [after lhs test] */
4376        *status|=DEC_Division_undefined;/* 0/0 will become NaN */
4377        }
4378       else {
4379        decNumberZero(res);
4380        if (op&(REMAINDER|REMNEAR)) *status|=DEC_Invalid_operation;
4381         else {
4382          *status|=DEC_Division_by_zero; /* x/0 */
4383          res->bits=bits|DECINF;         /* .. is +/- Infinity */
4384          }
4385        }
4386      break;}
4387
4388    /* handle 0 lhs (0/x) */
4389    if (ISZERO(lhs)) {                  /* 0/x [x!=0] */
4390      #if DECSUBSET
4391      if (!set->extended) decNumberZero(res);
4392       else {
4393      #endif
4394        if (op&DIVIDE) {
4395          residue=0;
4396          exponent=lhs->exponent-rhs->exponent; /* ideal exponent */
4397          decNumberCopy(res, lhs);      /* [zeros always fit] */
4398          res->bits=bits;               /* sign as computed */
4399          res->exponent=exponent;       /* exponent, too */
4400          decFinalize(res, set, &residue, status);   /* check exponent */
4401          }
4402         else if (op&DIVIDEINT) {
4403          decNumberZero(res);           /* integer 0 */
4404          res->bits=bits;               /* sign as computed */
4405          }
4406         else {                         /* a remainder */
4407          exponent=rhs->exponent;       /* [save in case overwrite] */
4408          decNumberCopy(res, lhs);      /* [zeros always fit] */
4409          if (exponent<res->exponent) res->exponent=exponent; /* use lower */
4410          }
4411      #if DECSUBSET
4412        }
4413      #endif
4414      break;}
4415
4416    /* Precalculate exponent.  This starts off adjusted (and hence fits */
4417    /* in 31 bits) and becomes the usual unadjusted exponent as the */
4418    /* division proceeds.  The order of evaluation is important, here, */
4419    /* to avoid wrap. */
4420    exponent=(lhs->exponent+lhs->digits)-(rhs->exponent+rhs->digits);
4421
4422    /* If the working exponent is -ve, then some quick exits are */
4423    /* possible because the quotient is known to be <1 */
4424    /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
4425    if (exponent<0 && !(op==DIVIDE)) {
4426      if (op&DIVIDEINT) {
4427        decNumberZero(res);                  /* integer part is 0 */
4428        #if DECSUBSET
4429        if (set->extended)
4430        #endif
4431          res->bits=bits;                    /* set +/- zero */
4432        break;}
4433      /* fastpath remainders so long as the lhs has the smaller */
4434      /* (or equal) exponent */
4435      if (lhs->exponent<=rhs->exponent) {
4436        if (op&REMAINDER || exponent<-1) {
4437          /* It is REMAINDER or safe REMNEAR; result is [finished */
4438          /* clone of] lhs  (r = x - 0*y) */
4439          residue=0;
4440          decCopyFit(res, lhs, set, &residue, status);
4441          decFinish(res, set, &residue, status);
4442          break;
4443          }
4444        /* [unsafe REMNEAR drops through] */
4445        }
4446      } /* fastpaths */
4447
4448    /* Long (slow) division is needed; roll up the sleeves... */
4449
4450    /* The accumulator will hold the quotient of the division. */
4451    /* If it needs to be too long for stack storage, then allocate. */
4452    acclength=D2U(reqdigits+DECDPUN);   /* in Units */
4453    if (acclength*sizeof(Unit)>sizeof(accbuff)) {
4454      /* printf("malloc dvacc %ld units\n", acclength); */
4455      allocacc=(Unit *)malloc(acclength*sizeof(Unit));
4456      if (allocacc==NULL) {             /* hopeless -- abandon */
4457        *status|=DEC_Insufficient_storage;
4458        break;}
4459      acc=allocacc;                     /* use the allocated space */
4460      }
4461
4462    /* var1 is the padded LHS ready for subtractions. */
4463    /* If it needs to be too long for stack storage, then allocate. */
4464    /* The maximum units needed for var1 (long subtraction) is: */
4465    /* Enough for */
4466    /*     (rhs->digits+reqdigits-1) -- to allow full slide to right */
4467    /* or  (lhs->digits)             -- to allow for long lhs */
4468    /* whichever is larger */
4469    /*   +1                -- for rounding of slide to right */
4470    /*   +1                -- for leading 0s */
4471    /*   +1                -- for pre-adjust if a remainder or DIVIDEINT */
4472    /* [Note: unused units do not participate in decUnitAddSub data] */
4473    maxdigits=rhs->digits+reqdigits-1;
4474    if (lhs->digits>maxdigits) maxdigits=lhs->digits;
4475    var1units=D2U(maxdigits)+2;
4476    /* allocate a guard unit above msu1 for REMAINDERNEAR */
4477    if (!(op&DIVIDE)) var1units++;
4478    if ((var1units+1)*sizeof(Unit)>sizeof(varbuff)) {
4479      /* printf("malloc dvvar %ld units\n", var1units+1); */
4480      varalloc=(Unit *)malloc((var1units+1)*sizeof(Unit));
4481      if (varalloc==NULL) {             /* hopeless -- abandon */
4482        *status|=DEC_Insufficient_storage;
4483        break;}
4484      var1=varalloc;                    /* use the allocated space */
4485      }
4486
4487    /* Extend the lhs and rhs to full long subtraction length.  The lhs */
4488    /* is truly extended into the var1 buffer, with 0 padding, so a */
4489    /* subtract in place is always possible.  The rhs (var2) has */
4490    /* virtual padding (implemented by decUnitAddSub). */
4491    /* One guard unit was allocated above msu1 for rem=rem+rem in */
4492    /* REMAINDERNEAR. */
4493    msu1=var1+var1units-1;              /* msu of var1 */
4494    source=lhs->lsu+D2U(lhs->digits)-1; /* msu of input array */
4495    for (target=msu1; source>=lhs->lsu; source--, target--) *target=*source;
4496    for (; target>=var1; target--) *target=0;
4497
4498    /* rhs (var2) is left-aligned with var1 at the start */
4499    var2ulen=var1units;                 /* rhs logical length (units) */
4500    var2units=D2U(rhs->digits);         /* rhs actual length (units) */
4501    var2=rhs->lsu;                      /* -> rhs array */
4502    msu2=var2+var2units-1;              /* -> msu of var2 [never changes] */
4503    /* now set up the variables which will be used for estimating the */
4504    /* multiplication factor.  If these variables are not exact, add */
4505    /* 1 to make sure that the multiplier is never overestimated. */
4506    msu2plus=*msu2;                     /* it's value .. */
4507    if (var2units>1) msu2plus++;        /* .. +1 if any more */
4508    msu2pair=(eInt)*msu2*(DECDPUNMAX+1);/* top two pair .. */
4509    if (var2units>1) {                  /* .. [else treat 2nd as 0] */
4510      msu2pair+=*(msu2-1);              /* .. */
4511      if (var2units>2) msu2pair++;      /* .. +1 if any more */
4512      }
4513
4514    /* The calculation is working in units, which may have leading zeros, */
4515    /* but the exponent was calculated on the assumption that they are */
4516    /* both left-aligned.  Adjust the exponent to compensate: add the */
4517    /* number of leading zeros in var1 msu and subtract those in var2 msu. */
4518    /* [This is actually done by counting the digits and negating, as */
4519    /* lead1=DECDPUN-digits1, and similarly for lead2.] */
4520    for (pow=&powers[1]; *msu1>=*pow; pow++) exponent--;
4521    for (pow=&powers[1]; *msu2>=*pow; pow++) exponent++;
4522
4523    /* Now, if doing an integer divide or remainder, ensure that */
4524    /* the result will be Unit-aligned.  To do this, shift the var1 */
4525    /* accumulator towards least if need be.  (It's much easier to */
4526    /* do this now than to reassemble the residue afterwards, if */
4527    /* doing a remainder.)  Also ensure the exponent is not negative. */
4528    if (!(op&DIVIDE)) {
4529      Unit *u;                          /* work */
4530      /* save the initial 'false' padding of var1, in digits */
4531      var1initpad=(var1units-D2U(lhs->digits))*DECDPUN;
4532      /* Determine the shift to do. */
4533      if (exponent<0) cut=-exponent;
4534       else cut=DECDPUN-exponent%DECDPUN;
4535      decShiftToLeast(var1, var1units, cut);
4536      exponent+=cut;                    /* maintain numerical value */
4537      var1initpad-=cut;                 /* .. and reduce padding */
4538      /* clean any most-significant units which were just emptied */
4539      for (u=msu1; cut>=DECDPUN; cut-=DECDPUN, u--) *u=0;
4540      } /* align */
4541     else { /* is DIVIDE */
4542      maxexponent=lhs->exponent-rhs->exponent;    /* save */
4543      /* optimization: if the first iteration will just produce 0, */
4544      /* preadjust to skip it [valid for DIVIDE only] */
4545      if (*msu1<*msu2) {
4546        var2ulen--;                     /* shift down */
4547        exponent-=DECDPUN;              /* update the exponent */
4548        }
4549      }
4550
4551    /* ---- start the long-division loops ------------------------------ */
4552    accunits=0;                         /* no units accumulated yet */
4553    accdigits=0;                        /* .. or digits */
4554    accnext=acc+acclength-1;            /* -> msu of acc [NB: allows digits+1] */
4555    for (;;) {                          /* outer forever loop */
4556      thisunit=0;                       /* current unit assumed 0 */
4557      /* find the next unit */
4558      for (;;) {                        /* inner forever loop */
4559        /* strip leading zero units [from either pre-adjust or from */
4560        /* subtract last time around].  Leave at least one unit. */
4561        for (; *msu1==0 && msu1>var1; msu1--) var1units--;
4562
4563        if (var1units<var2ulen) break;       /* var1 too low for subtract */
4564        if (var1units==var2ulen) {           /* unit-by-unit compare needed */
4565          /* compare the two numbers, from msu */
4566          const Unit *pv1, *pv2;
4567          Unit v2;                           /* units to compare */
4568          pv2=msu2;                          /* -> msu */
4569          for (pv1=msu1; ; pv1--, pv2--) {
4570            /* v1=*pv1 -- always OK */
4571            v2=0;                            /* assume in padding */
4572            if (pv2>=var2) v2=*pv2;          /* in range */
4573            if (*pv1!=v2) break;             /* no longer the same */
4574            if (pv1==var1) break;            /* done; leave pv1 as is */
4575            }
4576          /* here when all inspected or a difference seen */
4577          if (*pv1<v2) break;                /* var1 too low to subtract */
4578          if (*pv1==v2) {                    /* var1 == var2 */
4579            /* reach here if var1 and var2 are identical; subtraction */
4580            /* would increase digit by one, and the residue will be 0 so */
4581            /* the calculation is done; leave the loop with residue=0. */
4582            thisunit++;                      /* as though subtracted */
4583            *var1=0;                         /* set var1 to 0 */
4584            var1units=1;                     /* .. */
4585            break;  /* from inner */
4586            } /* var1 == var2 */
4587          /* *pv1>v2.  Prepare for real subtraction; the lengths are equal */
4588          /* Estimate the multiplier (there's always a msu1-1)... */
4589          /* Bring in two units of var2 to provide a good estimate. */
4590          mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2pair);
4591          } /* lengths the same */
4592         else { /* var1units > var2ulen, so subtraction is safe */
4593          /* The var2 msu is one unit towards the lsu of the var1 msu, */
4594          /* so only one unit for var2 can be used. */
4595          mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2plus);
4596          }
4597        if (mult==0) mult=1;                 /* must always be at least 1 */
4598        /* subtraction needed; var1 is > var2 */
4599        thisunit=(Unit)(thisunit+mult);      /* accumulate */
4600        /* subtract var1-var2, into var1; only the overlap needs */
4601        /* processing, as this is an in-place calculation */
4602        shift=var2ulen-var2units;
4603        #if DECTRACE
4604          decDumpAr('1', &var1[shift], var1units-shift);
4605          decDumpAr('2', var2, var2units);
4606          printf("m=%ld\n", -mult);
4607        #endif
4608        decUnitAddSub(&var1[shift], var1units-shift,
4609                      var2, var2units, 0,
4610                      &var1[shift], -mult);
4611        #if DECTRACE
4612          decDumpAr('#', &var1[shift], var1units-shift);
4613        #endif
4614        /* var1 now probably has leading zeros; these are removed at the */
4615        /* top of the inner loop. */
4616        } /* inner loop */
4617
4618      /* The next unit has been calculated in full; unless it's a */
4619      /* leading zero, add to acc */
4620      if (accunits!=0 || thisunit!=0) {      /* is first or non-zero */
4621        *accnext=thisunit;                   /* store in accumulator */
4622        /* account exactly for the new digits */
4623        if (accunits==0) {
4624          accdigits++;                       /* at least one */
4625          for (pow=&powers[1]; thisunit>=*pow; pow++) accdigits++;
4626          }
4627         else accdigits+=DECDPUN;
4628        accunits++;                          /* update count */
4629        accnext--;                           /* ready for next */
4630        if (accdigits>reqdigits) break;      /* have enough digits */
4631        }
4632
4633      /* if the residue is zero, the operation is done (unless divide */
4634      /* or divideInteger and still not enough digits yet) */
4635      if (*var1==0 && var1units==1) {        /* residue is 0 */
4636        if (op&(REMAINDER|REMNEAR)) break;
4637        if ((op&DIVIDE) && (exponent<=maxexponent)) break;
4638        /* [drop through if divideInteger] */
4639        }
4640      /* also done enough if calculating remainder or integer */
4641      /* divide and just did the last ('units') unit */
4642      if (exponent==0 && !(op&DIVIDE)) break;
4643
4644      /* to get here, var1 is less than var2, so divide var2 by the per- */
4645      /* Unit power of ten and go for the next digit */
4646      var2ulen--;                            /* shift down */
4647      exponent-=DECDPUN;                     /* update the exponent */
4648      } /* outer loop */
4649
4650    /* ---- division is complete --------------------------------------- */
4651    /* here: acc      has at least reqdigits+1 of good results (or fewer */
4652    /*                if early stop), starting at accnext+1 (its lsu) */
4653    /*       var1     has any residue at the stopping point */
4654    /*       accunits is the number of digits collected in acc */
4655    if (accunits==0) {             /* acc is 0 */
4656      accunits=1;                  /* show have a unit .. */
4657      accdigits=1;                 /* .. */
4658      *accnext=0;                  /* .. whose value is 0 */
4659      }
4660     else accnext++;               /* back to last placed */
4661    /* accnext now -> lowest unit of result */
4662
4663    residue=0;                     /* assume no residue */
4664    if (op&DIVIDE) {
4665      /* record the presence of any residue, for rounding */
4666      if (*var1!=0 || var1units>1) residue=1;
4667       else { /* no residue */
4668        /* Had an exact division; clean up spurious trailing 0s. */
4669        /* There will be at most DECDPUN-1, from the final multiply, */
4670        /* and then only if the result is non-0 (and even) and the */
4671        /* exponent is 'loose'. */
4672        #if DECDPUN>1
4673        Unit lsu=*accnext;
4674        if (!(lsu&0x01) && (lsu!=0)) {
4675          /* count the trailing zeros */
4676          Int drop=0;
4677          for (;; drop++) {    /* [will terminate because lsu!=0] */
4678            if (exponent>=maxexponent) break;     /* don't chop real 0s */
4679            #if DECDPUN<=4
4680              if ((lsu-QUOT10(lsu, drop+1)
4681                  *powers[drop+1])!=0) break;     /* found non-0 digit */
4682            #else
4683              if (lsu%powers[drop+1]!=0) break;   /* found non-0 digit */
4684            #endif
4685            exponent++;
4686            }
4687          if (drop>0) {
4688            accunits=decShiftToLeast(accnext, accunits, drop);
4689            accdigits=decGetDigits(accnext, accunits);
4690            accunits=D2U(accdigits);
4691            /* [exponent was adjusted in the loop] */
4692            }
4693          } /* neither odd nor 0 */
4694        #endif
4695        } /* exact divide */
4696      } /* divide */
4697     else /* op!=DIVIDE */ {
4698      /* check for coefficient overflow */
4699      if (accdigits+exponent>reqdigits) {
4700        *status|=DEC_Division_impossible;
4701        break;
4702        }
4703      if (op & (REMAINDER|REMNEAR)) {
4704        /* [Here, the exponent will be 0, because var1 was adjusted */
4705        /* appropriately.] */
4706        Int postshift;                       /* work */
4707        Flag wasodd=0;                       /* integer was odd */
4708        Unit *quotlsu;                       /* for save */
4709        Int  quotdigits;                     /* .. */
4710
4711        bits=lhs->bits;                      /* remainder sign is always as lhs */
4712
4713        /* Fastpath when residue is truly 0 is worthwhile [and */
4714        /* simplifies the code below] */
4715        if (*var1==0 && var1units==1) {      /* residue is 0 */
4716          Int exp=lhs->exponent;             /* save min(exponents) */
4717          if (rhs->exponent<exp) exp=rhs->exponent;
4718          decNumberZero(res);                /* 0 coefficient */
4719          #if DECSUBSET
4720          if (set->extended)
4721          #endif
4722          res->exponent=exp;                 /* .. with proper exponent */
4723          res->bits=(uByte)(bits&DECNEG);          /* [cleaned] */
4724          decFinish(res, set, &residue, status);   /* might clamp */
4725          break;
4726          }
4727        /* note if the quotient was odd */
4728        if (*accnext & 0x01) wasodd=1;       /* acc is odd */
4729        quotlsu=accnext;                     /* save in case need to reinspect */
4730        quotdigits=accdigits;                /* .. */
4731
4732        /* treat the residue, in var1, as the value to return, via acc */
4733        /* calculate the unused zero digits.  This is the smaller of: */
4734        /*   var1 initial padding (saved above) */
4735        /*   var2 residual padding, which happens to be given by: */
4736        postshift=var1initpad+exponent-lhs->exponent+rhs->exponent;
4737        /* [the 'exponent' term accounts for the shifts during divide] */
4738        if (var1initpad<postshift) postshift=var1initpad;
4739
4740        /* shift var1 the requested amount, and adjust its digits */
4741        var1units=decShiftToLeast(var1, var1units, postshift);
4742        accnext=var1;
4743        accdigits=decGetDigits(var1, var1units);
4744        accunits=D2U(accdigits);
4745
4746        exponent=lhs->exponent;         /* exponent is smaller of lhs & rhs */
4747        if (rhs->exponent<exponent) exponent=rhs->exponent;
4748
4749        /* Now correct the result if doing remainderNear; if it */
4750        /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
4751        /* the integer was odd then the result should be rem-rhs. */
4752        if (op&REMNEAR) {
4753          Int compare, tarunits;        /* work */
4754          Unit *up;                     /* .. */
4755          /* calculate remainder*2 into the var1 buffer (which has */
4756          /* 'headroom' of an extra unit and hence enough space) */
4757          /* [a dedicated 'double' loop would be faster, here] */
4758          tarunits=decUnitAddSub(accnext, accunits, accnext, accunits,
4759                                 0, accnext, 1);
4760          /* decDumpAr('r', accnext, tarunits); */
4761
4762          /* Here, accnext (var1) holds tarunits Units with twice the */
4763          /* remainder's coefficient, which must now be compared to the */
4764          /* RHS.  The remainder's exponent may be smaller than the RHS's. */
4765          compare=decUnitCompare(accnext, tarunits, rhs->lsu, D2U(rhs->digits),
4766                                 rhs->exponent-exponent);
4767          if (compare==BADINT) {             /* deep trouble */
4768            *status|=DEC_Insufficient_storage;
4769            break;}
4770
4771          /* now restore the remainder by dividing by two; the lsu */
4772          /* is known to be even. */
4773          for (up=accnext; up<accnext+tarunits; up++) {
4774            Int half;              /* half to add to lower unit */
4775            half=*up & 0x01;
4776            *up/=2;                /* [shift] */
4777            if (!half) continue;
4778            *(up-1)+=DIV_ROUND_UP(DECDPUNMAX, 2);
4779            }
4780          /* [accunits still describes the original remainder length] */
4781
4782          if (compare>0 || (compare==0 && wasodd)) { /* adjustment needed */
4783            Int exp, expunits, exprem;       /* work */
4784            /* This is effectively causing round-up of the quotient, */
4785            /* so if it was the rare case where it was full and all */
4786            /* nines, it would overflow and hence division-impossible */
4787            /* should be raised */
4788            Flag allnines=0;                 /* 1 if quotient all nines */
4789            if (quotdigits==reqdigits) {     /* could be borderline */
4790              for (up=quotlsu; ; up++) {
4791                if (quotdigits>DECDPUN) {
4792                  if (*up!=DECDPUNMAX) break;/* non-nines */
4793                  }
4794                 else {                      /* this is the last Unit */
4795                  if (*up==powers[quotdigits]-1) allnines=1;
4796                  break;
4797                  }
4798                quotdigits-=DECDPUN;         /* checked those digits */
4799                } /* up */
4800              } /* borderline check */
4801            if (allnines) {
4802              *status|=DEC_Division_impossible;
4803              break;}
4804
4805            /* rem-rhs is needed; the sign will invert.  Again, var1 */
4806            /* can safely be used for the working Units array. */
4807            exp=rhs->exponent-exponent;      /* RHS padding needed */
4808            /* Calculate units and remainder from exponent. */
4809            expunits=exp/DECDPUN;
4810            exprem=exp%DECDPUN;
4811            /* subtract [A+B*(-m)]; the result will always be negative */
4812            accunits=-decUnitAddSub(accnext, accunits,
4813                                    rhs->lsu, D2U(rhs->digits),
4814                                    expunits, accnext, -(Int)powers[exprem]);
4815            accdigits=decGetDigits(accnext, accunits); /* count digits exactly */
4816            accunits=D2U(accdigits);    /* and recalculate the units for copy */
4817            /* [exponent is as for original remainder] */
4818            bits^=DECNEG;               /* flip the sign */
4819            }
4820          } /* REMNEAR */
4821        } /* REMAINDER or REMNEAR */
4822      } /* not DIVIDE */
4823
4824    /* Set exponent and bits */
4825    res->exponent=exponent;
4826    res->bits=(uByte)(bits&DECNEG);          /* [cleaned] */
4827
4828    /* Now the coefficient. */
4829    decSetCoeff(res, set, accnext, accdigits, &residue, status);
4830
4831    decFinish(res, set, &residue, status);   /* final cleanup */
4832
4833    #if DECSUBSET
4834    /* If a divide then strip trailing zeros if subset [after round] */
4835    if (!set->extended && (op==DIVIDE)) decTrim(res, set, 0, &dropped);
4836    #endif
4837    } while(0);                              /* end protected */
4838
4839  if (varalloc!=NULL) free(varalloc);   /* drop any storage used */
4840  if (allocacc!=NULL) free(allocacc);   /* .. */
4841  #if DECSUBSET
4842  if (allocrhs!=NULL) free(allocrhs);   /* .. */
4843  if (alloclhs!=NULL) free(alloclhs);   /* .. */
4844  #endif
4845  return res;
4846  } /* decDivideOp */
4847
4848/* ------------------------------------------------------------------ */
4849/* decMultiplyOp -- multiplication operation                          */
4850/*                                                                    */
4851/*  This routine performs the multiplication C=A x B.                 */
4852/*                                                                    */
4853/*   res is C, the result.  C may be A and/or B (e.g., X=X*X)         */
4854/*   lhs is A                                                         */
4855/*   rhs is B                                                         */
4856/*   set is the context                                               */
4857/*   status is the usual accumulator                                  */
4858/*                                                                    */
4859/* C must have space for set->digits digits.                          */
4860/*                                                                    */
4861/* ------------------------------------------------------------------ */
4862/* 'Classic' multiplication is used rather than Karatsuba, as the     */
4863/* latter would give only a minor improvement for the short numbers   */
4864/* expected to be handled most (and uses much more memory).           */
4865/*                                                                    */
4866/* There are two major paths here: the general-purpose ('old code')   */
4867/* path which handles all DECDPUN values, and a fastpath version      */
4868/* which is used if 64-bit ints are available, DECDPUN<=4, and more   */
4869/* than two calls to decUnitAddSub would be made.                     */
4870/*                                                                    */
4871/* The fastpath version lumps units together into 8-digit or 9-digit  */
4872/* chunks, and also uses a lazy carry strategy to minimise expensive  */
4873/* 64-bit divisions.  The chunks are then broken apart again into     */
4874/* units for continuing processing.  Despite this overhead, the       */
4875/* fastpath can speed up some 16-digit operations by 10x (and much    */
4876/* more for higher-precision calculations).                           */
4877/*                                                                    */
4878/* A buffer always has to be used for the accumulator; in the         */
4879/* fastpath, buffers are also always needed for the chunked copies of */
4880/* of the operand coefficients.                                       */
4881/* Static buffers are larger than needed just for multiply, to allow  */
4882/* for calls from other operations (notably exp).                     */
4883/* ------------------------------------------------------------------ */
4884#define FASTMUL (DECUSE64 && DECDPUN<5)
4885static decNumber * decMultiplyOp(decNumber *res, const decNumber *lhs,
4886                                 const decNumber *rhs, decContext *set,
4887                                 uInt *status) {
4888  Int    accunits;                 /* Units of accumulator in use */
4889  Int    exponent;                 /* work */
4890  Int    residue=0;                /* rounding residue */
4891  uByte  bits;                     /* result sign */
4892  Unit  *acc;                      /* -> accumulator Unit array */
4893  Int    needbytes;                /* size calculator */
4894  void  *allocacc=NULL;            /* -> allocated accumulator, iff allocated */
4895  Unit  accbuff[SD2U(DECBUFFER*4+1)]; /* buffer (+1 for DECBUFFER==0, */
4896                                   /* *4 for calls from other operations) */
4897  const Unit *mer, *mermsup;       /* work */
4898  Int   madlength;                 /* Units in multiplicand */
4899  Int   shift;                     /* Units to shift multiplicand by */
4900
4901  #if FASTMUL
4902    /* if DECDPUN is 1 or 3 work in base 10**9, otherwise */
4903    /* (DECDPUN is 2 or 4) then work in base 10**8 */
4904    #if DECDPUN & 1                /* odd */
4905      #define FASTBASE 1000000000  /* base */
4906      #define FASTDIGS          9  /* digits in base */
4907      #define FASTLAZY         18  /* carry resolution point [1->18] */
4908    #else
4909      #define FASTBASE  100000000
4910      #define FASTDIGS          8
4911      #define FASTLAZY       1844  /* carry resolution point [1->1844] */
4912    #endif
4913    /* three buffers are used, two for chunked copies of the operands */
4914    /* (base 10**8 or base 10**9) and one base 2**64 accumulator with */
4915    /* lazy carry evaluation */
4916    uInt   zlhibuff[(DECBUFFER*2+1)/8+1]; /* buffer (+1 for DECBUFFER==0) */
4917    uInt  *zlhi=zlhibuff;                 /* -> lhs array */
4918    uInt  *alloclhi=NULL;                 /* -> allocated buffer, iff allocated */
4919    uInt   zrhibuff[(DECBUFFER*2+1)/8+1]; /* buffer (+1 for DECBUFFER==0) */
4920    uInt  *zrhi=zrhibuff;                 /* -> rhs array */
4921    uInt  *allocrhi=NULL;                 /* -> allocated buffer, iff allocated */
4922    uLong  zaccbuff[(DECBUFFER*2+1)/4+2]; /* buffer (+1 for DECBUFFER==0) */
4923    /* [allocacc is shared for both paths, as only one will run] */
4924    uLong *zacc=zaccbuff;          /* -> accumulator array for exact result */
4925    #if DECDPUN==1
4926    Int    zoff;                   /* accumulator offset */
4927    #endif
4928    uInt  *lip, *rip;              /* item pointers */
4929    uInt  *lmsi, *rmsi;            /* most significant items */
4930    Int    ilhs, irhs, iacc;       /* item counts in the arrays */
4931    Int    lazy;                   /* lazy carry counter */
4932    uLong  lcarry;                 /* uLong carry */
4933    uInt   carry;                  /* carry (NB not uLong) */
4934    Int    count;                  /* work */
4935    const  Unit *cup;              /* .. */
4936    Unit  *up;                     /* .. */
4937    uLong *lp;                     /* .. */
4938    Int    p;                      /* .. */
4939  #endif
4940
4941  #if DECSUBSET
4942    decNumber *alloclhs=NULL;      /* -> allocated buffer, iff allocated */
4943    decNumber *allocrhs=NULL;      /* -> allocated buffer, iff allocated */
4944  #endif
4945
4946  #if DECCHECK
4947  if (decCheckOperands(