OSDN Git Service

* libdecnumber: Import decNumber sources from the dfp-branch.
authorbje <bje@138bc75d-0d04-0410-961f-82ee72b054a4>
Mon, 28 Nov 2005 22:30:30 +0000 (22:30 +0000)
committerbje <bje@138bc75d-0d04-0410-961f-82ee72b054a4>
Mon, 28 Nov 2005 22:30:30 +0000 (22:30 +0000)
git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@107629 138bc75d-0d04-0410-961f-82ee72b054a4

18 files changed:
ChangeLog
libdecnumber/ChangeLog [new file with mode: 0644]
libdecnumber/decContext.c [new file with mode: 0644]
libdecnumber/decContext.h [new file with mode: 0644]
libdecnumber/decDPD.h [new file with mode: 0644]
libdecnumber/decLibrary.c [new file with mode: 0644]
libdecnumber/decNumber.c [new file with mode: 0644]
libdecnumber/decNumber.h [new file with mode: 0644]
libdecnumber/decNumberLocal.h [new file with mode: 0644]
libdecnumber/decRound.c [new file with mode: 0644]
libdecnumber/decUtility.c [new file with mode: 0644]
libdecnumber/decUtility.h [new file with mode: 0644]
libdecnumber/decimal128.c [new file with mode: 0644]
libdecnumber/decimal128.h [new file with mode: 0644]
libdecnumber/decimal32.c [new file with mode: 0644]
libdecnumber/decimal32.h [new file with mode: 0644]
libdecnumber/decimal64.c [new file with mode: 0644]
libdecnumber/decimal64.h [new file with mode: 0644]

index 584ca73..00518f2 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2005-11-29  Ben Elliston  <bje@au.ibm.com>
+
+       * libdecnumber: Import decNumber sources from the dfp-branch.
+
 2005-11-21  Kean Johnston  <jkj@sco.com>
 
        * config.sub, config.guess: Sync from upstream sources.
diff --git a/libdecnumber/ChangeLog b/libdecnumber/ChangeLog
new file mode 100644 (file)
index 0000000..64730ae
--- /dev/null
@@ -0,0 +1,9 @@
+2005-11-29  Ben Elliston  <bje@au.ibm.com>
+
+       * decimal32.h, decimal64.h, decimal128.h: New.
+        * decimal32.c, decimal64.c, decimal128.c: Likewise.
+       * decContext.c, decContext.h: Likewise.
+       * decUtility.c, decUtility.h: Likewise.
+       * decNumber.c, decNumber.h, decNumberLocal.h: Likewise.
+       * decDPD.h: Likewise.
+       * decLibrary.c, decRound.c: Likewise.
diff --git a/libdecnumber/decContext.c b/libdecnumber/decContext.c
new file mode 100644 (file)
index 0000000..26a7f3b
--- /dev/null
@@ -0,0 +1,218 @@
+/* Decimal context module for the decNumber C Library.
+   Copyright (C) 2005 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.
+
+   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, 59 Temple Place - Suite 330, Boston, MA
+   02111-1307, USA.  */
+
+/*  This module compirises the routines for handling the arithmetic
+    context structures. */
+
+#include <string.h>            /* for strcmp */
+#include "decContext.h"                /* context and base types */
+#include "decNumberLocal.h"    /* decNumber local types, etc. */
+
+/* ------------------------------------------------------------------ */
+/* decContextDefault -- initialize a context structure                */
+/*                                                                    */
+/*  context is the structure to be initialized                        */
+/*  kind selects the required set of default values, one of:          */
+/*      DEC_INIT_BASE       -- select ANSI X3-274 defaults            */
+/*      DEC_INIT_DECIMAL32  -- select IEEE 754r defaults, 32-bit      */
+/*      DEC_INIT_DECIMAL64  -- select IEEE 754r defaults, 64-bit      */
+/*      DEC_INIT_DECIMAL128 -- select IEEE 754r defaults, 128-bit     */
+/*      For any other value a valid context is returned, but with     */
+/*      Invalid_operation set in the status field.                    */
+/*  returns a context structure with the appropriate initial values.  */
+/* ------------------------------------------------------------------ */
+decContext *
+decContextDefault (decContext * context, Int kind)
+{
+  /* set defaults... */
+  context->digits = 9;         /* 9 digits */
+  context->emax = DEC_MAX_EMAX;        /* 9-digit exponents */
+  context->emin = DEC_MIN_EMIN;        /* .. balanced */
+  context->round = DEC_ROUND_HALF_UP;  /* 0.5 rises */
+  context->traps = DEC_Errors; /* all but informational */
+  context->status = 0;         /* cleared */
+  context->clamp = 0;          /* no clamping */
+#if DECSUBSET
+  context->extended = 0;       /* cleared */
+#endif
+  switch (kind)
+    {
+    case DEC_INIT_BASE:
+      /* [use defaults] */
+      break;
+    case DEC_INIT_DECIMAL32:
+      context->digits = 7;     /* digits */
+      context->emax = 96;      /* Emax */
+      context->emin = -95;     /* Emin */
+      context->round = DEC_ROUND_HALF_EVEN;    /* 0.5 to nearest even */
+      context->traps = 0;      /* no traps set */
+      context->clamp = 1;      /* clamp exponents */
+#if DECSUBSET
+      context->extended = 1;   /* set */
+#endif
+      break;
+    case DEC_INIT_DECIMAL64:
+      context->digits = 16;    /* digits */
+      context->emax = 384;     /* Emax */
+      context->emin = -383;    /* Emin */
+      context->round = DEC_ROUND_HALF_EVEN;    /* 0.5 to nearest even */
+      context->traps = 0;      /* no traps set */
+      context->clamp = 1;      /* clamp exponents */
+#if DECSUBSET
+      context->extended = 1;   /* set */
+#endif
+      break;
+    case DEC_INIT_DECIMAL128:
+      context->digits = 34;    /* digits */
+      context->emax = 6144;    /* Emax */
+      context->emin = -6143;   /* Emin */
+      context->round = DEC_ROUND_HALF_EVEN;    /* 0.5 to nearest even */
+      context->traps = 0;      /* no traps set */
+      context->clamp = 1;      /* clamp exponents */
+#if DECSUBSET
+      context->extended = 1;   /* set */
+#endif
+      break;
+
+    default:                   /* invalid Kind */
+      /* use defaults, and .. */
+      decContextSetStatus (context, DEC_Invalid_operation);    /* trap */
+    }
+  return context;
+}                              /* decContextDefault */
+
+/* ------------------------------------------------------------------ */
+/* decContextStatusToString -- convert status flags to a string       */
+/*                                                                    */
+/*  context is a context with valid status field                      */
+/*                                                                    */
+/*  returns a constant string describing the condition.  If multiple  */
+/*    (or no) flags are set, a generic constant message is returned.  */
+/* ------------------------------------------------------------------ */
+const char *
+decContextStatusToString (decContext * context)
+{
+  Int status = context->status;
+  if (status == DEC_Conversion_syntax)
+    return DEC_Condition_CS;
+  if (status == DEC_Division_by_zero)
+    return DEC_Condition_DZ;
+  if (status == DEC_Division_impossible)
+    return DEC_Condition_DI;
+  if (status == DEC_Division_undefined)
+    return DEC_Condition_DU;
+  if (status == DEC_Inexact)
+    return DEC_Condition_IE;
+  if (status == DEC_Insufficient_storage)
+    return DEC_Condition_IS;
+  if (status == DEC_Invalid_context)
+    return DEC_Condition_IC;
+  if (status == DEC_Invalid_operation)
+    return DEC_Condition_IO;
+#if DECSUBSET
+  if (status == DEC_Lost_digits)
+    return DEC_Condition_LD;
+#endif
+  if (status == DEC_Overflow)
+    return DEC_Condition_OV;
+  if (status == DEC_Clamped)
+    return DEC_Condition_PA;
+  if (status == DEC_Rounded)
+    return DEC_Condition_RO;
+  if (status == DEC_Subnormal)
+    return DEC_Condition_SU;
+  if (status == DEC_Underflow)
+    return DEC_Condition_UN;
+  if (status == 0)
+    return DEC_Condition_ZE;
+  return DEC_Condition_MU;     /* Multiple errors */
+}                              /* decContextStatusToString */
+
+/* ------------------------------------------------------------------ */
+/* decContextSetStatusFromString -- set status from a string          */
+/*                                                                    */
+/*  context is the controlling context                                */
+/*  string is a string exactly equal to one that might be returned    */
+/*            by decContextStatusToString                             */
+/*                                                                    */
+/*  The status bit corresponding to the string is set, and a trap     */
+/*  is raised if appropriate.                                         */
+/*                                                                    */
+/*  returns the context structure, unless the string is equal to      */
+/*    DEC_Condition_MU or is not recognized.  In these cases NULL is  */
+/*    returned.                                                       */
+/* ------------------------------------------------------------------ */
+decContext *
+decContextSetStatusFromString (decContext * context, char *string)
+{
+  if (strcmp (string, DEC_Condition_CS) == 0)
+    return decContextSetStatus (context, DEC_Conversion_syntax);
+  if (strcmp (string, DEC_Condition_DZ) == 0)
+    return decContextSetStatus (context, DEC_Division_by_zero);
+  if (strcmp (string, DEC_Condition_DI) == 0)
+    return decContextSetStatus (context, DEC_Division_impossible);
+  if (strcmp (string, DEC_Condition_DU) == 0)
+    return decContextSetStatus (context, DEC_Division_undefined);
+  if (strcmp (string, DEC_Condition_IE) == 0)
+    return decContextSetStatus (context, DEC_Inexact);
+  if (strcmp (string, DEC_Condition_IS) == 0)
+    return decContextSetStatus (context, DEC_Insufficient_storage);
+  if (strcmp (string, DEC_Condition_IC) == 0)
+    return decContextSetStatus (context, DEC_Invalid_context);
+  if (strcmp (string, DEC_Condition_IO) == 0)
+    return decContextSetStatus (context, DEC_Invalid_operation);
+#if DECSUBSET
+  if (strcmp (string, DEC_Condition_LD) == 0)
+    return decContextSetStatus (context, DEC_Lost_digits);
+#endif
+  if (strcmp (string, DEC_Condition_OV) == 0)
+    return decContextSetStatus (context, DEC_Overflow);
+  if (strcmp (string, DEC_Condition_PA) == 0)
+    return decContextSetStatus (context, DEC_Clamped);
+  if (strcmp (string, DEC_Condition_RO) == 0)
+    return decContextSetStatus (context, DEC_Rounded);
+  if (strcmp (string, DEC_Condition_SU) == 0)
+    return decContextSetStatus (context, DEC_Subnormal);
+  if (strcmp (string, DEC_Condition_UN) == 0)
+    return decContextSetStatus (context, DEC_Underflow);
+  if (strcmp (string, DEC_Condition_ZE) == 0)
+    return context;
+  return NULL;                 /* Multiple status, or unknown */
+}                              /* decContextSetStatusFromString */
+
+/* ------------------------------------------------------------------ */
+/* decContextSetStatus -- set status and raise trap if appropriate    */
+/*                                                                    */
+/*  context is the controlling context                                */
+/*  status  is the DEC_ exception code                                */
+/*  returns the context structure                                     */
+/*                                                                    */
+/* Control may never return from this routine, if there is a signal   */
+/* handler and it takes a long jump.                                  */
+/* ------------------------------------------------------------------ */
+decContext *
+decContextSetStatus (decContext * context, uInt status)
+{
+  context->status |= status;
+  if (status & context->traps)
+    raise (SIGFPE);
+  return context;
+}                              /* decContextSetStatus */
diff --git a/libdecnumber/decContext.h b/libdecnumber/decContext.h
new file mode 100644 (file)
index 0000000..d011f4f
--- /dev/null
@@ -0,0 +1,178 @@
+/* Decimal Context module header for the decNumber C Library
+   Copyright (C) 2005 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.
+
+   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, 59 Temple Place - Suite 330, Boston, MA
+   02111-1307, USA.  */
+
+/* ------------------------------------------------------------------ */
+/*                                                                    */
+/* Context must always be set correctly:                              */
+/*                                                                    */
+/*  digits   -- must be in the range 1 through 999999999              */
+/*  emax     -- must be in the range 0 through 999999999              */
+/*  emin     -- must be in the range 0 through -999999999             */
+/*  round    -- must be one of the enumerated rounding modes          */
+/*  traps    -- only defined bits may be set                          */
+/*  status   -- [any bits may be cleared, but not set, by user]       */
+/*  clamp    -- must be either 0 or 1                                 */
+/*  extended -- must be either 0 or 1 [present only if DECSUBSET]     */
+/*                                                                    */
+/* ------------------------------------------------------------------ */
+
+#if !defined(DECCONTEXT)
+#define DECCONTEXT
+#define DECCNAME     "decContext"      /* Short name */
+#define DECCFULLNAME "Decimal Context Descriptor"      /* Verbose name */
+#define DECCAUTHOR   "Mike Cowlishaw"  /* Who to blame */
+
+#include <stdint.h>            /* C99 standard integers */
+#include <signal.h>            /* for traps */
+
+
+  /* Conditional code flag -- set this to 0 for best performance */
+#define DECSUBSET 0            /* 1 to enable subset arithmetic */
+
+  /* Context for operations, with associated constants */
+enum rounding
+{
+  DEC_ROUND_CEILING,           /* round towards +infinity */
+  DEC_ROUND_UP,                        /* round away from 0 */
+  DEC_ROUND_HALF_UP,           /* 0.5 rounds up */
+  DEC_ROUND_HALF_EVEN,         /* 0.5 rounds to nearest even */
+  DEC_ROUND_HALF_DOWN,         /* 0.5 rounds down */
+  DEC_ROUND_DOWN,              /* round towards 0 (truncate) */
+  DEC_ROUND_FLOOR,             /* round towards -infinity */
+  DEC_ROUND_MAX                        /* enum must be less than this */
+};
+
+typedef struct
+{
+  int32_t digits;              /* working precision */
+  int32_t emax;                        /* maximum positive exponent */
+  int32_t emin;                        /* minimum negative exponent */
+  enum rounding round;         /* rounding mode */
+  uint32_t traps;              /* trap-enabler flags */
+  uint32_t status;             /* status flags */
+  uint8_t clamp;               /* flag: apply IEEE exponent clamp */
+#if DECSUBSET
+  uint8_t extended;            /* flag: special-values allowed */
+#endif
+} decContext;
+
+  /* Maxima and Minima */
+#define DEC_MAX_DIGITS 999999999
+#define DEC_MIN_DIGITS         1
+#define DEC_MAX_EMAX   999999999
+#define DEC_MIN_EMAX           0
+#define DEC_MAX_EMIN           0
+#define DEC_MIN_EMIN  -999999999
+
+  /* Trap-enabler and Status flags (exceptional conditions), and their names */
+  /* Top byte is reserved for internal use */
+#define DEC_Conversion_syntax    0x00000001
+#define DEC_Division_by_zero     0x00000002
+#define DEC_Division_impossible  0x00000004
+#define DEC_Division_undefined   0x00000008
+#define DEC_Insufficient_storage 0x00000010    /* [used if malloc fails] */
+#define DEC_Inexact              0x00000020
+#define DEC_Invalid_context      0x00000040
+#define DEC_Invalid_operation    0x00000080
+#if DECSUBSET
+#define DEC_Lost_digits          0x00000100
+#endif
+#define DEC_Overflow             0x00000200
+#define DEC_Clamped              0x00000400
+#define DEC_Rounded              0x00000800
+#define DEC_Subnormal            0x00001000
+#define DEC_Underflow            0x00002000
+
+  /* IEEE 854 groupings for the flags */
+  /* [DEC_Clamped, DEC_Lost_digits, DEC_Rounded, and DEC_Subnormal are */
+  /* not in IEEE 854] */
+#define DEC_IEEE_854_Division_by_zero  (DEC_Division_by_zero)
+#if DECSUBSET
+#define DEC_IEEE_854_Inexact           (DEC_Inexact | DEC_Lost_digits)
+#else
+#define DEC_IEEE_854_Inexact           (DEC_Inexact)
+#endif
+#define DEC_IEEE_854_Invalid_operation (DEC_Conversion_syntax |     \
+                                          DEC_Division_impossible |   \
+                                          DEC_Division_undefined |    \
+                                          DEC_Insufficient_storage |  \
+                                          DEC_Invalid_context |       \
+                                          DEC_Invalid_operation)
+#define DEC_IEEE_854_Overflow          (DEC_Overflow)
+#define DEC_IEEE_854_Underflow         (DEC_Underflow)
+
+  /* flags which are normally errors (results are qNaN, infinite, or 0) */
+#define DEC_Errors (DEC_IEEE_854_Division_by_zero |                 \
+                      DEC_IEEE_854_Invalid_operation |                \
+                      DEC_IEEE_854_Overflow | DEC_IEEE_854_Underflow)
+  /* flags which cause a result to become qNaN */
+#define DEC_NaNs    DEC_IEEE_854_Invalid_operation
+
+  /* flags which are normally for information only (have finite results) */
+#if DECSUBSET
+#define DEC_Information (DEC_Clamped | DEC_Rounded | DEC_Inexact     \
+                          | DEC_Lost_digits)
+#else
+#define DEC_Information (DEC_Clamped | DEC_Rounded | DEC_Inexact)
+#endif
+
+  /* name strings for the exceptional conditions */
+
+#define DEC_Condition_CS "Conversion syntax"
+#define DEC_Condition_DZ "Division by zero"
+#define DEC_Condition_DI "Division impossible"
+#define DEC_Condition_DU "Division undefined"
+#define DEC_Condition_IE "Inexact"
+#define DEC_Condition_IS "Insufficient storage"
+#define DEC_Condition_IC "Invalid context"
+#define DEC_Condition_IO "Invalid operation"
+#if DECSUBSET
+#define DEC_Condition_LD "Lost digits"
+#endif
+#define DEC_Condition_OV "Overflow"
+#define DEC_Condition_PA "Clamped"
+#define DEC_Condition_RO "Rounded"
+#define DEC_Condition_SU "Subnormal"
+#define DEC_Condition_UN "Underflow"
+#define DEC_Condition_ZE "No status"
+#define DEC_Condition_MU "Multiple status"
+#define DEC_Condition_Length 21        /* length of the longest string, */
+                                  /* including terminator */
+
+  /* Initialization descriptors, used by decContextDefault */
+#define DEC_INIT_BASE         0
+#define DEC_INIT_DECIMAL32   32
+#define DEC_INIT_DECIMAL64   64
+#define DEC_INIT_DECIMAL128 128
+
+  /* decContext routines */
+#ifdef IN_LIBGCC2
+#define decContextDefault __decContextDefault
+#define decContextSetStatus __decContextSetStatus
+#define decContextStatusToString __decContextStatusToString
+#define decContextSetStatusFromString __decContextSetStatusFromString
+#endif
+decContext *decContextDefault (decContext *, int32_t);
+decContext *decContextSetStatus (decContext *, uint32_t);
+const char *decContextStatusToString (decContext *);
+decContext *decContextSetStatusFromString (decContext *, char *);
+
+#endif
diff --git a/libdecnumber/decDPD.h b/libdecnumber/decDPD.h
new file mode 100644 (file)
index 0000000..827b675
--- /dev/null
@@ -0,0 +1,525 @@
+/* Binary Coded Decimal <--> Densely Packed Decimal lookup tables.
+   Copyright (C) 2005 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.
+
+   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, 59 Temple Place - Suite 330, Boston, MA
+   02111-1307, USA.  */
+
+/* ------------------------------------------------------------------------ */
+/* For details, see: http://www2.hursley.ibm.com/decimal/DPDecimal.html     */
+/*                                                                          */
+/* This include file defines conversion tables for DPD, as follows.         */
+/*                                                                          */
+/*   uint16_t BCD2DPD[2458];     // BCD -> DPD (0x999 => 2457)              */
+/*   uint16_t DPD2BCD[1024];     // DPD -> BCD (0x3FF => 0x999)             */
+/*   uint16_t BIN2DPD[1000];     // BIN -> DPD (999 => 2457)                */
+/*   uint16_t DPD2BIN[1024];     // DPD -> BIN (0x3FF => 999)               */
+/*                                                                          */
+/* In all cases the result (10 bits or 12 bits, or binary) is right-aligned */
+/* in the table entry.                                                      */
+/*                                                                          */
+/* To use a table, its name, prefixed with DEC_, must be defined with a     */
+/* value of 1 before this header file is included.  For example:            */
+/*    #define DEC_BCD2DPD 1                                                 */
+/* ------------------------------------------------------------------------ */
+
+#if DEC_BCD2DPD==1
+
+const uint16_t BCD2DPD[2458] = { 0, 1, 2, 3, 4, 5, 6, 7,
+  8, 9, 0, 0, 0, 0, 0, 0, 16, 17, 18, 19, 20,
+  21, 22, 23, 24, 25, 0, 0, 0, 0, 0, 0, 32, 33,
+  34, 35, 36, 37, 38, 39, 40, 41, 0, 0, 0, 0, 0,
+  0, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 0, 0,
+  0, 0, 0, 0, 64, 65, 66, 67, 68, 69, 70, 71, 72,
+  73, 0, 0, 0, 0, 0, 0, 80, 81, 82, 83, 84, 85,
+  86, 87, 88, 89, 0, 0, 0, 0, 0, 0, 96, 97, 98,
+  99, 100, 101, 102, 103, 104, 105, 0, 0, 0, 0, 0, 0,
+  112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 0, 0, 0,
+  0, 0, 0, 10, 11, 42, 43, 74, 75, 106, 107, 78, 79,
+  0, 0, 0, 0, 0, 0, 26, 27, 58, 59, 90, 91, 122,
+  123, 94, 95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 10, 11, 42, 43, 74,
+  75, 106, 107, 78, 79, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 0, 0,
+  0, 0, 0, 0, 144, 145, 146, 147, 148, 149, 150, 151, 152,
+  153, 0, 0, 0, 0, 0, 0, 160, 161, 162, 163, 164, 165,
+  166, 167, 168, 169, 0, 0, 0, 0, 0, 0, 176, 177, 178,
+  179, 180, 181, 182, 183, 184, 185, 0, 0, 0, 0, 0, 0,
+  192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 0, 0, 0,
+  0, 0, 0, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217,
+  0, 0, 0, 0, 0, 0, 224, 225, 226, 227, 228, 229, 230,
+  231, 232, 233, 0, 0, 0, 0, 0, 0, 240, 241, 242, 243,
+  244, 245, 246, 247, 248, 249, 0, 0, 0, 0, 0, 0, 138,
+  139, 170, 171, 202, 203, 234, 235, 206, 207, 0, 0, 0, 0,
+  0, 0, 154, 155, 186, 187, 218, 219, 250, 251, 222, 223, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 138, 139, 170, 171, 202, 203, 234, 235, 206,
+  207, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 256, 257, 258,
+  259, 260, 261, 262, 263, 264, 265, 0, 0, 0, 0, 0, 0,
+  272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 0, 0, 0,
+  0, 0, 0, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297,
+  0, 0, 0, 0, 0, 0, 304, 305, 306, 307, 308, 309, 310,
+  311, 312, 313, 0, 0, 0, 0, 0, 0, 320, 321, 322, 323,
+  324, 325, 326, 327, 328, 329, 0, 0, 0, 0, 0, 0, 336,
+  337, 338, 339, 340, 341, 342, 343, 344, 345, 0, 0, 0, 0,
+  0, 0, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 0,
+  0, 0, 0, 0, 0, 368, 369, 370, 371, 372, 373, 374, 375,
+  376, 377, 0, 0, 0, 0, 0, 0, 266, 267, 298, 299, 330,
+  331, 362, 363, 334, 335, 0, 0, 0, 0, 0, 0, 282, 283,
+  314, 315, 346, 347, 378, 379, 350, 351, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  266, 267, 298, 299, 330, 331, 362, 363, 334, 335, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 384, 385, 386, 387, 388, 389, 390,
+  391, 392, 393, 0, 0, 0, 0, 0, 0, 400, 401, 402, 403,
+  404, 405, 406, 407, 408, 409, 0, 0, 0, 0, 0, 0, 416,
+  417, 418, 419, 420, 421, 422, 423, 424, 425, 0, 0, 0, 0,
+  0, 0, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 0,
+  0, 0, 0, 0, 0, 448, 449, 450, 451, 452, 453, 454, 455,
+  456, 457, 0, 0, 0, 0, 0, 0, 464, 465, 466, 467, 468,
+  469, 470, 471, 472, 473, 0, 0, 0, 0, 0, 0, 480, 481,
+  482, 483, 484, 485, 486, 487, 488, 489, 0, 0, 0, 0, 0,
+  0, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 0, 0,
+  0, 0, 0, 0, 394, 395, 426, 427, 458, 459, 490, 491, 462,
+  463, 0, 0, 0, 0, 0, 0, 410, 411, 442, 443, 474, 475,
+  506, 507, 478, 479, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 394, 395, 426, 427,
+  458, 459, 490, 491, 462, 463, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 0,
+  0, 0, 0, 0, 0, 528, 529, 530, 531, 532, 533, 534, 535,
+  536, 537, 0, 0, 0, 0, 0, 0, 544, 545, 546, 547, 548,
+  549, 550, 551, 552, 553, 0, 0, 0, 0, 0, 0, 560, 561,
+  562, 563, 564, 565, 566, 567, 568, 569, 0, 0, 0, 0, 0,
+  0, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 0, 0,
+  0, 0, 0, 0, 592, 593, 594, 595, 596, 597, 598, 599, 600,
+  601, 0, 0, 0, 0, 0, 0, 608, 609, 610, 611, 612, 613,
+  614, 615, 616, 617, 0, 0, 0, 0, 0, 0, 624, 625, 626,
+  627, 628, 629, 630, 631, 632, 633, 0, 0, 0, 0, 0, 0,
+  522, 523, 554, 555, 586, 587, 618, 619, 590, 591, 0, 0, 0,
+  0, 0, 0, 538, 539, 570, 571, 602, 603, 634, 635, 606, 607,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 522, 523, 554, 555, 586, 587, 618, 619,
+  590, 591, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 640, 641,
+  642, 643, 644, 645, 646, 647, 648, 649, 0, 0, 0, 0, 0,
+  0, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 0, 0,
+  0, 0, 0, 0, 672, 673, 674, 675, 676, 677, 678, 679, 680,
+  681, 0, 0, 0, 0, 0, 0, 688, 689, 690, 691, 692, 693,
+  694, 695, 696, 697, 0, 0, 0, 0, 0, 0, 704, 705, 706,
+  707, 708, 709, 710, 711, 712, 713, 0, 0, 0, 0, 0, 0,
+  720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 0, 0, 0,
+  0, 0, 0, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745,
+  0, 0, 0, 0, 0, 0, 752, 753, 754, 755, 756, 757, 758,
+  759, 760, 761, 0, 0, 0, 0, 0, 0, 650, 651, 682, 683,
+  714, 715, 746, 747, 718, 719, 0, 0, 0, 0, 0, 0, 666,
+  667, 698, 699, 730, 731, 762, 763, 734, 735, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 650, 651, 682, 683, 714, 715, 746, 747, 718, 719, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 768, 769, 770, 771, 772, 773,
+  774, 775, 776, 777, 0, 0, 0, 0, 0, 0, 784, 785, 786,
+  787, 788, 789, 790, 791, 792, 793, 0, 0, 0, 0, 0, 0,
+  800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 0, 0, 0,
+  0, 0, 0, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825,
+  0, 0, 0, 0, 0, 0, 832, 833, 834, 835, 836, 837, 838,
+  839, 840, 841, 0, 0, 0, 0, 0, 0, 848, 849, 850, 851,
+  852, 853, 854, 855, 856, 857, 0, 0, 0, 0, 0, 0, 864,
+  865, 866, 867, 868, 869, 870, 871, 872, 873, 0, 0, 0, 0,
+  0, 0, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 0,
+  0, 0, 0, 0, 0, 778, 779, 810, 811, 842, 843, 874, 875,
+  846, 847, 0, 0, 0, 0, 0, 0, 794, 795, 826, 827, 858,
+  859, 890, 891, 862, 863, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 778, 779, 810,
+  811, 842, 843, 874, 875, 846, 847, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905,
+  0, 0, 0, 0, 0, 0, 912, 913, 914, 915, 916, 917, 918,
+  919, 920, 921, 0, 0, 0, 0, 0, 0, 928, 929, 930, 931,
+  932, 933, 934, 935, 936, 937, 0, 0, 0, 0, 0, 0, 944,
+  945, 946, 947, 948, 949, 950, 951, 952, 953, 0, 0, 0, 0,
+  0, 0, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 0,
+  0, 0, 0, 0, 0, 976, 977, 978, 979, 980, 981, 982, 983,
+  984, 985, 0, 0, 0, 0, 0, 0, 992, 993, 994, 995, 996,
+  997, 998, 999, 1000, 1001, 0, 0, 0, 0, 0, 0, 1008, 1009,
+  1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 0, 0, 0, 0, 0,
+  0, 906, 907, 938, 939, 970, 971, 1002, 1003, 974, 975, 0, 0,
+  0, 0, 0, 0, 922, 923, 954, 955, 986, 987, 1018, 1019, 990,
+  991, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 906, 907, 938, 939, 970, 971, 1002,
+  1003, 974, 975, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12,
+  13, 268, 269, 524, 525, 780, 781, 46, 47, 0, 0, 0, 0,
+  0, 0, 28, 29, 284, 285, 540, 541, 796, 797, 62, 63, 0,
+  0, 0, 0, 0, 0, 44, 45, 300, 301, 556, 557, 812, 813,
+  302, 303, 0, 0, 0, 0, 0, 0, 60, 61, 316, 317, 572,
+  573, 828, 829, 318, 319, 0, 0, 0, 0, 0, 0, 76, 77,
+  332, 333, 588, 589, 844, 845, 558, 559, 0, 0, 0, 0, 0,
+  0, 92, 93, 348, 349, 604, 605, 860, 861, 574, 575, 0, 0,
+  0, 0, 0, 0, 108, 109, 364, 365, 620, 621, 876, 877, 814,
+  815, 0, 0, 0, 0, 0, 0, 124, 125, 380, 381, 636, 637,
+  892, 893, 830, 831, 0, 0, 0, 0, 0, 0, 14, 15, 270,
+  271, 526, 527, 782, 783, 110, 111, 0, 0, 0, 0, 0, 0,
+  30, 31, 286, 287, 542, 543, 798, 799, 126, 127, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 14, 15, 270, 271, 526, 527, 782, 783, 110, 111, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 140, 141, 396, 397, 652,
+  653, 908, 909, 174, 175, 0, 0, 0, 0, 0, 0, 156, 157,
+  412, 413, 668, 669, 924, 925, 190, 191, 0, 0, 0, 0, 0,
+  0, 172, 173, 428, 429, 684, 685, 940, 941, 430, 431, 0, 0,
+  0, 0, 0, 0, 188, 189, 444, 445, 700, 701, 956, 957, 446,
+  447, 0, 0, 0, 0, 0, 0, 204, 205, 460, 461, 716, 717,
+  972, 973, 686, 687, 0, 0, 0, 0, 0, 0, 220, 221, 476,
+  477, 732, 733, 988, 989, 702, 703, 0, 0, 0, 0, 0, 0,
+  236, 237, 492, 493, 748, 749, 1004, 1005, 942, 943, 0, 0, 0,
+  0, 0, 0, 252, 253, 508, 509, 764, 765, 1020, 1021, 958, 959,
+  0, 0, 0, 0, 0, 0, 142, 143, 398, 399, 654, 655, 910,
+  911, 238, 239, 0, 0, 0, 0, 0, 0, 158, 159, 414, 415,
+  670, 671, 926, 927, 254, 255
+};
+#endif
+
+#if DEC_DPD2BCD==1
+
+const uint16_t DPD2BCD[1024] = { 0, 1, 2, 3, 4, 5, 6, 7,
+  8, 9, 128, 129, 2048, 2049, 2176, 2177, 16, 17, 18, 19, 20,
+  21, 22, 23, 24, 25, 144, 145, 2064, 2065, 2192, 2193, 32, 33,
+  34, 35, 36, 37, 38, 39, 40, 41, 130, 131, 2080, 2081, 2056,
+  2057, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 146, 147,
+  2096, 2097, 2072, 2073, 64, 65, 66, 67, 68, 69, 70, 71, 72,
+  73, 132, 133, 2112, 2113, 136, 137, 80, 81, 82, 83, 84, 85,
+  86, 87, 88, 89, 148, 149, 2128, 2129, 152, 153, 96, 97, 98,
+  99, 100, 101, 102, 103, 104, 105, 134, 135, 2144, 2145, 2184, 2185,
+  112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 150, 151, 2160,
+  2161, 2200, 2201, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265,
+  384, 385, 2304, 2305, 2432, 2433, 272, 273, 274, 275, 276, 277, 278,
+  279, 280, 281, 400, 401, 2320, 2321, 2448, 2449, 288, 289, 290, 291,
+  292, 293, 294, 295, 296, 297, 386, 387, 2336, 2337, 2312, 2313, 304,
+  305, 306, 307, 308, 309, 310, 311, 312, 313, 402, 403, 2352, 2353,
+  2328, 2329, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 388,
+  389, 2368, 2369, 392, 393, 336, 337, 338, 339, 340, 341, 342, 343,
+  344, 345, 404, 405, 2384, 2385, 408, 409, 352, 353, 354, 355, 356,
+  357, 358, 359, 360, 361, 390, 391, 2400, 2401, 2440, 2441, 368, 369,
+  370, 371, 372, 373, 374, 375, 376, 377, 406, 407, 2416, 2417, 2456,
+  2457, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 640, 641,
+  2050, 2051, 2178, 2179, 528, 529, 530, 531, 532, 533, 534, 535, 536,
+  537, 656, 657, 2066, 2067, 2194, 2195, 544, 545, 546, 547, 548, 549,
+  550, 551, 552, 553, 642, 643, 2082, 2083, 2088, 2089, 560, 561, 562,
+  563, 564, 565, 566, 567, 568, 569, 658, 659, 2098, 2099, 2104, 2105,
+  576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 644, 645, 2114,
+  2115, 648, 649, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601,
+  660, 661, 2130, 2131, 664, 665, 608, 609, 610, 611, 612, 613, 614,
+  615, 616, 617, 646, 647, 2146, 2147, 2184, 2185, 624, 625, 626, 627,
+  628, 629, 630, 631, 632, 633, 662, 663, 2162, 2163, 2200, 2201, 768,
+  769, 770, 771, 772, 773, 774, 775, 776, 777, 896, 897, 2306, 2307,
+  2434, 2435, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 912,
+  913, 2322, 2323, 2450, 2451, 800, 801, 802, 803, 804, 805, 806, 807,
+  808, 809, 898, 899, 2338, 2339, 2344, 2345, 816, 817, 818, 819, 820,
+  821, 822, 823, 824, 825, 914, 915, 2354, 2355, 2360, 2361, 832, 833,
+  834, 835, 836, 837, 838, 839, 840, 841, 900, 901, 2370, 2371, 904,
+  905, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 916, 917,
+  2386, 2387, 920, 921, 864, 865, 866, 867, 868, 869, 870, 871, 872,
+  873, 902, 903, 2402, 2403, 2440, 2441, 880, 881, 882, 883, 884, 885,
+  886, 887, 888, 889, 918, 919, 2418, 2419, 2456, 2457, 1024, 1025, 1026,
+  1027, 1028, 1029, 1030, 1031, 1032, 1033, 1152, 1153, 2052, 2053, 2180,
+    2181,
+  1040, 1041, 1042, 1043, 1044, 1045, 1046, 1047, 1048, 1049, 1168, 1169,
+    2068,
+  2069, 2196, 2197, 1056, 1057, 1058, 1059, 1060, 1061, 1062, 1063, 1064,
+    1065,
+  1154, 1155, 2084, 2085, 2120, 2121, 1072, 1073, 1074, 1075, 1076, 1077,
+    1078,
+  1079, 1080, 1081, 1170, 1171, 2100, 2101, 2136, 2137, 1088, 1089, 1090,
+    1091,
+  1092, 1093, 1094, 1095, 1096, 1097, 1156, 1157, 2116, 2117, 1160, 1161,
+    1104,
+  1105, 1106, 1107, 1108, 1109, 1110, 1111, 1112, 1113, 1172, 1173, 2132,
+    2133,
+  1176, 1177, 1120, 1121, 1122, 1123, 1124, 1125, 1126, 1127, 1128, 1129,
+    1158,
+  1159, 2148, 2149, 2184, 2185, 1136, 1137, 1138, 1139, 1140, 1141, 1142,
+    1143,
+  1144, 1145, 1174, 1175, 2164, 2165, 2200, 2201, 1280, 1281, 1282, 1283,
+    1284,
+  1285, 1286, 1287, 1288, 1289, 1408, 1409, 2308, 2309, 2436, 2437, 1296,
+    1297,
+  1298, 1299, 1300, 1301, 1302, 1303, 1304, 1305, 1424, 1425, 2324, 2325,
+    2452,
+  2453, 1312, 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320, 1321, 1410,
+    1411,
+  2340, 2341, 2376, 2377, 1328, 1329, 1330, 1331, 1332, 1333, 1334, 1335,
+    1336,
+  1337, 1426, 1427, 2356, 2357, 2392, 2393, 1344, 1345, 1346, 1347, 1348,
+    1349,
+  1350, 1351, 1352, 1353, 1412, 1413, 2372, 2373, 1416, 1417, 1360, 1361,
+    1362,
+  1363, 1364, 1365, 1366, 1367, 1368, 1369, 1428, 1429, 2388, 2389, 1432,
+    1433,
+  1376, 1377, 1378, 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1414, 1415,
+    2404,
+  2405, 2440, 2441, 1392, 1393, 1394, 1395, 1396, 1397, 1398, 1399, 1400,
+    1401,
+  1430, 1431, 2420, 2421, 2456, 2457, 1536, 1537, 1538, 1539, 1540, 1541,
+    1542,
+  1543, 1544, 1545, 1664, 1665, 2054, 2055, 2182, 2183, 1552, 1553, 1554,
+    1555,
+  1556, 1557, 1558, 1559, 1560, 1561, 1680, 1681, 2070, 2071, 2198, 2199,
+    1568,
+  1569, 1570, 1571, 1572, 1573, 1574, 1575, 1576, 1577, 1666, 1667, 2086,
+    2087,
+  2152, 2153, 1584, 1585, 1586, 1587, 1588, 1589, 1590, 1591, 1592, 1593,
+    1682,
+  1683, 2102, 2103, 2168, 2169, 1600, 1601, 1602, 1603, 1604, 1605, 1606,
+    1607,
+  1608, 1609, 1668, 1669, 2118, 2119, 1672, 1673, 1616, 1617, 1618, 1619,
+    1620,
+  1621, 1622, 1623, 1624, 1625, 1684, 1685, 2134, 2135, 1688, 1689, 1632,
+    1633,
+  1634, 1635, 1636, 1637, 1638, 1639, 1640, 1641, 1670, 1671, 2150, 2151,
+    2184,
+  2185, 1648, 1649, 1650, 1651, 1652, 1653, 1654, 1655, 1656, 1657, 1686,
+    1687,
+  2166, 2167, 2200, 2201, 1792, 1793, 1794, 1795, 1796, 1797, 1798, 1799,
+    1800,
+  1801, 1920, 1921, 2310, 2311, 2438, 2439, 1808, 1809, 1810, 1811, 1812,
+    1813,
+  1814, 1815, 1816, 1817, 1936, 1937, 2326, 2327, 2454, 2455, 1824, 1825,
+    1826,
+  1827, 1828, 1829, 1830, 1831, 1832, 1833, 1922, 1923, 2342, 2343, 2408,
+    2409,
+  1840, 1841, 1842, 1843, 1844, 1845, 1846, 1847, 1848, 1849, 1938, 1939,
+    2358,
+  2359, 2424, 2425, 1856, 1857, 1858, 1859, 1860, 1861, 1862, 1863, 1864,
+    1865,
+  1924, 1925, 2374, 2375, 1928, 1929, 1872, 1873, 1874, 1875, 1876, 1877,
+    1878,
+  1879, 1880, 1881, 1940, 1941, 2390, 2391, 1944, 1945, 1888, 1889, 1890,
+    1891,
+  1892, 1893, 1894, 1895, 1896, 1897, 1926, 1927, 2406, 2407, 2440, 2441,
+    1904,
+  1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1942, 1943, 2422,
+    2423,
+  2456, 2457
+};
+#endif
+
+#if DEC_BIN2DPD==1
+
+const uint16_t BIN2DPD[1000] = { 0, 1, 2, 3, 4, 5, 6, 7,
+  8, 9, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 32,
+  33, 34, 35, 36, 37, 38, 39, 40, 41, 48, 49, 50, 51,
+  52, 53, 54, 55, 56, 57, 64, 65, 66, 67, 68, 69, 70,
+  71, 72, 73, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
+  96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 112, 113, 114,
+  115, 116, 117, 118, 119, 120, 121, 10, 11, 42, 43, 74, 75,
+  106, 107, 78, 79, 26, 27, 58, 59, 90, 91, 122, 123, 94,
+  95, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 144, 145,
+  146, 147, 148, 149, 150, 151, 152, 153, 160, 161, 162, 163, 164,
+  165, 166, 167, 168, 169, 176, 177, 178, 179, 180, 181, 182, 183,
+  184, 185, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 208,
+  209, 210, 211, 212, 213, 214, 215, 216, 217, 224, 225, 226, 227,
+  228, 229, 230, 231, 232, 233, 240, 241, 242, 243, 244, 245, 246,
+  247, 248, 249, 138, 139, 170, 171, 202, 203, 234, 235, 206, 207,
+  154, 155, 186, 187, 218, 219, 250, 251, 222, 223, 256, 257, 258,
+  259, 260, 261, 262, 263, 264, 265, 272, 273, 274, 275, 276, 277,
+  278, 279, 280, 281, 288, 289, 290, 291, 292, 293, 294, 295, 296,
+  297, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 320, 321,
+  322, 323, 324, 325, 326, 327, 328, 329, 336, 337, 338, 339, 340,
+  341, 342, 343, 344, 345, 352, 353, 354, 355, 356, 357, 358, 359,
+  360, 361, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 266,
+  267, 298, 299, 330, 331, 362, 363, 334, 335, 282, 283, 314, 315,
+  346, 347, 378, 379, 350, 351, 384, 385, 386, 387, 388, 389, 390,
+  391, 392, 393, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409,
+  416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 432, 433, 434,
+  435, 436, 437, 438, 439, 440, 441, 448, 449, 450, 451, 452, 453,
+  454, 455, 456, 457, 464, 465, 466, 467, 468, 469, 470, 471, 472,
+  473, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 496, 497,
+  498, 499, 500, 501, 502, 503, 504, 505, 394, 395, 426, 427, 458,
+  459, 490, 491, 462, 463, 410, 411, 442, 443, 474, 475, 506, 507,
+  478, 479, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 528,
+  529, 530, 531, 532, 533, 534, 535, 536, 537, 544, 545, 546, 547,
+  548, 549, 550, 551, 552, 553, 560, 561, 562, 563, 564, 565, 566,
+  567, 568, 569, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585,
+  592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 608, 609, 610,
+  611, 612, 613, 614, 615, 616, 617, 624, 625, 626, 627, 628, 629,
+  630, 631, 632, 633, 522, 523, 554, 555, 586, 587, 618, 619, 590,
+  591, 538, 539, 570, 571, 602, 603, 634, 635, 606, 607, 640, 641,
+  642, 643, 644, 645, 646, 647, 648, 649, 656, 657, 658, 659, 660,
+  661, 662, 663, 664, 665, 672, 673, 674, 675, 676, 677, 678, 679,
+  680, 681, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 704,
+  705, 706, 707, 708, 709, 710, 711, 712, 713, 720, 721, 722, 723,
+  724, 725, 726, 727, 728, 729, 736, 737, 738, 739, 740, 741, 742,
+  743, 744, 745, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761,
+  650, 651, 682, 683, 714, 715, 746, 747, 718, 719, 666, 667, 698,
+  699, 730, 731, 762, 763, 734, 735, 768, 769, 770, 771, 772, 773,
+  774, 775, 776, 777, 784, 785, 786, 787, 788, 789, 790, 791, 792,
+  793, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 816, 817,
+  818, 819, 820, 821, 822, 823, 824, 825, 832, 833, 834, 835, 836,
+  837, 838, 839, 840, 841, 848, 849, 850, 851, 852, 853, 854, 855,
+  856, 857, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 880,
+  881, 882, 883, 884, 885, 886, 887, 888, 889, 778, 779, 810, 811,
+  842, 843, 874, 875, 846, 847, 794, 795, 826, 827, 858, 859, 890,
+  891, 862, 863, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905,
+  912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 928, 929, 930,
+  931, 932, 933, 934, 935, 936, 937, 944, 945, 946, 947, 948, 949,
+  950, 951, 952, 953, 960, 961, 962, 963, 964, 965, 966, 967, 968,
+  969, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 992, 993,
+  994, 995, 996, 997, 998, 999, 1000, 1001, 1008, 1009, 1010, 1011, 1012,
+  1013, 1014, 1015, 1016, 1017, 906, 907, 938, 939, 970, 971, 1002, 1003,
+  974, 975, 922, 923, 954, 955, 986, 987, 1018, 1019, 990, 991, 12,
+  13, 268, 269, 524, 525, 780, 781, 46, 47, 28, 29, 284, 285,
+  540, 541, 796, 797, 62, 63, 44, 45, 300, 301, 556, 557, 812,
+  813, 302, 303, 60, 61, 316, 317, 572, 573, 828, 829, 318, 319,
+  76, 77, 332, 333, 588, 589, 844, 845, 558, 559, 92, 93, 348,
+  349, 604, 605, 860, 861, 574, 575, 108, 109, 364, 365, 620, 621,
+  876, 877, 814, 815, 124, 125, 380, 381, 636, 637, 892, 893, 830,
+  831, 14, 15, 270, 271, 526, 527, 782, 783, 110, 111, 30, 31,
+  286, 287, 542, 543, 798, 799, 126, 127, 140, 141, 396, 397, 652,
+  653, 908, 909, 174, 175, 156, 157, 412, 413, 668, 669, 924, 925,
+  190, 191, 172, 173, 428, 429, 684, 685, 940, 941, 430, 431, 188,
+  189, 444, 445, 700, 701, 956, 957, 446, 447, 204, 205, 460, 461,
+  716, 717, 972, 973, 686, 687, 220, 221, 476, 477, 732, 733, 988,
+  989, 702, 703, 236, 237, 492, 493, 748, 749, 1004, 1005, 942, 943,
+  252, 253, 508, 509, 764, 765, 1020, 1021, 958, 959, 142, 143, 398,
+  399, 654, 655, 910, 911, 238, 239, 158, 159, 414, 415, 670, 671,
+  926, 927, 254, 255
+};
+#endif
+
+#if DEC_DPD2BIN==1
+
+const uint16_t DPD2BIN[1024] = { 0, 1, 2, 3, 4, 5, 6, 7,
+  8, 9, 80, 81, 800, 801, 880, 881, 10, 11, 12, 13, 14,
+  15, 16, 17, 18, 19, 90, 91, 810, 811, 890, 891, 20, 21,
+  22, 23, 24, 25, 26, 27, 28, 29, 82, 83, 820, 821, 808,
+  809, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 92, 93,
+  830, 831, 818, 819, 40, 41, 42, 43, 44, 45, 46, 47, 48,
+  49, 84, 85, 840, 841, 88, 89, 50, 51, 52, 53, 54, 55,
+  56, 57, 58, 59, 94, 95, 850, 851, 98, 99, 60, 61, 62,
+  63, 64, 65, 66, 67, 68, 69, 86, 87, 860, 861, 888, 889,
+  70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 96, 97, 870,
+  871, 898, 899, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109,
+  180, 181, 900, 901, 980, 981, 110, 111, 112, 113, 114, 115, 116,
+  117, 118, 119, 190, 191, 910, 911, 990, 991, 120, 121, 122, 123,
+  124, 125, 126, 127, 128, 129, 182, 183, 920, 921, 908, 909, 130,
+  131, 132, 133, 134, 135, 136, 137, 138, 139, 192, 193, 930, 931,
+  918, 919, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 184,
+  185, 940, 941, 188, 189, 150, 151, 152, 153, 154, 155, 156, 157,
+  158, 159, 194, 195, 950, 951, 198, 199, 160, 161, 162, 163, 164,
+  165, 166, 167, 168, 169, 186, 187, 960, 961, 988, 989, 170, 171,
+  172, 173, 174, 175, 176, 177, 178, 179, 196, 197, 970, 971, 998,
+  999, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 280, 281,
+  802, 803, 882, 883, 210, 211, 212, 213, 214, 215, 216, 217, 218,
+  219, 290, 291, 812, 813, 892, 893, 220, 221, 222, 223, 224, 225,
+  226, 227, 228, 229, 282, 283, 822, 823, 828, 829, 230, 231, 232,
+  233, 234, 235, 236, 237, 238, 239, 292, 293, 832, 833, 838, 839,
+  240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 284, 285, 842,
+  843, 288, 289, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259,
+  294, 295, 852, 853, 298, 299, 260, 261, 262, 263, 264, 265, 266,
+  267, 268, 269, 286, 287, 862, 863, 888, 889, 270, 271, 272, 273,
+  274, 275, 276, 277, 278, 279, 296, 297, 872, 873, 898, 899, 300,
+  301, 302, 303, 304, 305, 306, 307, 308, 309, 380, 381, 902, 903,
+  982, 983, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 390,
+  391, 912, 913, 992, 993, 320, 321, 322, 323, 324, 325, 326, 327,
+  328, 329, 382, 383, 922, 923, 928, 929, 330, 331, 332, 333, 334,
+  335, 336, 337, 338, 339, 392, 393, 932, 933, 938, 939, 340, 341,
+  342, 343, 344, 345, 346, 347, 348, 349, 384, 385, 942, 943, 388,
+  389, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 394, 395,
+  952, 953, 398, 399, 360, 361, 362, 363, 364, 365, 366, 367, 368,
+  369, 386, 387, 962, 963, 988, 989, 370, 371, 372, 373, 374, 375,
+  376, 377, 378, 379, 396, 397, 972, 973, 998, 999, 400, 401, 402,
+  403, 404, 405, 406, 407, 408, 409, 480, 481, 804, 805, 884, 885,
+  410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 490, 491, 814,
+  815, 894, 895, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429,
+  482, 483, 824, 825, 848, 849, 430, 431, 432, 433, 434, 435, 436,
+  437, 438, 439, 492, 493, 834, 835, 858, 859, 440, 441, 442, 443,
+  444, 445, 446, 447, 448, 449, 484, 485, 844, 845, 488, 489, 450,
+  451, 452, 453, 454, 455, 456, 457, 458, 459, 494, 495, 854, 855,
+  498, 499, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 486,
+  487, 864, 865, 888, 889, 470, 471, 472, 473, 474, 475, 476, 477,
+  478, 479, 496, 497, 874, 875, 898, 899, 500, 501, 502, 503, 504,
+  505, 506, 507, 508, 509, 580, 581, 904, 905, 984, 985, 510, 511,
+  512, 513, 514, 515, 516, 517, 518, 519, 590, 591, 914, 915, 994,
+  995, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 582, 583,
+  924, 925, 948, 949, 530, 531, 532, 533, 534, 535, 536, 537, 538,
+  539, 592, 593, 934, 935, 958, 959, 540, 541, 542, 543, 544, 545,
+  546, 547, 548, 549, 584, 585, 944, 945, 588, 589, 550, 551, 552,
+  553, 554, 555, 556, 557, 558, 559, 594, 595, 954, 955, 598, 599,
+  560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 586, 587, 964,
+  965, 988, 989, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579,
+  596, 597, 974, 975, 998, 999, 600, 601, 602, 603, 604, 605, 606,
+  607, 608, 609, 680, 681, 806, 807, 886, 887, 610, 611, 612, 613,
+  614, 615, 616, 617, 618, 619, 690, 691, 816, 817, 896, 897, 620,
+  621, 622, 623, 624, 625, 626, 627, 628, 629, 682, 683, 826, 827,
+  868, 869, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 692,
+  693, 836, 837, 878, 879, 640, 641, 642, 643, 644, 645, 646, 647,
+  648, 649, 684, 685, 846, 847, 688, 689, 650, 651, 652, 653, 654,
+  655, 656, 657, 658, 659, 694, 695, 856, 857, 698, 699, 660, 661,
+  662, 663, 664, 665, 666, 667, 668, 669, 686, 687, 866, 867, 888,
+  889, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 696, 697,
+  876, 877, 898, 899, 700, 701, 702, 703, 704, 705, 706, 707, 708,
+  709, 780, 781, 906, 907, 986, 987, 710, 711, 712, 713, 714, 715,
+  716, 717, 718, 719, 790, 791, 916, 917, 996, 997, 720, 721, 722,
+  723, 724, 725, 726, 727, 728, 729, 782, 783, 926, 927, 968, 969,
+  730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 792, 793, 936,
+  937, 978, 979, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749,
+  784, 785, 946, 947, 788, 789, 750, 751, 752, 753, 754, 755, 756,
+  757, 758, 759, 794, 795, 956, 957, 798, 799, 760, 761, 762, 763,
+  764, 765, 766, 767, 768, 769, 786, 787, 966, 967, 988, 989, 770,
+  771, 772, 773, 774, 775, 776, 777, 778, 779, 796, 797, 976, 977,
+  998, 999
+};
+#endif
diff --git a/libdecnumber/decLibrary.c b/libdecnumber/decLibrary.c
new file mode 100644 (file)
index 0000000..b7cbd10
--- /dev/null
@@ -0,0 +1,81 @@
+/* Temporary library support for decimal floating point.
+   Copyright (C) 2005 Free Software Foundation, Inc.
+
+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.
+
+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.  */
+
+#include "decContext.h"
+#include "decimal128.h"
+#include "decimal64.h"
+#include "decimal32.h"
+
+void __host_to_ieee_32 (_Decimal32, decimal32 *);
+void __host_to_ieee_64 (_Decimal64, decimal64 *);
+void __host_to_ieee_128 (_Decimal128, decimal128 *);
+
+extern int isinfd32 (_Decimal32);
+extern int isinfd64 (_Decimal64);
+extern int isinfd128 (_Decimal128);
+extern void __dfp_enable_traps (void);
+extern void __dfp_raise (int exception __attribute__ ((unused)));
+
+int
+isinfd32 (_Decimal32 arg)
+{
+  decNumber dn;
+  decimal32 d32;
+
+  __host_to_ieee_32 (arg, &d32);
+  decimal32ToNumber (&d32, &dn);
+  return (decNumberIsInfinite (&dn));
+}
+
+int
+isinfd64 (_Decimal64 arg)
+{
+  decNumber dn;
+  decimal64 d64;
+
+  __host_to_ieee_64 (arg, &d64);
+  decimal64ToNumber (&d64, &dn);
+  return (decNumberIsInfinite (&dn));
+}
+
+int
+isinfd128 (_Decimal128 arg)
+{
+  decNumber dn;
+  decimal128 d128;
+
+  __host_to_ieee_128 (arg, &d128);
+  decimal128ToNumber (&d128, &dn);
+  return (decNumberIsInfinite (&dn));
+}
+
+int __dfp_traps;
+
+void
+__dfp_enable_traps (void)
+{
+  __dfp_traps = 1;
+}
+
+void
+__dfp_raise (int exception __attribute__ ((unused)))
+{
+  raise (SIGFPE);
+}
diff --git a/libdecnumber/decNumber.c b/libdecnumber/decNumber.c
new file mode 100644 (file)
index 0000000..c2a8980
--- /dev/null
@@ -0,0 +1,5939 @@
+/* Decimal Number module for the decNumber C Library
+   Copyright (C) 2005 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.
+
+   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, 59 Temple Place - Suite 330, Boston, MA
+   02111-1307, USA.  */
+
+/* ------------------------------------------------------------------ */
+/* This module comprises the routines for Standard 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:                                   */
+/*                                                                    */
+/*    a) Line comments (double forward slash) are used.  (Most C      */
+/*       compilers accept these.  If yours does not, a simple script  */
+/*       can be used to convert them to ANSI C comments.)             */
+/*                                                                    */
+/*    b) Types from C99 stdint.h are used.  If you do not have this   */
+/*       header file, see the User's Guide section of the decNumber   */
+/*       documentation; this lists the necessary definitions.         */
+/*                                                                    */
+/*    c) If DECDPUN>4, non-ANSI 64-bit 'long long' types are used.    */
+/*       To avoid these, set 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).                                  */
+/*                                                                    */
+/* 3. 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 may not overlap.    */
+/*                                                                    */
+/* 4. 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).                                               */
+/*                                                                    */
+/* 5. The decNumber format is not an exchangeable concrete            */
+/*    representation as it comprises fields which may be machine-     */
+/*    dependent (big-endian or little-endian, for example).           */
+/*    Canonical conversions to and from strings are provided; other   */
+/*    conversions are available in separate modules.                  */
+/*                                                                    */
+/* 6. 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    */
+/*    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.     */
+/*                                                                    */
+/* 7. 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 from this.         */
+/*    Other routines are allowed to use the return statement inline.  */
+/*                                                                    */
+/*    Storage leak accounting can be enabled using DECALLOC.          */
+/*                                                                    */
+/* 2. All loops use the for(;;) construct.  Any do construct is for   */
+/*    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 we accumulate status and only finally apply it  */
+/*    by calling decContextSetStatus (via decStatus).                 */
+/*                                                                    */
+/*    Routines which allocate storage cannot, therefore, 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. We use the multiply-by-reciprocal 'trick' 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 we assumed floating-point multiply available). */
+/*                                                                    */
+/* 8. Unusual abbreviations possibly 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)                */
+/*      msu -- most significant Unit (of coefficient)                 */
+/*      rhs -- right hand side (operand, of an operation)             */
+/*      +ve -- positive                                               */
+/*      -ve -- negative                                               */
+/* ------------------------------------------------------------------ */
+
+#include <stdlib.h>            /* for malloc, free, etc. */
+#include <stdio.h>             /* for printf [if needed] */
+#include <string.h>            /* for strcpy */
+#include <ctype.h>             /* for lower */
+#include "decNumber.h"         /* base number library */
+#include "decNumberLocal.h"    /* decNumber local types, etc. */
+
+/* Constants */
+/* Public constant array: powers of ten (powers[n]==10**n) */
+const uInt powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
+  10000000, 100000000, 1000000000
+};
+
+/* 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 COMPNAN   0x04         /* .. [NaN processing] */
+
+#define DEC_sNaN 0x40000000    /* local status: sNaN signal */
+#define BADINT (Int)0x80000000 /* most-negative Int; error indicator */
+
+static Unit one[] = { 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 we currently use non-ANSI 64-bit types.  These could */
+  /* be replaced by subroutine calls later. */
+#ifdef long
+#undef long
+#endif
+typedef signed long long Long;
+typedef unsigned long long uLong;
+#define eInt  Long             /* extended integer */
+#define ueInt uLong            /* unsigned extended integer */
+#endif
+
+/* Local routines */
+static decNumber *decAddOp (decNumber *, decNumber *, decNumber *,
+                           decContext *, uByte, uInt *);
+static void decApplyRound (decNumber *, decContext *, Int, uInt *);
+static Int decCompare (decNumber * lhs, decNumber * rhs);
+static decNumber *decCompareOp (decNumber *, decNumber *, decNumber *,
+                               decContext *, Flag, uInt *);
+static void decCopyFit (decNumber *, decNumber *, decContext *,
+                       Int *, uInt *);
+static decNumber *decDivideOp (decNumber *, decNumber *, decNumber *,
+                              decContext *, Flag, uInt *);
+static void decFinalize (decNumber *, decContext *, Int *, uInt *);
+static Int decGetDigits (Unit *, Int);
+#if DECSUBSET
+static Int decGetInt (decNumber *, decContext *);
+#else
+static Int decGetInt (decNumber *);
+#endif
+static decNumber *decMultiplyOp (decNumber *, decNumber *, decNumber *,
+                                decContext *, uInt *);
+static decNumber *decNaNs (decNumber *, decNumber *, decNumber *, uInt *);
+static decNumber *decQuantizeOp (decNumber *, decNumber *, decNumber *,
+                                decContext *, Flag, uInt *);
+static void decSetCoeff (decNumber *, decContext *, Unit *,
+                        Int, Int *, uInt *);
+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 Flag decStrEq (const char *, const char *);
+static void decToString (decNumber *, char[], Flag);
+static decNumber *decTrim (decNumber *, Flag, Int *);
+static Int decUnitAddSub (Unit *, Int, Unit *, Int, Int, Unit *, Int);
+static Int decUnitCompare (Unit *, Int, 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 (decNumber *, decContext *, uInt *);
+#endif
+
+/* 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 should 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 operand 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'. */
+#define DECUNUSED (void *)(0xffffffff)
+static Flag decCheckOperands (decNumber *, decNumber *, decNumber *,
+                             decContext *);
+static Flag decCheckNumber (decNumber *, decContext *);
+#endif
+
+#if DECTRACE || DECCHECK
+/* Optional trace/debugging routines. */
+void decNumberShow (decNumber *);      /* displays the components of a number */
+static void decDumpAr (char, Unit *, Int);
+#endif
+
+/* ================================================================== */
+/* Conversions                                                        */
+/* ================================================================== */
+
+/* ------------------------------------------------------------------ */
+/* 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 (decNumber * dn, char *string)
+{
+  decToString (dn, string, 0);
+  return string;
+}
+
+char *
+decNumberToEngString (decNumber * dn, char *string)
+{
+  decToString (dn, string, 1);
+  return string;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, 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[D2U (DECBUFFER + 1)];   /* local buffer in case need temporary */
+  Unit *allocres = NULL;       /* -> allocated result, iff allocated */
+  Int need;                    /* units needed for result */
+  Int d = 0;                   /* count of digits found in decimal part */
+  char *dotchar = NULL;                /* where dot was found */
+  char *cfirst;                        /* -> first character of decimal part */
+  char *last = NULL;           /* -> last digit of decimal part */
+  char *firstexp;              /* -> first significant exponent digit */
+  char *c;                     /* work */
+  Unit *up;                    /* .. */
+#if DECDPUN>1
+  Int i;                       /* .. */
+#endif
+  Int residue = 0;             /* rounding residue */
+  uInt status = 0;             /* error code */
+
+#if DECCHECK
+  if (decCheckOperands (DECUNUSED, DECUNUSED, DECUNUSED, set))
+    return decNumberZero (dn);
+#endif
+
+  do
+    {                          /* status & malloc protection */
+      c = chars;               /* -> input character */
+      if (*c == '-')
+       {                       /* handle leading '-' */
+         bits = DECNEG;
+         c++;
+       }
+      else if (*c == '+')
+       c++;                    /* step over leading '+' */
+      /* We're at the start of the number [we think] */
+      cfirst = c;              /* save */
+      for (;; c++)
+       {
+         if (*c >= '0' && *c <= '9')
+           {                   /* test for Arabic digit */
+             last = c;
+             d++;              /* count of real digits */
+             continue;         /* still in decimal part */
+           }
+         if (*c != '.')
+           break;              /* done with decimal part */
+         /* dot: record, check, and ignore */
+         if (dotchar != NULL)
+           {                   /* two dots */
+             last = NULL;      /* indicate bad */
+             break;
+           }                   /* .. and go report */
+         dotchar = c;          /* offset into decimal part */
+       }                       /* c */
+
+      if (last == NULL)
+       {                       /* no decimal digits, or >1 . */
+#if DECSUBSET
+         /* If subset then infinities and NaNs are not allowed */
+         if (!set->extended)
+           {
+             status = DEC_Conversion_syntax;
+             break;            /* all done */
+           }
+         else
+           {
+#endif
+             /* Infinities and NaNs are possible, here */
+             decNumberZero (dn);       /* be optimistic */
+             if (decStrEq (c, "Infinity") || decStrEq (c, "Inf"))
+               {
+                 dn->bits = bits | DECINF;
+                 break;        /* all done */
+               }
+             else
+               {               /* a NaN expected */
+                 /* 2003.09.10 NaNs are now permitted to have a sign */
+                 status = DEC_Conversion_syntax;       /* assume the worst */
+                 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 nothing, or nnnn, 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)
+                   break;      /* too many digits */
+                 /* good; drop through and convert the integer */
+                 status = 0;
+                 bits = dn->bits;      /* for copy-back */
+               }               /* NaN expected */
+#if DECSUBSET
+           }
+#endif
+       }                       /* last==NULL */
+
+      if (*c != '\0')
+       {                       /* more there; exponent expected... */
+         Flag nege = 0;        /* 1=negative exponent */
+         if (*c != 'e' && *c != 'E')
+           {
+             status = DEC_Conversion_syntax;
+             break;
+           }
+
+         /* Found 'e' or 'E' -- now process explicit exponent */
+         /* 1998.07.11: sign no longer required */
+         c++;                  /* to (expected) sign */
+         if (*c == '-')
+           {
+             nege = 1;
+             c++;
+           }
+         else if (*c == '+')
+           c++;
+         if (*c == '\0')
+           {
+             status = DEC_Conversion_syntax;
+             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 we didn't end on '\0' must not be a digit */
+         if (*c != '\0')
+           {
+             status = DEC_Conversion_syntax;
+             break;
+           }
+
+         /* (this next test must be after the syntax check) */
+         /* 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 */
+       }                       /* had exponent */
+      /* Here when all inspected; syntax is good */
+
+      /* Handle decimal point... */
+      if (dotchar != NULL && dotchar < last)   /* embedded . found, so */
+       exponent = exponent - (last - dotchar); /* .. adjust exponent */
+      /* [we can now ignore the .] */
+
+      /* strip leading zeros/dot (leave final if all 0's) */
+      for (c = cfirst; c < last; c++)
+       {
+         if (*c == '0')
+           d--;                /* 0 stripped */
+         else if (*c != '.')
+           break;
+         cfirst++;             /* step past leader */
+       }                       /* c */
+
+#if DECSUBSET
+      /* We can now make a rapid exit for zeros if !extended */
+      if (*cfirst == '0' && !set->extended)
+       {
+         decNumberZero (dn);   /* clean result */
+         break;                /* [could be return] */
+       }
+#endif
+
+      /* OK, the digits string is good.  Copy to the decNumber, or to
+         a temporary decNumber if rounding is needed */
+      if (d <= set->digits)
+       res = dn->lsu;          /* fits into given decNumber */
+      else
+       {                       /* rounding needed */
+         need = D2U (d);       /* units needed */
+         res = resbuff;        /* assume use local buffer */
+         if (need * sizeof (Unit) > sizeof (resbuff))
+           {                   /* too big for local */
+             allocres = (Unit *) malloc (need * sizeof (Unit));
+             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 */
+#if DECDPUN>1
+      i = d % DECDPUN;         /* digits in top unit */
+      if (i == 0)
+       i = DECDPUN;
+      up = res + D2U (d) - 1;  /* -> msu */
+      *up = 0;
+      for (c = cfirst;; c++)
+       {                       /* along the digits */
+         if (*c == '.')
+           {                   /* ignore . [don't decrement i] */
+             if (c != last)
+               continue;
+             break;
+           }
+         *up = (Unit) (X10 (*up) + (Int) * c - (Int) '0');
+         i--;
+         if (i > 0)
+           continue;           /* more for this unit */
+         if (up == res)
+           break;              /* just filled the last unit */
+         i = DECDPUN;
+         up--;
+         *up = 0;
+       }                       /* c */
+#else
+      /* DECDPUN==1 */
+      up = res;                        /* -> lsu */
+      for (c = last; c >= cfirst; c--)
+       {                       /* over each character, from least */
+         if (*c == '.')
+           continue;           /* ignore . [don't step b] */
+         *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)
+       decSetCoeff (dn, set, res, d, &residue, &status);
+
+      /* Finally check for overflow or subnormal and round as needed */
+      decFinalize (dn, set, &residue, &status);
+      /* decNumberShow(dn); */
+    }
+  while (0);                   /* [for break] */
+
+  if (allocres != NULL)
+    free (allocres);           /* drop any storage we used */
+  if (status != 0)
+    decStatus (dn, status, set);
+  return dn;
+}
+
+/* ================================================================== */
+/* 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                                               */
+/*                                                                    */
+/* 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, 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);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs, decNumber * rhs,
+             decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decAddOp (res, lhs, rhs, set, 0, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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.                                   */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberCompare (decNumber * res, decNumber * lhs, decNumber * rhs,
+                 decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decCompareOp (res, lhs, rhs, set, COMPARE, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decDivideOp (res, lhs, rhs, set, DIVIDE, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                       decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decDivideOp (res, lhs, rhs, set, DIVIDEINT, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberMax -- compare two Numbers and return the maximum         */
+/*                                                                    */
+/*   This computes C = A ? B, returning the maximum or A if equal     */
+/*                                                                    */
+/*   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, decNumber * lhs, decNumber * rhs,
+             decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decCompareOp (res, lhs, rhs, set, COMPMAX, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberMin -- compare two Numbers and return the minimum         */
+/*                                                                    */
+/*   This computes C = A ? B, returning the minimum or A if equal     */
+/*                                                                    */
+/*   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, decNumber * lhs, decNumber * rhs,
+             decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decCompareOp (res, lhs, rhs, set, COMPMIN, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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                                               */
+/*                                                                    */
+/* C must have space for set->digits digits.                          */
+/* ------------------------------------------------------------------ */
+/* We simply use AddOp for the subtract, which will do the necessary. */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberMinus (decNumber * res, 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);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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                                               */
+/*                                                                    */
+/* C must have space for set->digits digits.                          */
+/* ------------------------------------------------------------------ */
+/* We simply use 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, 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);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                  decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decMultiplyOp (res, lhs, rhs, set, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberNormalize -- 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.                          */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberNormalize (decNumber * res, decNumber * rhs, decContext * set)
+{
+  decNumber *allocrhs = NULL;  /* non-NULL if rounded rhs allocated */
+  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] */
+
+      /* specials copy through, except NaNs need care */
+      if (decNumberIsNaN (rhs))
+       {
+         decNaNs (res, rhs, NULL, &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, 1, &dropped);      /* normalize in place */
+    }
+  while (0);                   /* end protected */
+
+  if (allocrhs != NULL)
+    free (allocrhs);           /* .. */
+  if (status != 0)
+    decStatus (res, status, set);      /* then report status */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberPower -- raise a number to an integer 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.                          */
+/*                                                                    */
+/* Specification restriction: abs(n) must be <=999999999              */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberPower (decNumber * res, decNumber * lhs,
+               decNumber * rhs, decContext * set)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  decNumber *allocdac = NULL;  /* -> allocated acc buffer, iff used */
+  decNumber *inrhs = rhs;      /* save original rhs */
+  Int reqdigits = set->digits; /* requested DIGITS */
+  Int n;                       /* RHS in binary */
+  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;             /* accumulator */
+  uByte bits = 0;              /* result sign if errors */
+  decContext workset;          /* working context */
+  decNumber dnOne;             /* work value 1... */
+  /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
+  uByte dacbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
+  /* same again for possible 1/lhs calculation */
+  uByte lhsbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
+  decNumber *dac = (decNumber *) dacbuff;      /* -> result accumulator */
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits status, as needed */
+         if (lhs->digits > reqdigits)
+           {
+             alloclhs = decRoundOperand (lhs, set, &status);
+             if (alloclhs == NULL)
+               break;
+             lhs = alloclhs;
+           }
+         /* rounding won't affect the result, but we might signal lostDigits */
+         /* as well as the error for non-integer [x**y would need this too] */
+         if (rhs->digits > reqdigits)
+           {
+             allocrhs = decRoundOperand (rhs, set, &status);
+             if (allocrhs == NULL)
+               break;
+             rhs = allocrhs;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* handle rhs Infinity */
+      if (decNumberIsInfinite (rhs))
+       {
+         status |= DEC_Invalid_operation;      /* bad */
+         break;
+       }
+      /* handle NaNs */
+      if ((lhs->bits | rhs->bits) & (DECNAN | DECSNAN))
+       {
+         decNaNs (res, lhs, rhs, &status);
+         break;
+       }
+
+      /* Original rhs must be an integer that fits and is in range */
+#if DECSUBSET
+      n = decGetInt (inrhs, set);
+#else
+      n = decGetInt (inrhs);
+#endif
+      if (n == BADINT || n > 999999999 || n < -999999999)
+       {
+         status |= DEC_Invalid_operation;
+         break;
+       }
+      if (n < 0)
+       {                       /* negative */
+         n = -n;               /* use the absolute value */
+       }
+      if (decNumberIsNegative (lhs)    /* -x .. */
+         && (n & 0x00000001))
+       bits = DECNEG;          /* .. to an odd power */
+
+      /* handle LHS infinity */
+      if (decNumberIsInfinite (lhs))
+       {                       /* [NaNs already handled] */
+         uByte rbits = rhs->bits;      /* save */
+         decNumberZero (res);
+         if (n == 0)
+           *res->lsu = 1;      /* [-]Inf**0 => 1 */
+         else
+           {
+             if (!(rbits & DECNEG))
+               bits |= DECINF; /* was not a **-n */
+             /* [otherwise will be 0 or -0] */
+             res->bits = bits;
+           }
+         break;
+       }
+
+      /* clone the context */
+      workset = *set;          /* copy all fields */
+      /* calculate the working DIGITS */
+      workset.digits = reqdigits + (inrhs->digits + inrhs->exponent) + 1;
+      /* it's an error if this is more than we can handle */
+      if (workset.digits > DECNUMMAXP)
+       {
+         status |= DEC_Invalid_operation;
+         break;
+       }
+
+      /* workset.digits is the count of digits for the accumulator we need */
+      /* if accumulator is too long for local storage, then allocate */
+      needbytes =
+       sizeof (decNumber) + (D2U (workset.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 */
+       }
+      decNumberZero (dac);     /* acc=1 */
+      *dac->lsu = 1;           /* .. */
+
+      if (n == 0)
+       {                       /* x**0 is usually 1 */
+         /* 0**0 is bad unless subset, when it becomes 1 */
+         if (ISZERO (lhs)
+#if DECSUBSET
+             && set->extended
+#endif
+           )
+           status |= DEC_Invalid_operation;
+         else
+           decNumberCopy (res, dac);   /* copy the 1 */
+         break;
+       }
+
+      /* if a negative power we'll need the constant 1, and if not subset */
+      /* we'll invert the lhs now rather than inverting the result later */
+      if (decNumberIsNegative (rhs))
+       {                       /* was a **-n [hence digits>0] */
+         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, &workset, DIVIDE, &status);
+             if (alloclhs != NULL)
+               {
+                 free (alloclhs);      /* done with intermediate */
+                 alloclhs = NULL;      /* indicate freed */
+               }
+             /* now locate or allocate space for the inverted lhs */
+             if (needbytes > sizeof (lhsbuff))
+               {
+                 alloclhs = (decNumber *) malloc (needbytes);
+                 if (alloclhs == NULL)
+                   {           /* hopeless -- abandon */
+                     status |= DEC_Insufficient_storage;
+                     break;
+                   }
+                 lhs = alloclhs;       /* use the allocated space */
+               }
+             else
+               lhs = (decNumber *) lhsbuff;    /* use stack storage */
+             /* [lhs now points to buffer or allocated storage] */
+             decNumberCopy (lhs, dac); /* copy the 1/lhs */
+             decNumberCopy (dac, &dnOne);      /* restore acc=1 */
+#if DECSUBSET
+           }
+#endif
+       }
+
+      /* Raise-to-the-power loop... */
+      seenbit = 0;             /* set once we've seen a 1-bit */
+      for (i = 1;; i++)
+       {                       /* for each bit [top bit ignored] */
+         /* abandon if we have 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, we're off */
+             decMultiplyOp (dac, dac, lhs, &workset, &status); /* dac=dac*x */
+           }
+         if (i == 31)
+           break;              /* that was the last bit */
+         if (!seenbit)
+           continue;           /* we don't have to square 1 */
+         decMultiplyOp (dac, dac, dac, &workset, &status);     /* dac=dac*dac [square] */
+       }                       /*i *//* 32 bits */
+
+      /* complete internal overflow or underflow processing */
+      if (status & (DEC_Overflow | DEC_Subnormal))
+       {
+#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 workset.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, &workset, DIVIDE, &status);
+       }
+#endif
+
+      /* 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, 0, &dropped);     /* trailing zeros */
+#endif
+    }
+  while (0);                   /* end protected */
+
+  if (allocdac != NULL)
+    free (allocdac);           /* drop any storage we used */
+  if (allocrhs != NULL)
+    free (allocrhs);           /* .. */
+  if (alloclhs != NULL)
+    free (alloclhs);           /* .. */
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                  decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decQuantizeOp (res, lhs, rhs, set, 1, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                 decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decQuantizeOp (res, lhs, rhs, set, 0, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                   decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decDivideOp (res, lhs, rhs, set, REMAINDER, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs,
+                       decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+  decDivideOp (res, lhs, rhs, set, REMNEAR, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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, decNumber * lhs, decNumber * rhs)
+{
+  uByte merged;                        /* merged flags */
+  Unit ret = 0;                        /* return value */
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, DECUNUSED))
+    return res;
+#endif
+
+  merged = (lhs->bits | rhs->bits) & DECSPECIAL;
+  if (merged)
+    {
+      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 */
+  *res->lsu = ret;
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* 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.                                  */
+/*                                                                    */
+/* % [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, decNumber * rhs, decContext * set)
+{
+  decContext workset, approxset;       /* work contexts */
+  decNumber dzero;             /* used for constant zero */
+  Int maxp = set->digits + 2;  /* largest working precision */
+  Int residue = 0;             /* rounding residue */
+  uInt status = 0, ignore = 0; /* status accumulators */
+  Int exp;                     /* working exponent */
+  Int ideal;                   /* ideal (preferred) exponent */
+  uInt needbytes;              /* work */
+  Int dropped;                 /* .. */
+
+  decNumber *allocrhs = NULL;  /* non-NULL if rounded rhs allocated */
+  /* buffer for f [needs +1 in case DECBUFFER 0] */
+  uByte buff[sizeof (decNumber) + (D2U (DECBUFFER + 1) - 1) * sizeof (Unit)];
+  /* buffer for a [needs +2 to match maxp] */
+  uByte bufa[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
+  /* buffer for temporary, b [must be same size as a] */
+  uByte bufb[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
+  decNumber *allocbuff = NULL; /* -> allocated buff, iff allocated */
+  decNumber *allocbufa = NULL; /* -> allocated bufa, iff allocated */
+  decNumber *allocbufb = NULL; /* -> allocated bufb, iff allocated */
+  decNumber *f = (decNumber *) buff;   /* reduced fraction */
+  decNumber *a = (decNumber *) bufa;   /* approximation to result */
+  decNumber *b = (decNumber *) bufb;   /* intermediate result */
+  /* buffer for temporary variable, up to 3 digits */
+  uByte buft[sizeof (decNumber) + (D2U (3) - 1) * sizeof (Unit)];
+  decNumber *t = (decNumber *) 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 we keep them separate for clarity.] */
+             rhs = allocrhs;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* handle infinities and NaNs */
+      if (rhs->bits & DECSPECIAL)
+       {
+         if (decNumberIsInfinite (rhs))
+           {                   /* an infinity */
+             if (decNumberIsNegative (rhs))
+               status |= DEC_Invalid_operation;
+             else
+               decNumberCopy (res, rhs);       /* +Infinity */
+           }
+         else
+           decNaNs (res, rhs, NULL, &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] */
+         break;
+       }
+
+      /* any other -x is an oops */
+      if (decNumberIsNegative (rhs))
+       {
+         status |= DEC_Invalid_operation;
+         break;
+       }
+
+      /* we need space for three working variables */
+      /*   f -- the same precision as the RHS, reduced to 0.01->0.99... */
+      /*   a -- Hull's approx -- precision, when assigned, is */
+      /*        currentprecision (we allow +2 for use as temporary) */
+      /*   b -- intermediate temporary result */
+      /* if any is too long for local storage, then allocate */
+      needbytes =
+       sizeof (decNumber) + (D2U (rhs->digits) - 1) * sizeof (Unit);
+      if (needbytes > 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 > 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 space */
+         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 contexts (the second is used for Numerical */
+      /* Turing assignment) */
+      decContextDefault (&workset, DEC_INIT_DECIMAL64);
+      decContextDefault (&approxset, DEC_INIT_DECIMAL64);
+      approxset.digits = set->digits;  /* approx's length */
+
+      /* [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 = set->digits;    /* 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 */
+         /* assign to approx [round to length] */
+         decAddOp (a, &dzero, a, &approxset, 0, &ignore);
+         if (workset.digits == maxp)
+           break;              /* just did final */
+       }                       /* loop */
+
+      /* a is now at currentprecision 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 */
+      workset.digits--;                /* maxp-1 is OK now */
+      t->exponent = -set->digits - 1;  /* make 0.5 ulp */
+      decNumberCopy (b, a);
+      decAddOp (b, b, t, &workset, DECNEG, &ignore);   /* b = a - 0.5 ulp */
+      workset.round = DEC_ROUND_UP;
+      decMultiplyOp (b, b, b, &workset, &ignore);      /* b = mulru(b, b) */
+      decCompareOp (b, f, b, &workset, COMPARE, &ignore);      /* b ? f, reversed */
+      if (decNumberIsNegative (b))
+       {                       /* f < b [i.e., b > f] */
+         /* this is the more common adjustment, though both are rare */
+         t->exponent++;        /* make 1.0 ulp */
+         t->lsu[0] = 1;        /* .. */
+         decAddOp (a, a, t, &workset, DECNEG, &ignore);        /* a = a - 1 ulp */
+         /* assign to approx [round to length] */
+         decAddOp (a, &dzero, a, &approxset, 0, &ignore);
+       }
+      else
+       {
+         decNumberCopy (b, a);
+         decAddOp (b, b, t, &workset, 0, &ignore);     /* b = a + 0.5 ulp */
+         workset.round = DEC_ROUND_DOWN;
+         decMultiplyOp (b, b, b, &workset, &ignore);   /* b = mulrd(b, b) */
+         decCompareOp (b, b, f, &workset, COMPARE, &ignore);   /* b ? f */
+         if (decNumberIsNegative (b))
+           {                   /* b < f */
+             t->exponent++;    /* make 1.0 ulp */
+             t->lsu[0] = 1;    /* .. */
+             decAddOp (a, a, t, &workset, 0, &ignore); /* a = a + 1 ulp */
+             /* assign to approx [round to length] */
+             decAddOp (a, &dzero, a, &approxset, 0, &ignore);
+           }
+       }
+      /* [no errors are possible in the above, and rounding/inexact during */
+      /* estimation are irrelevant, so status was not accumulated] */
+
+      /* Here, 0.1 <= a < 1  [Hull] */
+      a->exponent += exp / 2;  /* set correct exponent */
+
+      /* Process Subnormals */
+      decFinalize (a, set, &residue, &status);
+
+      /* count dropable zeros [after any subnormal rounding] */
+      decNumberCopy (b, a);
+      decTrim (b, 1, &dropped);        /* [drops trailing zeros] */
+
+      /* Finally set Inexact and Rounded.  The answer can only be exact if */
+      /* it is short enough so that squaring it could fit in set->digits, */
+      /* so this is the only (relatively rare) time we have to check */
+      /* carefully */
+      if (b->digits * 2 - 1 > set->digits)
+       {                       /* cannot fit */
+         status |= DEC_Inexact | DEC_Rounded;
+       }
+      else
+       {                       /* could be exact/unrounded */
+         uInt mstatus = 0;     /* local status */
+         decMultiplyOp (b, b, b, &workset, &mstatus);  /* try the multiply */
+         if (mstatus != 0)
+           {                   /* result won't fit */
+             status |= DEC_Inexact | DEC_Rounded;
+           }
+         else
+           {                   /* plausible */
+             decCompareOp (t, b, rhs, &workset, COMPARE, &mstatus);    /* b ? rhs */
+             if (!ISZERO (t))
+               {
+                 status |= DEC_Inexact | DEC_Rounded;
+               }
+             else
+               {               /* is Exact */
+                 /* here, dropped is the count of trailing zeros in 'a' */
+                 /* use closest exponent to ideal... */
+                 Int todrop = ideal - a->exponent;     /* most we can drop */
+
+                 if (todrop < 0)
+                   {           /* ideally would add 0s */
+                     status |= DEC_Rounded;
+                   }
+                 else
+                   {           /* unrounded */
+                     if (dropped < todrop)
+                       todrop = dropped;       /* clamp to those available */
+                     if (todrop > 0)
+                       {       /* OK, some to drop */
+                         decShiftToLeast (a->lsu, D2U (a->digits), todrop);
+                         a->exponent += todrop;        /* maintain numerical value */
+                         a->digits -= todrop;  /* new length */
+                       }
+                   }
+               }
+           }
+       }
+      decNumberCopy (res, a);  /* assume this is the result */
+    }
+  while (0);                   /* end protected */
+
+  if (allocbuff != NULL)
+    free (allocbuff);          /* drop any storage we used */
+  if (allocbufa != NULL)
+    free (allocbufa);          /* .. */
+  if (allocbufb != NULL)
+    free (allocbufb);          /* .. */
+  if (allocrhs != NULL)
+    free (allocrhs);           /* .. */
+  if (status != 0)
+    decStatus (res, status, set);      /* then report status */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberSubtract -- subtract 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.                          */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberSubtract (decNumber * res, decNumber * lhs,
+                  decNumber * rhs, decContext * set)
+{
+  uInt status = 0;             /* accumulator */
+
+  decAddOp (res, lhs, rhs, set, DECNEG, &status);
+  if (status != 0)
+    decStatus (res, status, set);
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberToIntegralValue -- round-to-integral-value                */
+/*                                                                    */
+/*   res is the result                                                */
+/*   rhs is input number                                              */
+/*   set is the context                                               */
+/*                                                                    */
+/* res must have space for any value of rhs.                          */
+/*                                                                    */
+/* This implements the IEEE special operator and therefore treats     */
+/* special values as valid, and also never sets Inexact.  For finite  */
+/* numbers it returns rescale(rhs, 0) if rhs->exponent is <0.         */
+/* Otherwise the result is rhs (so no error is possible).             */
+/*                                                                    */
+/* The context is used for rounding mode and status after sNaN, but   */
+/* the digits setting is ignored.                                     */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberToIntegralValue (decNumber * res, decNumber * rhs, decContext * set)
+{
+  decNumber dn;
+  decContext workset;          /* working context */
+
+#if DECCHECK
+  if (decCheckOperands (res, DECUNUSED, rhs, set))
+    return res;
+#endif
+
+  /* handle infinities and NaNs */
+  if (rhs->bits & DECSPECIAL)
+    {
+      uInt status = 0;
+      if (decNumberIsInfinite (rhs))
+       decNumberCopy (res, rhs);       /* an Infinity */
+      else
+       decNaNs (res, rhs, NULL, &status);      /* a NaN */
+      if (status != 0)
+       decStatus (res, status, set);
+      return res;
+    }
+
+  /* we have a finite number; no error possible */
+  if (rhs->exponent >= 0)
+    return decNumberCopy (res, rhs);
+  /* that was easy, but if negative exponent we have work to do... */
+  workset = *set;              /* clone rounding, etc. */
+  workset.digits = rhs->digits;        /* no length rounding */
+  workset.traps = 0;           /* no traps */
+  decNumberZero (&dn);         /* make a number with exponent 0 */
+  return decNumberQuantize (res, rhs, &dn, &workset);
+}
+
+/* ================================================================== */
+/* Utility routines                                                   */
+/* ================================================================== */
+
+/* ------------------------------------------------------------------ */
+/* decNumberCopy -- copy a number                                     */
+/*                                                                    */
+/*   dest is the target decNumber                                     */
+/*   src  is the source decNumber                                     */
+/*   returns dest                                                     */
+/*                                                                    */
+/* (dest==src is allowed and is a no-op)                              */
+/* All fields are updated as required.  This is a utility operation,  */
+/* so special values are unchanged and no error is possible.          */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberCopy (decNumber * dest, decNumber * src)
+{
+
+#if DECCHECK
+  if (src == NULL)
+    return decNumberZero (dest);
+#endif
+
+  if (dest == src)
+    return dest;               /* no copy required */
+
+  /* We use explicit assignments here as structure assignment can copy */
+  /* more than just the lsu (for small DECDPUN).  This would not affect */
+  /* the value of the results, but would disturb test harness spill */
+  /* checking. */
+  dest->bits = src->bits;
+  dest->exponent = src->exponent;
+  dest->digits = src->digits;
+  dest->lsu[0] = src->lsu[0];
+  if (src->digits > DECDPUN)
+    {                          /* more Units to come */
+      Unit *s, *d, *smsup;     /* work */
+      /* memcpy for the remaining Units would be safe as they cannot */
+      /* overlap.  However, this explicit loop is faster in short cases. */
+      d = dest->lsu + 1;       /* -> first destination */
+      smsup = src->lsu + D2U (src->digits);    /* -> source msu+1 */
+      for (s = src->lsu + 1; s < smsup; s++, d++)
+       *d = *s;
+    }
+  return dest;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberTrim -- remove insignificant zeros                        */
+/*                                                                    */
+/*   dn is the number to trim                                         */
+/*   returns dn                                                       */
+/*                                                                    */
+/* All fields are updated as required.  This is a utility operation,  */
+/* so special values are unchanged and no error is possible.          */
+/* ------------------------------------------------------------------ */
+decNumber *
+decNumberTrim (decNumber * dn)
+{
+  Int dropped;                 /* work */
+  return decTrim (dn, 0, &dropped);
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberVersion -- return the name and version of this module     */
+/*                                                                    */
+/* No error is possible.                                              */
+/* ------------------------------------------------------------------ */
+const char *
+decNumberVersion (void)
+{
+  return DECVERSION;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNumberZero -- set a number to 0                                 */
+/*                                                                    */
+/*   dn is the number to set, with space for one digit                */
+/*   returns dn                                                       */
+/*                                                                    */
+/* No error is possible.                                              */
+/* ------------------------------------------------------------------ */
+/* Memset is not used as it is much slower in some environments. */
+decNumber *
+decNumberZero (decNumber * dn)
+{
+
+#if DECCHECK
+  if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
+    return dn;
+#endif
+
+  dn->bits = 0;
+  dn->exponent = 0;
+  dn->digits = 1;
+  dn->lsu[0] = 0;
+  return dn;
+}
+
+/* ================================================================== */
+/* Local routines                                                     */
+/* ================================================================== */
+
+/* ------------------------------------------------------------------ */
+/* decToString -- lay out a number into a string                      */
+/*                                                                    */
+/*   dn     is the number to lay out                                  */
+/*   string is where to lay out the number                            */
+/*   eng    is 1 if Engineering, 0 if Scientific                      */
+/*                                                                    */
+/* str must be at least dn->digits+14 characters long                 */
+/* No error is possible.                                              */
+/*                                                                    */
+/* Note that this routine can generate a -0 or 0.000.  These are      */
+/* never generated in subset to-number or arithmetic, but can occur   */
+/* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234).              */
+/* ------------------------------------------------------------------ */
+/* If DECCHECK is enabled the string "?" is returned if a number is */
+/* invalid. */
+
+/* TODIGIT -- macro to remove the leading digit from the unsigned */
+/* integer u at column cut (counting from the right, LSD=0) and place */
+/* it as an ASCII character into the character pointed to by c.  Note */
+/* that cut must be <= 9, and the maximum value for u is 2,000,000,000 */
+/* (as is needed for negative exponents of subnormals).  The unsigned */
+/* integer pow is used as a temporary variable. */
+#define TODIGIT(u, cut, c) {            \
+  *(c)='0';                             \
+  pow=powers[cut]*2;                    \
+  if ((u)>pow) {                        \
+    pow*=4;                             \
+    if ((u)>=pow) {(u)-=pow; *(c)+=8;}  \
+    pow/=2;                             \
+    if ((u)>=pow) {(u)-=pow; *(c)+=4;}  \
+    pow/=2;                             \
+    }                                   \
+  if ((u)>=pow) {(u)-=pow; *(c)+=2;}    \
+  pow/=2;                               \
+  if ((u)>=pow) {(u)-=pow; *(c)+=1;}    \
+  }
+
+static void
+decToString (decNumber * dn, char *string, Flag eng)
+{
+  Int exp = dn->exponent;      /* local copy */
+  Int e;                       /* E-part value */
+  Int pre;                     /* digits before the '.' */
+  Int cut;                     /* for counting digits in a Unit */
+  char *c = string;            /* work [output pointer] */
+  Unit *up = dn->lsu + D2U (dn->digits) - 1;   /* -> msu [input pointer] */
+  uInt u, pow;                 /* work */
+
+#if DECCHECK
+  if (decCheckOperands (DECUNUSED, dn, DECUNUSED, DECUNUSED))
+    {
+      strcpy (string, "?");
+      return;
+    }
+#endif
+
+  if (decNumberIsNegative (dn))
+    {                          /* Negatives get a minus (except */
+      *c = '-';                        /* NaNs, which remove the '-' below) */
+      c++;
+    }
+  if (dn->bits & DECSPECIAL)
+    {                          /* Is a special value */
+      if (decNumberIsInfinite (dn))
+       {
+         strcpy (c, "Infinity");
+         return;
+       }
+      /* a NaN */
+      if (dn->bits & DECSNAN)
+       {                       /* signalling NaN */
+         *c = 's';
+         c++;
+       }
+      strcpy (c, "NaN");
+      c += 3;                  /* step past */
+      /* if not a clean non-zero coefficient, that's all we have in a */
+      /* NaN string */
+      if (exp != 0 || (*dn->lsu == 0 && dn->digits == 1))
+       return;
+      /* [drop through to add integer] */
+    }
+
+  /* calculate how many digits in msu, and hence first cut */
+  cut = dn->digits % DECDPUN;
+  if (cut == 0)
+    cut = DECDPUN;             /* msu is full */
+  cut--;                       /* power of ten for digit */
+
+  if (exp == 0)
+    {                          /* simple integer [common fastpath, */
+      /*   used for NaNs, too] */
+      for (; up >= dn->lsu; up--)
+       {                       /* each Unit from msu */
+         u = *up;              /* contains DECDPUN digits to lay out */
+         for (; cut >= 0; c++, cut--)
+           TODIGIT (u, cut, c);
+         cut = DECDPUN - 1;    /* next Unit has all digits */
+       }
+      *c = '\0';               /* terminate the string */
+      return;
+    }
+
+  /* non-0 exponent -- assume plain form */
+  pre = dn->digits + exp;      /* digits before '.' */
+  e = 0;                       /* no E */
+  if ((exp > 0) || (pre < -5))
+    {                          /* need exponential form */
+      e = exp + dn->digits - 1;        /* calculate E value */
+      pre = 1;                 /* assume one digit before '.' */
+      if (eng && (e != 0))
+       {                       /* may need to adjust */
+         Int adj;              /* adjustment */
+         /* The C remainder operator is undefined for negative numbers, so */
+         /* we must use positive remainder calculation here */
+         if (e < 0)
+           {
+             adj = (-e) % 3;
+             if (adj != 0)
+               adj = 3 - adj;
+           }
+         else
+           {                   /* e>0 */
+             adj = e % 3;
+           }
+         e = e - adj;
+         /* if we are dealing with zero we will use exponent which is a */
+         /* multiple of three, as expected, but there will only be the */
+         /* one zero before the E, still.  Otherwise note the padding. */
+         if (!ISZERO (dn))
+           pre += adj;
+         else
+           {                   /* is zero */
+             if (adj != 0)
+               {               /* 0.00Esnn needed */
+                 e = e + 3;
+                 pre = -(2 - adj);
+               }
+           }                   /* zero */
+       }                       /* eng */
+    }
+
+  /* lay out the digits of the coefficient, adding 0s and . as needed */
+  u = *up;
+  if (pre > 0)
+    {                          /* xxx.xxx or xx00 (engineering) form */
+      for (; pre > 0; pre--, c++, cut--)
+       {
+         if (cut < 0)
+           {                   /* need new Unit */
+             if (up == dn->lsu)
+               break;          /* out of input digits (pre>digits) */
+             up--;
+             cut = DECDPUN - 1;
+             u = *up;
+           }
+         TODIGIT (u, cut, c);
+       }
+      if (up > dn->lsu || (up == dn->lsu && cut >= 0))
+       {                       /* more to come, after '.' */
+         *c = '.';
+         c++;
+         for (;; c++, cut--)
+           {
+             if (cut < 0)
+               {               /* need new Unit */
+                 if (up == dn->lsu)
+                   break;      /* out of input digits */
+                 up--;
+                 cut = DECDPUN - 1;
+                 u = *up;
+               }
+             TODIGIT (u, cut, c);
+           }
+       }
+      else
+       for (; pre > 0; pre--, c++)
+         *c = '0';             /* 0 padding (for engineering) needed */
+    }
+  else
+    {                          /* 0.xxx or 0.000xxx form */
+      *c = '0';
+      c++;
+      *c = '.';
+      c++;
+      for (; pre < 0; pre++, c++)
+       *c = '0';               /* add any 0's after '.' */
+      for (;; c++, cut--)
+       {
+         if (cut < 0)
+           {                   /* need new Unit */
+             if (up == dn->lsu)
+               break;          /* out of input digits */
+             up--;
+             cut = DECDPUN - 1;
+             u = *up;
+           }
+         TODIGIT (u, cut, c);
+       }
+    }
+
+  /* Finally add the E-part, if needed.  It will never be 0, has a
+     base maximum and minimum of +999999999 through -999999999, but
+     could range down to -1999999998 for subnormal numbers */
+  if (e != 0)
+    {
+      Flag had = 0;            /* 1=had non-zero */
+      *c = 'E';
+      c++;
+      *c = '+';
+      c++;                     /* assume positive */
+      u = e;                   /* .. */
+      if (e < 0)
+       {
+         *(c - 1) = '-';       /* oops, need - */
+         u = -e;               /* uInt, please */
+       }
+      /* layout the exponent (_itoa is not ANSI C) */
+      for (cut = 9; cut >= 0; cut--)
+       {
+         TODIGIT (u, cut, c);
+         if (*c == '0' && !had)
+           continue;           /* skip leading zeros */
+         had = 1;              /* had non-0 */
+         c++;                  /* step for next */
+       }                       /* cut */
+    }
+  *c = '\0';                   /* terminate the string (all paths) */
+  return;
+}
+
+/* ------------------------------------------------------------------ */
+/* decAddOp -- add/subtract operation                                 */
+/*                                                                    */
+/*   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                                               */
+/*   negate is DECNEG if rhs should be negated, or 0 otherwise        */
+/*   status accumulates status for the caller                         */
+/*                                                                    */
+/* C must have space for set->digits digits.                          */
+/* ------------------------------------------------------------------ */
+/* If possible, we calculate the coefficient directly into C.         */
+/* However, if:                                                       */
+/*   -- we need a digits+1 calculation because numbers are unaligned  */
+/*      and span more than set->digits digits                         */
+/*   -- a carry to digits+1 digits looks possible                     */
+/*   -- C is the same as A or B, and the result would destructively   */
+/*      overlap the A or B coefficient                                */
+/* then we must calculate into a temporary buffer.  In this latter    */
+/* case we use the local (stack) buffer if possible, and only if too  */
+/* long for that do we resort to malloc.                              */
+/*                                                                    */
+/* Misalignment is handled as follows:                                */
+/*   Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp.    */
+/*   BPad: Apply the padding by a combination of shifting (whole      */
+/*         units) and multiplication (part units).                    */
+/*                                                                    */
+/* Addition, especially x=x+1, is speed-critical, so we take pains    */
+/* to make returning as fast as possible, by flagging any allocation. */
+/* ------------------------------------------------------------------ */
+static decNumber *
+decAddOp (decNumber * res, decNumber * lhs,
+         decNumber * rhs, decContext * set, uByte negate, uInt * status)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  Int rhsshift;                        /* working shift (in Units) */
+  Int maxdigits;               /* longest logical length */
+  Int mult;                    /* multiplier */
+  Int residue;                 /* rounding accumulator */
+  uByte bits;                  /* result bits */
+  Flag diffsign;               /* non-0 if arguments have different sign */
+  Unit *acc;                   /* accumulator for result */
+  Unit accbuff[D2U (DECBUFFER + 1)];   /* local buffer [+1 is for possible */
+  /* final carry digit or DECBUFFER=0] */
+  Unit *allocacc = NULL;       /* -> allocated acc buffer, iff allocated */
+  Flag alloced = 0;            /* set non-0 if any allocations */
+  Int reqdigits = set->digits; /* local copy; requested DIGITS */
+  uByte merged;                        /* merged flags */
+  Int padding;                 /* work */
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits status, as needed */
+         if (lhs->digits > reqdigits)
+           {
+             alloclhs = decRoundOperand (lhs, set, status);
+             if (alloclhs == NULL)
+               break;
+             lhs = alloclhs;
+             alloced = 1;
+           }
+         if (rhs->digits > reqdigits)
+           {
+             allocrhs = decRoundOperand (rhs, set, status);
+             if (allocrhs == NULL)
+               break;
+             rhs = allocrhs;
+             alloced = 1;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* note whether signs differ */
+      diffsign = (Flag) ((lhs->bits ^ rhs->bits ^ negate) & DECNEG);
+
+      /* handle infinities and NaNs */
+      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
+      if (merged)
+       {                       /* a special bit set */
+         if (merged & (DECSNAN | DECNAN))      /* a NaN */
+           decNaNs (res, lhs, rhs, status);
+         else
+           {                   /* one or two infinities */
+             if (decNumberIsInfinite (lhs))
+               {               /* LHS is infinity */
+                 /* two infinities with different signs is invalid */
+                 if (decNumberIsInfinite (rhs) && diffsign)
+                   {
+                     *status |= DEC_Invalid_operation;
+                     break;
+                   }
+                 bits = lhs->bits & DECNEG;    /* get sign from LHS */
+               }
+             else
+               bits = (rhs->bits ^ negate) & DECNEG;   /* RHS must be Infinity */
+             bits |= DECINF;
+             decNumberZero (res);
+             res->bits = bits; /* set +/- infinity */
+           }                   /* an infinity */
+         break;
+       }
+
+      /* Quick exit for add 0s; return the non-0, modified as need be */
+      if (ISZERO (lhs))
+       {
+         Int adjust;           /* work */
+         Int lexp = lhs->exponent;     /* save in case LHS==RES */
+         bits = lhs->bits;     /* .. */
+         residue = 0;          /* clear accumulator */
+         decCopyFit (res, rhs, set, &residue, status); /* copy (as needed) */
+         res->bits ^= negate;  /* flip if rhs was negated */
+#if DECSUBSET
+         if (set->extended)
+           {                   /* exponents on zeros count */
+#endif
+             /* exponent will be the lower of the two */
+             adjust = lexp - res->exponent;    /* adjustment needed [if -ve] */
+             if (ISZERO (res))
+               {               /* both 0: special IEEE 854 rules */
+                 if (adjust < 0)
+                   res->exponent = lexp;       /* set exponent */
+                 /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
+                 if (diffsign)
+                   {
+                     if (set->round != DEC_ROUND_FLOOR)
+                       res->bits = 0;
+                     else
+                       res->bits = DECNEG;     /* preserve 0 sign */
+                   }
+               }
+             else
+               {               /* non-0 res */
+                 if (adjust < 0)
+                   {           /* 0-padding needed */
+                     if ((res->digits - adjust) > set->digits)
+                       {
+                         adjust = res->digits - set->digits;   /* to fit exactly */
+                         *status |= DEC_Rounded;       /* [but exact] */
+                       }
+                     res->digits =
+                       decShiftToMost (res->lsu, res->digits, -adjust);
+                     res->exponent += adjust;  /* set the exponent. */
+                   }
+               }               /* non-0 res */
+#if DECSUBSET
+           }                   /* extended */
+#endif
+         decFinish (res, set, &residue, status);       /* clean and finalize */
+         break;
+       }
+
+      if (ISZERO (rhs))
+       {                       /* [lhs is non-zero] */
+         Int adjust;           /* work */
+         Int rexp = rhs->exponent;     /* save in case RHS==RES */
+         bits = rhs->bits;     /* be clean */
+         residue = 0;          /* clear accumulator */
+         decCopyFit (res, lhs, set, &residue, status); /* copy (as needed) */
+#if DECSUBSET
+         if (set->extended)
+           {                   /* exponents on zeros count */
+#endif
+             /* exponent will be the lower of the two */
+             /* [0-0 case handled above] */
+             adjust = rexp - res->exponent;    /* adjustment needed [if -ve] */
+             if (adjust < 0)
+               {               /* 0-padding needed */
+                 if ((res->digits - adjust) > set->digits)
+                   {
+                     adjust = res->digits - set->digits;       /* to fit exactly */
+                     *status |= DEC_Rounded;   /* [but exact] */
+                   }
+                 res->digits =
+                   decShiftToMost (res->lsu, res->digits, -adjust);
+                 res->exponent += adjust;      /* set the exponent. */
+               }
+#if DECSUBSET
+           }                   /* extended */
+#endif
+         decFinish (res, set, &residue, status);       /* clean and finalize */
+         break;
+       }
+      /* [both fastpath and mainpath code below assume these cases */
+      /* (notably 0-0) have already been handled] */
+
+      /* calculate the padding needed to align the operands */
+      padding = rhs->exponent - lhs->exponent;
+
+      /* Fastpath cases where the numbers are aligned and normal, the RHS */
+      /* is all in one unit, no operand rounding is needed, and no carry, */
+      /* lengthening, or borrow is needed */
+      if (rhs->digits <= DECDPUN && padding == 0 && rhs->exponent >= set->emin /* [some normals drop through] */
+         && rhs->digits <= reqdigits && lhs->digits <= reqdigits)
+       {
+         Int partial = *lhs->lsu;
+         if (!diffsign)
+           {                   /* adding */
+             Int maxv = DECDPUNMAX;    /* highest no-overflow */
+             if (lhs->digits < DECDPUN)
+               maxv = powers[lhs->digits] - 1;
+             partial += *rhs->lsu;
+             if (partial <= maxv)
+               {               /* no carry */
+                 if (res != lhs)
+                   decNumberCopy (res, lhs);   /* not in place */
+                 *res->lsu = (Unit) partial;   /* [copy could have overwritten RHS] */
+                 break;
+               }
+             /* else drop out for careful add */
+           }
+         else
+           {                   /* signs differ */
+             partial -= *rhs->lsu;
+             if (partial > 0)
+               {               /* no borrow needed, and non-0 result */
+                 if (res != lhs)
+                   decNumberCopy (res, lhs);   /* not in place */
+                 *res->lsu = (Unit) partial;
+                 /* this could have reduced digits [but result>0] */
+                 res->digits = decGetDigits (res->lsu, D2U (res->digits));
+                 break;
+               }
+             /* else drop out for careful subtract */
+           }
+       }
+
+      /* Now align (pad) the lhs or rhs so we can add or subtract them, as
+         necessary.  If one number is much larger than the other (that is,
+         if in plain form there is a least one digit between the lowest
+         digit or one and the highest of the other) we need to pad with up
+         to DIGITS-1 trailing zeros, and then apply rounding (as exotic
+         rounding modes may be affected by the residue).
+       */
+      rhsshift = 0;            /* rhs shift to left (padding) in Units */
+      bits = lhs->bits;                /* assume sign is that of LHS */
+      mult = 1;                        /* likely multiplier */
+
+      /* if padding==0 the operands are aligned; no padding needed */
+      if (padding != 0)
+       {
+         /* some padding needed */
+         /* We always pad the RHS, as we can then effect any required */
+         /* padding by a combination of shifts and a multiply */
+         Flag swapped = 0;
+         if (padding < 0)
+           {                   /* LHS needs the padding */
+             decNumber *t;
+             padding = -padding;       /* will be +ve */
+             bits = (uByte) (rhs->bits ^ negate);      /* assumed sign is now that of RHS */
+             t = lhs;
+             lhs = rhs;
+             rhs = t;
+             swapped = 1;
+           }
+
+         /* If, after pad, rhs would be longer than lhs by digits+1 or */
+         /* more then lhs cannot affect the answer, except as a residue, */
+         /* so we only need to pad up to a length of DIGITS+1. */
+         if (rhs->digits + padding > lhs->digits + reqdigits + 1)
+           {
+             /* The RHS is sufficient */
+             /* for residue we use the relative sign indication... */
+             Int shift = reqdigits - rhs->digits;      /* left shift needed */
+             residue = 1;      /* residue for rounding */
+             if (diffsign)
+               residue = -residue;     /* signs differ */
+             /* copy, shortening if necessary */
+             decCopyFit (res, rhs, set, &residue, status);
+             /* if it was already shorter, then need to pad with zeros */
+             if (shift > 0)
+               {
+                 res->digits = decShiftToMost (res->lsu, res->digits, shift);
+                 res->exponent -= shift;       /* adjust the exponent. */
+               }
+             /* flip the result sign if unswapped and rhs was negated */
+             if (!swapped)
+               res->bits ^= negate;
+             decFinish (res, set, &residue, status);   /* done */
+             break;
+           }
+
+         /* LHS digits may affect result */
+         rhsshift = D2U (padding + 1) - 1;     /* this much by Unit shift .. */
+         mult = powers[padding - (rhsshift * DECDPUN)];        /* .. this by multiplication */
+       }                       /* padding needed */
+
+      if (diffsign)
+       mult = -mult;           /* signs differ */
+
+      /* determine the longer operand */
+      maxdigits = rhs->digits + padding;       /* virtual length of RHS */
+      if (lhs->digits > maxdigits)
+       maxdigits = lhs->digits;
+
+      /* Decide on the result buffer to use; if possible place directly */
+      /* into result. */
+      acc = res->lsu;          /* assume build direct */
+      /* If destructive overlap, or the number is too long, or a carry or */
+      /* borrow to DIGITS+1 might be possible we must use a buffer. */
+      /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
+      if ((maxdigits >= reqdigits)     /* is, or could be, too large */
+         || (res == rhs && rhsshift > 0))
+       {                       /* destructive overlap */
+         /* buffer needed; choose it */
+         /* we'll need units for maxdigits digits, +1 Unit for carry or borrow */
+         Int need = D2U (maxdigits) + 1;
+         acc = accbuff;        /* assume use local buffer */
+         if (need * sizeof (Unit) > sizeof (accbuff))
+           {
+             allocacc = (Unit *) malloc (need * sizeof (Unit));
+             if (allocacc == NULL)
+               {               /* hopeless -- abandon */
+                 *status |= DEC_Insufficient_storage;
+                 break;
+               }
+             acc = allocacc;
+             alloced = 1;
+           }
+       }
+
+      res->bits = (uByte) (bits & DECNEG);     /* it's now safe to overwrite.. */
+      res->exponent = lhs->exponent;   /* .. operands (even if aliased) */
+
+#if DECTRACE
+      decDumpAr ('A', lhs->lsu, D2U (lhs->digits));
+      decDumpAr ('B', rhs->lsu, D2U (rhs->digits));
+      printf ("  :h: %d %d\n", rhsshift, mult);
+#endif
+
+      /* add [A+B*m] or subtract [A+B*(-m)] */
+      res->digits = decUnitAddSub (lhs->lsu, D2U (lhs->digits), rhs->lsu, D2U (rhs->digits), rhsshift, acc, mult) * DECDPUN;   /* [units -> digits] */
+      if (res->digits < 0)
+       {                       /* we borrowed */
+         res->digits = -res->digits;
+         res->bits ^= DECNEG;  /* flip the sign */
+       }
+#if DECTRACE
+      decDumpAr ('+', acc, D2U (res->digits));
+#endif
+
+      /* If we used a buffer we need to copy back, possibly shortening */
+      /* (If we didn't use buffer it must have fit, so can't need rounding */
+      /* and residue must be 0.) */
+      residue = 0;             /* clear accumulator */
+      if (acc != res->lsu)
+       {
+#if DECSUBSET
+         if (set->extended)
+           {                   /* round from first significant digit */
+#endif
+             /* remove leading zeros that we added due to rounding up to */
+             /* integral Units -- before the test for rounding. */
+             if (res->digits > reqdigits)
+               res->digits = decGetDigits (acc, D2U (res->digits));
+             decSetCoeff (res, set, acc, res->digits, &residue, status);
+#if DECSUBSET
+           }
+         else
+           {                   /* subset arithmetic rounds from original significant digit */
+             /* We may have an underestimate.  This only occurs when both */
+             /* numbers fit in DECDPUN digits and we are padding with a */
+             /* negative multiple (-10, -100...) and the top digit(s) become */
+             /* 0.  (This only matters if we are using X3.274 rules where the */
+             /* leading zero could be included in the rounding.) */
+             if (res->digits < maxdigits)
+               {
+                 *(acc + D2U (res->digits)) = 0;       /* ensure leading 0 is there */
+                 res->digits = maxdigits;
+               }
+             else
+               {
+                 /* remove leading zeros that we added due to rounding up to */
+                 /* integral Units (but only those in excess of the original */
+                 /* maxdigits length, unless extended) before test for rounding. */
+                 if (res->digits > reqdigits)
+                   {
+                     res->digits = decGetDigits (acc, D2U (res->digits));
+                     if (res->digits < maxdigits)
+                       res->digits = maxdigits;
+                   }
+               }
+             decSetCoeff (res, set, acc, res->digits, &residue, status);
+             /* Now apply rounding if needed before removing leading zeros. */
+             /* This is safe because subnormals are not a possibility */
+             if (residue != 0)
+               {
+                 decApplyRound (res, set, residue, status);
+                 residue = 0;  /* we did what we had to do */
+               }
+           }                   /* subset */
+#endif
+       }                       /* used buffer */
+
+      /* strip leading zeros [these were left on in case of subset subtract] */
+      res->digits = decGetDigits (res->lsu, D2U (res->digits));
+
+      /* apply checks and rounding */
+      decFinish (res, set, &residue, status);
+
+      /* "When the sum of two operands with opposite signs is exactly */
+      /* zero, the sign of that sum shall be '+' in all rounding modes */
+      /* except round toward -Infinity, in which mode that sign shall be */
+      /* '-'."  [Subset zeros also never have '-', set by decFinish.] */
+      if (ISZERO (res) && diffsign
+#if DECSUBSET
+         && set->extended
+#endif
+         && (*status & DEC_Inexact) == 0)
+       {
+         if (set->round == DEC_ROUND_FLOOR)
+           res->bits |= DECNEG;        /* sign - */
+         else
+           res->bits &= ~DECNEG;       /* sign + */
+       }
+    }
+  while (0);                   /* end protected */
+
+  if (alloced)
+    {
+      if (allocacc != NULL)
+       free (allocacc);        /* drop any storage we used */
+      if (allocrhs != NULL)
+       free (allocrhs);        /* .. */
+      if (alloclhs != NULL)
+       free (alloclhs);        /* .. */
+    }
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decDivideOp -- division operation                                  */
+/*                                                                    */
+/*  This routine performs the calculations for all four division      */
+/*  operators (divide, divideInteger, remainder, remainderNear).      */
+/*                                                                    */
+/*  C=A op 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                                               */
+/*   op  is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively.    */
+/*   status is the usual accumulator                                  */
+/*                                                                    */
+/* C must have space for set->digits digits.                          */
+/*                                                                    */
+/* ------------------------------------------------------------------ */
+/*   The underlying algorithm of this routine is the same as in the   */
+/*   1981 S/370 implementation, that is, non-restoring long division  */
+/*   with bi-unit (rather than bi-digit) estimation for each unit     */
+/*   multiplier.  In this pseudocode overview, complications for the  */
+/*   Remainder operators and division residues for exact rounding are */
+/*   omitted for clarity.                                             */
+/*                                                                    */
+/*     Prepare operands and handle special values                     */
+/*     Test for x/0 and then 0/x                                      */
+/*     Exp =Exp1 - Exp2                                               */
+/*     Exp =Exp +len(var1) -len(var2)                                 */
+/*     Sign=Sign1 * Sign2                                             */
+/*     Pad accumulator (Var1) to double-length with 0's (pad1)        */
+/*     Pad Var2 to same length as Var1                                */
+/*     msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round  */
+/*     have=0                                                         */
+/*     Do until (have=digits+1 OR residue=0)                          */
+/*       if exp<0 then if integer divide/residue then leave           */
+/*       this_unit=0                                                  */
+/*       Do forever                                                   */
+/*          compare numbers                                           */
+/*          if <0 then leave inner_loop                               */
+/*          if =0 then (* quick exit without subtract *) do           */
+/*             this_unit=this_unit+1; output this_unit                */
+/*             leave outer_loop; end                                  */
+/*          Compare lengths of numbers (mantissae):                   */
+/*          If same then tops2=msu2pair -- {units 1&2 of var2}        */
+/*                  else tops2=msu2plus -- {0, unit 1 of var2}        */
+/*          tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
+/*          mult=tops1/tops2  -- Good and safe guess at divisor       */
+/*          if mult=0 then mult=1                                     */
+/*          this_unit=this_unit+mult                                  */
+/*          subtract                                                  */
+/*          end inner_loop                                            */
+/*        if have\=0 | this_unit\=0 then do                           */
+/*          output this_unit                                          */
+/*          have=have+1; end                                          */
+/*        var2=var2/10                                                */
+/*        exp=exp-1                                                   */
+/*        end outer_loop                                              */
+/*     exp=exp+1   -- set the proper exponent                         */
+/*     if have=0 then generate answer=0                               */
+/*     Return (Result is defined by Var1)                             */
+/*                                                                    */
+/* ------------------------------------------------------------------ */
+/* We need two working buffers during the long division; one (digits+ */
+/* 1) to accumulate the result, and the other (up to 2*digits+1) for  */
+/* long subtractions.  These are acc and var1 respectively.           */
+/* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
+/* ------------------------------------------------------------------ */
+static decNumber *
+decDivideOp (decNumber * res,
+            decNumber * lhs, decNumber * rhs,
+            decContext * set, Flag op, uInt * status)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  Unit accbuff[D2U (DECBUFFER + DECDPUN)];     /* local buffer */
+  Unit *acc = accbuff;         /* -> accumulator array for result */
+  Unit *allocacc = NULL;       /* -> allocated buffer, iff allocated */
+  Unit *accnext;               /* -> where next digit will go */
+  Int acclength;               /* length of acc needed [Units] */
+  Int accunits;                        /* count of units accumulated */
+  Int accdigits;               /* count of digits accumulated */
+
+  Unit varbuff[D2U (DECBUFFER * 2 + DECDPUN) * sizeof (Unit)]; /* buffer for var1 */
+  Unit *var1 = varbuff;                /* -> var1 array for long subtraction */
+  Unit *varalloc = NULL;       /* -> allocated buffer, iff used */
+
+  Unit *var2;                  /* -> var2 array */
+
+  Int var1units, var2units;    /* actual lengths */
+  Int var2ulen;                        /* logical length (units) */
+  Int var1initpad = 0;         /* var1 initial padding (digits) */
+  Unit *msu1, *msu2;           /* -> msu of each var */
+  Int msu2plus;                        /* msu2 plus one [does not vary] */
+  eInt msu2pair;               /* msu2 pair plus one [does not vary] */
+  Int maxdigits;               /* longest LHS or required acc length */
+  Int mult;                    /* multiplier for subtraction */
+  Unit thisunit;               /* current unit being accumulated */
+  Int residue;                 /* for rounding */
+  Int reqdigits = set->digits; /* requested DIGITS */
+  Int exponent;                        /* working exponent */
+  Int maxexponent = 0;         /* DIVIDE maximum exponent if unrounded */
+  uByte bits;                  /* working sign */
+  uByte merged;                        /* merged flags */
+  Unit *target, *source;       /* work */
+  uInt const *pow;             /* .. */
+  Int shift, cut;              /* .. */
+#if DECSUBSET
+  Int dropped;                 /* work */
+#endif
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits 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] */
+
+      bits = (lhs->bits ^ rhs->bits) & DECNEG; /* assumed sign for divisions */
+
+      /* handle infinities and NaNs */
+      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
+      if (merged)
+       {                       /* a special bit set */
+         if (merged & (DECSNAN | DECNAN))
+           {                   /* one or two NaNs */
+             decNaNs (res, lhs, rhs, status);
+             break;
+           }
+         /* one or two infinities */
+         if (decNumberIsInfinite (lhs))
+           {                   /* LHS (dividend) is infinite */
+             if (decNumberIsInfinite (rhs) ||  /* two infinities are invalid .. */
+                 op & (REMAINDER | REMNEAR))
+               {               /* as is remainder of infinity */
+                 *status |= DEC_Invalid_operation;
+                 break;
+               }
+             /* [Note that infinity/0 raises no exceptions] */
+             decNumberZero (res);
+             res->bits = bits | DECINF;        /* set +/- infinity */
+             break;
+           }
+         else
+           {                   /* RHS (divisor) is infinite */
+             residue = 0;
+             if (op & (REMAINDER | REMNEAR))
+               {
+                 /* result is [finished clone of] lhs */
+                 decCopyFit (res, lhs, set, &residue, status);
+               }
+             else
+               {               /* a division */
+                 decNumberZero (res);
+                 res->bits = bits;     /* set +/- zero */
+                 /* for DIVIDEINT the exponent is always 0.  For DIVIDE, result */
+                 /* is a 0 with infinitely negative exponent, clamped to minimum */
+                 if (op & DIVIDE)
+                   {
+                     res->exponent = set->emin - set->digits + 1;
+                     *status |= DEC_Clamped;
+                   }
+               }
+             decFinish (res, set, &residue, status);
+             break;
+           }
+       }
+
+      /* handle 0 rhs (x/0) */
+      if (ISZERO (rhs))
+       {                       /* x/0 is always exceptional */
+         if (ISZERO (lhs))
+           {
+             decNumberZero (res);      /* [after lhs test] */
+             *status |= DEC_Division_undefined;        /* 0/0 will become NaN */
+           }
+         else
+           {
+             decNumberZero (res);
+             if (op & (REMAINDER | REMNEAR))
+               *status |= DEC_Invalid_operation;
+             else
+               {
+                 *status |= DEC_Division_by_zero;      /* x/0 */
+                 res->bits = bits | DECINF;    /* .. is +/- Infinity */
+               }
+           }
+         break;
+       }
+
+      /* handle 0 lhs (0/x) */
+      if (ISZERO (lhs))
+       {                       /* 0/x [x!=0] */
+#if DECSUBSET
+         if (!set->extended)
+           decNumberZero (res);
+         else
+           {
+#endif
+             if (op & DIVIDE)
+               {
+                 residue = 0;
+                 exponent = lhs->exponent - rhs->exponent;     /* ideal exponent */
+                 decNumberCopy (res, lhs);     /* [zeros always fit] */
+                 res->bits = bits;     /* sign as computed */
+                 res->exponent = exponent;     /* exponent, too */
+                 decFinalize (res, set, &residue, status);     /* check exponent */
+               }
+             else if (op & DIVIDEINT)
+               {
+                 decNumberZero (res);  /* integer 0 */
+                 res->bits = bits;     /* sign as computed */
+               }
+             else
+               {               /* a remainder */
+                 exponent = rhs->exponent;     /* [save in case overwrite] */
+                 decNumberCopy (res, lhs);     /* [zeros always fit] */
+                 if (exponent < res->exponent)
+                   res->exponent = exponent;   /* use lower */
+               }
+#if DECSUBSET
+           }
+#endif
+         break;
+       }
+
+      /* Precalculate exponent.  This starts off adjusted (and hence fits */
+      /* in 31 bits) and becomes the usual unadjusted exponent as the */
+      /* division proceeds.  The order of evaluation is important, here, */
+      /* to avoid wrap. */
+      exponent =
+       (lhs->exponent + lhs->digits) - (rhs->exponent + rhs->digits);
+
+      /* If the working exponent is -ve, then some quick exits are */
+      /* possible because the quotient is known to be <1 */
+      /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
+      if (exponent < 0 && !(op == DIVIDE))
+       {
+         if (op & DIVIDEINT)
+           {
+             decNumberZero (res);      /* integer part is 0 */
+#if DECSUBSET
+             if (set->extended)
+#endif
+               res->bits = bits;       /* set +/- zero */
+             break;
+           }
+         /* we can fastpath remainders so long as the lhs has the */
+         /* smaller (or equal) exponent */
+         if (lhs->exponent <= rhs->exponent)
+           {
+             if (op & REMAINDER || exponent < -1)
+               {
+                 /* It is REMAINDER or safe REMNEAR; result is [finished */
+                 /* clone of] lhs  (r = x - 0*y) */
+                 residue = 0;
+                 decCopyFit (res, lhs, set, &residue, status);
+                 decFinish (res, set, &residue, status);
+                 break;
+               }
+             /* [unsafe REMNEAR drops through] */
+           }
+       }                       /* fastpaths */
+
+      /* We need long (slow) division; roll up the sleeves... */
+
+      /* The accumulator will hold the quotient of the division. */
+      /* If it needs to be too long for stack storage, then allocate. */
+      acclength = D2U (reqdigits + DECDPUN);   /* in Units */
+      if (acclength * sizeof (Unit) > sizeof (accbuff))
+       {
+         allocacc = (Unit *) malloc (acclength * sizeof (Unit));
+         if (allocacc == NULL)
+           {                   /* hopeless -- abandon */
+             *status |= DEC_Insufficient_storage;
+             break;
+           }
+         acc = allocacc;       /* use the allocated space */
+       }
+
+      /* var1 is the padded LHS ready for subtractions. */
+      /* If it needs to be too long for stack storage, then allocate. */
+      /* The maximum units we need for var1 (long subtraction) is: */
+      /* Enough for */
+      /*     (rhs->digits+reqdigits-1) -- to allow full slide to right */
+      /* or  (lhs->digits)             -- to allow for long lhs */
+      /* whichever is larger */
+      /*   +1                -- for rounding of slide to right */
+      /*   +1                -- for leading 0s */
+      /*   +1                -- for pre-adjust if a remainder or DIVIDEINT */
+      /* [Note: unused units do not participate in decUnitAddSub data] */
+      maxdigits = rhs->digits + reqdigits - 1;
+      if (lhs->digits > maxdigits)
+       maxdigits = lhs->digits;
+      var1units = D2U (maxdigits) + 2;
+      /* allocate a guard unit above msu1 for REMAINDERNEAR */
+      if (!(op & DIVIDE))
+       var1units++;
+      if ((var1units + 1) * sizeof (Unit) > sizeof (varbuff))
+       {
+         varalloc = (Unit *) malloc ((var1units + 1) * sizeof (Unit));
+         if (varalloc == NULL)
+           {                   /* hopeless -- abandon */
+             *status |= DEC_Insufficient_storage;
+             break;
+           }
+         var1 = varalloc;      /* use the allocated space */
+       }
+
+      /* Extend the lhs and rhs to full long subtraction length.  The lhs */
+      /* is truly extended into the var1 buffer, with 0 padding, so we can */
+      /* subtract in place.  The rhs (var2) has virtual padding */
+      /* (implemented by decUnitAddSub). */
+      /* We allocated one guard unit above msu1 for rem=rem+rem in REMAINDERNEAR */
+      msu1 = var1 + var1units - 1;     /* msu of var1 */
+      source = lhs->lsu + D2U (lhs->digits) - 1;       /* msu of input array */
+      for (target = msu1; source >= lhs->lsu; source--, target--)
+       *target = *source;
+      for (; target >= var1; target--)
+       *target = 0;
+
+      /* rhs (var2) is left-aligned with var1 at the start */
+      var2ulen = var1units;    /* rhs logical length (units) */
+      var2units = D2U (rhs->digits);   /* rhs actual length (units) */
+      var2 = rhs->lsu;         /* -> rhs array */
+      msu2 = var2 + var2units - 1;     /* -> msu of var2 [never changes] */
+      /* now set up the variables which we'll use for estimating the */
+      /* multiplication factor.  If these variables are not exact, we add */
+      /* 1 to make sure that we never overestimate the multiplier. */
+      msu2plus = *msu2;                /* it's value .. */
+      if (var2units > 1)
+       msu2plus++;             /* .. +1 if any more */
+      msu2pair = (eInt) * msu2 * (DECDPUNMAX + 1);     /* top two pair .. */
+      if (var2units > 1)
+       {                       /* .. [else treat 2nd as 0] */
+         msu2pair += *(msu2 - 1);      /* .. */
+         if (var2units > 2)
+           msu2pair++;         /* .. +1 if any more */
+       }
+
+      /* Since we are working in units, the units may have leading zeros, */
+      /* but we calculated the exponent on the assumption that they are */
+      /* both left-aligned.  Adjust the exponent to compensate: add the */
+      /* number of leading zeros in var1 msu and subtract those in var2 msu. */
+      /* [We actually do this by counting the digits and negating, as */
+      /* lead1=DECDPUN-digits1, and similarly for lead2.] */
+      for (pow = &powers[1]; *msu1 >= *pow; pow++)
+       exponent--;
+      for (pow = &powers[1]; *msu2 >= *pow; pow++)
+       exponent++;
+
+      /* Now, if doing an integer divide or remainder, we want to ensure */
+      /* that the result will be Unit-aligned.  To do this, we shift the */
+      /* var1 accumulator towards least if need be.  (It's much easier to */
+      /* do this now than to reassemble the residue afterwards, if we are */
+      /* doing a remainder.)  Also ensure the exponent is not negative. */
+      if (!(op & DIVIDE))
+       {
+         Unit *u;
+         /* save the initial 'false' padding of var1, in digits */
+         var1initpad = (var1units - D2U (lhs->digits)) * DECDPUN;
+         /* Determine the shift to do. */
+         if (exponent < 0)
+           cut = -exponent;
+         else
+           cut = DECDPUN - exponent % DECDPUN;
+         decShiftToLeast (var1, var1units, cut);
+         exponent += cut;      /* maintain numerical value */
+         var1initpad -= cut;   /* .. and reduce padding */
+         /* clean any most-significant units we just emptied */
+         for (u = msu1; cut >= DECDPUN; cut -= DECDPUN, u--)
+           *u = 0;
+       }                       /* align */
+      else
+       {                       /* is DIVIDE */
+         maxexponent = lhs->exponent - rhs->exponent;  /* save */
+         /* optimization: if the first iteration will just produce 0, */
+         /* preadjust to skip it [valid for DIVIDE only] */
+         if (*msu1 < *msu2)
+           {
+             var2ulen--;       /* shift down */
+             exponent -= DECDPUN;      /* update the exponent */
+           }
+       }
+
+      /* ---- start the long-division loops ------------------------------ */
+      accunits = 0;            /* no units accumulated yet */
+      accdigits = 0;           /* .. or digits */
+      accnext = acc + acclength - 1;   /* -> msu of acc [NB: allows digits+1] */
+      for (;;)
+       {                       /* outer forever loop */
+         thisunit = 0;         /* current unit assumed 0 */
+         /* find the next unit */
+         for (;;)
+           {                   /* inner forever loop */
+             /* strip leading zero units [from either pre-adjust or from */
+             /* subtract last time around].  Leave at least one unit. */
+             for (; *msu1 == 0 && msu1 > var1; msu1--)
+               var1units--;
+
+             if (var1units < var2ulen)
+               break;          /* var1 too low for subtract */
+             if (var1units == var2ulen)
+               {               /* unit-by-unit compare needed */
+                 /* compare the two numbers, from msu */
+                 Unit *pv1, *pv2, v2;  /* units to compare */
+                 pv2 = msu2;   /* -> msu */
+                 for (pv1 = msu1;; pv1--, pv2--)
+                   {
+                     /* v1=*pv1 -- always OK */
+                     v2 = 0;   /* assume in padding */
+                     if (pv2 >= var2)
+                       v2 = *pv2;      /* in range */
+                     if (*pv1 != v2)
+                       break;  /* no longer the same */
+                     if (pv1 == var1)
+                       break;  /* done; leave pv1 as is */
+                   }
+                 /* here when all inspected or a difference seen */
+                 if (*pv1 < v2)
+                   break;      /* var1 too low to subtract */
+                 if (*pv1 == v2)
+                   {           /* var1 == var2 */
+                     /* reach here if var1 and var2 are identical; subtraction */
+                     /* would increase digit by one, and the residue will be 0 so */
+                     /* we are done; leave the loop with residue set to 0. */
+                     thisunit++;       /* as though subtracted */
+                     *var1 = 0;        /* set var1 to 0 */
+                     var1units = 1;    /* .. */
+                     break;    /* from inner */
+                   }           /* var1 == var2 */
+                 /* *pv1>v2.  Prepare for real subtraction; the lengths are equal */
+                 /* Estimate the multiplier (there's always a msu1-1)... */
+                 /* Bring in two units of var2 to provide a good estimate. */
+                 mult =
+                   (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
+                           *(msu1 - 1)) / msu2pair);
+               }               /* lengths the same */
+             else
+               {               /* var1units > var2ulen, so subtraction is safe */
+                 /* The var2 msu is one unit towards the lsu of the var1 msu, */
+                 /* so we can only use one unit for var2. */
+                 mult =
+                   (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
+                           *(msu1 - 1)) / msu2plus);
+               }
+             if (mult == 0)
+               mult = 1;       /* must always be at least 1 */
+             /* subtraction needed; var1 is > var2 */
+             thisunit = (Unit) (thisunit + mult);      /* accumulate */
+             /* subtract var1-var2, into var1; only the overlap needs */
+             /* processing, as we are in place */
+             shift = var2ulen - var2units;
+#if DECTRACE
+             decDumpAr ('1', &var1[shift], var1units - shift);
+             decDumpAr ('2', var2, var2units);
+             printf ("m=%d\n", -mult);
+#endif
+             decUnitAddSub (&var1[shift], var1units - shift,
+                            var2, var2units, 0, &var1[shift], -mult);
+#if DECTRACE
+             decDumpAr ('#', &var1[shift], var1units - shift);
+#endif
+             /* var1 now probably has leading zeros; these are removed at the */
+             /* top of the inner loop. */
+           }                   /* inner loop */
+
+         /* We have the next unit; unless it's a leading zero, add to acc */
+         if (accunits != 0 || thisunit != 0)
+           {                   /* put the unit we got */
+             *accnext = thisunit;      /* store in accumulator */
+             /* account exactly for the digits we got */
+             if (accunits == 0)
+               {
+                 accdigits++;  /* at least one */
+                 for (pow = &powers[1]; thisunit >= *pow; pow++)
+                   accdigits++;
+               }
+             else
+               accdigits += DECDPUN;
+             accunits++;       /* update count */
+             accnext--;        /* ready for next */
+             if (accdigits > reqdigits)
+               break;          /* we have all we need */
+           }
+
+         /* if the residue is zero, we're done (unless divide or */
+         /* divideInteger and we haven't got enough digits yet) */
+         if (*var1 == 0 && var1units == 1)
+           {                   /* residue is 0 */
+             if (op & (REMAINDER | REMNEAR))
+               break;
+             if ((op & DIVIDE) && (exponent <= maxexponent))
+               break;
+             /* [drop through if divideInteger] */
+           }
+         /* we've also done enough if calculating remainder or integer */
+         /* divide and we just did the last ('units') unit */
+         if (exponent == 0 && !(op & DIVIDE))
+           break;
+
+         /* to get here, var1 is less than var2, so divide var2 by the per- */
+         /* Unit power of ten and go for the next digit */
+         var2ulen--;           /* shift down */
+         exponent -= DECDPUN;  /* update the exponent */
+       }                       /* outer loop */
+
+      /* ---- division is complete --------------------------------------- */
+      /* here: acc      has at least reqdigits+1 of good results (or fewer */
+      /*                if early stop), starting at accnext+1 (its lsu) */
+      /*       var1     has any residue at the stopping point */
+      /*       accunits is the number of digits we collected in acc */
+      if (accunits == 0)
+       {                       /* acc is 0 */
+         accunits = 1;         /* show we have one .. */
+         accdigits = 1;        /* .. */
+         *accnext = 0;         /* .. whose value is 0 */
+       }
+      else
+       accnext++;              /* back to last placed */
+      /* accnext now -> lowest unit of result */
+
+      residue = 0;             /* assume no residue */
+      if (op & DIVIDE)
+       {
+         /* record the presence of any residue, for rounding */
+         if (*var1 != 0 || var1units > 1)
+           residue = 1;
+         else
+           {                   /* no residue */
+             /* We had an exact division; clean up spurious trailing 0s. */
+             /* There will be at most DECDPUN-1, from the final multiply, */
+             /* and then only if the result is non-0 (and even) and the */
+             /* exponent is 'loose'. */
+#if DECDPUN>1
+             Unit lsu = *accnext;
+             if (!(lsu & 0x01) && (lsu != 0))
+               {
+                 /* count the trailing zeros */
+                 Int drop = 0;
+                 for (;; drop++)
+                   {           /* [will terminate because lsu!=0] */
+                     if (exponent >= maxexponent)
+                       break;  /* don't chop real 0s */
+#if DECDPUN<=4
+                     if ((lsu - QUOT10 (lsu, drop + 1)
+                          * powers[drop + 1]) != 0)
+                       break;  /* found non-0 digit */
+#else
+                     if (lsu % powers[drop + 1] != 0)
+                       break;  /* found non-0 digit */
+#endif
+                     exponent++;
+                   }
+                 if (drop > 0)
+                   {
+                     accunits = decShiftToLeast (accnext, accunits, drop);
+                     accdigits = decGetDigits (accnext, accunits);
+                     accunits = D2U (accdigits);
+                     /* [exponent was adjusted in the loop] */
+                   }
+               }               /* neither odd nor 0 */
+#endif
+           }                   /* exact divide */
+       }                       /* divide */
+      else                     /* op!=DIVIDE */
+       {
+         /* check for coefficient overflow */
+         if (accdigits + exponent > reqdigits)
+           {
+             *status |= DEC_Division_impossible;
+             break;
+           }
+         if (op & (REMAINDER | REMNEAR))
+           {
+             /* [Here, the exponent will be 0, because we adjusted var1 */
+             /* appropriately.] */
+             Int postshift;    /* work */
+             Flag wasodd = 0;  /* integer was odd */
+             Unit *quotlsu;    /* for save */
+             Int quotdigits;   /* .. */
+
+             /* Fastpath when residue is truly 0 is worthwhile [and */
+             /* simplifies the code below] */
+             if (*var1 == 0 && var1units == 1)
+               {               /* residue is 0 */
+                 Int exp = lhs->exponent;      /* save min(exponents) */
+                 if (rhs->exponent < exp)
+                   exp = rhs->exponent;
+                 decNumberZero (res);  /* 0 coefficient */
+#if DECSUBSET
+                 if (set->extended)
+#endif
+                   res->exponent = exp;        /* .. with proper exponent */
+                 break;
+               }
+             /* note if the quotient was odd */
+             if (*accnext & 0x01)
+               wasodd = 1;     /* acc is odd */
+             quotlsu = accnext;        /* save in case need to reinspect */
+             quotdigits = accdigits;   /* .. */
+
+             /* treat the residue, in var1, as the value to return, via acc */
+             /* calculate the unused zero digits.  This is the smaller of: */
+             /*   var1 initial padding (saved above) */
+             /*   var2 residual padding, which happens to be given by: */
+             postshift =
+               var1initpad + exponent - lhs->exponent + rhs->exponent;
+             /* [the 'exponent' term accounts for the shifts during divide] */
+             if (var1initpad < postshift)
+               postshift = var1initpad;
+
+             /* shift var1 the requested amount, and adjust its digits */
+             var1units = decShiftToLeast (var1, var1units, postshift);
+             accnext = var1;
+             accdigits = decGetDigits (var1, var1units);
+             accunits = D2U (accdigits);
+
+             exponent = lhs->exponent; /* exponent is smaller of lhs & rhs */
+             if (rhs->exponent < exponent)
+               exponent = rhs->exponent;
+             bits = lhs->bits; /* remainder sign is always as lhs */
+
+             /* Now correct the result if we are doing remainderNear; if it */
+             /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
+             /* the integer was odd then the result should be rem-rhs. */
+             if (op & REMNEAR)
+               {
+                 Int compare, tarunits;        /* work */
+                 Unit *up;     /* .. */
+
+
+                 /* calculate remainder*2 into the var1 buffer (which has */
+                 /* 'headroom' of an extra unit and hence enough space) */
+                 /* [a dedicated 'double' loop would be faster, here] */
+                 tarunits =
+                   decUnitAddSub (accnext, accunits, accnext, accunits, 0,
+                                  accnext, 1);
+                 /* decDumpAr('r', accnext, tarunits); */
+
+                 /* Here, accnext (var1) holds tarunits Units with twice the */
+                 /* remainder's coefficient, which we must now compare to the */
+                 /* RHS.  The remainder's exponent may be smaller than the RHS's. */
+                 compare =
+                   decUnitCompare (accnext, tarunits, rhs->lsu,
+                                   D2U (rhs->digits),
+                                   rhs->exponent - exponent);
+                 if (compare == BADINT)
+                   {           /* deep trouble */
+                     *status |= DEC_Insufficient_storage;
+                     break;
+                   }
+
+                 /* now restore the remainder by dividing by two; we know the */
+                 /* lsu is even. */
+                 for (up = accnext; up < accnext + tarunits; up++)
+                   {
+                     Int half; /* half to add to lower unit */
+                     half = *up & 0x01;
+                     *up /= 2; /* [shift] */
+                     if (!half)
+                       continue;
+                     *(up - 1) += (DECDPUNMAX + 1) / 2;
+                   }
+                 /* [accunits still describes the original remainder length] */
+
+                 if (compare > 0 || (compare == 0 && wasodd))
+                   {           /* adjustment needed */
+                     Int exp, expunits, exprem;        /* work */
+                     /* This is effectively causing round-up of the quotient, */
+                     /* so if it was the rare case where it was full and all */
+                     /* nines, it would overflow and hence division-impossible */
+                     /* should be raised */
+                     Flag allnines = 0;        /* 1 if quotient all nines */
+                     if (quotdigits == reqdigits)
+                       {       /* could be borderline */
+                         for (up = quotlsu;; up++)
+                           {
+                             if (quotdigits > DECDPUN)
+                               {
+                                 if (*up != DECDPUNMAX)
+                                   break;      /* non-nines */
+                               }
+                             else
+                               {       /* this is the last Unit */
+                                 if (*up == powers[quotdigits] - 1)
+                                   allnines = 1;
+                                 break;
+                               }
+                             quotdigits -= DECDPUN;    /* checked those digits */
+                           }   /* up */
+                       }       /* borderline check */
+                     if (allnines)
+                       {
+                         *status |= DEC_Division_impossible;
+                         break;
+                       }
+
+                     /* we need rem-rhs; the sign will invert.  Again we can */
+                     /* safely use var1 for the working Units array. */
+                     exp = rhs->exponent - exponent;   /* RHS padding needed */
+                     /* Calculate units and remainder from exponent. */
+                     expunits = exp / DECDPUN;
+                     exprem = exp % DECDPUN;
+                     /* subtract [A+B*(-m)]; the result will always be negative */
+                     accunits = -decUnitAddSub (accnext, accunits,
+                                                rhs->lsu, D2U (rhs->digits),
+                                                expunits, accnext,
+                                                -(Int) powers[exprem]);
+                     accdigits = decGetDigits (accnext, accunits);     /* count digits exactly */
+                     accunits = D2U (accdigits);       /* and recalculate the units for copy */
+                     /* [exponent is as for original remainder] */
+                     bits ^= DECNEG;   /* flip the sign */
+                   }
+               }               /* REMNEAR */
+           }                   /* REMAINDER or REMNEAR */
+       }                       /* not DIVIDE */
+
+      /* Set exponent and bits */
+      res->exponent = exponent;
+      res->bits = (uByte) (bits & DECNEG);     /* [cleaned] */
+
+      /* Now the coefficient. */
+      decSetCoeff (res, set, accnext, accdigits, &residue, status);
+
+      decFinish (res, set, &residue, status);  /* final cleanup */
+
+#if DECSUBSET
+      /* If a divide then strip trailing zeros if subset [after round] */
+      if (!set->extended && (op == DIVIDE))
+       decTrim (res, 0, &dropped);
+#endif
+    }
+  while (0);                   /* end protected */
+
+  if (varalloc != NULL)
+    free (varalloc);           /* drop any storage we used */
+  if (allocacc != NULL)
+    free (allocacc);           /* .. */
+  if (allocrhs != NULL)
+    free (allocrhs);           /* .. */
+  if (alloclhs != NULL)
+    free (alloclhs);           /* .. */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decMultiplyOp -- multiplication operation                          */
+/*                                                                    */
+/*  This routine performs the multiplication 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                                               */
+/*   status is the usual accumulator                                  */
+/*                                                                    */
+/* C must have space for set->digits digits.                          */
+/*                                                                    */
+/* ------------------------------------------------------------------ */
+/* Note: We use 'long' multiplication rather than Karatsuba, as the   */
+/* latter would give only a minor improvement for the short numbers   */
+/* we expect to handle most (and uses much more memory).              */
+/*                                                                    */
+/* We always have to use a buffer for the accumulator.                */
+/* ------------------------------------------------------------------ */
+static decNumber *
+decMultiplyOp (decNumber * res, decNumber * lhs,
+              decNumber * rhs, decContext * set, uInt * status)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  Unit accbuff[D2U (DECBUFFER * 2 + 1)];       /* local buffer (+1 in case DECBUFFER==0) */
+  Unit *acc = accbuff;         /* -> accumulator array for exact result */
+  Unit *allocacc = NULL;       /* -> allocated buffer, iff allocated */
+  Unit *mer, *mermsup;         /* work */
+  Int accunits;                        /* Units of accumulator in use */
+  Int madlength;               /* Units in multiplicand */
+  Int shift;                   /* Units to shift multiplicand by */
+  Int need;                    /* Accumulator units needed */
+  Int exponent;                        /* work */
+  Int residue = 0;             /* rounding residue */
+  uByte bits;                  /* result sign */
+  uByte merged;                        /* merged flags */
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits status, as needed */
+         if (lhs->digits > set->digits)
+           {
+             alloclhs = decRoundOperand (lhs, set, status);
+             if (alloclhs == NULL)
+               break;
+             lhs = alloclhs;
+           }
+         if (rhs->digits > set->digits)
+           {
+             allocrhs = decRoundOperand (rhs, set, status);
+             if (allocrhs == NULL)
+               break;
+             rhs = allocrhs;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* precalculate result sign */
+      bits = (uByte) ((lhs->bits ^ rhs->bits) & DECNEG);
+
+      /* handle infinities and NaNs */
+      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
+      if (merged)
+       {                       /* a special bit set */
+         if (merged & (DECSNAN | DECNAN))
+           {                   /* one or two NaNs */
+             decNaNs (res, lhs, rhs, status);
+             break;
+           }
+         /* one or two infinities. Infinity * 0 is invalid */
+         if (((lhs->bits & DECSPECIAL) == 0 && ISZERO (lhs))
+             || ((rhs->bits & DECSPECIAL) == 0 && ISZERO (rhs)))
+           {
+             *status |= DEC_Invalid_operation;
+             break;
+           }
+         decNumberZero (res);
+         res->bits = bits | DECINF;    /* infinity */
+         break;
+       }
+
+      /* For best speed, as in DMSRCN, we use the shorter number as the */
+      /* multiplier (rhs) and the longer as the multiplicand (lhs) */
+      if (lhs->digits < rhs->digits)
+       {                       /* swap... */
+         decNumber *hold = lhs;
+         lhs = rhs;
+         rhs = hold;
+       }
+
+      /* if accumulator is too long for local storage, then allocate */
+      need = D2U (lhs->digits) + D2U (rhs->digits);    /* maximum units in result */
+      if (need * sizeof (Unit) > sizeof (accbuff))
+       {
+         allocacc = (Unit *) malloc (need * sizeof (Unit));
+         if (allocacc == NULL)
+           {
+             *status |= DEC_Insufficient_storage;
+             break;
+           }
+         acc = allocacc;       /* use the allocated space */
+       }
+
+      /* Now the main long multiplication loop */
+      /* Unlike the equivalent in the IBM Java implementation, there */
+      /* is no advantage in calculating from msu to lsu.  So we do it */
+      /* by the book, as it were. */
+      /* Each iteration calculates ACC=ACC+MULTAND*MULT */
+      accunits = 1;            /* accumulator starts at '0' */
+      *acc = 0;                        /* .. (lsu=0) */
+      shift = 0;               /* no multiplicand shift at first */
+      madlength = D2U (lhs->digits);   /* we know this won't change */
+      mermsup = rhs->lsu + D2U (rhs->digits);  /* -> msu+1 of multiplier */
+
+      for (mer = rhs->lsu; mer < mermsup; mer++)
+       {
+         /* Here, *mer is the next Unit in the multiplier to use */
+         /* If non-zero [optimization] add it... */
+         if (*mer != 0)
+           {
+             accunits =
+               decUnitAddSub (&acc[shift], accunits - shift, lhs->lsu,
+                              madlength, 0, &acc[shift], *mer) + shift;
+           }
+         else
+           {                   /* extend acc with a 0; we'll use it shortly */
+             /* [this avoids length of <=0 later] */
+             *(acc + accunits) = 0;
+             accunits++;
+           }
+         /* multiply multiplicand by 10**DECDPUN for next Unit to left */
+         shift++;              /* add this for 'logical length' */
+       }                       /* n */
+#if DECTRACE
+      /* Show exact result */
+      decDumpAr ('*', acc, accunits);
+#endif
+
+      /* acc now contains the exact result of the multiplication */
+      /* Build a decNumber from it, noting if any residue */
+      res->bits = bits;                /* set sign */
+      res->digits = decGetDigits (acc, accunits);      /* count digits exactly */
+
+      /* We might have a 31-bit wrap in calculating the exponent. */
+      /* This can only happen if both input exponents are negative and */
+      /* both their magnitudes are large.  If we did wrap, we set a safe */
+      /* very negative exponent, from which decFinalize() will raise a */
+      /* hard underflow. */
+      exponent = lhs->exponent + rhs->exponent;        /* calculate exponent */
+      if (lhs->exponent < 0 && rhs->exponent < 0 && exponent > 0)
+       exponent = -2 * DECNUMMAXE;     /* force underflow */
+      res->exponent = exponent;        /* OK to overwrite now */
+
+      /* Set the coefficient.  If any rounding, residue records */
+      decSetCoeff (res, set, acc, res->digits, &residue, status);
+
+      decFinish (res, set, &residue, status);  /* final cleanup */
+    }
+  while (0);                   /* end protected */
+
+  if (allocacc != NULL)
+    free (allocacc);           /* drop any storage we used */
+  if (allocrhs != NULL)
+    free (allocrhs);           /* .. */
+  if (alloclhs != NULL)
+    free (alloclhs);           /* .. */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decQuantizeOp  -- 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 or matches the 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 requested exponent                                 */
+/*   set is the context                                               */
+/*   quant is 1 for quantize or 0 for rescale                         */
+/*   status is the status accumulator (this can be called without     */
+/*          risk of control loss)                                     */
+/*                                                                    */
+/* 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 that requested.            */
+/* ------------------------------------------------------------------ */
+static decNumber *
+decQuantizeOp (decNumber * res, decNumber * lhs,
+              decNumber * rhs, decContext * set, Flag quant, uInt * status)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  decNumber *inrhs = rhs;      /* save original rhs */
+  Int reqdigits = set->digits; /* requested DIGITS */
+  Int reqexp;                  /* requested exponent [-scale] */
+  Int residue = 0;             /* rounding residue */
+  uByte merged;                        /* merged flags */
+  Int etiny = set->emin - (set->digits - 1);
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits status, as needed */
+         if (lhs->digits > reqdigits)
+           {
+             alloclhs = decRoundOperand (lhs, set, status);
+             if (alloclhs == NULL)
+               break;
+             lhs = alloclhs;
+           }
+         if (rhs->digits > reqdigits)
+           {                   /* [this only checks lostDigits] */
+             allocrhs = decRoundOperand (rhs, set, status);
+             if (allocrhs == NULL)
+               break;
+             rhs = allocrhs;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* Handle special values */
+      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
+      if ((lhs->bits | rhs->bits) & DECSPECIAL)
+       {
+         /* NaNs get usual processing */
+         if (merged & (DECSNAN | DECNAN))
+           decNaNs (res, lhs, rhs, status);
+         /* one infinity but not both is bad */
+         else if ((lhs->bits ^ rhs->bits) & DECINF)
+           *status |= DEC_Invalid_operation;
+         /* both infinity: return lhs */
+         else
+           decNumberCopy (res, lhs);   /* [nop if in place] */
+         break;
+       }
+
+      /* set requested exponent */
+      if (quant)
+       reqexp = inrhs->exponent;       /* quantize -- match exponents */
+      else
+       {                       /* rescale -- use value of rhs */
+         /* Original rhs must be an integer that fits and is in range */
+#if DECSUBSET
+         reqexp = decGetInt (inrhs, set);
+#else
+         reqexp = decGetInt (inrhs);
+#endif
+       }
+
+#if DECSUBSET
+      if (!set->extended)
+       etiny = set->emin;      /* no subnormals */
+#endif
+
+      if (reqexp == BADINT     /* bad (rescale only) or .. */
+         || (reqexp < etiny)   /* < lowest */
+         || (reqexp > set->emax))
+       {                       /* > Emax */
+         *status |= DEC_Invalid_operation;
+         break;
+       }
+
+      /* we've processed the RHS, so we can overwrite it now if necessary */
+      if (ISZERO (lhs))
+       {                       /* zero coefficient unchanged */
+         decNumberCopy (res, lhs);     /* [nop if in place] */
+         res->exponent = reqexp;       /* .. just set exponent */
+#if DECSUBSET
+         if (!set->extended)
+           res->bits = 0;      /* subset specification; no -0 */
+#endif
+       }
+      else
+       {                       /* non-zero lhs */
+         Int adjust = reqexp - lhs->exponent;  /* digit adjustment needed */
+         /* if adjusted coefficient will not fit, give up now */
+         if ((lhs->digits - adjust) > reqdigits)
+           {
+             *status |= DEC_Invalid_operation;
+             break;
+           }
+
+         if (adjust > 0)
+           {                   /* increasing exponent */
+             /* this will decrease the length of the coefficient by adjust */
+             /* digits, and must round as it does so */
+             decContext workset;       /* work */
+             workset = *set;   /* clone rounding, etc. */
+             workset.digits = lhs->digits - adjust;    /* set requested length */
+             /* [note that the latter can be <1, here] */
+             decCopyFit (res, lhs, &workset, &residue, status);        /* fit to result */
+             decApplyRound (res, &workset, residue, status);   /* .. and round */
+             residue = 0;      /* [used] */
+             /* If we rounded a 999s case, exponent will be off by one; */
+             /* adjust back if so. */
+             if (res->exponent > reqexp)
+               {
+                 res->digits = decShiftToMost (res->lsu, res->digits, 1);      /* shift */
+                 res->exponent--;      /* (re)adjust the exponent. */
+               }
+#if DECSUBSET
+             if (ISZERO (res) && !set->extended)
+               res->bits = 0;  /* subset; no -0 */
+#endif
+           }                   /* increase */
+         else                  /* adjust<=0 */
+           {                   /* decreasing or = exponent */
+             /* this will increase the length of the coefficient by -adjust */
+             /* digits, by adding trailing zeros. */
+             decNumberCopy (res, lhs); /* [it will fit] */
+             /* if padding needed (adjust<0), add it now... */
+             if (adjust < 0)
+               {
+                 res->digits =
+                   decShiftToMost (res->lsu, res->digits, -adjust);
+                 res->exponent += adjust;      /* adjust the exponent */
+               }
+           }                   /* decrease */
+       }                       /* non-zero */
+
+      /* Check for overflow [do not use Finalize in this case, as an */
+      /* overflow here is a "don't fit" situation] */
+      if (res->exponent > set->emax - res->digits + 1)
+       {                       /* too big */
+         *status |= DEC_Invalid_operation;
+         break;
+       }
+      else
+       {
+         decFinalize (res, set, &residue, status);     /* set subnormal flags */
+         *status &= ~DEC_Underflow;    /* suppress Underflow [754r] */
+       }
+    }
+  while (0);                   /* end protected */
+
+  if (allocrhs != NULL)
+    free (allocrhs);           /* drop any storage we used */
+  if (alloclhs != NULL)
+    free (alloclhs);           /* .. */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decCompareOp -- compare, min, or max two Numbers                   */
+/*                                                                    */
+/*   This computes C = A ? B and returns the signum (as a Number)     */
+/*   for COMPARE or the maximum or minimum (for COMPMAX and COMPMIN). */
+/*                                                                    */
+/*   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                                               */
+/*   op  is the operation flag                                        */
+/*   status is the usual accumulator                                  */
+/*                                                                    */
+/* C must have space for one digit for COMPARE or set->digits for     */
+/* COMPMAX and COMPMIN.                                               */
+/* ------------------------------------------------------------------ */
+/* The emphasis here is on speed for common cases, and avoiding       */
+/* coefficient comparison if possible.                                */
+/* ------------------------------------------------------------------ */
+decNumber *
+decCompareOp (decNumber * res, decNumber * lhs, decNumber * rhs,
+             decContext * set, Flag op, uInt * status)
+{
+  decNumber *alloclhs = NULL;  /* non-NULL if rounded lhs allocated */
+  decNumber *allocrhs = NULL;  /* .., rhs */
+  Int result = 0;              /* default result value */
+  uByte merged;                        /* merged flags */
+  uByte bits = 0;              /* non-0 for NaN */
+
+#if DECCHECK
+  if (decCheckOperands (res, lhs, rhs, set))
+    return res;
+#endif
+
+  do
+    {                          /* protect allocated storage */
+#if DECSUBSET
+      if (!set->extended)
+       {
+         /* reduce operands and set lostDigits status, as needed */
+         if (lhs->digits > set->digits)
+           {
+             alloclhs = decRoundOperand (lhs, set, status);
+             if (alloclhs == NULL)
+               {
+                 result = BADINT;
+                 break;
+               }
+             lhs = alloclhs;
+           }
+         if (rhs->digits > set->digits)
+           {
+             allocrhs = decRoundOperand (rhs, set, status);
+             if (allocrhs == NULL)
+               {
+                 result = BADINT;
+                 break;
+               }
+             rhs = allocrhs;
+           }
+       }
+#endif
+      /* [following code does not require input rounding] */
+
+      /* handle NaNs now; let infinities drop through */
+      /* +++ review sNaN handling with 754r, for now assumes sNaN */
+      /* (even just one) leads to NaN. */
+      merged = (lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
+      if (merged)
+       {                       /* a NaN bit set */
+         if (op == COMPARE);
+         else if (merged & DECSNAN);
+         else
+           {                   /* 754r rules for MIN and MAX ignore single NaN */
+             /* here if MIN or MAX, and one or two quiet NaNs */
+             if (lhs->bits & rhs->bits & DECNAN);
+             else
+               {               /* just one quiet NaN */
+                 /* force choice to be the non-NaN operand */
+                 op = COMPMAX;
+                 if (lhs->bits & DECNAN)
+                   result = -1;        /* pick rhs */
+                 else
+                   result = +1;        /* pick lhs */
+                 break;
+               }
+           }
+         op = COMPNAN;         /* use special path */
+         decNaNs (res, lhs, rhs, status);
+         break;
+       }
+
+      result = decCompare (lhs, rhs);  /* we have numbers */
+    }
+  while (0);                   /* end protected */
+
+  if (result == BADINT)
+    *status |= DEC_Insufficient_storage;       /* rare */
+  else
+    {
+      if (op == COMPARE)
+       {                       /* return signum */
+         decNumberZero (res);  /* [always a valid result] */
+         if (result == 0)
+           res->bits = bits;   /* (maybe qNaN) */
+         else
+           {
+             *res->lsu = 1;
+             if (result < 0)
+               res->bits = DECNEG;
+           }
+       }
+      else if (op == COMPNAN); /* special, drop through */
+      else
+       {                       /* MAX or MIN, non-NaN result */
+         Int residue = 0;      /* rounding accumulator */
+         /* choose the operand for the result */
+         decNumber *choice;
+         if (result == 0)
+           {                   /* operands are numerically equal */
+             /* choose according to sign then exponent (see 754r) */
+             uByte slhs = (lhs->bits & DECNEG);
+             uByte srhs = (rhs->bits & DECNEG);
+#if DECSUBSET
+             if (!set->extended)
+               {               /* subset: force left-hand */
+                 op = COMPMAX;
+                 result = +1;
+               }
+             else
+#endif
+             if (slhs != srhs)
+               {               /* signs differ */
+                 if (slhs)
+                   result = -1;        /* rhs is max */
+                 else
+                   result = +1;        /* lhs is max */
+               }
+             else if (slhs && srhs)
+               {               /* both negative */
+                 if (lhs->exponent < rhs->exponent)
+                   result = +1;
+                 else
+                   result = -1;
+                 /* [if equal, we use lhs, technically identical] */
+               }
+             else
+               {               /* both positive */
+                 if (lhs->exponent > rhs->exponent)
+                   result = +1;
+                 else
+                   result = -1;
+                 /* [ditto] */
+               }
+           }                   /* numerically equal */
+         /* here result will be non-0 */
+         if (op == COMPMIN)
+           result = -result;   /* reverse if looking for MIN */
+         choice = (result > 0 ? lhs : rhs);    /* choose */
+         /* copy chosen to result, rounding if need be */
+         decCopyFit (res, choice, set, &residue, status);
+         decFinish (res, set, &residue, status);
+       }
+    }
+  if (allocrhs != NULL)
+    free (allocrhs);           /* free any storage we used */
+  if (alloclhs != NULL)
+    free (alloclhs);           /* .. */
+  return res;
+}
+
+/* ------------------------------------------------------------------ */
+/* decCompare -- compare two decNumbers by numerical value            */
+/*                                                                    */
+/*  This routine compares A ? B without altering them.                */
+/*                                                                    */
+/*  Arg1 is A, a decNumber which is not a NaN                         */
+/*  Arg2 is B, a decNumber which is not a NaN                         */
+/*                                                                    */
+/*  returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure   */
+/*  (the only possible failure is an allocation error)                */
+/* ------------------------------------------------------------------ */
+/* This could be merged into decCompareOp */
+static Int
+decCompare (decNumber * lhs, decNumber * rhs)
+{
+  Int result;                  /* result value */
+  Int sigr;                    /* rhs signum */
+  Int compare;                 /* work */
+  result = 1;                  /* assume signum(lhs) */
+  if (ISZERO (lhs))
+    result = 0;
+  else if (decNumberIsNegative (lhs))
+    result = -1;
+  sigr = 1;                    /* compute signum(rhs) */
+  if (ISZERO (rhs))
+    sigr = 0;
+  else if (decNumberIsNegative (rhs))
+    sigr = -1;
+  if (result > sigr)
+    return +1;                 /* L > R, return 1 */
+  if (result < sigr)
+    return -1;                 /* R < L, return -1 */
+
+  /* signums are the same */
+  if (result == 0)
+    return 0;                  /* both 0 */
+  /* Both non-zero */
+  if ((lhs->bits | rhs->bits) & DECINF)
+    {                          /* one or more infinities */
+      if (lhs->bits == rhs->bits)
+       result = 0;             /* both the same */
+      else if (decNumberIsInfinite (rhs))
+       result = -result;
+      return result;
+    }
+
+  /* we must compare the coefficients, allowing for exponents */
+  if (lhs->exponent > rhs->exponent)
+    {                          /* LHS exponent larger */
+      /* swap sides, and sign */
+      decNumber *temp = lhs;
+      lhs = rhs;
+      rhs = temp;
+      result = -result;
+    }
+
+  compare = decUnitCompare (lhs->lsu, D2U (lhs->digits),
+                           rhs->lsu, D2U (rhs->digits),
+                           rhs->exponent - lhs->exponent);
+
+  if (compare != BADINT)
+    compare *= result;         /* comparison succeeded */
+  return compare;              /* what we got */
+}
+
+/* ------------------------------------------------------------------ */
+/* decUnitCompare -- compare two >=0 integers in Unit arrays          */
+/*                                                                    */
+/*  This routine compares A ? B*10**E where A and B are unit arrays   */
+/*  A is a plain integer                                              */
+/*  B has an exponent of E (which must be non-negative)               */
+/*                                                                    */
+/*  Arg1 is A first Unit (lsu)                                        */
+/*  Arg2 is A length in Units                                         */
+/*  Arg3 is B first Unit (lsu)                                        */
+/*  Arg4 is B length in Units                                         */
+/*  Arg5 is E                                                         */
+/*                                                                    */
+/*  returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure   */
+/*  (the only possible failure is an allocation error)                */
+/* ------------------------------------------------------------------ */
+static Int
+decUnitCompare (Unit * a, Int alength, Unit * b, Int blength, Int exp)
+{
+  Unit *acc;                   /* accumulator for result */
+  Unit accbuff[D2U (DECBUFFER + 1)];   /* local buffer */
+  Unit *allocacc = NULL;       /* -> allocated acc buffer, iff allocated */
+  Int accunits, need;          /* units in use or needed for acc */
+  Unit *l, *r, *u;             /* work */
+  Int expunits, exprem, result;        /* .. */
+
+  if (exp == 0)
+    {                          /* aligned; fastpath */
+      if (alength > blength)
+       return 1;
+      if (alength < blength)
+       return -1;
+      /* same number of units in both -- need unit-by-unit compare */
+      l = a + alength - 1;
+      r = b + alength - 1;
+      for (; l >= a; l--, r--)
+       {
+         if (*l > *r)
+           return 1;
+         if (*l < *r)
+           return -1;
+       }
+      return 0;                        /* all units match */
+    }                          /* aligned */
+
+  /* Unaligned.  If one is >1 unit longer than the other, padded */
+  /* approximately, then we can return easily */
+  if (alength > blength + (Int) D2U (exp))
+    return 1;
+  if (alength + 1 < blength + (Int) D2U (exp))
+    return -1;
+
+  /* We need to do a real subtract.  For this, we need a result buffer */
+  /* even though we only are interested in the sign.  Its length needs */
+  /* to be the larger of alength and padded blength, +2 */
+  need = blength + D2U (exp);  /* maximum real length of B */
+  if (need < alength)
+    need = alength;
+  need += 2;
+  acc = accbuff;               /* assume use local buffer */
+  if (need * sizeof (Unit) > sizeof (accbuff))
+    {
+      allocacc = (Unit *) malloc (need * sizeof (Unit));
+      if (allocacc == NULL)
+       return BADINT;          /* hopeless -- abandon */
+      acc = allocacc;
+    }
+  /* Calculate units and remainder from exponent. */
+  expunits = exp / DECDPUN;
+  exprem = exp % DECDPUN;
+  /* subtract [A+B*(-m)] */
+  accunits = decUnitAddSub (a, alength, b, blength, expunits, acc,
+                           -(Int) powers[exprem]);
+  /* [UnitAddSub result may have leading zeros, even on zero] */
+  if (accunits < 0)
+    result = -1;               /* negative result */
+  else
+    {                          /* non-negative result */
+      /* check units of the result before freeing any storage */
+      for (u = acc; u < acc + accunits - 1 && *u == 0;)
+       u++;
+      result = (*u == 0 ? 0 : +1);
+    }
+  /* clean up and return the result */
+  if (allocacc != NULL)
+    free (allocacc);           /* drop any storage we used */
+  return result;
+}
+
+/* ------------------------------------------------------------------ */
+/* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays   */
+/*                                                                    */
+/*  This routine performs the calculation:                            */
+/*                                                                    */
+/*  C=A+(B*M)                                                         */
+/*                                                                    */
+/*  Where M is in the range -DECDPUNMAX through +DECDPUNMAX.          */
+/*                                                                    */
+/*  A may be shorter or longer than B.                                */
+/*                                                                    */
+/*  Leading zeros are not removed after a calculation.  The result is */
+/*  either the same length as the longer of A and B (adding any       */
+/*  shift), or one Unit longer than that (if a Unit carry occurred).  */
+/*                                                                    */
+/*  A and B content are not altered unless C is also A or B.          */
+/*  C may be the same array as A or B, but only if no zero padding is */
+/*  requested (that is, C may be B only if bshift==0).                */
+/*  C is filled from the lsu; only those units necessary to complete  */
+/*  the calculation are referenced.                                   */
+/*                                                                    */
+/*  Arg1 is A first Unit (lsu)                                        */
+/*  Arg2 is A length in Units                                         */
+/*  Arg3 is B first Unit (lsu)                                        */
+/*  Arg4 is B length in Units                                         */
+/*  Arg5 is B shift in Units  (>=0; pads with 0 units if positive)    */
+/*  Arg6 is C first Unit (lsu)                                        */
+/*  Arg7 is M, the multiplier                                         */
+/*                                                                    */
+/*  returns the count of Units written to C, which will be non-zero   */
+/*  and negated if the result is negative.  That is, the sign of the  */
+/*  returned Int is the sign of the result (positive for zero) and    */
+/*  the absolute value of the Int is the count of Units.              */
+/*                                                                    */
+/*  It is the caller's responsibility to make sure that C size is     */
+/*  safe, allowing space if necessary for a one-Unit carry.           */
+/*                                                                    */
+/*  This routine is severely performance-critical; *any* change here  */
+/*  must be measured (timed) to assure no performance degradation.    */
+/*  In particular, trickery here tends to be counter-productive, as   */
+/*  increased complexity of code hurts register optimizations on      */
+/*  register-poor architectures.  Avoiding divisions is nearly        */
+/*  always a Good Idea, however.                                      */
+/*                                                                    */
+/* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark  */
+/* (IBM Warwick, UK) for some of the ideas used in this routine.      */
+/* ------------------------------------------------------------------ */
+static Int
+decUnitAddSub (Unit * a, Int alength,
+              Unit * b, Int blength, Int bshift, Unit * c, Int m)
+{
+  Unit *alsu = a;              /* A lsu [need to remember it] */
+  Unit *clsu = c;              /* C ditto */
+  Unit *minC;                  /* low water mark for C */
+  Unit *maxC;                  /* high water mark for C */
+  eInt carry = 0;              /* carry integer (could be Long) */
+  Int add;                     /* work */
+#if DECDPUN==4                 /* myriadal */
+  Int est;                     /* estimated quotient */
+#endif
+
+#if DECTRACE
+  if (alength < 1 || blength < 1)
+    printf ("decUnitAddSub: alen blen m %d %d [%d]\n", alength, blength, m);
+#endif
+
+  maxC = c + alength;          /* A is usually the longer */
+  minC = c + blength;          /* .. and B the shorter */
+  if (bshift != 0)
+    {                          /* B is shifted; low As copy across */
+      minC += bshift;
+      /* if in place [common], skip copy unless there's a gap [rare] */
+      if (a == c && bshift <= alength)
+       {
+         c += bshift;
+         a += bshift;
+       }
+      else
+       for (; c < clsu + bshift; a++, c++)
+         {                     /* copy needed */
+           if (a < alsu + alength)
+             *c = *a;
+           else
+             *c = 0;
+         }
+    }
+  if (minC > maxC)
+    {                          /* swap */
+      Unit *hold = minC;
+      minC = maxC;
+      maxC = hold;
+    }
+
+  /* For speed, we do the addition as two loops; the first where both A */
+  /* and B contribute, and the second (if necessary) where only one or */
+  /* other of the numbers contribute. */
+  /* Carry handling is the same (i.e., duplicated) in each case. */
+  for (; c < minC; c++)
+    {
+      carry += *a;
+      a++;
+      carry += ((eInt) * b) * m;       /* [special-casing m=1/-1 */
+      b++;                     /* here is not a win] */
+      /* here carry is new Unit of digits; it could be +ve or -ve */
+      if ((ueInt) carry <= DECDPUNMAX)
+       {                       /* fastpath 0-DECDPUNMAX */
+         *c = (Unit) carry;
+         carry = 0;
+         continue;
+       }
+      /* remainder operator is undefined if negative, so we must test */
+#if DECDPUN==4                 /* use divide-by-multiply */
+      if (carry >= 0)
+       {
+         est = (((ueInt) carry >> 11) * 53687) >> 18;
+         *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
+         carry = est;          /* likely quotient [89%] */
+         if (*c < DECDPUNMAX + 1)
+           continue;           /* estimate was correct */
+         carry++;
+         *c -= DECDPUNMAX + 1;
+         continue;
+       }
+      /* negative case */
+      carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);      /* make positive */
+      est = (((ueInt) carry >> 11) * 53687) >> 18;
+      *c = (Unit) (carry - est * (DECDPUNMAX + 1));
+      carry = est - (DECDPUNMAX + 1);  /* correctly negative */
+      if (*c < DECDPUNMAX + 1)
+       continue;               /* was OK */
+      carry++;
+      *c -= DECDPUNMAX + 1;
+#else
+      if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
+       {                       /* fastpath carry +1 */
+         *c = (Unit) (carry - (DECDPUNMAX + 1));       /* [helps additions] */
+         carry = 1;
+         continue;
+       }
+      if (carry >= 0)
+       {
+         *c = (Unit) (carry % (DECDPUNMAX + 1));
+         carry = carry / (DECDPUNMAX + 1);
+         continue;
+       }
+      /* negative case */
+      carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);      /* make positive */
+      *c = (Unit) (carry % (DECDPUNMAX + 1));
+      carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
+#endif
+    }                          /* c */
+
+  /* we now may have one or other to complete */
+  /* [pretest to avoid loop setup/shutdown] */
+  if (c < maxC)
+    for (; c < maxC; c++)
+      {
+       if (a < alsu + alength)
+         {                     /* still in A */
+           carry += *a;
+           a++;
+         }
+       else
+         {                     /* inside B */
+           carry += ((eInt) * b) * m;
+           b++;
+         }
+       /* here carry is new Unit of digits; it could be +ve or -ve and */
+       /* magnitude up to DECDPUNMAX squared */
+       if ((ueInt) carry <= DECDPUNMAX)
+         {                     /* fastpath 0-DECDPUNMAX */
+           *c = (Unit) carry;
+           carry = 0;
+           continue;
+         }
+       /* result for this unit is negative or >DECDPUNMAX */
+#if DECDPUN==4                 /* use divide-by-multiply */
+       /* remainder is undefined if negative, so we must test */
+       if (carry >= 0)
+         {
+           est = (((ueInt) carry >> 11) * 53687) >> 18;
+           *c = (Unit) (carry - est * (DECDPUNMAX + 1));       /* remainder */
+           carry = est;        /* likely quotient [79.7%] */
+           if (*c < DECDPUNMAX + 1)
+             continue;         /* estimate was correct */
+           carry++;
+           *c -= DECDPUNMAX + 1;
+           continue;
+         }
+       /* negative case */
+       carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);     /* make positive */
+       est = (((ueInt) carry >> 11) * 53687) >> 18;
+       *c = (Unit) (carry - est * (DECDPUNMAX + 1));
+       carry = est - (DECDPUNMAX + 1); /* correctly negative */
+       if (*c < DECDPUNMAX + 1)
+         continue;             /* was OK */
+       carry++;
+       *c -= DECDPUNMAX + 1;
+#else
+       if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
+         {                     /* fastpath carry 1 */
+           *c = (Unit) (carry - (DECDPUNMAX + 1));
+           carry = 1;
+           continue;
+         }
+       /* remainder is undefined if negative, so we must test */
+       if (carry >= 0)
+         {
+           *c = (Unit) (carry % (DECDPUNMAX + 1));
+           carry = carry / (DECDPUNMAX + 1);
+           continue;
+         }
+       /* negative case */
+       carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);     /* make positive */
+       *c = (Unit) (carry % (DECDPUNMAX + 1));
+       carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
+#endif
+      }                                /* c */
+
+  /* OK, all A and B processed; might still have carry or borrow */
+  /* return number of Units in the result, negated if a borrow */
+  if (carry == 0)
+    return c - clsu;           /* no carry, we're done */
+  if (carry > 0)
+    {                          /* positive carry */
+      *c = (Unit) carry;       /* place as new unit */
+      c++;                     /* .. */
+      return c - clsu;
+    }
+  /* -ve carry: it's a borrow; complement needed */
+  add = 1;                     /* temporary carry... */
+  for (c = clsu; c < maxC; c++)
+    {
+      add = DECDPUNMAX + add - *c;
+      if (add <= DECDPUNMAX)
+       {
+         *c = (Unit) add;
+         add = 0;
+       }
+      else
+       {
+         *c = 0;
+         add = 1;
+       }
+    }
+  /* add an extra unit iff it would be non-zero */
+#if DECTRACE
+  printf ("UAS borrow: add %d, carry %d\n", add, carry);
+#endif
+  if ((add - carry - 1) != 0)
+    {
+      *c = (Unit) (add - carry - 1);
+      c++;                     /* interesting, include it */
+    }
+  return clsu - c;             /* -ve result indicates borrowed */
+}
+
+/* ------------------------------------------------------------------ */
+/* decTrim -- trim trailing zeros or normalize                        */
+/*                                                                    */
+/*   dn is the number to trim or normalize                            */
+/*   all is 1 to remove all trailing zeros, 0 for just fraction ones  */
+/*   dropped returns the number of discarded trailing zeros           */
+/*   returns dn                                                       */
+/*                                                                    */
+/* All fields are updated as required.  This is a utility operation,  */
+/* so special values are unchanged and no error is possible.          */
+/* ------------------------------------------------------------------ */
+static decNumber *
+decTrim (decNumber * dn, Flag all, Int * dropped)
+{
+  Int d, exp;                  /* work */
+  uInt cut;                    /* .. */
+  Unit *up;                    /* -> current Unit */
+
+#if DECCHECK
+  if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
+    return dn;
+#endif
+
+  *dropped = 0;                        /* assume no zeros dropped */
+  if ((dn->bits & DECSPECIAL)  /* fast exit if special .. */
+      || (*dn->lsu & 0x01))
+    return dn;                 /* .. or odd */
+  if (ISZERO (dn))
+    {                          /* .. or 0 */
+      dn->exponent = 0;                /* (sign is preserved) */
+      return dn;
+    }
+
+  /* we have a finite number which is even */
+  exp = dn->exponent;
+  cut = 1;                     /* digit (1-DECDPUN) in Unit */
+  up = dn->lsu;                        /* -> current Unit */
+  for (d = 0; d < dn->digits - 1; d++)
+    {                          /* [don't strip the final digit] */
+      /* slice by powers */
+#if DECDPUN<=4
+      uInt quot = QUOT10 (*up, cut);
+      if ((*up - quot * powers[cut]) != 0)
+       break;                  /* found non-0 digit */
+#else
+      if (*up % powers[cut] != 0)
+       break;                  /* found non-0 digit */
+#endif
+      /* have a trailing 0 */
+      if (!all)
+       {                       /* trimming */
+         /* [if exp>0 then all trailing 0s are significant for trim] */
+         if (exp <= 0)
+           {                   /* if digit might be significant */
+             if (exp == 0)
+               break;          /* then quit */
+             exp++;            /* next digit might be significant */
+           }
+       }
+      cut++;                   /* next power */
+      if (cut > DECDPUN)
+       {                       /* need new Unit */
+         up++;
+         cut = 1;
+       }
+    }                          /* d */
+  if (d == 0)
+    return dn;                 /* none dropped */
+
+  /* effect the drop */
+  decShiftToLeast (dn->lsu, D2U (dn->digits), d);
+  dn->exponent += d;           /* maintain numerical value */
+  dn->digits -= d;             /* new length */
+  *dropped = d;                        /* report the count */
+  return dn;
+}
+
+/* ------------------------------------------------------------------ */
+/* decShiftToMost -- shift digits in array towards most significant   */
+/*                                                                    */
+/*   uar    is the array                                              */
+/*   digits is the count of digits in use in the array                */
+/*   shift  is the number of zeros to pad with (least significant);   */
+/*     it must be zero or positive                                    */
+/*                                                                    */
+/*   returns the new length of the integer in the array, in digits    */
+/*                                                                    */
+/* No overflow is permitted (that is, the uar array must be known to  */
+/* be large enough to hold the result, after shifting).               */
+/* ------------------------------------------------------------------ */
+static Int
+decShiftToMost (Unit * uar, Int digits, Int shift)
+{
+  Unit *target, *source, *first;       /* work */
+  uInt rem;                    /* for division */
+  Int cut;                     /* odd 0's to add */
+  uInt next;                   /* work */
+
+  if (shift == 0)
+    return digits;             /* [fastpath] nothing to do */
+  if ((digits + shift) <= DECDPUN)
+    {                          /* [fastpath] single-unit case */
+      *uar = (Unit) (*uar * powers[shift]);
+      return digits + shift;
+    }
+
+  cut = (DECDPUN - shift % DECDPUN) % DECDPUN;
+  source = uar + D2U (digits) - 1;     /* where msu comes from */
+  first = uar + D2U (digits + shift) - 1;      /* where msu of source will end up */
+  target = source + D2U (shift);       /* where upper part of first cut goes */
+  next = 0;
+
+  for (; source >= uar; source--, target--)
+    {
+      /* split the source Unit and accumulate remainder for next */
+#if DECDPUN<=4
+      uInt quot = QUOT10 (*source, cut);
+      rem = *source - quot * powers[cut];
+      next += quot;
+#else
+      rem = *source % powers[cut];
+      next += *source / powers[cut];
+#endif
+      if (target <= first)
+       *target = (Unit) next;  /* write to target iff valid */
+      next = rem * powers[DECDPUN - cut];      /* save remainder for next Unit */
+    }
+  /* propagate to one below and clear the rest */
+  for (; target >= uar; target--)
+    {
+      *target = (Unit) next;
+      next = 0;
+    }
+  return digits + shift;
+}
+
+/* ------------------------------------------------------------------ */
+/* decShiftToLeast -- shift digits in array towards least significant */
+/*                                                                    */
+/*   uar   is the array                                               */
+/*   units is length of the array, in units                           */
+/*   shift is the number of digits to remove from the lsu end; it     */
+/*     must be zero or positive and less than units*DECDPUN.          */
+/*                                                                    */
+/*   returns the new length of the integer in the array, in units     */
+/*                                                                    */
+/* Removed digits are discarded (lost).  Units not required to hold   */
+/* the final result are unchanged.                                    */
+/* ------------------------------------------------------------------ */
+static Int
+decShiftToLeast (Unit * uar, Int units, Int shift)
+{
+  Unit *target, *up;           /* work */
+  Int cut, count;              /* work */
+  Int quot, rem;               /* for division */
+
+  if (shift == 0)
+    return units;              /* [fastpath] nothing to do */
+
+  up = uar + shift / DECDPUN;  /* source; allow for whole Units */
+  cut = shift % DECDPUN;       /* odd 0's to drop */
+  target = uar;                        /* both paths */
+  if (cut == 0)
+    {                          /* whole units shift */
+      for (; up < uar + units; target++, up++)
+       *target = *up;
+      return target - uar;
+    }
+  /* messier */
+  count = units * DECDPUN - shift;     /* the maximum new length */
+#if DECDPUN<=4
+  quot = QUOT10 (*up, cut);
+#else
+  quot = *up / powers[cut];
+#endif
+  for (;; target++)
+    {
+      *target = (Unit) quot;
+      count -= (DECDPUN - cut);
+      if (count <= 0)
+       break;
+      up++;
+      quot = *up;
+#if DECDPUN<=4
+      quot = QUOT10 (quot, cut);
+      rem = *up - quot * powers[cut];
+#else
+      rem = quot % powers[cut];
+      quot = quot / powers[cut];
+#endif
+      *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
+      count -= cut;
+      if (count <= 0)
+       break;
+    }
+  return target - uar + 1;
+}
+
+#if DECSUBSET
+/* ------------------------------------------------------------------ */
+/* decRoundOperand -- round an operand  [used for subset only]        */
+/*                                                                    */
+/*   dn is the number to round (dn->digits is > set->digits)          */
+/*   set is the relevant context                                      */
+/*   status is the status accumulator                                 */
+/*                                                                    */
+/*   returns an allocated decNumber with the rounded result.          */
+/*                                                                    */
+/* lostDigits and other status may be set by this.                    */
+/*                                                                    */
+/* Since the input is an operand, we are not permitted to modify it.  */
+/* We therefore return an allocated decNumber, rounded as required.   */
+/* It is the caller's responsibility to free the allocated storage.   */
+/*                                                                    */
+/* If no storage is available then the result cannot be used, so NULL */
+/* is returned.                                                       */
+/* ------------------------------------------------------------------ */
+static decNumber *
+decRoundOperand (decNumber * dn, decContext * set, uInt * status)
+{
+  decNumber *res;              /* result structure */
+  uInt newstatus = 0;          /* status from round */
+  Int residue = 0;             /* rounding accumulator */
+
+  /* Allocate storage for the returned decNumber, big enough for the */
+  /* length specified by the context */
+  res = (decNumber *) malloc (sizeof (decNumber)
+                             + (D2U (set->digits) - 1) * sizeof (Unit));
+  if (res == NULL)
+    {
+      *status |= DEC_Insufficient_storage;
+      return NULL;
+    }
+  decCopyFit (res, dn, set, &residue, &newstatus);
+  decApplyRound (res, set, residue, &newstatus);
+
+  /* If that set Inexact then we "lost digits" */
+  if (newstatus & DEC_Inexact)
+    newstatus |= DEC_Lost_digits;
+  *status |= newstatus;
+  return res;
+}
+#endif
+
+/* ------------------------------------------------------------------ */
+/* decCopyFit -- copy a number, shortening the coefficient if needed  */
+/*                                                                    */
+/*   dest is the target decNumber                                     */
+/*   src  is the source decNumber                                     */
+/*   set is the context [used for length (digits) and rounding mode]  */
+/*   residue is the residue accumulator                               */
+/*   status contains the current status to be updated                 */
+/*                                                                    */
+/* (dest==src is allowed and will be a no-op if fits)                 */
+/* All fields are updated as required.                                */
+/* ------------------------------------------------------------------ */
+static void
+decCopyFit (decNumber * dest, decNumber * src, decContext * set,
+           Int * residue, uInt * status)
+{
+  dest->bits = src->bits;
+  dest->exponent = src->exponent;
+  decSetCoeff (dest, set, src->lsu, src->digits, residue, status);
+}
+
+/* ------------------------------------------------------------------ */
+/* decSetCoeff -- set the coefficient of a number                     */
+/*                                                                    */
+/*   dn    is the number whose coefficient array is to be set.        */
+/*         It must have space for set->digits digits                  */
+/*   set   is the context [for size]                                  */
+/*   lsu   -> lsu of the source coefficient [may be dn->lsu]          */
+/*   len   is digits in the source coefficient [may be dn->digits]    */
+/*   residue is the residue accumulator.  This has values as in       */
+/*         decApplyRound, and will be unchanged unless the            */
+/*         target size is less than len.  In this case, the           */
+/*         coefficient is truncated and the residue is updated to     */
+/*         reflect the previous residue and the dropped digits.       */
+/*   status is the status accumulator, as usual                       */
+/*                                                                    */
+/* The coefficient may already be in the number, or it can be an      */
+/* external intermediate array.  If it is in the number, lsu must ==  */
+/* dn->lsu and len must == dn->digits.                                */
+/*                                                                    */
+/* Note that the coefficient length (len) may be < set->digits, and   */
+/* in this case this merely copies the coefficient (or is a no-op     */
+/* if dn->lsu==lsu).                                                  */
+/*                                                                    */
+/* Note also that (only internally, from decNumberRescale and         */
+/* decSetSubnormal) the value of set->digits may be less than one,    */
+/* indicating a round to left.                                        */
+/* This routine handles that case correctly; caller ensures space.    */
+/*                                                                    */
+/* dn->digits, dn->lsu (and as required), and dn->exponent are        */
+/* updated as necessary.   dn->bits (sign) is unchanged.              */
+/*                                                                    */
+/* DEC_Rounded status is set if any digits are discarded.             */
+/* DEC_Inexact status is set if any non-zero digits are discarded, or */
+/*                       incoming residue was non-0 (implies rounded) */
+/* ------------------------------------------------------------------ */
+/* mapping array: maps 0-9 to canonical residues, so that we can */
+/* adjust by a residue in range [-1, +1] and achieve correct rounding */
+/*                             0  1  2  3  4  5  6  7  8  9 */
+static const uByte resmap[10] = { 0, 3, 3, 3, 3, 5, 7, 7, 7, 7 };
+static void
+decSetCoeff (decNumber * dn, decContext * set, Unit * lsu,
+            Int len, Int * residue, uInt * status)
+{
+  Int discard;                 /* number of digits to discard */
+  uInt discard1;               /* first discarded digit */
+  uInt cut;                    /* cut point in Unit */
+  uInt quot, rem;              /* for divisions */
+  Unit *up, *target;           /* work */
+  Int count;                   /* .. */
+#if DECDPUN<=4
+  uInt temp;                   /* .. */
+#endif
+
+  discard = len - set->digits; /* digits to discard */
+  if (discard <= 0)
+    {                          /* no digits are being discarded */
+      if (dn->lsu != lsu)
+       {                       /* copy needed */
+         /* copy the coefficient array to the result number; no shift needed */
+         up = lsu;
+         for (target = dn->lsu; target < dn->lsu + D2U (len); target++, up++)
+           {
+             *target = *up;
+           }
+         dn->digits = len;     /* set the new length */
+       }
+      /* dn->exponent and residue are unchanged */
+      if (*residue != 0)
+       *status |= (DEC_Inexact | DEC_Rounded); /* record inexactitude */
+      return;
+    }
+
+  /* we have to discard some digits */
+  *status |= DEC_Rounded;      /* accumulate Rounded status */
+  if (*residue > 1)
+    *residue = 1;              /* previous residue now to right, so -1 to +1 */
+
+  if (discard > len)
+    {                          /* everything, +1, is being discarded */
+      /* guard digit is 0 */
+      /* residue is all the number [NB could be all 0s] */
+      if (*residue <= 0)
+       for (up = lsu + D2U (len) - 1; up >= lsu; up--)
+         {
+           if (*up != 0)
+             {                 /* found a non-0 */
+               *residue = 1;
+               break;          /* no need to check any others */
+             }
+         }
+      if (*residue != 0)
+       *status |= DEC_Inexact; /* record inexactitude */
+      *dn->lsu = 0;            /* coefficient will now be 0 */
+      dn->digits = 1;          /* .. */
+      dn->exponent += discard; /* maintain numerical value */
+      return;
+    }                          /* total discard */
+
+  /* partial discard [most common case] */
+  /* here, at least the first (most significant) discarded digit exists */
+
+  /* spin up the number, noting residue as we pass, until we get to */
+  /* the Unit with the first discarded digit.  When we get there, */
+  /* extract it and remember where we're at */
+  count = 0;
+  for (up = lsu;; up++)
+    {
+      count += DECDPUN;
+      if (count >= discard)
+       break;                  /* full ones all checked */
+      if (*up != 0)
+       *residue = 1;
+    }                          /* up */
+
+  /* here up -> Unit with discarded digit */
+  cut = discard - (count - DECDPUN) - 1;
+  if (cut == DECDPUN - 1)
+    {                          /* discard digit is at top */
+#if DECDPUN<=4
+      discard1 = QUOT10 (*up, DECDPUN - 1);
+      rem = *up - discard1 * powers[DECDPUN - 1];
+#else
+      rem = *up % powers[DECDPUN - 1];
+      discard1 = *up / powers[DECDPUN - 1];
+#endif
+      if (rem != 0)
+       *residue = 1;
+      up++;                    /* move to next */
+      cut = 0;                 /* bottom digit of result */
+      quot = 0;                        /* keep a certain compiler happy */
+    }
+  else
+    {
+      /* discard digit is in low digit(s), not top digit */
+      if (cut == 0)
+       quot = *up;
+      else                     /* cut>0 */
+       {                       /* it's not at bottom of Unit */
+#if DECDPUN<=4
+         quot = QUOT10 (*up, cut);
+         rem = *up - quot * powers[cut];
+#else
+         rem = *up % powers[cut];
+         quot = *up / powers[cut];
+#endif
+         if (rem != 0)
+           *residue = 1;
+       }
+      /* discard digit is now at bottom of quot */
+#if DECDPUN<=4
+      temp = (quot * 6554) >> 16;      /* fast /10 */
+      /* Vowels algorithm here not a win (9 instructions) */
+      discard1 = quot - X10 (temp);
+      quot = temp;
+#else
+      discard1 = quot % 10;
+      quot = quot / 10;
+#endif
+      cut++;                   /* update cut */
+    }
+
+  /* here: up -> Unit of the array with discarded digit */
+  /*       cut is the division point for each Unit */
+  /*       quot holds the uncut high-order digits for the current */
+  /*            Unit, unless cut==0 in which case it's still in *up */
+  /* copy the coefficient array to the result number, shifting as we go */
+  count = set->digits;         /* digits to end up with */
+  if (count <= 0)
+    {                          /* special for Rescale/Subnormal :-( */
+      *dn->lsu = 0;            /* .. result is 0 */
+      dn->digits = 1;          /* .. */
+    }
+  else
+    {                          /* shift to least */
+      /* [this is similar to decShiftToLeast code, with copy] */
+      dn->digits = count;      /* set the new length */
+      if (cut == 0)
+       {
+         /* on unit boundary, so simple shift down copy loop suffices */
+         for (target = dn->lsu; target < dn->lsu + D2U (count);
+              target++, up++)
+           {
+             *target = *up;
+           }
+       }
+      else
+       for (target = dn->lsu;; target++)
+         {
+           *target = (Unit) quot;
+           count -= (DECDPUN - cut);
+           if (count <= 0)
+             break;
+           up++;
+           quot = *up;
+#if DECDPUN<=4
+           quot = QUOT10 (quot, cut);
+           rem = *up - quot * powers[cut];
+#else
+           rem = quot % powers[cut];
+           quot = quot / powers[cut];
+#endif
+           *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
+           count -= cut;
+           if (count <= 0)
+             break;
+         }
+    }                          /* shift to least needed */
+  dn->exponent += discard;     /* maintain numerical value */
+
+  /* here, discard1 is the guard digit, and residue is everything else */
+  /* [use mapping to accumulate residue safely] */
+  *residue += resmap[discard1];
+
+  if (*residue != 0)
+    *status |= DEC_Inexact;    /* record inexactitude */
+  return;
+}
+
+/* ------------------------------------------------------------------ */
+/* decApplyRound -- apply pending rounding to a number                */
+/*                                                                    */
+/*   dn    is the number, with space for set->digits digits           */
+/*   set   is the context [for size and rounding mode]                */
+/*   residue indicates pending rounding, being any accumulated        */
+/*         guard and sticky information.  It may be:                  */
+/*         6-9: rounding digit is >5                                  */
+/*         5:   rounding digit is exactly half-way                    */
+/*         1-4: rounding digit is <5 and >0                           */
+/*         0:   the coefficient is exact                              */
+/*        -1:   as 1, but the hidden digits are subtractive, that     */
+/*              is, of the opposite sign to dn.  In this case the     */
+/*              coefficient must be non-0.                            */
+/*   status is the status accumulator, as usual                       */
+/*                                                                    */
+/* This routine applies rounding while keeping the length of the      */
+/* coefficient constant.  The exponent and status are unchanged       */
+/* except if:                                                         */
+/*                                                                    */
+/*   -- the coefficient was increased and is all nines (in which      */
+/*      case Overflow could occur, and is handled directly here so    */
+/*      the caller does not need to re-test for overflow)             */
+/*                                                                    */
+/*   -- the coefficient was decreased and becomes all nines (in which */
+/*      case Underflow could occur, and is also handled directly).    */
+/*                                                                    */
+/* All fields in dn are updated as required.                          */
+/*                                                                    */
+/* ------------------------------------------------------------------ */
+static void
+decApplyRound (decNumber * dn, decContext * set, Int residue, uInt * status)
+{
+  Int bump;                    /* 1 if coefficient needs to be incremented */
+  /* -1 if coefficient needs to be decremented */
+
+  if (residue == 0)
+    return;                    /* nothing to apply */
+
+  bump = 0;                    /* assume a smooth ride */
+
+  /* now decide whether, and how, to round, depending on mode */
+  switch (set->round)
+    {
+    case DEC_ROUND_DOWN:
+      {
+       /* no change, except if negative residue */
+       if (residue < 0)
+         bump = -1;
+       break;
+      }                                /* r-d */
+
+    case DEC_ROUND_HALF_DOWN:
+      {
+       if (residue > 5)
+         bump = 1;
+       break;
+      }                                /* r-h-d */
+
+    case DEC_ROUND_HALF_EVEN:
+      {
+       if (residue > 5)
+         bump = 1;             /* >0.5 goes up */
+       else if (residue == 5)
+         {                     /* exactly 0.5000... */
+           /* 0.5 goes up iff [new] lsd is odd */
+           if (*dn->lsu & 0x01)
+             bump = 1;
+         }
+       break;
+      }                                /* r-h-e */
+
+    case DEC_ROUND_HALF_UP:
+      {
+       if (residue >= 5)
+         bump = 1;
+       break;
+      }                                /* r-h-u */
+
+    case DEC_ROUND_UP:
+      {
+       if (residue > 0)
+         bump = 1;
+       break;
+      }                                /* r-u */
+
+    case DEC_ROUND_CEILING:
+      {
+       /* same as _UP for positive numbers, and as _DOWN for negatives */
+       /* [negative residue cannot occur on 0] */
+       if (decNumberIsNegative (dn))
+         {
+           if (residue < 0)
+             bump = -1;
+         }
+       else
+         {
+           if (residue > 0)
+             bump = 1;
+         }
+       break;
+      }                                /* r-c */
+
+    case DEC_ROUND_FLOOR:
+      {
+       /* same as _UP for negative numbers, and as _DOWN for positive */
+       /* [negative residue cannot occur on 0] */
+       if (!decNumberIsNegative (dn))
+         {
+           if (residue < 0)
+             bump = -1;
+         }
+       else
+         {
+           if (residue > 0)
+             bump = 1;
+         }
+       break;
+      }                                /* r-f */
+
+    default:
+      {                                /* e.g., DEC_ROUND_MAX */
+       *status |= DEC_Invalid_context;
+#if DECTRACE
+       printf ("Unknown rounding mode: %d\n", set->round);
+#endif
+       break;
+      }
+    }                          /* switch */
+
+  /* now bump the number, up or down, if need be */
+  if (bump == 0)
+    return;                    /* no action required */
+
+  /* Simply use decUnitAddSub unless we are bumping up and the number */
+  /* is all nines.  In this special case we set to 1000... and adjust */
+  /* the exponent by one (as otherwise we could overflow the array) */
+  /* Similarly handle all-nines result if bumping down. */
+  if (bump > 0)
+    {
+      Unit *up;                        /* work */
+      uInt count = dn->digits; /* digits to be checked */
+      for (up = dn->lsu;; up++)
+       {
+         if (count <= DECDPUN)
+           {
+             /* this is the last Unit (the msu) */
+             if (*up != powers[count] - 1)
+               break;          /* not still 9s */
+             /* here if it, too, is all nines */
+             *up = (Unit) powers[count - 1];   /* here 999 -> 100 etc. */
+             for (up = up - 1; up >= dn->lsu; up--)
+               *up = 0;        /* others all to 0 */
+             dn->exponent++;   /* and bump exponent */
+             /* [which, very rarely, could cause Overflow...] */
+             if ((dn->exponent + dn->digits) > set->emax + 1)
+               {
+                 decSetOverflow (dn, set, status);
+               }
+             return;           /* done */
+           }
+         /* a full unit to check, with more to come */
+         if (*up != DECDPUNMAX)
+           break;              /* not still 9s */
+         count -= DECDPUN;
+       }                       /* up */
+    }                          /* bump>0 */
+  else
+    {                          /* -1 */
+      /* here we are lookng for a pre-bump of 1000... (leading 1, */
+      /* all other digits zero) */
+      Unit *up, *sup;          /* work */
+      uInt count = dn->digits; /* digits to be checked */
+      for (up = dn->lsu;; up++)
+       {
+         if (count <= DECDPUN)
+           {
+             /* this is the last Unit (the msu) */
+             if (*up != powers[count - 1])
+               break;          /* not 100.. */
+             /* here if we have the 1000... case */
+             sup = up;         /* save msu pointer */
+             *up = (Unit) powers[count] - 1;   /* here 100 in msu -> 999 */
+             /* others all to all-nines, too */
+             for (up = up - 1; up >= dn->lsu; up--)
+               *up = (Unit) powers[DECDPUN] - 1;
+             dn->exponent--;   /* and bump exponent */
+
+             /* iff the number was at the subnormal boundary (exponent=etiny) */
+             /* then the exponent is now out of range, so it will in fact get */
+             /* clamped to etiny and the final 9 dropped. */
+             /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
+             /*        dn->exponent, set->digits); */
+             if (dn->exponent + 1 == set->emin - set->digits + 1)
+               {
+                 if (count == 1 && dn->digits == 1)
+                   *sup = 0;   /* here 9 -> 0[.9] */
+                 else
+                   {
+                     *sup = (Unit) powers[count - 1] - 1;      /* here 999.. in msu -> 99.. */
+                     dn->digits--;
+                   }
+                 dn->exponent++;
+                 *status |=
+                   DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
+               }
+             return;           /* done */
+           }
+
+         /* a full unit to check, with more to come */
+         if (*up != 0)
+           break;              /* not still 0s */
+         count -= DECDPUN;
+       }                       /* up */
+
+    }                          /* bump<0 */
+
+  /* Actual bump needed.  Do it. */
+  decUnitAddSub (dn->lsu, D2U (dn->digits), one, 1, 0, dn->lsu, bump);
+}
+
+#if DECSUBSET
+/* ------------------------------------------------------------------ */
+/* decFinish -- finish processing a number                            */
+/*                                                                    */
+/*   dn is the number                                                 */
+/*   set is the context                                               */
+/*   residue is the rounding accumulator (as in decApplyRound)        */
+/*   status is the accumulator                                        */
+/*                                                                    */
+/* This finishes off the current number by:                           */
+/*    1. If not extended:                                             */
+/*       a. Converting a zero result to clean '0'                     */
+/*       b. Reducing positive exponents to 0, if would fit in digits  */
+/*    2. Checking for overflow and subnormals (always)                */
+/* Note this is just Finalize when no subset arithmetic.              */
+/* All fields are updated as required.                                */
+/* ------------------------------------------------------------------ */
+static void
+decFinish (decNumber * dn, decContext * set, Int * residue, uInt * status)
+{
+  if (!set->extended)
+    {
+      if ISZERO
+       (dn)
+       {                       /* value is zero */
+         dn->exponent = 0;     /* clean exponent .. */
+         dn->bits = 0;         /* .. and sign */
+         return;               /* no error possible */
+       }
+      if (dn->exponent >= 0)
+       {                       /* non-negative exponent */
+         /* >0; reduce to integer if possible */
+         if (set->digits >= (dn->exponent + dn->digits))
+           {
+             dn->digits = decShiftToMost (dn->lsu, dn->digits, dn->exponent);
+             dn->exponent = 0;
+           }
+       }
+    }                          /* !extended */
+
+  decFinalize (dn, set, residue, status);
+}
+#endif
+
+/* ------------------------------------------------------------------ */
+/* decFinalize -- final check, clamp, and round of a number           */
+/*                                                                    */
+/*   dn is the number                                                 */
+/*   set is the context                                               */
+/*   residue is the rounding accumulator (as in decApplyRound)        */
+/*   status is the status accumulator                                 */
+/*                                                                    */
+/* This finishes off the current number by checking for subnormal     */
+/* results, applying any pending rounding, checking for overflow,     */
+/* and applying any clamping.                                         */
+/* Underflow and overflow conditions are raised as appropriate.       */
+/* All fields are updated as required.                                */
+/* ------------------------------------------------------------------ */
+static void
+decFinalize (decNumber * dn, decContext * set, Int * residue, uInt * status)
+{
+  Int shift;                   /* shift needed if clamping */
+
+  /* We have to be careful when checking the exponent as the adjusted */
+  /* exponent could overflow 31 bits [because it may already be up */
+  /* to twice the expected]. */
+
+  /* First test for subnormal.  This must be done before any final */
+  /* round as the result could be rounded to Nmin or 0. */
+  if (dn->exponent < 0         /* negative exponent */
+      && (dn->exponent < set->emin - dn->digits + 1))
+    {
+      /* Go handle subnormals; this will apply round if needed. */
+      decSetSubnormal (dn, set, residue, status);
+      return;
+    }
+
+  /* now apply any pending round (this could raise overflow). */
+  if (*residue != 0)
+    decApplyRound (dn, set, *residue, status);
+
+  /* Check for overflow [redundant in the 'rare' case] or clamp */
+  if (dn->exponent <= set->emax - set->digits + 1)
+    return;                    /* neither needed */
+
+  /* here when we might have an overflow or clamp to do */
+  if (dn->exponent > set->emax - dn->digits + 1)
+    {                          /* too big */
+      decSetOverflow (dn, set, status);
+      return;
+    }
+  /* here when the result is normal but in clamp range */
+  if (!set->clamp)
+    return;
+
+  /* here when we need to apply the IEEE exponent clamp (fold-down) */
+  shift = dn->exponent - (set->emax - set->digits + 1);
+
+  /* shift coefficient (if non-zero) */
+  if (!ISZERO (dn))
+    {
+      dn->digits = decShiftToMost (dn->lsu, dn->digits, shift);
+    }
+  dn->exponent -= shift;       /* adjust the exponent to match */
+  *status |= DEC_Clamped;      /* and record the dirty deed */
+  return;
+}
+
+/* ------------------------------------------------------------------ */
+/* decSetOverflow -- set number to proper overflow value              */
+/*                                                                    */
+/*   dn is the number (used for sign [only] and result)               */
+/*   set is the context [used for the rounding mode]                  */
+/*   status contains the current status to be updated                 */
+/*                                                                    */
+/* This sets the sign of a number and sets its value to either        */
+/* Infinity or the maximum finite value, depending on the sign of     */
+/* dn and therounding mode, following IEEE 854 rules.                 */
+/* ------------------------------------------------------------------ */
+static void
+decSetOverflow (decNumber * dn, decContext * set, uInt * status)
+{
+  Flag needmax = 0;            /* result is maximum finite value */
+  uByte sign = dn->bits & DECNEG;      /* clean and save sign bit */
+
+  if (ISZERO (dn))
+    {                          /* zero does not overflow magnitude */
+      Int emax = set->emax;    /* limit value */
+      if (set->clamp)
+       emax -= set->digits - 1;        /* lower if clamping */
+      if (dn->exponent > emax)
+       {                       /* clamp required */
+         dn->exponent = emax;
+         *status |= DEC_Clamped;
+       }
+      return;
+    }
+
+  decNumberZero (dn);
+  switch (set->round)
+    {
+    case DEC_ROUND_DOWN:
+      {
+       needmax = 1;            /* never Infinity */
+       break;
+      }                                /* r-d */
+    case DEC_ROUND_CEILING:
+      {
+       if (sign)
+         needmax = 1;          /* Infinity if non-negative */
+       break;
+      }                                /* r-c */
+    case DEC_ROUND_FLOOR:
+      {
+       if (!sign)
+         needmax = 1;          /* Infinity if negative */
+       break;
+      }                                /* r-f */
+    default:
+      break;                   /* Infinity in all other cases */
+    }
+  if (needmax)
+    {
+      Unit *up;                        /* work */
+      Int count = set->digits; /* nines to add */
+      dn->digits = count;
+      /* fill in all nines to set maximum value */
+      for (up = dn->lsu;; up++)
+       {
+         if (count > DECDPUN)
+           *up = DECDPUNMAX;   /* unit full o'nines */
+         else
+           {                   /* this is the msu */
+             *up = (Unit) (powers[count] - 1);
+             break;
+           }
+         count -= DECDPUN;     /* we filled those digits */
+       }                       /* up */
+      dn->bits = sign;         /* sign */
+      dn->exponent = set->emax - set->digits + 1;
+    }
+  else
+    dn->bits = sign | DECINF;  /* Value is +/-Infinity */
+  *status |= DEC_Overflow | DEC_Inexact | DEC_Rounded;
+}
+
+/* ------------------------------------------------------------------ */
+/* decSetSubnormal -- process value whose exponent is <Emin           */
+/*                                                                    */
+/*   dn is the number (used as input as well as output; it may have   */
+/*         an allowed subnormal value, which may need to be rounded)  */
+/*   set is the context [used for the rounding mode]                  */
+/*   residue is any pending residue                                   */
+/*   status contains the current status to be updated                 */
+/*                                                                    */
+/* If subset mode, set result to zero and set Underflow flags.        */
+/*                                                                    */
+/* Value may be zero with a low exponent; this does not set Subnormal */
+/* but the exponent will be clamped to Etiny.                         */
+/*                                                                    */
+/* Otherwise ensure exponent is not out of range, and round as        */
+/* necessary.  Underflow is set if the result is Inexact.             */
+/* ------------------------------------------------------------------ */
+static void
+decSetSubnormal (decNumber * dn, decContext * set, Int * residue,
+                uInt * status)
+{
+  decContext workset;          /* work */
+  Int etiny, adjust;           /* .. */
+
+#if DECSUBSET
+  /* simple set to zero and 'hard underflow' for subset */
+  if (!set->extended)
+    {
+      decNumberZero (dn);
+      /* always full overflow */
+      *status |= DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
+      return;
+    }
+#endif
+
+  /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
+  /* (Etiny) if needed */
+  etiny = set->emin - (set->digits - 1);       /* smallest allowed exponent */
+
+  if ISZERO
+    (dn)
+    {                          /* value is zero */
+      /* residue can never be non-zero here */
+#if DECCHECK
+      if (*residue != 0)
+       {
+         printf ("++ Subnormal 0 residue %d\n", *residue);
+         *status |= DEC_Invalid_operation;
+       }
+#endif
+      if (dn->exponent < etiny)
+       {                       /* clamp required */
+         dn->exponent = etiny;
+         *status |= DEC_Clamped;
+       }
+      return;
+    }
+
+  *status |= DEC_Subnormal;    /* we have a non-zero subnormal */
+
+  adjust = etiny - dn->exponent;       /* calculate digits to remove */
+  if (adjust <= 0)
+    {                          /* not out of range; unrounded */
+      /* residue can never be non-zero here, so fast-path out */
+#if DECCHECK
+      if (*residue != 0)
+       {
+         printf ("++ Subnormal no-adjust residue %d\n", *residue);
+         *status |= DEC_Invalid_operation;
+       }
+#endif
+      /* it may already be inexact (from setting the coefficient) */
+      if (*status & DEC_Inexact)
+       *status |= DEC_Underflow;
+      return;
+    }
+
+  /* adjust>0.  we need to rescale the result so exponent becomes Etiny */
+  /* [this code is similar to that in rescale] */
+  workset = *set;              /* clone rounding, etc. */
+  workset.digits = dn->digits - adjust;        /* set requested length */
+  workset.emin -= adjust;      /* and adjust emin to match */
+  /* [note that the latter can be <1, here, similar to Rescale case] */
+  decSetCoeff (dn, &workset, dn->lsu, dn->digits, residue, status);
+  decApplyRound (dn, &workset, *residue, status);
+
+  /* Use 754R/854 default rule: Underflow is set iff Inexact */
+  /* [independent of whether trapped] */
+  if (*status & DEC_Inexact)
+    *status |= DEC_Underflow;
+
+  /* if we rounded up a 999s case, exponent will be off by one; adjust */
+  /* back if so [it will fit, because we shortened] */
+  if (dn->exponent > etiny)
+    {
+      dn->digits = decShiftToMost (dn->lsu, dn->digits, 1);
+      dn->exponent--;          /* (re)adjust the exponent. */
+    }
+}
+
+/* ------------------------------------------------------------------ */
+/* decGetInt -- get integer from a number                             */
+/*                                                                    */
+/*   dn is the number [which will not be altered]                     */
+/*   set is the context [requested digits], subset only               */
+/*   returns the converted integer, or BADINT if error                */
+/*                                                                    */
+/* This checks and gets a whole number from the input decNumber.      */
+/* The magnitude of the integer must be <2^31.                        */
+/* Any discarded fractional part must be 0.                           */
+/* If subset it must also fit in set->digits                          */
+/* ------------------------------------------------------------------ */
+#if DECSUBSET
+static Int
+decGetInt (decNumber * dn, decContext * set)
+{
+#else
+static Int
+decGetInt (decNumber * dn)
+{
+#endif
+  Int theInt;                  /* result accumulator */
+  Unit *up;                    /* work */
+  Int got;                     /* digits (real or not) processed */
+  Int ilength = dn->digits + dn->exponent;     /* integral length */
+
+  /* The number must be an integer that fits in 10 digits */
+  /* Assert, here, that 10 is enough for any rescale Etiny */
+#if DEC_MAX_EMAX > 999999999
+#error GetInt may need updating [for Emax]
+#endif
+#if DEC_MIN_EMIN < -999999999
+#error GetInt may need updating [for Emin]
+#endif
+  if (ISZERO (dn))
+    return 0;                  /* zeros are OK, with any exponent */
+  if (ilength > 10)
+    return BADINT;             /* always too big */
+#if DECSUBSET
+  if (!set->extended && ilength > set->digits)
+    return BADINT;
+#endif
+
+  up = dn->lsu;                        /* ready for lsu */
+  theInt = 0;                  /* ready to accumulate */
+  if (dn->exponent >= 0)
+    {                          /* relatively easy */
+      /* no fractional part [usual]; allow for positive exponent */
+      got = dn->exponent;
+    }
+  else
+    {                          /* -ve exponent; some fractional part to check and discard */
+      Int count = -dn->exponent;       /* digits to discard */
+      /* spin up whole units until we get to the Unit with the unit digit */
+      for (; count >= DECDPUN; up++)
+       {
+         if (*up != 0)
+           return BADINT;      /* non-zero Unit to discard */
+         count -= DECDPUN;
+       }
+      if (count == 0)
+       got = 0;                /* [a multiple of DECDPUN] */
+      else
+       {                       /* [not multiple of DECDPUN] */
+         Int rem;              /* work */
+         /* slice off fraction digits and check for non-zero */
+#if DECDPUN<=4
+         theInt = QUOT10 (*up, count);
+         rem = *up - theInt * powers[count];
+#else
+         rem = *up % powers[count];    /* slice off discards */
+         theInt = *up / powers[count];
+#endif
+         if (rem != 0)
+           return BADINT;      /* non-zero fraction */
+         /* OK, we're good */
+         got = DECDPUN - count;        /* number of digits so far */
+         up++;                 /* ready for next */
+       }
+    }
+  /* collect the rest */
+  for (; got < ilength; up++)
+    {
+      theInt += *up * powers[got];
+      got += DECDPUN;
+    }
+  if ((ilength == 10)          /* check no wrap */
+      && (theInt / (Int) powers[got - DECDPUN] != *(up - 1)))
+    return BADINT;
+  /* [that test also disallows the BADINT result case] */
+
+  /* apply any sign and return */
+  if (decNumberIsNegative (dn))
+    theInt = -theInt;
+  return theInt;
+}
+
+/* ------------------------------------------------------------------ */
+/* decStrEq -- caseless comparison of strings                         */
+/*                                                                    */
+/*   str1 is one of the strings to compare                            */
+/*   str2 is the other                                                */
+/*                                                                    */
+/*   returns 1 if strings caseless-compare equal, 0 otherwise         */
+/*                                                                    */
+/* Note that the strings must be the same length if they are to       */
+/* compare equal; there is no padding.                                */
+/* ------------------------------------------------------------------ */
+/* [strcmpi is not in ANSI C] */
+static Flag
+decStrEq (const char *str1, const char *str2)
+{
+  for (;; str1++, str2++)
+    {
+      if (*str1 == *str2)
+       {
+         if (*str1 == '\0')
+           break;
+       }
+      else
+       {
+         if (tolower (*str1) != tolower (*str2))
+           return 0;
+       }
+    }                          /* stepping */
+  return 1;
+}
+
+/* ------------------------------------------------------------------ */
+/* decNaNs -- handle NaN operand or operands     &n