Merge remote-tracking branch 'remotes/mcayland/tags/qemu-openbios-20200707' into...

[qemu.git] / libdecnumber / decNumber.c

/* Decimal number arithmetic module for the decNumber C Library.
   Copyright (C) 2005, 2007 Free Software Foundation, Inc.
   Contributed by IBM Corporation.  Author Mike Cowlishaw.

   This file is part of GCC.

   GCC is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free
   Software Foundation; either version 2, or (at your option) any later
   version.

   In addition to the permissions in the GNU General Public License,
   the Free Software Foundation gives you unlimited permission to link
   the compiled version of this file into combinations with other
   programs, and to distribute those combinations without any
   restriction coming from the use of this file.  (The General Public
   License restrictions do apply in other respects; for example, they
   cover modification of the file, and distribution when not linked
   into a combine executable.)

   GCC is distributed in the hope that it will be useful, but WITHOUT ANY
   WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public License
   along with GCC; see the file COPYING.  If not, write to the Free
   Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301, USA.  */

/* ------------------------------------------------------------------ */
/* Decimal Number arithmetic module                                   */
/* ------------------------------------------------------------------ */
/* This module comprises the routines for General Decimal Arithmetic  */
/* as defined in the specification which may be found on the          */
/* http://www2.hursley.ibm.com/decimal web pages.  It implements both */
/* the full ('extended') arithmetic and the simpler ('subset')        */
/* arithmetic.                                                        */
/*                                                                    */
/* Usage notes:                                                       */
/*                                                                    */
/* 1. This code is ANSI C89 except:                                   */
/*                                                                    */
/*       If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and       */
/*       uint64_t types may be used.  To avoid these, set DECUSE64=0  */
/*       and DECDPUN<=4 (see documentation).                          */
/*                                                                    */
/* 2. The decNumber format which this library uses is optimized for   */
/*    efficient processing of relatively short numbers; in particular */
/*    it allows the use of fixed sized structures and minimizes copy  */
/*    and move operations.  It does, however, support arbitrary       */
/*    precision (up to 999,999,999 digits) and arbitrary exponent     */
/*    range (Emax in the range 0 through 999,999,999 and Emin in the  */
/*    range -999,999,999 through 0).  Mathematical functions (for     */
/*    example decNumberExp) as identified below are restricted more   */
/*    tightly: digits, emax, and -emin in the context must be <=      */
/*    DEC_MAX_MATH (999999), and their operand(s) must be within      */
/*    these bounds.                                                   */
/*                                                                    */
/* 3. Logical functions are further restricted; their operands must   */
/*    be finite, positive, have an exponent of zero, and all digits   */
/*    must be either 0 or 1.  The result will only contain digits     */
/*    which are 0 or 1 (and will have exponent=0 and a sign of 0).    */
/*                                                                    */
/* 4. Operands to operator functions are never modified unless they   */
/*    are also specified to be the result number (which is always     */
/*    permitted).  Other than that case, operands must not overlap.   */
/*                                                                    */
/* 5. Error handling: the type of the error is ORed into the status   */
/*    flags in the current context (decContext structure).  The       */
/*    SIGFPE signal is then raised if the corresponding trap-enabler  */
/*    flag in the decContext is set (is 1).                           */
/*                                                                    */
/*    It is the responsibility of the caller to clear the status      */
/*    flags as required.                                              */
/*                                                                    */
/*    The result of any routine which returns a number will always    */
/*    be a valid number (which may be a special value, such as an     */
/*    Infinity or NaN).                                               */
/*                                                                    */
/* 6. The decNumber format is not an exchangeable concrete            */
/*    representation as it comprises fields which may be machine-     */
/*    dependent (packed or unpacked, or special length, for example). */
/*    Canonical conversions to and from strings are provided; other   */
/*    conversions are available in separate modules.                  */
/*                                                                    */
/* 7. Normally, input operands are assumed to be valid.  Set DECCHECK */
/*    to 1 for extended operand checking (including NULL operands).   */
/*    Results are undefined if a badly-formed structure (or a NULL    */
/*    pointer to a structure) is provided, though with DECCHECK       */
/*    enabled the operator routines are protected against exceptions. */
/*    (Except if the result pointer is NULL, which is unrecoverable.) */
/*                                                                    */
/*    However, the routines will never cause exceptions if they are   */
/*    given well-formed operands, even if the value of the operands   */
/*    is inappropriate for the operation and DECCHECK is not set.     */
/*    (Except for SIGFPE, as and where documented.)                   */
/*                                                                    */
/* 8. Subset arithmetic is available only if DECSUBSET is set to 1.   */
/* ------------------------------------------------------------------ */
/* Implementation notes for maintenance of this module:               */
/*                                                                    */
/* 1. Storage leak protection:  Routines which use malloc are not     */
/*    permitted to use return for fastpath or error exits (i.e.,      */
/*    they follow strict structured programming conventions).         */
/*    Instead they have a do{}while(0); construct surrounding the     */
/*    code which is protected -- break may be used to exit this.      */
/*    Other routines can safely use the return statement inline.      */
/*                                                                    */
/*    Storage leak accounting can be enabled using DECALLOC.          */
/*                                                                    */
/* 2. All loops use the for(;;) construct.  Any do construct does     */
/*    not loop; it is for allocation protection as just described.    */
/*                                                                    */
/* 3. Setting status in the context must always be the very last      */
/*    action in a routine, as non-0 status may raise a trap and hence */
/*    the call to set status may not return (if the handler uses long */
/*    jump).  Therefore all cleanup must be done first.  In general,  */
/*    to achieve this status is accumulated and is only applied just  */
/*    before return by calling decContextSetStatus (via decStatus).   */
/*                                                                    */
/*    Routines which allocate storage cannot, in general, use the     */
/*    'top level' routines which could cause a non-returning          */
/*    transfer of control.  The decXxxxOp routines are safe (do not   */
/*    call decStatus even if traps are set in the context) and should */
/*    be used instead (they are also a little faster).                */
/*                                                                    */
/* 4. Exponent checking is minimized by allowing the exponent to      */
/*    grow outside its limits during calculations, provided that      */
/*    the decFinalize function is called later.  Multiplication and   */
/*    division, and intermediate calculations in exponentiation,      */
/*    require more careful checks because of the risk of 31-bit       */
/*    overflow (the most negative valid exponent is -1999999997, for  */
/*    a 999999999-digit number with adjusted exponent of -999999999). */
/*                                                                    */
/* 5. Rounding is deferred until finalization of results, with any    */
/*    'off to the right' data being represented as a single digit     */
/*    residue (in the range -1 through 9).  This avoids any double-   */
/*    rounding when more than one shortening takes place (for         */
/*    example, when a result is subnormal).                           */
/*                                                                    */
/* 6. The digits count is allowed to rise to a multiple of DECDPUN    */
/*    during many operations, so whole Units are handled and exact    */
/*    accounting of digits is not needed.  The correct digits value   */
/*    is found by decGetDigits, which accounts for leading zeros.     */
/*    This must be called before any rounding if the number of digits */
/*    is not known exactly.                                           */
/*                                                                    */
/* 7. The multiply-by-reciprocal 'trick' is used for partitioning     */
/*    numbers up to four digits, using appropriate constants.  This   */
/*    is not useful for longer numbers because overflow of 32 bits    */
/*    would lead to 4 multiplies, which is almost as expensive as     */
/*    a divide (unless a floating-point or 64-bit multiply is         */
/*    assumed to be available).                                       */
/*                                                                    */
/* 8. Unusual abbreviations that may be used in the commentary:       */
/*      lhs -- left hand side (operand, of an operation)              */
/*      lsd -- least significant digit (of coefficient)               */
/*      lsu -- least significant Unit (of coefficient)                */
/*      msd -- most significant digit (of coefficient)                */
/*      msi -- most significant item (in an array)                    */
/*      msu -- most significant Unit (of coefficient)                 */
/*      rhs -- right hand side (operand, of an operation)             */
/*      +ve -- positive                                               */
/*      -ve -- negative                                               */
/*      **  -- raise to the power                                     */
/* ------------------------------------------------------------------ */

#include "qemu/osdep.h"
#include "libdecnumber/dconfig.h"
#include "libdecnumber/decNumber.h"
#include "libdecnumber/decNumberLocal.h"

/* Constants */
/* Public lookup table used by the D2U macro */
const uByte d2utable[DECMAXD2U+1]=D2UTABLE;

#define DECVERB     1              /* set to 1 for verbose DECCHECK */
#define powers      DECPOWERS      /* old internal name */

/* Local constants */
#define DIVIDE      0x80           /* Divide operators */
#define REMAINDER   0x40           /* .. */
#define DIVIDEINT   0x20           /* .. */
#define REMNEAR     0x10           /* .. */
#define COMPARE     0x01           /* Compare operators */
#define COMPMAX     0x02           /* .. */
#define COMPMIN     0x03           /* .. */
#define COMPTOTAL   0x04           /* .. */
#define COMPNAN     0x05           /* .. [NaN processing] */
#define COMPSIG     0x06           /* .. [signaling COMPARE] */
#define COMPMAXMAG  0x07           /* .. */
#define COMPMINMAG  0x08           /* .. */

#define DEC_sNaN     0x40000000    /* local status: sNaN signal */
#define BADINT  (Int)0x80000000    /* most-negative Int; error indicator */
/* Next two indicate an integer >= 10**6, and its parity (bottom bit) */
#define BIGEVEN (Int)0x80000002
#define BIGODD  (Int)0x80000003

static Unit uarrone[1]={1};   /* Unit array of 1, used for incrementing */

/* Granularity-dependent code */
#if DECDPUN<=4
  #define eInt  Int           /* extended integer */
  #define ueInt uInt          /* unsigned extended integer */
  /* Constant multipliers for divide-by-power-of five using reciprocal */
  /* multiply, after removing powers of 2 by shifting, and final shift */
  /* of 17 [we only need up to **4] */
  static const uInt multies[]={131073, 26215, 5243, 1049, 210};
  /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
  #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
#else
  /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
  #if !DECUSE64
    #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4
  #endif
  #define eInt  Long          /* extended integer */
  #define ueInt uLong         /* unsigned extended integer */
#endif

/* Local routines */
static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *,
                              decContext *, uByte, uInt *);
static Flag        decBiStr(const char *, const char *, const char *);
static uInt        decCheckMath(const decNumber *, decContext *, uInt *);
static void        decApplyRound(decNumber *, decContext *, Int, uInt *);
static Int         decCompare(const decNumber *lhs, const decNumber *rhs, Flag);
static decNumber * decCompareOp(decNumber *, const decNumber *,
                              const decNumber *, decContext *,
                              Flag, uInt *);
static void        decCopyFit(decNumber *, const decNumber *, decContext *,
                              Int *, uInt *);
static decNumber * decDecap(decNumber *, Int);
static decNumber * decDivideOp(decNumber *, const decNumber *,
                              const decNumber *, decContext *, Flag, uInt *);
static decNumber * decExpOp(decNumber *, const decNumber *,
                              decContext *, uInt *);
static void        decFinalize(decNumber *, decContext *, Int *, uInt *);
static Int         decGetDigits(Unit *, Int);
static Int         decGetInt(const decNumber *);
static decNumber * decLnOp(decNumber *, const decNumber *,
                              decContext *, uInt *);
static decNumber * decMultiplyOp(decNumber *, const decNumber *,
                              const decNumber *, decContext *,
                              uInt *);
static decNumber * decNaNs(decNumber *, const decNumber *,
                              const decNumber *, decContext *, uInt *);
static decNumber * decQuantizeOp(decNumber *, const decNumber *,
                              const decNumber *, decContext *, Flag,
                              uInt *);
static void        decReverse(Unit *, Unit *);
static void        decSetCoeff(decNumber *, decContext *, const Unit *,
                              Int, Int *, uInt *);
static void        decSetMaxValue(decNumber *, decContext *);
static void        decSetOverflow(decNumber *, decContext *, uInt *);
static void        decSetSubnormal(decNumber *, decContext *, Int *, uInt *);
static Int         decShiftToLeast(Unit *, Int, Int);
static Int         decShiftToMost(Unit *, Int, Int);
static void        decStatus(decNumber *, uInt, decContext *);
static void        decToString(const decNumber *, char[], Flag);
static decNumber * decTrim(decNumber *, decContext *, Flag, Int *);
static Int         decUnitAddSub(const Unit *, Int, const Unit *, Int, Int,
                              Unit *, Int);
static Int         decUnitCompare(const Unit *, Int, const Unit *, Int, Int);

#if !DECSUBSET
/* decFinish == decFinalize when no subset arithmetic needed */
#define decFinish(a,b,c,d) decFinalize(a,b,c,d)
#else
static void        decFinish(decNumber *, decContext *, Int *, uInt *);
static decNumber * decRoundOperand(const decNumber *, decContext *, uInt *);
#endif

/* Local macros */
/* masked special-values bits */
#define SPECIALARG  (rhs->bits & DECSPECIAL)
#define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)

/* Diagnostic macros, etc. */
#if DECALLOC
/* Handle malloc/free accounting.  If enabled, our accountable routines */
/* are used; otherwise the code just goes straight to the system malloc */
/* and free routines. */
#define malloc(a) decMalloc(a)
#define free(a) decFree(a)
#define DECFENCE 0x5a              /* corruption detector */
/* 'Our' malloc and free: */
static void *decMalloc(size_t);
static void  decFree(void *);
uInt decAllocBytes=0;              /* count of bytes allocated */
/* Note that DECALLOC code only checks for storage buffer overflow. */
/* To check for memory leaks, the decAllocBytes variable must be */
/* checked to be 0 at appropriate times (e.g., after the test */
/* harness completes a set of tests).  This checking may be unreliable */
/* if the testing is done in a multi-thread environment. */
#endif

#if DECCHECK
/* Optional checking routines.  Enabling these means that decNumber */
/* and decContext operands to operator routines are checked for */
/* correctness.  This roughly doubles the execution time of the */
/* fastest routines (and adds 600+ bytes), so should not normally be */
/* used in 'production'. */
/* decCheckInexact is used to check that inexact results have a full */
/* complement of digits (where appropriate -- this is not the case */
/* for Quantize, for example) */
#define DECUNRESU ((decNumber *)(void *)0xffffffff)
#define DECUNUSED ((const decNumber *)(void *)0xffffffff)
#define DECUNCONT ((decContext *)(void *)(0xffffffff))
static Flag decCheckOperands(decNumber *, const decNumber *,
                             const decNumber *, decContext *);
static Flag decCheckNumber(const decNumber *);
static void decCheckInexact(const decNumber *, decContext *);
#endif

#if DECTRACE || DECCHECK
/* Optional trace/debugging routines (may or may not be used) */
void decNumberShow(const decNumber *);  /* displays the components of a number */
static void decDumpAr(char, const Unit *, Int);
#endif

/* ================================================================== */
/* Conversions                                                        */
/* ================================================================== */

/* ------------------------------------------------------------------ */
/* from-int32 -- conversion from Int or uInt                          */
/*                                                                    */
/*  dn is the decNumber to receive the integer                        */
/*  in or uin is the integer to be converted                          */
/*  returns dn                                                        */
/*                                                                    */
/* No error is possible.                                              */
/* ------------------------------------------------------------------ */
decNumber * decNumberFromInt32(decNumber *dn, Int in) {
  uInt unsig;
  if (in>=0) unsig=in;
   else {                               /* negative (possibly BADINT) */
    if (in==BADINT) unsig=(uInt)1073741824*2; /* special case */
     else unsig=-in;                    /* invert */
    }
  /* in is now positive */
  decNumberFromUInt32(dn, unsig);
  if (in<0) dn->bits=DECNEG;            /* sign needed */
  return dn;
  } /* decNumberFromInt32 */

decNumber * decNumberFromUInt32(decNumber *dn, uInt uin) {
  Unit *up;                             /* work pointer */
  decNumberZero(dn);                    /* clean */
  if (uin==0) return dn;                /* [or decGetDigits bad call] */
  for (up=dn->lsu; uin>0; up++) {
    *up=(Unit)(uin%(DECDPUNMAX+1));
    uin=uin/(DECDPUNMAX+1);
    }
  dn->digits=decGetDigits(dn->lsu, up-dn->lsu);
  return dn;
  } /* decNumberFromUInt32 */

/* ------------------------------------------------------------------ */
/* to-int32 -- conversion to Int or uInt                              */
/*                                                                    */
/*  dn is the decNumber to convert                                    */
/*  set is the context for reporting errors                           */
/*  returns the converted decNumber, or 0 if Invalid is set           */
/*                                                                    */
/* Invalid is set if the decNumber does not have exponent==0 or if    */
/* it is a NaN, Infinite, or out-of-range.                            */
/* ------------------------------------------------------------------ */
Int decNumberToInt32(const decNumber *dn, decContext *set) {
  #if DECCHECK
  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
  #endif

  /* special or too many digits, or bad exponent */
  if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0) ; /* bad */
   else { /* is a finite integer with 10 or fewer digits */
    Int d;                         /* work */
    const Unit *up;                /* .. */
    uInt hi=0, lo;                 /* .. */
    up=dn->lsu;                    /* -> lsu */
    lo=*up;                        /* get 1 to 9 digits */
    #if DECDPUN>1                  /* split to higher */
      hi=lo/10;
      lo=lo%10;
    #endif
    up++;
    /* collect remaining Units, if any, into hi */
    for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];
    /* now low has the lsd, hi the remainder */
    if (hi>214748364 || (hi==214748364 && lo>7)) { /* out of range? */
      /* most-negative is a reprieve */
      if (dn->bits&DECNEG && hi==214748364 && lo==8) return 0x80000000;
      /* bad -- drop through */
      }
     else { /* in-range always */
      Int i=X10(hi)+lo;
      if (dn->bits&DECNEG) return -i;
      return i;
      }
    } /* integer */
  decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
  return 0;
  } /* decNumberToInt32 */

uInt decNumberToUInt32(const decNumber *dn, decContext *set) {
  #if DECCHECK
  if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
  #endif
  /* special or too many digits, or bad exponent, or negative (<0) */
  if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0
    || (dn->bits&DECNEG && !ISZERO(dn)));                   /* bad */
   else { /* is a finite integer with 10 or fewer digits */
    Int d;                         /* work */
    const Unit *up;                /* .. */
    uInt hi=0, lo;                 /* .. */
    up=dn->lsu;                    /* -> lsu */
    lo=*up;                        /* get 1 to 9 digits */
    #if DECDPUN>1                  /* split to higher */
      hi=lo/10;
      lo=lo%10;
    #endif
    up++;
    /* collect remaining Units, if any, into hi */
    for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];

    /* now low has the lsd, hi the remainder */
    if (hi>429496729 || (hi==429496729 && lo>5)) ; /* no reprieve possible */
     else return X10(hi)+lo;
    } /* integer */
  decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
  return 0;
  } /* decNumberToUInt32 */

decNumber *decNumberFromInt64(decNumber *dn, int64_t in)
{
    uint64_t unsig = in;
    if (in < 0) {
        unsig = -unsig;
    }

    decNumberFromUInt64(dn, unsig);
    if (in < 0) {
        dn->bits = DECNEG;        /* sign needed */
    }
    return dn;
} /* decNumberFromInt64 */

decNumber *decNumberFromUInt64(decNumber *dn, uint64_t uin)
{
    Unit *up;                             /* work pointer */
    decNumberZero(dn);                    /* clean */
    if (uin == 0) {
        return dn;                /* [or decGetDigits bad call] */
    }
    for (up = dn->lsu; uin > 0; up++) {
        *up = (Unit)(uin % (DECDPUNMAX + 1));
        uin = uin / (DECDPUNMAX + 1);
    }
    dn->digits = decGetDigits(dn->lsu, up-dn->lsu);
    return dn;
} /* decNumberFromUInt64 */

/* ------------------------------------------------------------------ */
/* to-int64 -- conversion to int64                                    */
/*                                                                    */
/*  dn is the decNumber to convert.  dn is assumed to have been       */
/*    rounded to a floating point integer value.                      */
/*  set is the context for reporting errors                           */
/*  returns the converted decNumber, or 0 if Invalid is set           */
/*                                                                    */
/* Invalid is set if the decNumber is a NaN, Infinite or is out of    */
/* range for a signed 64 bit integer.                                 */
/* ------------------------------------------------------------------ */

int64_t decNumberIntegralToInt64(const decNumber *dn, decContext *set)
{
    if (decNumberIsSpecial(dn) || (dn->exponent < 0) ||
       (dn->digits + dn->exponent > 19)) {
        goto Invalid;
    } else {
        int64_t d;        /* work */
        const Unit *up;   /* .. */
        uint64_t hi = 0;
        up = dn->lsu;     /* -> lsu */

        for (d = 1; d <= dn->digits; up++, d += DECDPUN) {
            uint64_t prev = hi;
            hi += *up * powers[d-1];
            if ((hi < prev) || (hi > INT64_MAX)) {
                goto Invalid;
            }
        }

        uint64_t prev = hi;
        hi *= (uint64_t)powers[dn->exponent];
        if ((hi < prev) || (hi > INT64_MAX)) {
            goto Invalid;
        }
        return (decNumberIsNegative(dn)) ? -((int64_t)hi) : (int64_t)hi;
    }

Invalid:
    decContextSetStatus(set, DEC_Invalid_operation);
    return 0;
} /* decNumberIntegralToInt64 */


/* ------------------------------------------------------------------ */
/* to-scientific-string -- conversion to numeric string               */
/* to-engineering-string -- conversion to numeric string              */
/*                                                                    */
/*   decNumberToString(dn, string);                                   */
/*   decNumberToEngString(dn, string);                                */
/*                                                                    */
/*  dn is the decNumber to convert                                    */
/*  string is the string where the result will be laid out            */
/*                                                                    */
/*  string must be at least dn->digits+14 characters long             */
/*                                                                    */
/*  No error is possible, and no status can be set.                   */
/* ------------------------------------------------------------------ */
char * decNumberToString(const decNumber *dn, char *string){
  decToString(dn, string, 0);
  return string;
  } /* DecNumberToString */

char * decNumberToEngString(const decNumber *dn, char *string){
  decToString(dn, string, 1);
  return string;
  } /* DecNumberToEngString */

/* ------------------------------------------------------------------ */
/* to-number -- conversion from numeric string                        */
/*                                                                    */
/* decNumberFromString -- convert string to decNumber                 */
/*   dn        -- the number structure to fill                        */
/*   chars[]   -- the string to convert ('\0' terminated)             */
/*   set       -- the context used for processing any error,          */
/*                determining the maximum precision available         */
/*                (set.digits), determining the maximum and minimum   */
/*                exponent (set.emax and set.emin), determining if    */
/*                extended values are allowed, and checking the       */
/*                rounding mode if overflow occurs or rounding is     */
/*                needed.                                             */
/*                                                                    */
/* The length of the coefficient and the size of the exponent are     */
/* checked by this routine, so the correct error (Underflow or        */
/* Overflow) can be reported or rounding applied, as necessary.       */
/*                                                                    */
/* If bad syntax is detected, the result will be a quiet NaN.         */
/* ------------------------------------------------------------------ */
decNumber * decNumberFromString(decNumber *dn, const char chars[],
                                decContext *set) {
  Int   exponent=0;                /* working exponent [assume 0] */
  uByte bits=0;                    /* working flags [assume +ve] */
  Unit  *res;                      /* where result will be built */
  Unit  resbuff[SD2U(DECBUFFER+9)];/* local buffer in case need temporary */
                                   /* [+9 allows for ln() constants] */
  Unit  *allocres=NULL;            /* -> allocated result, iff allocated */
  Int   d=0;                       /* count of digits found in decimal part */
  const char *dotchar=NULL;        /* where dot was found */
  const char *cfirst=chars;        /* -> first character of decimal part */
  const char *last=NULL;           /* -> last digit of decimal part */
  const char *c;                   /* work */
  Unit  *up;                       /* .. */
  #if DECDPUN>1
  Int   cut, out;                  /* .. */
  #endif
  Int   residue;                   /* rounding residue */
  uInt  status=0;                  /* error code */

  #if DECCHECK
  if (decCheckOperands(DECUNRESU, DECUNUSED, DECUNUSED, set))
    return decNumberZero(dn);
  #endif

  do {                             /* status & malloc protection */
    for (c=chars;; c++) {          /* -> input character */
      if (*c>='0' && *c<='9') {    /* test for Arabic digit */
        last=c;
        d++;                       /* count of real digits */
        continue;                  /* still in decimal part */
        }
      if (*c=='.' && dotchar==NULL) { /* first '.' */
        dotchar=c;                 /* record offset into decimal part */
        if (c==cfirst) cfirst++;   /* first digit must follow */
        continue;}
      if (c==chars) {              /* first in string... */
        if (*c=='-') {             /* valid - sign */
          cfirst++;
          bits=DECNEG;
          continue;}
        if (*c=='+') {             /* valid + sign */
          cfirst++;
          continue;}
        }
      /* *c is not a digit, or a valid +, -, or '.' */
      break;
      } /* c */

    if (last==NULL) {              /* no digits yet */
      status=DEC_Conversion_syntax;/* assume the worst */
      if (*c=='\0') break;         /* and no more to come... */
      #if DECSUBSET
      /* if subset then infinities and NaNs are not allowed */
      if (!set->extended) break;   /* hopeless */
      #endif
      /* Infinities and NaNs are possible, here */
      if (dotchar!=NULL) break;    /* .. unless had a dot */
      decNumberZero(dn);           /* be optimistic */
      if (decBiStr(c, "infinity", "INFINITY")
       || decBiStr(c, "inf", "INF")) {
        dn->bits=bits | DECINF;
        status=0;                  /* is OK */
        break; /* all done */
        }
      /* a NaN expected */
      /* 2003.09.10 NaNs are now permitted to have a sign */
      dn->bits=bits | DECNAN;      /* assume simple NaN */
      if (*c=='s' || *c=='S') {    /* looks like an sNaN */
        c++;
        dn->bits=bits | DECSNAN;
        }
      if (*c!='n' && *c!='N') break;    /* check caseless "NaN" */
      c++;
      if (*c!='a' && *c!='A') break;    /* .. */
      c++;
      if (*c!='n' && *c!='N') break;    /* .. */
      c++;
      /* now either nothing, or nnnn payload, expected */
      /* -> start of integer and skip leading 0s [including plain 0] */
      for (cfirst=c; *cfirst=='0';) cfirst++;
      if (*cfirst=='\0') {         /* "NaN" or "sNaN", maybe with all 0s */
        status=0;                  /* it's good */
        break;                     /* .. */
        }
      /* something other than 0s; setup last and d as usual [no dots] */
      for (c=cfirst;; c++, d++) {
        if (*c<'0' || *c>'9') break; /* test for Arabic digit */
        last=c;
        }
      if (*c!='\0') break;         /* not all digits */
      if (d>set->digits-1) {
        /* [NB: payload in a decNumber can be full length unless */
        /* clamped, in which case can only be digits-1] */
        if (set->clamp) break;
        if (d>set->digits) break;
        } /* too many digits? */
      /* good; drop through to convert the integer to coefficient */
      status=0;                    /* syntax is OK */
      bits=dn->bits;               /* for copy-back */
      } /* last==NULL */

     else if (*c!='\0') {          /* more to process... */
      /* had some digits; exponent is only valid sequence now */
      Flag nege;                   /* 1=negative exponent */
      const char *firstexp;        /* -> first significant exponent digit */
      status=DEC_Conversion_syntax;/* assume the worst */
      if (*c!='e' && *c!='E') break;
      /* Found 'e' or 'E' -- now process explicit exponent */
      /* 1998.07.11: sign no longer required */
      nege=0;
      c++;                         /* to (possible) sign */
      if (*c=='-') {nege=1; c++;}
       else if (*c=='+') c++;
      if (*c=='\0') break;

      for (; *c=='0' && *(c+1)!='\0';) c++;  /* strip insignificant zeros */
      firstexp=c;                            /* save exponent digit place */
      for (; ;c++) {
        if (*c<'0' || *c>'9') break;         /* not a digit */
        exponent=X10(exponent)+(Int)*c-(Int)'0';
        } /* c */
      /* if not now on a '\0', *c must not be a digit */
      if (*c!='\0') break;

      /* (this next test must be after the syntax checks) */
      /* if it was too long the exponent may have wrapped, so check */
      /* carefully and set it to a certain overflow if wrap possible */
      if (c>=firstexp+9+1) {
        if (c>firstexp+9+1 || *firstexp>'1') exponent=DECNUMMAXE*2;
        /* [up to 1999999999 is OK, for example 1E-1000000998] */
        }
      if (nege) exponent=-exponent;     /* was negative */
      status=0;                         /* is OK */
      } /* stuff after digits */

    /* Here when whole string has been inspected; syntax is good */
    /* cfirst->first digit (never dot), last->last digit (ditto) */

    /* strip leading zeros/dot [leave final 0 if all 0's] */
    if (*cfirst=='0') {                 /* [cfirst has stepped over .] */
      for (c=cfirst; c<last; c++, cfirst++) {
        if (*c=='.') continue;          /* ignore dots */
        if (*c!='0') break;             /* non-zero found */
        d--;                            /* 0 stripped */
        } /* c */
      #if DECSUBSET
      /* make a rapid exit for easy zeros if !extended */
      if (*cfirst=='0' && !set->extended) {
        decNumberZero(dn);              /* clean result */
        break;                          /* [could be return] */
        }
      #endif
      } /* at least one leading 0 */

    /* Handle decimal point... */
    if (dotchar!=NULL && dotchar<last)  /* non-trailing '.' found? */
      exponent-=(last-dotchar);         /* adjust exponent */
    /* [we can now ignore the .] */

    /* OK, the digits string is good.  Assemble in the decNumber, or in */
    /* a temporary units array if rounding is needed */
    if (d<=set->digits) res=dn->lsu;    /* fits into supplied decNumber */
     else {                             /* rounding needed */
      Int needbytes=D2U(d)*sizeof(Unit);/* bytes needed */
      res=resbuff;                      /* assume use local buffer */
      if (needbytes>(Int)sizeof(resbuff)) { /* too big for local */
        allocres=(Unit *)malloc(needbytes);
        if (allocres==NULL) {status|=DEC_Insufficient_storage; break;}
        res=allocres;
        }
      }
    /* res now -> number lsu, buffer, or allocated storage for Unit array */

    /* Place the coefficient into the selected Unit array */
    /* [this is often 70% of the cost of this function when DECDPUN>1] */
    #if DECDPUN>1
    out=0;                         /* accumulator */
    up=res+D2U(d)-1;               /* -> msu */
    cut=d-(up-res)*DECDPUN;        /* digits in top unit */
    for (c=cfirst;; c++) {         /* along the digits */
      if (*c=='.') continue;       /* ignore '.' [don't decrement cut] */
      out=X10(out)+(Int)*c-(Int)'0';
      if (c==last) break;          /* done [never get to trailing '.'] */
      cut--;
      if (cut>0) continue;         /* more for this unit */
      *up=(Unit)out;               /* write unit */
      up--;                        /* prepare for unit below.. */
      cut=DECDPUN;                 /* .. */
      out=0;                       /* .. */
      } /* c */
    *up=(Unit)out;                 /* write lsu */

    #else
    /* DECDPUN==1 */
    up=res;                        /* -> lsu */
    for (c=last; c>=cfirst; c--) { /* over each character, from least */
      if (*c=='.') continue;       /* ignore . [don't step up] */
      *up=(Unit)((Int)*c-(Int)'0');
      up++;
      } /* c */
    #endif

    dn->bits=bits;
    dn->exponent=exponent;
    dn->digits=d;

    /* if not in number (too long) shorten into the number */
    if (d>set->digits) {
      residue=0;
      decSetCoeff(dn, set, res, d, &residue, &status);
      /* always check for overflow or subnormal and round as needed */
      decFinalize(dn, set, &residue, &status);
      }
     else { /* no rounding, but may still have overflow or subnormal */
      /* [these tests are just for performance; finalize repeats them] */
      if ((dn->exponent-1<set->emin-dn->digits)
       || (dn->exponent-1>set->emax-set->digits)) {
        residue=0;
        decFinalize(dn, set, &residue, &status);
        }
      }
    /* decNumberShow(dn); */
    } while(0);                         /* [for break] */

  if (allocres!=NULL) free(allocres);   /* drop any storage used */
  if (status!=0) decStatus(dn, status, set);
  return dn;
  } /* decNumberFromString */

/* ================================================================== */
/* Operators                                                          */
/* ================================================================== */

/* ------------------------------------------------------------------ */
/* decNumberAbs -- absolute value operator                            */
/*                                                                    */
/*   This computes C = abs(A)                                         */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* See also decNumberCopyAbs for a quiet bitwise version of this.     */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* This has the same effect as decNumberPlus unless A is negative,    */
/* in which case it has the same effect as decNumberMinus.            */
/* ------------------------------------------------------------------ */
decNumber * decNumberAbs(decNumber *res, const decNumber *rhs,
                         decContext *set) {
  decNumber dzero;                      /* for 0 */
  uInt status=0;                        /* accumulator */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  decNumberZero(&dzero);                /* set 0 */
  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
  decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberAbs */

/* ------------------------------------------------------------------ */
/* decNumberAdd -- add two Numbers                                    */
/*                                                                    */
/*   This computes C = A + B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* This just calls the routine shared with Subtract                   */
decNumber * decNumberAdd(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decAddOp(res, lhs, rhs, set, 0, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberAdd */

/* ------------------------------------------------------------------ */
/* decNumberAnd -- AND two Numbers, digitwise                         */
/*                                                                    */
/*   This computes C = A & B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X&X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context (used for result length and error report)     */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Logical function restrictions apply (see above); a NaN is          */
/* returned with Invalid_operation if a restriction is violated.      */
/* ------------------------------------------------------------------ */
decNumber * decNumberAnd(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  const Unit *ua, *ub;                  /* -> operands */
  const Unit *msua, *msub;              /* -> operand msus */
  Unit *uc,  *msuc;                     /* -> result and its msu */
  Int   msudigs;                        /* digits in res msu */
  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
   || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
    decStatus(res, DEC_Invalid_operation, set);
    return res;
    }

  /* operands are valid */
  ua=lhs->lsu;                          /* bottom-up */
  ub=rhs->lsu;                          /* .. */
  uc=res->lsu;                          /* .. */
  msua=ua+D2U(lhs->digits)-1;           /* -> msu of lhs */
  msub=ub+D2U(rhs->digits)-1;           /* -> msu of rhs */
  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
  for (; uc<=msuc; ua++, ub++, uc++) {  /* Unit loop */
    Unit a, b;                          /* extract units */
    if (ua>msua) a=0;
     else a=*ua;
    if (ub>msub) b=0;
     else b=*ub;
    *uc=0;                              /* can now write back */
    if (a|b) {                          /* maybe 1 bits to examine */
      Int i, j;
      *uc=0;                            /* can now write back */
      /* This loop could be unrolled and/or use BIN2BCD tables */
      for (i=0; i<DECDPUN; i++) {
        if (a&b&1) *uc=*uc+(Unit)powers[i];  /* effect AND */
        j=a%10;
        a=a/10;
        j|=b%10;
        b=b/10;
        if (j>1) {
          decStatus(res, DEC_Invalid_operation, set);
          return res;
          }
        if (uc==msuc && i==msudigs-1) break; /* just did final digit */
        } /* each digit */
      } /* both OK */
    } /* each unit */
  /* [here uc-1 is the msu of the result] */
  res->digits=decGetDigits(res->lsu, uc-res->lsu);
  res->exponent=0;                      /* integer */
  res->bits=0;                          /* sign=0 */
  return res;  /* [no status to set] */
  } /* decNumberAnd */

/* ------------------------------------------------------------------ */
/* decNumberCompare -- compare two Numbers                            */
/*                                                                    */
/*   This computes C = A ? B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for one digit (or NaN).                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberCompare(decNumber *res, const decNumber *lhs,
                             const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPARE, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberCompare */

/* ------------------------------------------------------------------ */
/* decNumberCompareSignal -- compare, signalling on all NaNs          */
/*                                                                    */
/*   This computes C = A ? B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for one digit (or NaN).                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberCompareSignal(decNumber *res, const decNumber *lhs,
                                   const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPSIG, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberCompareSignal */

/* ------------------------------------------------------------------ */
/* decNumberCompareTotal -- compare two Numbers, using total ordering */
/*                                                                    */
/*   This computes C = A ? B, under total ordering                    */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for one digit; the result will always be one of  */
/* -1, 0, or 1.                                                       */
/* ------------------------------------------------------------------ */
decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs,
                                  const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberCompareTotal */

/* ------------------------------------------------------------------ */
/* decNumberCompareTotalMag -- compare, total ordering of magnitudes  */
/*                                                                    */
/*   This computes C = |A| ? |B|, under total ordering                */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for one digit; the result will always be one of  */
/* -1, 0, or 1.                                                       */
/* ------------------------------------------------------------------ */
decNumber * decNumberCompareTotalMag(decNumber *res, const decNumber *lhs,
                                     const decNumber *rhs, decContext *set) {
  uInt status=0;                   /* accumulator */
  uInt needbytes;                  /* for space calculations */
  decNumber bufa[D2N(DECBUFFER+1)];/* +1 in case DECBUFFER=0 */
  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
  decNumber bufb[D2N(DECBUFFER+1)];
  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
  decNumber *a, *b;                /* temporary pointers */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  do {                                  /* protect allocated storage */
    /* if either is negative, take a copy and absolute */
    if (decNumberIsNegative(lhs)) {     /* lhs<0 */
      a=bufa;
      needbytes=sizeof(decNumber)+(D2U(lhs->digits)-1)*sizeof(Unit);
      if (needbytes>sizeof(bufa)) {     /* need malloc space */
        allocbufa=(decNumber *)malloc(needbytes);
        if (allocbufa==NULL) {          /* hopeless -- abandon */
          status|=DEC_Insufficient_storage;
          break;}
        a=allocbufa;                    /* use the allocated space */
        }
      decNumberCopy(a, lhs);            /* copy content */
      a->bits&=~DECNEG;                 /* .. and clear the sign */
      lhs=a;                            /* use copy from here on */
      }
    if (decNumberIsNegative(rhs)) {     /* rhs<0 */
      b=bufb;
      needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
      if (needbytes>sizeof(bufb)) {     /* need malloc space */
        allocbufb=(decNumber *)malloc(needbytes);
        if (allocbufb==NULL) {          /* hopeless -- abandon */
          status|=DEC_Insufficient_storage;
          break;}
        b=allocbufb;                    /* use the allocated space */
        }
      decNumberCopy(b, rhs);            /* copy content */
      b->bits&=~DECNEG;                 /* .. and clear the sign */
      rhs=b;                            /* use copy from here on */
      }
    decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
    } while(0);                         /* end protected */

  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
  if (allocbufb!=NULL) free(allocbufb); /* .. */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberCompareTotalMag */

/* ------------------------------------------------------------------ */
/* decNumberDivide -- divide one number by another                    */
/*                                                                    */
/*   This computes C = A / B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberDivide(decNumber *res, const decNumber *lhs,
                            const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decDivideOp(res, lhs, rhs, set, DIVIDE, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberDivide */

/* ------------------------------------------------------------------ */
/* decNumberDivideInteger -- divide and return integer quotient       */
/*                                                                    */
/*   This computes C = A # B, where # is the integer divide operator  */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X#X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs,
                                   const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decDivideOp(res, lhs, rhs, set, DIVIDEINT, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberDivideInteger */

/* ------------------------------------------------------------------ */
/* decNumberExp -- exponentiation                                     */
/*                                                                    */
/*   This computes C = exp(A)                                         */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context; note that rounding mode has no effect        */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Mathematical function restrictions apply (see above); a NaN is     */
/* returned with Invalid_operation if a restriction is violated.      */
/*                                                                    */
/* Finite results will always be full precision and Inexact, except   */
/* when A is a zero or -Infinity (giving 1 or 0 respectively).        */
/*                                                                    */
/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
/* almost always be correctly rounded, but may be up to 1 ulp in      */
/* error in rare cases.                                               */
/* ------------------------------------------------------------------ */
/* This is a wrapper for decExpOp which can handle the slightly wider */
/* (double) range needed by Ln (which has to be able to calculate     */
/* exp(-a) where a can be the tiniest number (Ntiny).                 */
/* ------------------------------------------------------------------ */
decNumber * decNumberExp(decNumber *res, const decNumber *rhs,
                         decContext *set) {
  uInt status=0;                        /* accumulator */
  #if DECSUBSET
  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
  #endif

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* Check restrictions; these restrictions ensure that if h=8 (see */
  /* decExpOp) then the result will either overflow or underflow to 0. */
  /* Other math functions restrict the input range, too, for inverses. */
  /* If not violated then carry out the operation. */
  if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
    #if DECSUBSET
    if (!set->extended) {
      /* reduce operand and set lostDigits status, as needed */
      if (rhs->digits>set->digits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        rhs=allocrhs;
        }
      }
    #endif
    decExpOp(res, rhs, set, &status);
    } while(0);                         /* end protected */

  #if DECSUBSET
  if (allocrhs !=NULL) free(allocrhs);  /* drop any storage used */
  #endif
  /* apply significant status */
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberExp */

/* ------------------------------------------------------------------ */
/* decNumberFMA -- fused multiply add                                 */
/*                                                                    */
/*   This computes D = (A * B) + C with only one rounding             */
/*                                                                    */
/*   res is D, the result.  D may be A or B or C (e.g., X=FMA(X,X,X)) */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   fhs is C [far hand side]                                         */
/*   set is the context                                               */
/*                                                                    */
/* Mathematical function restrictions apply (see above); a NaN is     */
/* returned with Invalid_operation if a restriction is violated.      */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberFMA(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, const decNumber *fhs,
                         decContext *set) {
  uInt status=0;                   /* accumulator */
  decContext dcmul;                /* context for the multiplication */
  uInt needbytes;                  /* for space calculations */
  decNumber bufa[D2N(DECBUFFER*2+1)];
  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
  decNumber *acc;                  /* accumulator pointer */
  decNumber dzero;                 /* work */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  if (decCheckOperands(res, fhs, DECUNUSED, set)) return res;
  #endif

  do {                                  /* protect allocated storage */
    #if DECSUBSET
    if (!set->extended) {               /* [undefined if subset] */
      status|=DEC_Invalid_operation;
      break;}
    #endif
    /* Check math restrictions [these ensure no overflow or underflow] */
    if ((!decNumberIsSpecial(lhs) && decCheckMath(lhs, set, &status))
     || (!decNumberIsSpecial(rhs) && decCheckMath(rhs, set, &status))
     || (!decNumberIsSpecial(fhs) && decCheckMath(fhs, set, &status))) break;
    /* set up context for multiply */
    dcmul=*set;
    dcmul.digits=lhs->digits+rhs->digits; /* just enough */
    /* [The above may be an over-estimate for subset arithmetic, but that's OK] */
    dcmul.emax=DEC_MAX_EMAX;            /* effectively unbounded .. */
    dcmul.emin=DEC_MIN_EMIN;            /* [thanks to Math restrictions] */
    /* set up decNumber space to receive the result of the multiply */
    acc=bufa;                           /* may fit */
    needbytes=sizeof(decNumber)+(D2U(dcmul.digits)-1)*sizeof(Unit);
    if (needbytes>sizeof(bufa)) {       /* need malloc space */
      allocbufa=(decNumber *)malloc(needbytes);
      if (allocbufa==NULL) {            /* hopeless -- abandon */
        status|=DEC_Insufficient_storage;
        break;}
      acc=allocbufa;                    /* use the allocated space */
      }
    /* multiply with extended range and necessary precision */
    /*printf("emin=%ld\n", dcmul.emin); */
    decMultiplyOp(acc, lhs, rhs, &dcmul, &status);
    /* Only Invalid operation (from sNaN or Inf * 0) is possible in */
    /* status; if either is seen than ignore fhs (in case it is */
    /* another sNaN) and set acc to NaN unless we had an sNaN */
    /* [decMultiplyOp leaves that to caller] */
    /* Note sNaN has to go through addOp to shorten payload if */
    /* necessary */
    if ((status&DEC_Invalid_operation)!=0) {
      if (!(status&DEC_sNaN)) {         /* but be true invalid */
        decNumberZero(res);             /* acc not yet set */
        res->bits=DECNAN;
        break;
        }
      decNumberZero(&dzero);            /* make 0 (any non-NaN would do) */
      fhs=&dzero;                       /* use that */
      }
    #if DECCHECK
     else { /* multiply was OK */
      if (status!=0) printf("Status=%08lx after FMA multiply\n", status);
      }
    #endif
    /* add the third operand and result -> res, and all is done */
    decAddOp(res, acc, fhs, set, 0, &status);
    } while(0);                         /* end protected */

  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberFMA */

/* ------------------------------------------------------------------ */
/* decNumberInvert -- invert a Number, digitwise                      */
/*                                                                    */
/*   This computes C = ~A                                             */
/*                                                                    */
/*   res is C, the result.  C may be A (e.g., X=~X)                   */
/*   rhs is A                                                         */
/*   set is the context (used for result length and error report)     */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Logical function restrictions apply (see above); a NaN is          */
/* returned with Invalid_operation if a restriction is violated.      */
/* ------------------------------------------------------------------ */
decNumber * decNumberInvert(decNumber *res, const decNumber *rhs,
                            decContext *set) {
  const Unit *ua, *msua;                /* -> operand and its msu */
  Unit  *uc, *msuc;                     /* -> result and its msu */
  Int   msudigs;                        /* digits in res msu */
  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  if (rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
    decStatus(res, DEC_Invalid_operation, set);
    return res;
    }
  /* operand is valid */
  ua=rhs->lsu;                          /* bottom-up */
  uc=res->lsu;                          /* .. */
  msua=ua+D2U(rhs->digits)-1;           /* -> msu of rhs */
  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
  for (; uc<=msuc; ua++, uc++) {        /* Unit loop */
    Unit a;                             /* extract unit */
    Int  i, j;                          /* work */
    if (ua>msua) a=0;
     else a=*ua;
    *uc=0;                              /* can now write back */
    /* always need to examine all bits in rhs */
    /* This loop could be unrolled and/or use BIN2BCD tables */
    for (i=0; i<DECDPUN; i++) {
      if ((~a)&1) *uc=*uc+(Unit)powers[i];   /* effect INVERT */
      j=a%10;
      a=a/10;
      if (j>1) {
        decStatus(res, DEC_Invalid_operation, set);
        return res;
        }
      if (uc==msuc && i==msudigs-1) break;   /* just did final digit */
      } /* each digit */
    } /* each unit */
  /* [here uc-1 is the msu of the result] */
  res->digits=decGetDigits(res->lsu, uc-res->lsu);
  res->exponent=0;                      /* integer */
  res->bits=0;                          /* sign=0 */
  return res;  /* [no status to set] */
  } /* decNumberInvert */

/* ------------------------------------------------------------------ */
/* decNumberLn -- natural logarithm                                   */
/*                                                                    */
/*   This computes C = ln(A)                                          */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context; note that rounding mode has no effect        */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Notable cases:                                                     */
/*   A<0 -> Invalid                                                   */
/*   A=0 -> -Infinity (Exact)                                         */
/*   A=+Infinity -> +Infinity (Exact)                                 */
/*   A=1 exactly -> 0 (Exact)                                         */
/*                                                                    */
/* Mathematical function restrictions apply (see above); a NaN is     */
/* returned with Invalid_operation if a restriction is violated.      */
/*                                                                    */
/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
/* almost always be correctly rounded, but may be up to 1 ulp in      */
/* error in rare cases.                                               */
/* ------------------------------------------------------------------ */
/* This is a wrapper for decLnOp which can handle the slightly wider  */
/* (+11) range needed by Ln, Log10, etc. (which may have to be able   */
/* to calculate at p+e+2).                                            */
/* ------------------------------------------------------------------ */
decNumber * decNumberLn(decNumber *res, const decNumber *rhs,
                        decContext *set) {
  uInt status=0;                   /* accumulator */
  #if DECSUBSET
  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
  #endif

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* Check restrictions; this is a math function; if not violated */
  /* then carry out the operation. */
  if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
    #if DECSUBSET
    if (!set->extended) {
      /* reduce operand and set lostDigits status, as needed */
      if (rhs->digits>set->digits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        rhs=allocrhs;
        }
      /* special check in subset for rhs=0 */
      if (ISZERO(rhs)) {                /* +/- zeros -> error */
        status|=DEC_Invalid_operation;
        break;}
      } /* extended=0 */
    #endif
    decLnOp(res, rhs, set, &status);
    } while(0);                         /* end protected */

  #if DECSUBSET
  if (allocrhs !=NULL) free(allocrhs);  /* drop any storage used */
  #endif
  /* apply significant status */
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberLn */

/* ------------------------------------------------------------------ */
/* decNumberLogB - get adjusted exponent, by 754r rules               */
/*                                                                    */
/*   This computes C = adjustedexponent(A)                            */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context, used only for digits and status              */
/*                                                                    */
/* C must have space for 10 digits (A might have 10**9 digits and     */
/* an exponent of +999999999, or one digit and an exponent of         */
/* -1999999999).                                                      */
/*                                                                    */
/* This returns the adjusted exponent of A after (in theory) padding  */
/* with zeros on the right to set->digits digits while keeping the    */
/* same value.  The exponent is not limited by emin/emax.             */
/*                                                                    */
/* Notable cases:                                                     */
/*   A<0 -> Use |A|                                                   */
/*   A=0 -> -Infinity (Division by zero)                              */
/*   A=Infinite -> +Infinity (Exact)                                  */
/*   A=1 exactly -> 0 (Exact)                                         */
/*   NaNs are propagated as usual                                     */
/* ------------------------------------------------------------------ */
decNumber * decNumberLogB(decNumber *res, const decNumber *rhs,
                          decContext *set) {
  uInt status=0;                   /* accumulator */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* NaNs as usual; Infinities return +Infinity; 0->oops */
  if (decNumberIsNaN(rhs)) decNaNs(res, rhs, NULL, set, &status);
   else if (decNumberIsInfinite(rhs)) decNumberCopyAbs(res, rhs);
   else if (decNumberIsZero(rhs)) {
    decNumberZero(res);                 /* prepare for Infinity */
    res->bits=DECNEG|DECINF;            /* -Infinity */
    status|=DEC_Division_by_zero;       /* as per 754r */
    }
   else { /* finite non-zero */
    Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */
    decNumberFromInt32(res, ae);        /* lay it out */
    }

  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberLogB */

/* ------------------------------------------------------------------ */
/* decNumberLog10 -- logarithm in base 10                             */
/*                                                                    */
/*   This computes C = log10(A)                                       */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context; note that rounding mode has no effect        */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Notable cases:                                                     */
/*   A<0 -> Invalid                                                   */
/*   A=0 -> -Infinity (Exact)                                         */
/*   A=+Infinity -> +Infinity (Exact)                                 */
/*   A=10**n (if n is an integer) -> n (Exact)                        */
/*                                                                    */
/* Mathematical function restrictions apply (see above); a NaN is     */
/* returned with Invalid_operation if a restriction is violated.      */
/*                                                                    */
/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
/* almost always be correctly rounded, but may be up to 1 ulp in      */
/* error in rare cases.                                               */
/* ------------------------------------------------------------------ */
/* This calculates ln(A)/ln(10) using appropriate precision.  For     */
/* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the      */
/* requested digits and t is the number of digits in the exponent     */
/* (maximum 6).  For ln(10) it is p + 3; this is often handled by the */
/* fastpath in decLnOp.  The final division is done to the requested  */
/* precision.                                                         */
/* ------------------------------------------------------------------ */
decNumber * decNumberLog10(decNumber *res, const decNumber *rhs,
                          decContext *set) {
  uInt status=0, ignore=0;         /* status accumulators */
  uInt needbytes;                  /* for space calculations */
  Int p;                           /* working precision */
  Int t;                           /* digits in exponent of A */

  /* buffers for a and b working decimals */
  /* (adjustment calculator, same size) */
  decNumber bufa[D2N(DECBUFFER+2)];
  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
  decNumber *a=bufa;               /* temporary a */
  decNumber bufb[D2N(DECBUFFER+2)];
  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
  decNumber *b=bufb;               /* temporary b */
  decNumber bufw[D2N(10)];         /* working 2-10 digit number */
  decNumber *w=bufw;               /* .. */
  #if DECSUBSET
  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
  #endif

  decContext aset;                 /* working context */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* Check restrictions; this is a math function; if not violated */
  /* then carry out the operation. */
  if (!decCheckMath(rhs, set, &status)) do { /* protect malloc */
    #if DECSUBSET
    if (!set->extended) {
      /* reduce operand and set lostDigits status, as needed */
      if (rhs->digits>set->digits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        rhs=allocrhs;
        }
      /* special check in subset for rhs=0 */
      if (ISZERO(rhs)) {                /* +/- zeros -> error */
        status|=DEC_Invalid_operation;
        break;}
      } /* extended=0 */
    #endif

    decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */

    /* handle exact powers of 10; only check if +ve finite */
    if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) {
      Int residue=0;               /* (no residue) */
      uInt copystat=0;             /* clean status */

      /* round to a single digit... */
      aset.digits=1;
      decCopyFit(w, rhs, &aset, &residue, &copystat); /* copy & shorten */
      /* if exact and the digit is 1, rhs is a power of 10 */
      if (!(copystat&DEC_Inexact) && w->lsu[0]==1) {
        /* the exponent, conveniently, is the power of 10; making */
        /* this the result needs a little care as it might not fit, */
        /* so first convert it into the working number, and then move */
        /* to res */
        decNumberFromInt32(w, w->exponent);
        residue=0;
        decCopyFit(res, w, set, &residue, &status); /* copy & round */
        decFinish(res, set, &residue, &status);     /* cleanup/set flags */
        break;
        } /* not a power of 10 */
      } /* not a candidate for exact */

    /* simplify the information-content calculation to use 'total */
    /* number of digits in a, including exponent' as compared to the */
    /* requested digits, as increasing this will only rarely cost an */
    /* iteration in ln(a) anyway */
    t=6;                                /* it can never be >6 */

    /* allocate space when needed... */
    p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3;
    needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
    if (needbytes>sizeof(bufa)) {       /* need malloc space */
      allocbufa=(decNumber *)malloc(needbytes);
      if (allocbufa==NULL) {            /* hopeless -- abandon */
        status|=DEC_Insufficient_storage;
        break;}
      a=allocbufa;                      /* use the allocated space */
      }
    aset.digits=p;                      /* as calculated */
    aset.emax=DEC_MAX_MATH;             /* usual bounds */
    aset.emin=-DEC_MAX_MATH;            /* .. */
    aset.clamp=0;                       /* and no concrete format */
    decLnOp(a, rhs, &aset, &status);    /* a=ln(rhs) */

    /* skip the division if the result so far is infinite, NaN, or */
    /* zero, or there was an error; note NaN from sNaN needs copy */
    if (status&DEC_NaNs && !(status&DEC_sNaN)) break;
    if (a->bits&DECSPECIAL || ISZERO(a)) {
      decNumberCopy(res, a);            /* [will fit] */
      break;}

    /* for ln(10) an extra 3 digits of precision are needed */
    p=set->digits+3;
    needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
    if (needbytes>sizeof(bufb)) {       /* need malloc space */
      allocbufb=(decNumber *)malloc(needbytes);
      if (allocbufb==NULL) {            /* hopeless -- abandon */
        status|=DEC_Insufficient_storage;
        break;}
      b=allocbufb;                      /* use the allocated space */
      }
    decNumberZero(w);                   /* set up 10... */
    #if DECDPUN==1
    w->lsu[1]=1; w->lsu[0]=0;           /* .. */
    #else
    w->lsu[0]=10;                       /* .. */
    #endif
    w->digits=2;                        /* .. */

    aset.digits=p;
    decLnOp(b, w, &aset, &ignore);      /* b=ln(10) */

    aset.digits=set->digits;            /* for final divide */
    decDivideOp(res, a, b, &aset, DIVIDE, &status); /* into result */
    } while(0);                         /* [for break] */

  if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
  if (allocbufb!=NULL) free(allocbufb); /* .. */
  #if DECSUBSET
  if (allocrhs !=NULL) free(allocrhs);  /* .. */
  #endif
  /* apply significant status */
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberLog10 */

/* ------------------------------------------------------------------ */
/* decNumberMax -- compare two Numbers and return the maximum         */
/*                                                                    */
/*   This computes C = A ? B, returning the maximum by 754R rules     */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberMax(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPMAX, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMax */

/* ------------------------------------------------------------------ */
/* decNumberMaxMag -- compare and return the maximum by magnitude     */
/*                                                                    */
/*   This computes C = A ? B, returning the maximum by 754R rules     */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberMaxMag(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPMAXMAG, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMaxMag */

/* ------------------------------------------------------------------ */
/* decNumberMin -- compare two Numbers and return the minimum         */
/*                                                                    */
/*   This computes C = A ? B, returning the minimum by 754R rules     */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberMin(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPMIN, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMin */

/* ------------------------------------------------------------------ */
/* decNumberMinMag -- compare and return the minimum by magnitude     */
/*                                                                    */
/*   This computes C = A ? B, returning the minimum by 754R rules     */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberMinMag(decNumber *res, const decNumber *lhs,
                         const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decCompareOp(res, lhs, rhs, set, COMPMINMAG, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMinMag */

/* ------------------------------------------------------------------ */
/* decNumberMinus -- prefix minus operator                            */
/*                                                                    */
/*   This computes C = 0 - A                                          */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* See also decNumberCopyNegate for a quiet bitwise version of this.  */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* Simply use AddOp for the subtract, which will do the necessary.    */
/* ------------------------------------------------------------------ */
decNumber * decNumberMinus(decNumber *res, const decNumber *rhs,
                           decContext *set) {
  decNumber dzero;
  uInt status=0;                        /* accumulator */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  decNumberZero(&dzero);                /* make 0 */
  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
  decAddOp(res, &dzero, rhs, set, DECNEG, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMinus */

/* ------------------------------------------------------------------ */
/* decNumberNextMinus -- next towards -Infinity                       */
/*                                                                    */
/*   This computes C = A - infinitesimal, rounded towards -Infinity   */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* This is a generalization of 754r NextDown.                         */
/* ------------------------------------------------------------------ */
decNumber * decNumberNextMinus(decNumber *res, const decNumber *rhs,
                               decContext *set) {
  decNumber dtiny;                           /* constant */
  decContext workset=*set;                   /* work */
  uInt status=0;                             /* accumulator */
  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* +Infinity is the special case */
  if ((rhs->bits&(DECINF|DECNEG))==DECINF) {
    decSetMaxValue(res, set);                /* is +ve */
    /* there is no status to set */
    return res;
    }
  decNumberZero(&dtiny);                     /* start with 0 */
  dtiny.lsu[0]=1;                            /* make number that is .. */
  dtiny.exponent=DEC_MIN_EMIN-1;             /* .. smaller than tiniest */
  workset.round=DEC_ROUND_FLOOR;
  decAddOp(res, rhs, &dtiny, &workset, DECNEG, &status);
  status&=DEC_Invalid_operation|DEC_sNaN;    /* only sNaN Invalid please */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberNextMinus */

/* ------------------------------------------------------------------ */
/* decNumberNextPlus -- next towards +Infinity                        */
/*                                                                    */
/*   This computes C = A + infinitesimal, rounded towards +Infinity   */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* This is a generalization of 754r NextUp.                           */
/* ------------------------------------------------------------------ */
decNumber * decNumberNextPlus(decNumber *res, const decNumber *rhs,
                              decContext *set) {
  decNumber dtiny;                           /* constant */
  decContext workset=*set;                   /* work */
  uInt status=0;                             /* accumulator */
  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  /* -Infinity is the special case */
  if ((rhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
    decSetMaxValue(res, set);
    res->bits=DECNEG;                        /* negative */
    /* there is no status to set */
    return res;
    }
  decNumberZero(&dtiny);                     /* start with 0 */
  dtiny.lsu[0]=1;                            /* make number that is .. */
  dtiny.exponent=DEC_MIN_EMIN-1;             /* .. smaller than tiniest */
  workset.round=DEC_ROUND_CEILING;
  decAddOp(res, rhs, &dtiny, &workset, 0, &status);
  status&=DEC_Invalid_operation|DEC_sNaN;    /* only sNaN Invalid please */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberNextPlus */

/* ------------------------------------------------------------------ */
/* decNumberNextToward -- next towards rhs                            */
/*                                                                    */
/*   This computes C = A +/- infinitesimal, rounded towards           */
/*   +/-Infinity in the direction of B, as per 754r nextafter rules   */
/*                                                                    */
/*   res is C, the result.  C may be A or B.                          */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* This is a generalization of 754r NextAfter.                        */
/* ------------------------------------------------------------------ */
decNumber * decNumberNextToward(decNumber *res, const decNumber *lhs,
                                const decNumber *rhs, decContext *set) {
  decNumber dtiny;                           /* constant */
  decContext workset=*set;                   /* work */
  Int result;                                /* .. */
  uInt status=0;                             /* accumulator */
  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) {
    decNaNs(res, lhs, rhs, set, &status);
    }
   else { /* Is numeric, so no chance of sNaN Invalid, etc. */
    result=decCompare(lhs, rhs, 0);     /* sign matters */
    if (result==BADINT) status|=DEC_Insufficient_storage; /* rare */
     else { /* valid compare */
      if (result==0) decNumberCopySign(res, lhs, rhs); /* easy */
       else { /* differ: need NextPlus or NextMinus */
        uByte sub;                      /* add or subtract */
        if (result<0) {                 /* lhs<rhs, do nextplus */
          /* -Infinity is the special case */
          if ((lhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
            decSetMaxValue(res, set);
            res->bits=DECNEG;           /* negative */
            return res;                 /* there is no status to set */
            }
          workset.round=DEC_ROUND_CEILING;
          sub=0;                        /* add, please */
          } /* plus */
         else {                         /* lhs>rhs, do nextminus */
          /* +Infinity is the special case */
          if ((lhs->bits&(DECINF|DECNEG))==DECINF) {
            decSetMaxValue(res, set);
            return res;                 /* there is no status to set */
            }
          workset.round=DEC_ROUND_FLOOR;
          sub=DECNEG;                   /* subtract, please */
          } /* minus */
        decNumberZero(&dtiny);          /* start with 0 */
        dtiny.lsu[0]=1;                 /* make number that is .. */
        dtiny.exponent=DEC_MIN_EMIN-1;  /* .. smaller than tiniest */
        decAddOp(res, lhs, &dtiny, &workset, sub, &status); /* + or - */
        /* turn off exceptions if the result is a normal number */
        /* (including Nmin), otherwise let all status through */
        if (decNumberIsNormal(res, set)) status=0;
        } /* unequal */
      } /* compare OK */
    } /* numeric */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberNextToward */

/* ------------------------------------------------------------------ */
/* decNumberOr -- OR two Numbers, digitwise                           */
/*                                                                    */
/*   This computes C = A | B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X|X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context (used for result length and error report)     */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Logical function restrictions apply (see above); a NaN is          */
/* returned with Invalid_operation if a restriction is violated.      */
/* ------------------------------------------------------------------ */
decNumber * decNumberOr(decNumber *res, const decNumber *lhs,
                        const decNumber *rhs, decContext *set) {
  const Unit *ua, *ub;                  /* -> operands */
  const Unit *msua, *msub;              /* -> operand msus */
  Unit  *uc, *msuc;                     /* -> result and its msu */
  Int   msudigs;                        /* digits in res msu */
  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
   || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
    decStatus(res, DEC_Invalid_operation, set);
    return res;
    }
  /* operands are valid */
  ua=lhs->lsu;                          /* bottom-up */
  ub=rhs->lsu;                          /* .. */
  uc=res->lsu;                          /* .. */
  msua=ua+D2U(lhs->digits)-1;           /* -> msu of lhs */
  msub=ub+D2U(rhs->digits)-1;           /* -> msu of rhs */
  msuc=uc+D2U(set->digits)-1;           /* -> msu of result */
  msudigs=MSUDIGITS(set->digits);       /* [faster than remainder] */
  for (; uc<=msuc; ua++, ub++, uc++) {  /* Unit loop */
    Unit a, b;                          /* extract units */
    if (ua>msua) a=0;
     else a=*ua;
    if (ub>msub) b=0;
     else b=*ub;
    *uc=0;                              /* can now write back */
    if (a|b) {                          /* maybe 1 bits to examine */
      Int i, j;
      /* This loop could be unrolled and/or use BIN2BCD tables */
      for (i=0; i<DECDPUN; i++) {
        if ((a|b)&1) *uc=*uc+(Unit)powers[i];     /* effect OR */
        j=a%10;
        a=a/10;
        j|=b%10;
        b=b/10;
        if (j>1) {
          decStatus(res, DEC_Invalid_operation, set);
          return res;
          }
        if (uc==msuc && i==msudigs-1) break;      /* just did final digit */
        } /* each digit */
      } /* non-zero */
    } /* each unit */
  /* [here uc-1 is the msu of the result] */
  res->digits=decGetDigits(res->lsu, uc-res->lsu);
  res->exponent=0;                      /* integer */
  res->bits=0;                          /* sign=0 */
  return res;  /* [no status to set] */
  } /* decNumberOr */

/* ------------------------------------------------------------------ */
/* decNumberPlus -- prefix plus operator                              */
/*                                                                    */
/*   This computes C = 0 + A                                          */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* See also decNumberCopy for a quiet bitwise version of this.        */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* This simply uses AddOp; Add will take fast path after preparing A. */
/* Performance is a concern here, as this routine is often used to    */
/* check operands and apply rounding and overflow/underflow testing.  */
/* ------------------------------------------------------------------ */
decNumber * decNumberPlus(decNumber *res, const decNumber *rhs,
                          decContext *set) {
  decNumber dzero;
  uInt status=0;                        /* accumulator */
  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  decNumberZero(&dzero);                /* make 0 */
  dzero.exponent=rhs->exponent;         /* [no coefficient expansion] */
  decAddOp(res, &dzero, rhs, set, 0, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberPlus */

/* ------------------------------------------------------------------ */
/* decNumberMultiply -- multiply two Numbers                          */
/*                                                                    */
/*   This computes C = A x B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs,
                              const decNumber *rhs, decContext *set) {
  uInt status=0;                   /* accumulator */
  decMultiplyOp(res, lhs, rhs, set, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberMultiply */

/* ------------------------------------------------------------------ */
/* decNumberPower -- raise a number to a power                        */
/*                                                                    */
/*   This computes C = A ** B                                         */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X**X)        */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Mathematical function restrictions apply (see above); a NaN is     */
/* returned with Invalid_operation if a restriction is violated.      */
/*                                                                    */
/* However, if 1999999997<=B<=999999999 and B is an integer then the  */
/* restrictions on A and the context are relaxed to the usual bounds, */
/* for compatibility with the earlier (integer power only) version    */
/* of this function.                                                  */
/*                                                                    */
/* When B is an integer, the result may be exact, even if rounded.    */
/*                                                                    */
/* The final result is rounded according to the context; it will      */
/* almost always be correctly rounded, but may be up to 1 ulp in      */
/* error in rare cases.                                               */
/* ------------------------------------------------------------------ */
decNumber * decNumberPower(decNumber *res, const decNumber *lhs,
                           const decNumber *rhs, decContext *set) {
  #if DECSUBSET
  decNumber *alloclhs=NULL;        /* non-NULL if rounded lhs allocated */
  decNumber *allocrhs=NULL;        /* .., rhs */
  #endif
  decNumber *allocdac=NULL;        /* -> allocated acc buffer, iff used */
  decNumber *allocinv=NULL;        /* -> allocated 1/x buffer, iff used */
  Int   reqdigits=set->digits;     /* requested DIGITS */
  Int   n;                         /* rhs in binary */
  Flag  rhsint=0;                  /* 1 if rhs is an integer */
  Flag  useint=0;                  /* 1 if can use integer calculation */
  Flag  isoddint=0;                /* 1 if rhs is an integer and odd */
  Int   i;                         /* work */
  #if DECSUBSET
  Int   dropped;                   /* .. */
  #endif
  uInt  needbytes;                 /* buffer size needed */
  Flag  seenbit;                   /* seen a bit while powering */
  Int   residue=0;                 /* rounding residue */
  uInt  status=0;                  /* accumulators */
  uByte bits=0;                    /* result sign if errors */
  decContext aset;                 /* working context */
  decNumber dnOne;                 /* work value 1... */
  /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
  decNumber dacbuff[D2N(DECBUFFER+9)];
  decNumber *dac=dacbuff;          /* -> result accumulator */
  /* same again for possible 1/lhs calculation */
  decNumber invbuff[D2N(DECBUFFER+9)];

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  do {                             /* protect allocated storage */
    #if DECSUBSET
    if (!set->extended) { /* reduce operands and set status, as needed */
      if (lhs->digits>reqdigits) {
        alloclhs=decRoundOperand(lhs, set, &status);
        if (alloclhs==NULL) break;
        lhs=alloclhs;
        }
      if (rhs->digits>reqdigits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        rhs=allocrhs;
        }
      }
    #endif
    /* [following code does not require input rounding] */

    /* handle NaNs and rhs Infinity (lhs infinity is harder) */
    if (SPECIALARGS) {
      if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { /* NaNs */
        decNaNs(res, lhs, rhs, set, &status);
        break;}
      if (decNumberIsInfinite(rhs)) {   /* rhs Infinity */
        Flag rhsneg=rhs->bits&DECNEG;   /* save rhs sign */
        if (decNumberIsNegative(lhs)    /* lhs<0 */
         && !decNumberIsZero(lhs))      /* .. */
          status|=DEC_Invalid_operation;
         else {                         /* lhs >=0 */
          decNumberZero(&dnOne);        /* set up 1 */
          dnOne.lsu[0]=1;
          decNumberCompare(dac, lhs, &dnOne, set); /* lhs ? 1 */
          decNumberZero(res);           /* prepare for 0/1/Infinity */
          if (decNumberIsNegative(dac)) {    /* lhs<1 */
            if (rhsneg) res->bits|=DECINF;   /* +Infinity [else is +0] */
            }
           else if (dac->lsu[0]==0) {        /* lhs=1 */
            /* 1**Infinity is inexact, so return fully-padded 1.0000 */
            Int shift=set->digits-1;
            *res->lsu=1;                     /* was 0, make int 1 */
            res->digits=decShiftToMost(res->lsu, 1, shift);
            res->exponent=-shift;            /* make 1.0000... */
            status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */
            }
           else {                            /* lhs>1 */
            if (!rhsneg) res->bits|=DECINF;  /* +Infinity [else is +0] */
            }
          } /* lhs>=0 */
        break;}
      /* [lhs infinity drops through] */
      } /* specials */

    /* Original rhs may be an integer that fits and is in range */
    n=decGetInt(rhs);
    if (n!=BADINT) {                    /* it is an integer */
      rhsint=1;                         /* record the fact for 1**n */
      isoddint=(Flag)n&1;               /* [works even if big] */
      if (n!=BIGEVEN && n!=BIGODD)      /* can use integer path? */
        useint=1;                       /* looks good */
      }

    if (decNumberIsNegative(lhs)        /* -x .. */
      && isoddint) bits=DECNEG;         /* .. to an odd power */

    /* handle LHS infinity */
    if (decNumberIsInfinite(lhs)) {     /* [NaNs already handled] */
      uByte rbits=rhs->bits;            /* save */
      decNumberZero(res);               /* prepare */
      if (n==0) *res->lsu=1;            /* [-]Inf**0 => 1 */
       else {
        /* -Inf**nonint -> error */
        if (!rhsint && decNumberIsNegative(lhs)) {
          status|=DEC_Invalid_operation;     /* -Inf**nonint is error */
          break;}
        if (!(rbits & DECNEG)) bits|=DECINF; /* was not a **-n */
        /* [otherwise will be 0 or -0] */
        res->bits=bits;
        }
      break;}

    /* similarly handle LHS zero */
    if (decNumberIsZero(lhs)) {
      if (n==0) {                            /* 0**0 => Error */
        #if DECSUBSET
        if (!set->extended) {                /* [unless subset] */
          decNumberZero(res);
          *res->lsu=1;                       /* return 1 */
          break;}
        #endif
        status|=DEC_Invalid_operation;
        }
       else {                                /* 0**x */
        uByte rbits=rhs->bits;               /* save */
        if (rbits & DECNEG) {                /* was a 0**(-n) */
          #if DECSUBSET
          if (!set->extended) {              /* [bad if subset] */
            status|=DEC_Invalid_operation;
            break;}
          #endif
          bits|=DECINF;
          }
        decNumberZero(res);                  /* prepare */
        /* [otherwise will be 0 or -0] */
        res->bits=bits;
        }
      break;}

    /* here both lhs and rhs are finite; rhs==0 is handled in the */
    /* integer path.  Next handle the non-integer cases */
    if (!useint) {                      /* non-integral rhs */
      /* any -ve lhs is bad, as is either operand or context out of */
      /* bounds */
      if (decNumberIsNegative(lhs)) {
        status|=DEC_Invalid_operation;
        break;}
      if (decCheckMath(lhs, set, &status)
       || decCheckMath(rhs, set, &status)) break; /* variable status */

      decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */
      aset.emax=DEC_MAX_MATH;           /* usual bounds */
      aset.emin=-DEC_MAX_MATH;          /* .. */
      aset.clamp=0;                     /* and no concrete format */

      /* calculate the result using exp(ln(lhs)*rhs), which can */
      /* all be done into the accumulator, dac.  The precision needed */
      /* is enough to contain the full information in the lhs (which */
      /* is the total digits, including exponent), or the requested */
      /* precision, if larger, + 4; 6 is used for the exponent */
      /* maximum length, and this is also used when it is shorter */
      /* than the requested digits as it greatly reduces the >0.5 ulp */
      /* cases at little cost (because Ln doubles digits each */
      /* iteration so a few extra digits rarely causes an extra */
      /* iteration) */
      aset.digits=MAXI(lhs->digits, set->digits)+6+4;
      } /* non-integer rhs */

     else { /* rhs is in-range integer */
      if (n==0) {                       /* x**0 = 1 */
        /* (0**0 was handled above) */
        decNumberZero(res);             /* result=1 */
        *res->lsu=1;                    /* .. */
        break;}
      /* rhs is a non-zero integer */
      if (n<0) n=-n;                    /* use abs(n) */

      aset=*set;                        /* clone the context */
      aset.round=DEC_ROUND_HALF_EVEN;   /* internally use balanced */
      /* calculate the working DIGITS */
      aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2;
      #if DECSUBSET
      if (!set->extended) aset.digits--;     /* use classic precision */
      #endif
      /* it's an error if this is more than can be handled */
      if (aset.digits>DECNUMMAXP) {status|=DEC_Invalid_operation; break;}
      } /* integer path */

    /* aset.digits is the count of digits for the accumulator needed */
    /* if accumulator is too long for local storage, then allocate */
    needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit);
    /* [needbytes also used below if 1/lhs needed] */
    if (needbytes>sizeof(dacbuff)) {
      allocdac=(decNumber *)malloc(needbytes);
      if (allocdac==NULL) {   /* hopeless -- abandon */
        status|=DEC_Insufficient_storage;
        break;}
      dac=allocdac;           /* use the allocated space */
      }
    /* here, aset is set up and accumulator is ready for use */

    if (!useint) {                           /* non-integral rhs */
      /* x ** y; special-case x=1 here as it will otherwise always */
      /* reduce to integer 1; decLnOp has a fastpath which detects */
      /* the case of x=1 */
      decLnOp(dac, lhs, &aset, &status);     /* dac=ln(lhs) */
      /* [no error possible, as lhs 0 already handled] */
      if (ISZERO(dac)) {                     /* x==1, 1.0, etc. */
        /* need to return fully-padded 1.0000 etc., but rhsint->1 */
        *dac->lsu=1;                         /* was 0, make int 1 */
        if (!rhsint) {                       /* add padding */
          Int shift=set->digits-1;
          dac->digits=decShiftToMost(dac->lsu, 1, shift);
          dac->exponent=-shift;              /* make 1.0000... */
          status|=DEC_Inexact|DEC_Rounded;   /* deemed inexact */
          }
        }
       else {
        decMultiplyOp(dac, dac, rhs, &aset, &status);  /* dac=dac*rhs */
        decExpOp(dac, dac, &aset, &status);            /* dac=exp(dac) */
        }
      /* and drop through for final rounding */
      } /* non-integer rhs */

     else {                             /* carry on with integer */
      decNumberZero(dac);               /* acc=1 */
      *dac->lsu=1;                      /* .. */

      /* if a negative power the constant 1 is needed, and if not subset */
      /* invert the lhs now rather than inverting the result later */
      if (decNumberIsNegative(rhs)) {   /* was a **-n [hence digits>0] */
        decNumber *inv=invbuff;         /* assume use fixed buffer */
        decNumberCopy(&dnOne, dac);     /* dnOne=1;  [needed now or later] */
        #if DECSUBSET
        if (set->extended) {            /* need to calculate 1/lhs */
        #endif
          /* divide lhs into 1, putting result in dac [dac=1/dac] */
          decDivideOp(dac, &dnOne, lhs, &aset, DIVIDE, &status);
          /* now locate or allocate space for the inverted lhs */
          if (needbytes>sizeof(invbuff)) {
            allocinv=(decNumber *)malloc(needbytes);
            if (allocinv==NULL) {       /* hopeless -- abandon */
              status|=DEC_Insufficient_storage;
              break;}
            inv=allocinv;               /* use the allocated space */
            }
          /* [inv now points to big-enough buffer or allocated storage] */
          decNumberCopy(inv, dac);      /* copy the 1/lhs */
          decNumberCopy(dac, &dnOne);   /* restore acc=1 */
          lhs=inv;                      /* .. and go forward with new lhs */
        #if DECSUBSET
          }
        #endif
        }

      /* Raise-to-the-power loop... */
      seenbit=0;                   /* set once a 1-bit is encountered */
      for (i=1;;i++){              /* for each bit [top bit ignored] */
        /* abandon if had overflow or terminal underflow */
        if (status & (DEC_Overflow|DEC_Underflow)) { /* interesting? */
          if (status&DEC_Overflow || ISZERO(dac)) break;
          }
        /* [the following two lines revealed an optimizer bug in a C++ */
        /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
        n=n<<1;                    /* move next bit to testable position */
        if (n<0) {                 /* top bit is set */
          seenbit=1;               /* OK, significant bit seen */
          decMultiplyOp(dac, dac, lhs, &aset, &status); /* dac=dac*x */
          }
        if (i==31) break;          /* that was the last bit */
        if (!seenbit) continue;    /* no need to square 1 */
        decMultiplyOp(dac, dac, dac, &aset, &status); /* dac=dac*dac [square] */
        } /*i*/ /* 32 bits */

      /* complete internal overflow or underflow processing */
      if (status & (DEC_Overflow|DEC_Underflow)) {
        #if DECSUBSET
        /* If subset, and power was negative, reverse the kind of -erflow */
        /* [1/x not yet done] */
        if (!set->extended && decNumberIsNegative(rhs)) {
          if (status & DEC_Overflow)
            status^=DEC_Overflow | DEC_Underflow | DEC_Subnormal;
           else { /* trickier -- Underflow may or may not be set */
            status&=~(DEC_Underflow | DEC_Subnormal); /* [one or both] */
            status|=DEC_Overflow;
            }
          }
        #endif
        dac->bits=(dac->bits & ~DECNEG) | bits; /* force correct sign */
        /* round subnormals [to set.digits rather than aset.digits] */
        /* or set overflow result similarly as required */
        decFinalize(dac, set, &residue, &status);
        decNumberCopy(res, dac);   /* copy to result (is now OK length) */
        break;
        }

      #if DECSUBSET
      if (!set->extended &&                  /* subset math */
          decNumberIsNegative(rhs)) {        /* was a **-n [hence digits>0] */
        /* so divide result into 1 [dac=1/dac] */
        decDivideOp(dac, &dnOne, dac, &aset, DIVIDE, &status);
        }
      #endif
      } /* rhs integer path */

    /* reduce result to the requested length and copy to result */
    decCopyFit(res, dac, set, &residue, &status);
    decFinish(res, set, &residue, &status);  /* final cleanup */
    #if DECSUBSET
    if (!set->extended) decTrim(res, set, 0, &dropped); /* trailing zeros */
    #endif
    } while(0);                         /* end protected */

  if (allocdac!=NULL) free(allocdac);   /* drop any storage used */
  if (allocinv!=NULL) free(allocinv);   /* .. */
  #if DECSUBSET
  if (alloclhs!=NULL) free(alloclhs);   /* .. */
  if (allocrhs!=NULL) free(allocrhs);   /* .. */
  #endif
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberPower */

/* ------------------------------------------------------------------ */
/* decNumberQuantize -- force exponent to requested value             */
/*                                                                    */
/*   This computes C = op(A, B), where op adjusts the coefficient     */
/*   of C (by rounding or shifting) such that the exponent (-scale)   */
/*   of C has exponent of B.  The numerical value of C will equal A,  */
/*   except for the effects of any rounding that occurred.            */
/*                                                                    */
/*   res is C, the result.  C may be A or B                           */
/*   lhs is A, the number to adjust                                   */
/*   rhs is B, the number with exponent to match                      */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Unless there is an error or the result is infinite, the exponent   */
/* after the operation is guaranteed to be equal to that of B.        */
/* ------------------------------------------------------------------ */
decNumber * decNumberQuantize(decNumber *res, const decNumber *lhs,
                              const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decQuantizeOp(res, lhs, rhs, set, 1, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberQuantize */

/* ------------------------------------------------------------------ */
/* decNumberReduce -- remove trailing zeros                           */
/*                                                                    */
/*   This computes C = 0 + A, and normalizes the result               */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* Previously known as Normalize */
decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs,
                               decContext *set) {
  return decNumberReduce(res, rhs, set);
  } /* decNumberNormalize */

decNumber * decNumberReduce(decNumber *res, const decNumber *rhs,
                            decContext *set) {
  #if DECSUBSET
  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
  #endif
  uInt status=0;                   /* as usual */
  Int  residue=0;                  /* as usual */
  Int  dropped;                    /* work */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  do {                             /* protect allocated storage */
    #if DECSUBSET
    if (!set->extended) {
      /* reduce operand and set lostDigits status, as needed */
      if (rhs->digits>set->digits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        rhs=allocrhs;
        }
      }
    #endif
    /* [following code does not require input rounding] */

    /* Infinities copy through; NaNs need usual treatment */
    if (decNumberIsNaN(rhs)) {
      decNaNs(res, rhs, NULL, set, &status);
      break;
      }

    /* reduce result to the requested length and copy to result */
    decCopyFit(res, rhs, set, &residue, &status); /* copy & round */
    decFinish(res, set, &residue, &status);       /* cleanup/set flags */
    decTrim(res, set, 1, &dropped);               /* normalize in place */
    } while(0);                              /* end protected */

  #if DECSUBSET
  if (allocrhs !=NULL) free(allocrhs);       /* .. */
  #endif
  if (status!=0) decStatus(res, status, set);/* then report status */
  return res;
  } /* decNumberReduce */

/* ------------------------------------------------------------------ */
/* decNumberRescale -- force exponent to requested value              */
/*                                                                    */
/*   This computes C = op(A, B), where op adjusts the coefficient     */
/*   of C (by rounding or shifting) such that the exponent (-scale)   */
/*   of C has the value B.  The numerical value of C will equal A,    */
/*   except for the effects of any rounding that occurred.            */
/*                                                                    */
/*   res is C, the result.  C may be A or B                           */
/*   lhs is A, the number to adjust                                   */
/*   rhs is B, the requested exponent                                 */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* Unless there is an error or the result is infinite, the exponent   */
/* after the operation is guaranteed to be equal to B.                */
/* ------------------------------------------------------------------ */
decNumber * decNumberRescale(decNumber *res, const decNumber *lhs,
                             const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decQuantizeOp(res, lhs, rhs, set, 0, &status);
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberRescale */

/* ------------------------------------------------------------------ */
/* decNumberRemainder -- divide and return remainder                  */
/*                                                                    */
/*   This computes C = A % B                                          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs,
                               const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decDivideOp(res, lhs, rhs, set, REMAINDER, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberRemainder */

/* ------------------------------------------------------------------ */
/* decNumberRemainderNear -- divide and return remainder from nearest */
/*                                                                    */
/*   This computes C = A % B, where % is the IEEE remainder operator  */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
/*   lhs is A                                                         */
/*   rhs is B                                                         */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs,
                                   const decNumber *rhs, decContext *set) {
  uInt status=0;                        /* accumulator */
  decDivideOp(res, lhs, rhs, set, REMNEAR, &status);
  if (status!=0) decStatus(res, status, set);
  #if DECCHECK
  decCheckInexact(res, set);
  #endif
  return res;
  } /* decNumberRemainderNear */

/* ------------------------------------------------------------------ */
/* decNumberRotate -- rotate the coefficient of a Number left/right   */
/*                                                                    */
/*   This computes C = A rot B  (in base ten and rotating set->digits */
/*   digits).                                                         */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=XrotX)       */
/*   lhs is A                                                         */
/*   rhs is B, the number of digits to rotate (-ve to right)          */
/*   set is the context                                               */
/*                                                                    */
/* The digits of the coefficient of A are rotated to the left (if B   */
/* is positive) or to the right (if B is negative) without adjusting  */
/* the exponent or the sign of A.  If lhs->digits is less than        */
/* set->digits the coefficient is padded with zeros on the left       */
/* before the rotate.  Any leading zeros in the result are removed    */
/* as usual.                                                          */
/*                                                                    */
/* B must be an integer (q=0) and in the range -set->digits through   */
/* +set->digits.                                                      */
/* C must have space for set->digits digits.                          */
/* NaNs are propagated as usual.  Infinities are unaffected (but      */
/* B must be valid).  No status is set unless B is invalid or an      */
/* operand is an sNaN.                                                */
/* ------------------------------------------------------------------ */
decNumber * decNumberRotate(decNumber *res, const decNumber *lhs,
                           const decNumber *rhs, decContext *set) {
  uInt status=0;              /* accumulator */
  Int  rotate;                /* rhs as an Int */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  /* NaNs propagate as normal */
  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
    decNaNs(res, lhs, rhs, set, &status);
   /* rhs must be an integer */
   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
    status=DEC_Invalid_operation;
   else { /* both numeric, rhs is an integer */
    rotate=decGetInt(rhs);                   /* [cannot fail] */
    if (rotate==BADINT                       /* something bad .. */
     || rotate==BIGODD || rotate==BIGEVEN    /* .. very big .. */
     || abs(rotate)>set->digits)             /* .. or out of range */
      status=DEC_Invalid_operation;
     else {                                  /* rhs is OK */
      decNumberCopy(res, lhs);
      /* convert -ve rotate to equivalent positive rotation */
      if (rotate<0) rotate=set->digits+rotate;
      if (rotate!=0 && rotate!=set->digits   /* zero or full rotation */
       && !decNumberIsInfinite(res)) {       /* lhs was infinite */
        /* left-rotate to do; 0 < rotate < set->digits */
        uInt units, shift;                   /* work */
        uInt msudigits;                      /* digits in result msu */
        Unit *msu=res->lsu+D2U(res->digits)-1;    /* current msu */
        Unit *msumax=res->lsu+D2U(set->digits)-1; /* rotation msu */
        for (msu++; msu<=msumax; msu++) *msu=0;   /* ensure high units=0 */
        res->digits=set->digits;                  /* now full-length */
        msudigits=MSUDIGITS(res->digits);         /* actual digits in msu */

        /* rotation here is done in-place, in three steps */
        /* 1. shift all to least up to one unit to unit-align final */
        /*    lsd [any digits shifted out are rotated to the left, */
        /*    abutted to the original msd (which may require split)] */
        /* */
        /*    [if there are no whole units left to rotate, the */
        /*    rotation is now complete] */
        /* */
        /* 2. shift to least, from below the split point only, so that */
        /*    the final msd is in the right place in its Unit [any */
        /*    digits shifted out will fit exactly in the current msu, */
        /*    left aligned, no split required] */
        /* */
        /* 3. rotate all the units by reversing left part, right */
        /*    part, and then whole */
        /* */
        /* example: rotate right 8 digits (2 units + 2), DECDPUN=3. */
        /* */
        /*   start: 00a bcd efg hij klm npq */
        /* */
        /*      1a  000 0ab cde fgh|ijk lmn [pq saved] */
        /*      1b  00p qab cde fgh|ijk lmn */
        /* */
        /*      2a  00p qab cde fgh|00i jkl [mn saved] */
        /*      2b  mnp qab cde fgh|00i jkl */
        /* */
        /*      3a  fgh cde qab mnp|00i jkl */
        /*      3b  fgh cde qab mnp|jkl 00i */
        /*      3c  00i jkl mnp qab cde fgh */

        /* Step 1: amount to shift is the partial right-rotate count */
        rotate=set->digits-rotate;      /* make it right-rotate */
        units=rotate/DECDPUN;           /* whole units to rotate */
        shift=rotate%DECDPUN;           /* left-over digits count */
        if (shift>0) {                  /* not an exact number of units */
          uInt save=res->lsu[0]%powers[shift];    /* save low digit(s) */
          decShiftToLeast(res->lsu, D2U(res->digits), shift);
          if (shift>msudigits) {        /* msumax-1 needs >0 digits */
            uInt rem=save%powers[shift-msudigits];/* split save */
            *msumax=(Unit)(save/powers[shift-msudigits]); /* and insert */
            *(msumax-1)=*(msumax-1)
                       +(Unit)(rem*powers[DECDPUN-(shift-msudigits)]); /* .. */
            }
           else { /* all fits in msumax */
            *msumax=*msumax+(Unit)(save*powers[msudigits-shift]); /* [maybe *1] */
            }
          } /* digits shift needed */

        /* If whole units to rotate... */
        if (units>0) {                  /* some to do */
          /* Step 2: the units to touch are the whole ones in rotate, */
          /*   if any, and the shift is DECDPUN-msudigits (which may be */
          /*   0, again) */
          shift=DECDPUN-msudigits;
          if (shift>0) {                /* not an exact number of units */
            uInt save=res->lsu[0]%powers[shift];  /* save low digit(s) */
            decShiftToLeast(res->lsu, units, shift);
            *msumax=*msumax+(Unit)(save*powers[msudigits]);
            } /* partial shift needed */

          /* Step 3: rotate the units array using triple reverse */
          /* (reversing is easy and fast) */
          decReverse(res->lsu+units, msumax);     /* left part */
          decReverse(res->lsu, res->lsu+units-1); /* right part */
          decReverse(res->lsu, msumax);           /* whole */
          } /* whole units to rotate */
        /* the rotation may have left an undetermined number of zeros */
        /* on the left, so true length needs to be calculated */
        res->digits=decGetDigits(res->lsu, msumax-res->lsu+1);
        } /* rotate needed */
      } /* rhs OK */
    } /* numerics */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberRotate */

/* ------------------------------------------------------------------ */
/* decNumberSameQuantum -- test for equal exponents                   */
/*                                                                    */
/*   res is the result number, which will contain either 0 or 1       */
/*   lhs is a number to test                                          */
/*   rhs is the second (usually a pattern)                            */
/*                                                                    */
/* No errors are possible and no context is needed.                   */
/* ------------------------------------------------------------------ */
decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs,
                                 const decNumber *rhs) {
  Unit ret=0;                      /* return value */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, DECUNCONT)) return res;
  #endif

  if (SPECIALARGS) {
    if (decNumberIsNaN(lhs) && decNumberIsNaN(rhs)) ret=1;
     else if (decNumberIsInfinite(lhs) && decNumberIsInfinite(rhs)) ret=1;
     /* [anything else with a special gives 0] */
    }
   else if (lhs->exponent==rhs->exponent) ret=1;

  decNumberZero(res);              /* OK to overwrite an operand now */
  *res->lsu=ret;
  return res;
  } /* decNumberSameQuantum */

/* ------------------------------------------------------------------ */
/* decNumberScaleB -- multiply by a power of 10                       */
/*                                                                    */
/* This computes C = A x 10**B where B is an integer (q=0) with       */
/* maximum magnitude 2*(emax+digits)                                  */
/*                                                                    */
/*   res is C, the result.  C may be A or B                           */
/*   lhs is A, the number to adjust                                   */
/*   rhs is B, the requested power of ten to use                      */
/*   set is the context                                               */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/*                                                                    */
/* The result may underflow or overflow.                              */
/* ------------------------------------------------------------------ */
decNumber * decNumberScaleB(decNumber *res, const decNumber *lhs,
                            const decNumber *rhs, decContext *set) {
  Int  reqexp;                /* requested exponent change [B] */
  uInt status=0;              /* accumulator */
  Int  residue;               /* work */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  /* Handle special values except lhs infinite */
  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
    decNaNs(res, lhs, rhs, set, &status);
    /* rhs must be an integer */
   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
    status=DEC_Invalid_operation;
   else {
    /* lhs is a number; rhs is a finite with q==0 */
    reqexp=decGetInt(rhs);                   /* [cannot fail] */
    if (reqexp==BADINT                       /* something bad .. */
     || reqexp==BIGODD || reqexp==BIGEVEN    /* .. very big .. */
     || abs(reqexp)>(2*(set->digits+set->emax))) /* .. or out of range */
      status=DEC_Invalid_operation;
     else {                                  /* rhs is OK */
      decNumberCopy(res, lhs);               /* all done if infinite lhs */
      if (!decNumberIsInfinite(res)) {       /* prepare to scale */
        res->exponent+=reqexp;               /* adjust the exponent */
        residue=0;
        decFinalize(res, set, &residue, &status); /* .. and check */
        } /* finite LHS */
      } /* rhs OK */
    } /* rhs finite */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberScaleB */

/* ------------------------------------------------------------------ */
/* decNumberShift -- shift the coefficient of a Number left or right  */
/*                                                                    */
/*   This computes C = A << B or C = A >> -B  (in base ten).          */
/*                                                                    */
/*   res is C, the result.  C may be A and/or B (e.g., X=X<<X)        */
/*   lhs is A                                                         */
/*   rhs is B, the number of digits to shift (-ve to right)           */
/*   set is the context                                               */
/*                                                                    */
/* The digits of the coefficient of A are shifted to the left (if B   */
/* is positive) or to the right (if B is negative) without adjusting  */
/* the exponent or the sign of A.                                     */
/*                                                                    */
/* B must be an integer (q=0) and in the range -set->digits through   */
/* +set->digits.                                                      */
/* C must have space for set->digits digits.                          */
/* NaNs are propagated as usual.  Infinities are unaffected (but      */
/* B must be valid).  No status is set unless B is invalid or an      */
/* operand is an sNaN.                                                */
/* ------------------------------------------------------------------ */
decNumber * decNumberShift(decNumber *res, const decNumber *lhs,
                           const decNumber *rhs, decContext *set) {
  uInt status=0;              /* accumulator */
  Int  shift;                 /* rhs as an Int */

  #if DECCHECK
  if (decCheckOperands(res, lhs, rhs, set)) return res;
  #endif

  /* NaNs propagate as normal */
  if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
    decNaNs(res, lhs, rhs, set, &status);
   /* rhs must be an integer */
   else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
    status=DEC_Invalid_operation;
   else { /* both numeric, rhs is an integer */
    shift=decGetInt(rhs);                    /* [cannot fail] */
    if (shift==BADINT                        /* something bad .. */
     || shift==BIGODD || shift==BIGEVEN      /* .. very big .. */
     || abs(shift)>set->digits)              /* .. or out of range */
      status=DEC_Invalid_operation;
     else {                                  /* rhs is OK */
      decNumberCopy(res, lhs);
      if (shift!=0 && !decNumberIsInfinite(res)) { /* something to do */
        if (shift>0) {                       /* to left */
          if (shift==set->digits) {          /* removing all */
            *res->lsu=0;                     /* so place 0 */
            res->digits=1;                   /* .. */
            }
           else {                            /* */
            /* first remove leading digits if necessary */
            if (res->digits+shift>set->digits) {
              decDecap(res, res->digits+shift-set->digits);
              /* that updated res->digits; may have gone to 1 (for a */
              /* single digit or for zero */
              }
            if (res->digits>1 || *res->lsu)  /* if non-zero.. */
              res->digits=decShiftToMost(res->lsu, res->digits, shift);
            } /* partial left */
          } /* left */
         else { /* to right */
          if (-shift>=res->digits) {         /* discarding all */
            *res->lsu=0;                     /* so place 0 */
            res->digits=1;                   /* .. */
            }
           else {
            decShiftToLeast(res->lsu, D2U(res->digits), -shift);
            res->digits-=(-shift);
            }
          } /* to right */
        } /* non-0 non-Inf shift */
      } /* rhs OK */
    } /* numerics */
  if (status!=0) decStatus(res, status, set);
  return res;
  } /* decNumberShift */

/* ------------------------------------------------------------------ */
/* decNumberSquareRoot -- square root operator                        */
/*                                                                    */
/*   This computes C = squareroot(A)                                  */
/*                                                                    */
/*   res is C, the result.  C may be A                                */
/*   rhs is A                                                         */
/*   set is the context; note that rounding mode has no effect        */
/*                                                                    */
/* C must have space for set->digits digits.                          */
/* ------------------------------------------------------------------ */
/* This uses the following varying-precision algorithm in:            */
/*                                                                    */
/*   Properly Rounded Variable Precision Square Root, T. E. Hull and  */
/*   A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
/*   pp229-237, ACM, September 1985.                                  */
/*                                                                    */
/* The square-root is calculated using Newton's method, after which   */
/* a check is made to ensure the result is correctly rounded.         */
/*                                                                    */
/* % [Reformatted original Numerical Turing source code follows.]     */
/* function sqrt(x : real) : real                                     */
/* % sqrt(x) returns the properly rounded approximation to the square */
/* % root of x, in the precision of the calling environment, or it    */
/* % fails if x < 0.                                                  */
/* % t e hull and a abrham, august, 1984                              */
/* if x <= 0 then                                                     */
/*   if x < 0 then                                                    */
/*     assert false                                                   */
/*   else                                                             */
/*     result 0                                                       */
/*   end if                                                           */
/* end if                                                             */
/* var f := setexp(x, 0)  % fraction part of x   [0.1 <= x < 1]       */
/* var e := getexp(x)     % exponent part of x                        */
/* var approx : real                                                  */
/* if e mod 2 = 0  then                                               */
/*   approx := .259 + .819 * f   % approx to root of f                */
/* else                                                               */
/*   f := f/l0                   % adjustments                        */
/*   e := e + 1                  %   for odd                          */
/*   approx := .0819 + 2.59 * f  %   exponent                         */
/* end if                                                             */
/*                                                                    */
/* var p:= 3                                                          */
/* const maxp := currentprecision + 2                                 */
/* loop                                                               */
/*   p := min(2*p - 2, maxp)     % p = 4,6,10, . . . , maxp           */
/*   precision p                                                      */
/*   approx := .5 * (approx + f/approx)                               */
/*   exit when p = maxp                                               */
/* end loop                                                           */
/*                                                                    */
/* % approx is now within 1 ulp of the properly rounded square root   */
/* % of f; to ensure proper rounding, compare squares of (approx -    */
/* % l/2 ulp) and (approx + l/2 ulp) with f.                          */
/* p := currentprecision                                              */
/* begin                                                              */
/*   precision p + 2                                                  */
/*   const approxsubhalf := approx - setexp(.5, -p)                   */
/*   if mulru(approxsubhalf, approxsubhalf) > f then                  */
/*     approx := approx - setexp(.l, -p + 1)                          */
/*   else                                                             */
/*     const approxaddhalf := approx + setexp(.5, -p)                 */
/*     if mulrd(approxaddhalf, approxaddhalf) < f then                */
/*       approx := approx + setexp(.l, -p + 1)                        */
/*     end if                                                         */
/*   end if                                                           */
/* end                                                                */
/* result setexp(approx, e div 2)  % fix exponent                     */
/* end sqrt                                                           */
/* ------------------------------------------------------------------ */
decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs,
                                decContext *set) {
  decContext workset, approxset;   /* work contexts */
  decNumber dzero;                 /* used for constant zero */
  Int  maxp;                       /* largest working precision */
  Int  workp;                      /* working precision */
  Int  residue=0;                  /* rounding residue */
  uInt status=0, ignore=0;         /* status accumulators */
  uInt rstatus;                    /* .. */
  Int  exp;                        /* working exponent */
  Int  ideal;                      /* ideal (preferred) exponent */
  Int  needbytes;                  /* work */
  Int  dropped;                    /* .. */

  #if DECSUBSET
  decNumber *allocrhs=NULL;        /* non-NULL if rounded rhs allocated */
  #endif
  /* buffer for f [needs +1 in case DECBUFFER 0] */
  decNumber buff[D2N(DECBUFFER+1)];
  /* buffer for a [needs +2 to match likely maxp] */
  decNumber bufa[D2N(DECBUFFER+2)];
  /* buffer for temporary, b [must be same size as a] */
  decNumber bufb[D2N(DECBUFFER+2)];
  decNumber *allocbuff=NULL;       /* -> allocated buff, iff allocated */
  decNumber *allocbufa=NULL;       /* -> allocated bufa, iff allocated */
  decNumber *allocbufb=NULL;       /* -> allocated bufb, iff allocated */
  decNumber *f=buff;               /* reduced fraction */
  decNumber *a=bufa;               /* approximation to result */
  decNumber *b=bufb;               /* intermediate result */
  /* buffer for temporary variable, up to 3 digits */
  decNumber buft[D2N(3)];
  decNumber *t=buft;               /* up-to-3-digit constant or work */

  #if DECCHECK
  if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
  #endif

  do {                             /* protect allocated storage */
    #if DECSUBSET
    if (!set->extended) {
      /* reduce operand and set lostDigits status, as needed */
      if (rhs->digits>set->digits) {
        allocrhs=decRoundOperand(rhs, set, &status);
        if (allocrhs==NULL) break;
        /* [Note: 'f' allocation below could reuse this buffer if */
        /* used, but as this is rare they are kept separate for clarity.] */
        rhs=allocrhs;
        }
      }
    #endif
    /* [following code does not require input rounding] */

    /* handle infinities and NaNs */
    if (SPECIALARG) {
      if (decNumberIsInfinite(rhs)) {         /* an infinity */
        if (decNumberIsNegative(rhs)) status|=DEC_Invalid_operation;
         else decNumberCopy(res, rhs);        /* +Infinity */
        }
       else decNaNs(res, rhs, NULL, set, &status); /* a NaN */
      break;
      }

    /* calculate the ideal (preferred) exponent [floor(exp/2)] */
    /* [We would like to write: ideal=rhs->exponent>>1, but this */
    /* generates a compiler warning.  Generated code is the same.] */
    ideal=(rhs->exponent&~1)/2;         /* target */

    /* handle zeros */
    if (ISZERO(rhs)) {
      decNumberCopy(res, rhs);          /* could be 0 or -0 */
      res->exponent=ideal;              /* use the ideal [safe] */
      /* use decFinish to clamp any out-of-range exponent, etc. */
      decFinish(res, set, &residue, &status);
      break;
      }

    /* any other -x is an oops */
    if (decNumberIsNegative(rhs)) {
      status|=DEC_Invalid_operation;
      break;
      }

    /* space is needed for three working variables */
    /*   f -- the same precision as the RHS, reduced to 0.01->0.99... */
    /*   a -- Hull's approximation -- precision, when assigned, is */
    /*        currentprecision+1 or the input argument precision, */
    /*        whichever is larger (+2 for use as temporary) */
    /*   b -- intermediate temporary result (same size as a) */
    /* if any is too long for local storage, then allocate */
    workp=MAXI(set->digits+1, rhs->digits);  /* actual rounding precision */
    maxp=workp+2;                            /* largest working precision */

    needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
    if (needbytes>(Int)sizeof(buff)) {
      allocbuff=(decNumber *)malloc(needbytes);
      if (allocbuff==NULL) {  /* hopeless -- abandon */
        status|=DEC_Insufficient_storage;
        break;}
      f=allocbuff;            /* use the allocated space */
      }
    /* a and b both need to be able to hold a maxp-length number */
    needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit);
    if (needbytes>(Int)sizeof(bufa)) {            /* [same applies to b] */
      allocbufa=(decNumber *)malloc(needbytes);
      allocbufb=(decNumber *)malloc(needbytes);
      if (allocbufa==NULL || allocbufb==NULL) {   /* hopeless */
        status|=DEC_Insufficient_storage;
        break;}
      a=allocbufa;            /* use the allocated spaces */
      b=allocbufb;            /* .. */
      }

    /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
    decNumberCopy(f, rhs);
    exp=f->exponent+f->digits;               /* adjusted to Hull rules */
    f->exponent=-(f->digits);                /* to range */

    /* set up working context */
    decContextDefault(&workset, DEC_INIT_DECIMAL64);

    /* [Until further notice, no error is possible and status bits */
    /* (Rounded, etc.) should be ignored, not accumulated.] */

    /* Calculate initial approximation, and allow for odd exponent */
    workset.digits=workp;                    /* p for initial calculation */
    t->bits=0; t->digits=3;
    a->bits=0; a->digits=3;
    if ((exp & 1)==0) {                      /* even exponent */
      /* Set t=0.259, a=0.819 */
      t->exponent=-3;
      a->exponent=-3;
      #if DECDPUN>=3
        t->lsu[0]=259;
        a->lsu[0]=819;
      #elif DECDPUN==2
        t->lsu[0]=59; t->lsu[1]=2;
        a->lsu[0]=19; a->lsu[1]=8;
      #else
        t->lsu[0]=9; t->lsu[1]=5; t->lsu[2]=2;
        a->lsu[0]=9; a->lsu[1]=1; a->lsu[2]=8;
      #endif
      }
     else {                                  /* odd exponent */
      /* Set t=0.0819, a=2.59 */
      f->exponent--;                         /* f=f/10 */
      exp++;                                 /* e=e+1 */
      t->exponent=-4;
      a->exponent=-2;
      #if DECDPUN>=3
        t->lsu[0]=819;
        a->lsu[0]=259;
      #elif DECDPUN==2
        t->lsu[0]=19; t->lsu[1]=8;
        a->lsu[0]=59; a->lsu[1]=2;
      #else
        t->lsu[0]=9; t->lsu[1]=1; t->lsu[2]=8;
        a->lsu[0]=9; a->lsu[1]=5; a->lsu[2]=2;
      #endif
      }
    decMultiplyOp(a, a, f, &workset, &ignore);    /* a=a*f */
    decAddOp(a, a, t, &workset, 0, &ignore);      /* ..+t */
    /* [a is now the initial approximation for sqrt(f), calculated with */
    /* currentprecision, which is also a's precision.] */

    /* the main calculation loop */
    decNumberZero(&dzero);                   /* make 0 */
    decNumberZero(t);                        /* set t = 0.5 */
    t->lsu[0]=5;                             /* .. */
    t->exponent=-1;                          /* .. */
    workset.digits=3;                        /* initial p */
    for (;;) {
      /* set p to min(2*p - 2, maxp)  [hence 3; or: 4, 6, 10, ... , maxp] */
      workset.digits=workset.digits*2-2;
      if (workset.digits>maxp) workset.digits=maxp;
      /* a = 0.5 * (a + f/a) */
      /* [calculated at p then rounded to currentprecision] */
      decDivideOp(b, f, a, &workset, DIVIDE, &ignore); /* b=f/a */
      decAddOp(b, b, a, &workset, 0, &ignore);    /* b=b+a */
      decMultiplyOp(a, b, t, &workset, &ignore);  /* a=b*0.5 */
      if (a->digits==maxp) break;            /* have required digits */
      } /* loop */

    /* Here, 0.1 <= a < 1 [Hull], and a has maxp digits */
    /* now reduce to length, etc.; this needs to be done with a */
    /* having the correct exponent so as to handle subnormals */
    /* correctly */
    approxset=*set;                          /* get emin, emax, etc. */
    approxset.round=DEC_ROUND_HALF_EVEN;
    a->exponent+=exp/2;                      /* set correct exponent */

    rstatus=0;                               /* clear status */
    residue=0;                               /* .. and accumulator */
    decCopyFit(a, a, &approxset, &residue, &rstatus);  /* reduce (if needed) */
    decFinish(a, &approxset, &residue, &rstatus);      /* clean and finalize */

    /* Overflow was possible if the input exponent was out-of-range, */
    /* in which case quit */
    if (rstatus&DEC_Overflow) {
      status=rstatus;                        /* use the status as-is */
      decNumberCopy(res, a);                 /* copy to result */
      break;
      }

    /* Preserve status except Inexact/Rounded */
    status|=(rstatus & ~(DEC_Rounded|DEC_Inexact));

    /* Carry out the Hull correction */
    a->exponent-=exp/2;                      /* back to 0.1->1 */

    /* a is now at final precision and within 1 ulp of the properly */
    /* rounded square root of f; to ensure proper rounding, compare */
    /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
    /* Here workset.digits=maxp and t=0.5, and a->digits determines */
    /* the ulp */
    workset.digits--;                             /* maxp-1 is OK now */
    t->exponent=-a->digits-1;                     /* make 0.5 ulp */
    decAddOp(b, a, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */
    workset.round=DEC_ROUND_UP;
    decMultiplyOp(b, b, b, &work