/buildworker/marxinbox-gcc-clang-static-analyzer/build/libdecnumber/decNumber.c

Bug Summary

File:	build/libdecnumber/decNumber.c
Warning:	line 4581, column 16 The left operand of '==' is a garbage value
Annotated Source Code

Press '?' to see keyboard shortcuts
Show analyzer invocation
clang -cc1 -cc1 -triple x86_64-suse-linux -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name decNumber.c -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -analyzer-config-compatibility-mode=true -mrelocation-model pic -pic-level 2 -pic-is-pie -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fcoverage-compilation-dir=/buildworker/marxinbox-gcc-clang-static-analyzer/objdir/libdecnumber -resource-dir /usr/lib64/clang/15.0.7 -I /buildworker/marxinbox-gcc-clang-static-analyzer/build/libdecnumber -I . -I /buildworker/marxinbox-gcc-clang-static-analyzer/build/libdecnumber -I . -internal-isystem /usr/lib64/clang/15.0.7/include -internal-isystem /usr/local/include -internal-isystem /usr/bin/../lib64/gcc/x86_64-suse-linux/13/../../../../x86_64-suse-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -Wwrite-strings -Wno-long-long -fconst-strings -fdebug-compilation-dir=/buildworker/marxinbox-gcc-clang-static-analyzer/objdir/libdecnumber -ferror-limit 19 -fgnuc-version=4.2.1 -vectorize-loops -vectorize-slp -analyzer-output=plist-html -analyzer-config silence-checkers=core.NullDereference -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /buildworker/marxinbox-gcc-clang-static-analyzer/objdir/clang-static-analyzer/2023-03-27-141847-20772-1/report-GfSLy1.plist -x c /buildworker/marxinbox-gcc-clang-static-analyzer/build/libdecnumber/decNumber.c
1/* Decimal number arithmetic module for the decNumber C Library.
 Copyright (C) 2005-2023 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 3, 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.

17Under Section 7 of GPL version 3, you are granted additional
18permissions described in the GCC Runtime Library Exception, version
1, as published by the Free Software Foundation.

21You should have received a copy of the GNU General Public License and
22a copy of the GCC Runtime Library Exception along with this program;
23see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24<http://www.gnu.org/licenses/>.  */

26/* ------------------------------------------------------------------ */
27/* Decimal Number arithmetic module				      */
28/* ------------------------------------------------------------------ */
29/* This module comprises the routines for arbitrary-precision General */
30/* Decimal Arithmetic as defined in the specification which may be    */
31/* found on the General Decimal Arithmetic pages.  It implements both */
32/* the full ('extended') arithmetic and the simpler ('subset')	      */
33/* arithmetic.							      */
34/*								      */
35/* Usage notes: 						      */
36/*								      */
37/* 1. This code is ANSI C89 except:				      */
38/*								      */
39/*    a) C99 line comments (double forward slash) are used.  (Most C  */
40/*	 compilers accept these.  If yours does not, a simple script  */
41/*	 can be used to convert them to ANSI C comments.)	      */
42/*								      */
43/*    b) Types from C99 stdint.h are used.  If you do not have this   */
44/*	 header file, see the User's Guide section of the decNumber   */
45/*	 documentation; this lists the necessary definitions.	      */
46/*								      */
47/*    c) If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and       */
48/*	 uint64_t types may be used.  To avoid these, set DECUSE64=0  */
49/*	 and DECDPUN<=4 (see documentation).			      */
50/*								      */
51/*    The code also conforms to C99 restrictions; in particular,      */
52/*    strict aliasing rules are observed.			      */
53/*								      */
54/* 2. The decNumber format which this library uses is optimized for   */
55/*    efficient processing of relatively short numbers; in particular */
56/*    it allows the use of fixed sized structures and minimizes copy  */
57/*    and move operations.  It does, however, support arbitrary       */
58/*    precision (up to 999,999,999 digits) and arbitrary exponent     */
59/*    range (Emax in the range 0 through 999,999,999 and Emin in the  */
60/*    range -999,999,999 through 0).  Mathematical functions (for     */
61/*    example decNumberExp) as identified below are restricted more   */
62/*    tightly: digits, emax, and -emin in the context must be <=      */
63/*    DEC_MAX_MATH (999999), and their operand(s) must be within      */
64/*    these bounds.						      */
65/*								      */
66/* 3. Logical functions are further restricted; their operands must   */
67/*    be finite, positive, have an exponent of zero, and all digits   */
68/*    must be either 0 or 1.  The result will only contain digits     */
69/*    which are 0 or 1 (and will have exponent=0 and a sign of 0).    */
70/*								      */
71/* 4. Operands to operator functions are never modified unless they   */
72/*    are also specified to be the result number (which is always     */
73/*    permitted).  Other than that case, operands must not overlap.   */
74/*								      */
75/* 5. Error handling: the type of the error is ORed into the status   */
76/*    flags in the current context (decContext structure).  The       */
77/*    SIGFPE signal is then raised if the corresponding trap-enabler  */
78/*    flag in the decContext is set (is 1).			      */
79/*								      */
80/*    It is the responsibility of the caller to clear the status      */
81/*    flags as required.					      */
82/*								      */
83/*    The result of any routine which returns a number will always    */
84/*    be a valid number (which may be a special value, such as an     */
85/*    Infinity or NaN). 					      */
86/*								      */
87/* 6. The decNumber format is not an exchangeable concrete	      */
88/*    representation as it comprises fields which may be machine-     */
89/*    dependent (packed or unpacked, or special length, for example). */
90/*    Canonical conversions to and from strings are provided; other   */
91/*    conversions are available in separate modules.		      */
92/*								      */
93/* 7. Normally, input operands are assumed to be valid.  Set DECCHECK */
94/*    to 1 for extended operand checking (including NULL operands).   */
95/*    Results are undefined if a badly-formed structure (or a NULL    */
96/*    pointer to a structure) is provided, though with DECCHECK       */
97/*    enabled the operator routines are protected against exceptions. */
98/*    (Except if the result pointer is NULL, which is unrecoverable.) */
99/*								      */
100/*    However, the routines will never cause exceptions if they are   */
101/*    given well-formed operands, even if the value of the operands   */
102/*    is inappropriate for the operation and DECCHECK is not set.     */
103/*    (Except for SIGFPE, as and where documented.)		      */
104/*								      */
105/* 8. Subset arithmetic is available only if DECSUBSET is set to 1.   */
106/* ------------------------------------------------------------------ */
107/* Implementation notes for maintenance of this module: 	      */
108/*								      */
109/* 1. Storage leak protection:	Routines which use malloc are not     */
110/*    permitted to use return for fastpath or error exits (i.e.,      */
111/*    they follow strict structured programming conventions).	      */
112/*    Instead they have a do{}while(0); construct surrounding the     */
113/*    code which is protected -- break may be used to exit this.      */
114/*    Other routines can safely use the return statement inline.      */
115/*								      */
116/*    Storage leak accounting can be enabled using DECALLOC.	      */
117/*								      */
118/* 2. All loops use the for(;;) construct.  Any do construct does     */
119/*    not loop; it is for allocation protection as just described.    */
120/*								      */
121/* 3. Setting status in the context must always be the very last      */
122/*    action in a routine, as non-0 status may raise a trap and hence */
123/*    the call to set status may not return (if the handler uses long */
124/*    jump).  Therefore all cleanup must be done first.  In general,  */
125/*    to achieve this status is accumulated and is only applied just  */
126/*    before return by calling decContextSetStatus (via decStatus).   */
127/*								      */
128/*    Routines which allocate storage cannot, in general, use the     */
129/*    'top level' routines which could cause a non-returning	      */
130/*    transfer of control.  The decXxxxOp routines are safe (do not   */
131/*    call decStatus even if traps are set in the context) and should */
132/*    be used instead (they are also a little faster).		      */
133/*								      */
134/* 4. Exponent checking is minimized by allowing the exponent to      */
135/*    grow outside its limits during calculations, provided that      */
136/*    the decFinalize function is called later.  Multiplication and   */
137/*    division, and intermediate calculations in exponentiation,      */
138/*    require more careful checks because of the risk of 31-bit       */
139/*    overflow (the most negative valid exponent is -1999999997, for  */
140/*    a 999999999-digit number with adjusted exponent of -999999999). */
141/*								      */
142/* 5. Rounding is deferred until finalization of results, with any    */
143/*    'off to the right' data being represented as a single digit     */
144/*    residue (in the range -1 through 9).  This avoids any double-   */
145/*    rounding when more than one shortening takes place (for	      */
146/*    example, when a result is subnormal).			      */
147/*								      */
148/* 6. The digits count is allowed to rise to a multiple of DECDPUN    */
149/*    during many operations, so whole Units are handled and exact    */
150/*    accounting of digits is not needed.  The correct digits value   */
151/*    is found by decGetDigits, which accounts for leading zeros.     */
152/*    This must be called before any rounding if the number of digits */
153/*    is not known exactly.					      */
154/*								      */
155/* 7. The multiply-by-reciprocal 'trick' is used for partitioning     */
156/*    numbers up to four digits, using appropriate constants.  This   */
157/*    is not useful for longer numbers because overflow of 32 bits    */
158/*    would lead to 4 multiplies, which is almost as expensive as     */
159/*    a divide (unless a floating-point or 64-bit multiply is	      */
160/*    assumed to be available). 				      */
161/*								      */
162/* 8. Unusual abbreviations that may be used in the commentary:       */
163/*	lhs -- left hand side (operand, of an operation)	      */
164/*	lsd -- least significant digit (of coefficient) 	      */
165/*	lsu -- least significant Unit (of coefficient)		      */
166/*	msd -- most significant digit (of coefficient)		      */
167/*	msi -- most significant item (in an array)		      */
168/*	msu -- most significant Unit (of coefficient)		      */
169/*	rhs -- right hand side (operand, of an operation)	      */
170/*	+ve -- positive 					      */
171/*	-ve -- negative 					      */
172/*	**  -- raise to the power				      */
173/* ------------------------------------------------------------------ */

175#include <stdlib.h>		   /* for malloc, free, etc. */
176#include <stdio.h>		   /* for printf [if needed] */
177#include <string.h>		   /* for strcpy */
178#include <ctype.h>		   /* for lower */
179#include "dconfig.h"		   /* for GCC definitions */
180#include "decNumber.h"		   /* base number library */
181#include "decNumberLocal.h"	   /* decNumber local types, etc. */

183/* Constants */
184/* Public lookup table used by the D2U macro */
185const uByteuint8_t d2utable[DECMAXD2U49+1]=D2UTABLE{0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7, 8,8,8,9,9,9,10,
10,10,11,11,11,12,12,12,13,13, 13,14,14,14,15,15,15,16,16,16,
17};

187#define DECVERB1     1		   /* set to 1 for verbose DECCHECK */
188#define powersDECPOWERS	    DECPOWERS	   /* old internal name */

190/* Local constants */
191#define DIVIDE0x80	    0x80	   /* Divide operators */
192#define REMAINDER0x40   0x40	   /* .. */
193#define DIVIDEINT0x20   0x20	   /* .. */
194#define REMNEAR0x10     0x10	   /* .. */
195#define COMPARE0x01     0x01	   /* Compare operators */
196#define COMPMAX0x02     0x02	   /* .. */
197#define COMPMIN0x03     0x03	   /* .. */
198#define COMPTOTAL0x04   0x04	   /* .. */
199#define COMPNAN0x05     0x05	   /* .. [NaN processing] */
200#define COMPSIG0x06     0x06	   /* .. [signaling COMPARE] */
201#define COMPMAXMAG0x07  0x07	   /* .. */
202#define COMPMINMAG0x08  0x08	   /* .. */

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

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

212/* Granularity-dependent code */
213#if DECDPUN3<=4
#define eIntint32_t	Intint32_t	      /* extended integer */
#define ueIntuint32_t uIntuint32_t	      /* 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 uIntuint32_t multies[]={131073, 26215, 5243, 1049, 210};
/* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
#define QUOT10(u, n)((((uint32_t)(u)>>(n))*multies[n])>>17) ((((uIntuint32_t)(u)>>(n))*multies[n])>>17)
222#else
/* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
#if !DECUSE641
  #error decNumber.c: DECUSE641 must be 1 when DECDPUN3>4
#endif
#define eIntint32_t	Longint64_t	      /* extended integer */
#define ueIntuint32_t uLonguint64_t	      /* unsigned extended integer */
229#endif

231/* Local routines */
232static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *,
	      decContext *, uByteuint8_t, uIntuint32_t *);
234static Flaguint8_t	   decBiStr(const char *, const char *, const char *);
235static uIntuint32_t	   decCheckMath(const decNumber *, decContext *, uIntuint32_t *);
236static void	   decApplyRound(decNumber *, decContext *, Intint32_t, uIntuint32_t *);
237static Intint32_t	   decCompare(const decNumber *lhs, const decNumber *rhs, Flaguint8_t);
238static decNumber * decCompareOp(decNumber *, const decNumber *,
	      const decNumber *, decContext *,
	      Flaguint8_t, uIntuint32_t *);
241static void	   decCopyFit(decNumber *, const decNumber *, decContext *,
	      Intint32_t *, uIntuint32_t *);
243static decNumber * decDecap(decNumber *, Intint32_t);
244static decNumber * decDivideOp(decNumber *, const decNumber *,
	      const decNumber *, decContext *, Flaguint8_t, uIntuint32_t *);
246static decNumber * decExpOp(decNumber *, const decNumber *,
	      decContext *, uIntuint32_t *);
248static void	   decFinalize(decNumber *, decContext *, Intint32_t *, uIntuint32_t *);
249static Intint32_t	   decGetDigits(Unituint16_t *, Intint32_t);
250static Intint32_t	   decGetInt(const decNumber *);
251static decNumber * decLnOp(decNumber *, const decNumber *,
	      decContext *, uIntuint32_t *);
253static decNumber * decMultiplyOp(decNumber *, const decNumber *,
	      const decNumber *, decContext *,
	      uIntuint32_t *);
256static decNumber * decNaNs(decNumber *, const decNumber *,
	      const decNumber *, decContext *, uIntuint32_t *);
258static decNumber * decQuantizeOp(decNumber *, const decNumber *,
	      const decNumber *, decContext *, Flaguint8_t,
	      uIntuint32_t *);
261static void	   decReverse(Unituint16_t *, Unituint16_t *);
262static void	   decSetCoeff(decNumber *, decContext *, const Unituint16_t *,
	      Intint32_t, Intint32_t *, uIntuint32_t *);
264static void	   decSetMaxValue(decNumber *, decContext *);
265static void	   decSetOverflow(decNumber *, decContext *, uIntuint32_t *);
266static void	   decSetSubnormal(decNumber *, decContext *, Intint32_t *, uIntuint32_t *);
267static Intint32_t	   decShiftToLeast(Unituint16_t *, Intint32_t, Intint32_t);
268static Intint32_t	   decShiftToMost(Unituint16_t *, Intint32_t, Intint32_t);
269static void	   decStatus(decNumber *, uIntuint32_t, decContext *);
270static void	   decToString(const decNumber *, char[], Flaguint8_t);
271static decNumber * decTrim(decNumber *, decContext *, Flaguint8_t, Flaguint8_t, Intint32_t *);
272static Intint32_t	   decUnitAddSub(const Unituint16_t *, Intint32_t, const Unituint16_t *, Intint32_t, Intint32_t,
	      Unituint16_t *, Intint32_t);
274static Intint32_t	   decUnitCompare(const Unituint16_t *, Intint32_t, const Unituint16_t *, Intint32_t, Intint32_t);

276#if !DECSUBSET0
277/* decFinish == decFinalize when no subset arithmetic needed */
278#define decFinish(a,b,c,d)decFinalize(a,b,c,d) decFinalize(a,b,c,d)
279#else
280static void	   decFinish(decNumber *, decContext *, Int *, uInt *)decFinalize(decNumber *,decContext *,int32_t *,uint32_t *);
281static decNumber * decRoundOperand(const decNumber *, decContext *, uIntuint32_t *);
282#endif

284/* Local macros */
285/* masked special-values bits */
286#define SPECIALARG(rhs->bits & (0x40|0x20|0x10))  (rhs->bits & DECSPECIAL(0x40|0x20|0x10))
287#define SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10)) ((lhs->bits | rhs->bits) & DECSPECIAL(0x40|0x20|0x10))

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

308#if DECCHECK0
309/* Optional checking routines.	Enabling these means that decNumber */
310/* and decContext operands to operator routines are checked for */
311/* correctness.  This roughly doubles the execution time of the */
312/* fastest routines (and adds 600+ bytes), so should not normally be */
313/* used in 'production'. */
314/* decCheckInexact is used to check that inexact results have a full */
315/* complement of digits (where appropriate -- this is not the case */
316/* for Quantize, for example) */
317#define DECUNRESU ((decNumber *)(void *)0xffffffff)
318#define DECUNUSED ((const decNumber *)(void *)0xffffffff)
319#define DECUNCONT ((decContext *)(void *)(0xffffffff))
320static Flaguint8_t decCheckOperands(decNumber *, const decNumber *,
	     const decNumber *, decContext *);
322static Flaguint8_t decCheckNumber(const decNumber *);
323static void decCheckInexact(const decNumber *, decContext *);
324#endif

326#if DECTRACE0 || DECCHECK0
327/* Optional trace/debugging routines (may or may not be used) */
328void decNumberShow(const decNumber *);	/* displays the components of a number */
329static void decDumpAr(char, const Unituint16_t *, Intint32_t);
330#endif

332/* ================================================================== */
333/* Conversions							      */
334/* ================================================================== */

336/* ------------------------------------------------------------------ */
337/* from-int32 -- conversion from Int or uInt			      */
338/*								      */
339/*  dn is the decNumber to receive the integer			      */
340/*  in or uin is the integer to be converted			      */
341/*  returns dn							      */
342/*								      */
343/* No error is possible.					      */
344/* ------------------------------------------------------------------ */
345decNumber * decNumberFromInt32(decNumber *dn, Intint32_t in) {
uIntuint32_t unsig;
if (in>=0) unsig=in;
 else {				/* negative (possibly BADINT) */
  if (in==BADINT(int32_t)0x80000000) unsig=(uIntuint32_t)1073741824*2; /* special case */
   else unsig=-in;			/* invert */
  }
/* in is now positive */
decNumberFromUInt32(dn, unsig);
if (in<0) dn->bits=DECNEG0x80;		/* sign needed */
return dn;
} /* decNumberFromInt32 */

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

370/* ------------------------------------------------------------------ */
371/* to-int32 -- conversion to Int or uInt			      */
372/*								      */
373/*  dn is the decNumber to convert				      */
374/*  set is the context for reporting errors			      */
375/*  returns the converted decNumber, or 0 if Invalid is set	      */
376/*								      */
377/* Invalid is set if the decNumber does not have exponent==0 or if    */
378/* it is a NaN, Infinite, or out-of-range.			      */
379/* ------------------------------------------------------------------ */
380Intint32_t decNumberToInt32(const decNumber *dn, decContext *set) {
#if DECCHECK0
if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
#endif

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

416uIntuint32_t decNumberToUInt32(const decNumber *dn, decContext *set) {
#if DECCHECK0
if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
#endif
/* special or too many digits, or bad exponent, or negative (<0) */
if (dn->bits&DECSPECIAL(0x40|0x20|0x10) || dn->digits>10 || dn->exponent!=0
  || (dn->bits&DECNEG0x80 && !ISZERO(dn)(*(dn)->lsu==0 && (dn)->digits==1 && ((
(dn)->bits&(0x40|0x20|0x10))==0))));		    /* bad */
 else { /* is a finite integer with 10 or fewer digits */
  Intint32_t d;			   /* work */
  const Unituint16_t *up;		   /* .. */
  uIntuint32_t hi=0, lo;		   /* .. */
  up=dn->lsu; 		   /* -> lsu */
  lo=*up;			   /* get 1 to 9 digits */
  #if DECDPUN3>1		   /* split to higher */
    hi=lo/10;
    lo=lo%10;
  #endif
  up++;
  /* collect remaining Units, if any, into hi */
  for (d=DECDPUN3; d<dn->digits; up++, d+=DECDPUN3) hi+=*up*powersDECPOWERS[d-1];

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

445/* ------------------------------------------------------------------ */
446/* to-scientific-string -- conversion to numeric string 	      */
447/* to-engineering-string -- conversion to numeric string	      */
448/*								      */
449/*   decNumberToString(dn, string);				      */
450/*   decNumberToEngString(dn, string);				      */
451/*								      */
452/*  dn is the decNumber to convert				      */
453/*  string is the string where the result will be laid out	      */
454/*								      */
455/*  string must be at least dn->digits+14 characters long	      */
456/*								      */
457/*  No error is possible, and no status can be set.		      */
458/* ------------------------------------------------------------------ */
459char * decNumberToString(const decNumber *dn, char *string){
decToString(dn, string, 0);
return string;
} /* DecNumberToString */

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

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

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

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

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

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

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

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

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

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

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

  /* OK, the digits string is good.  Assemble in the decNumber, or in */
  /* a temporary units array if rounding is needed */
  if (d<=set->digits) res=dn->lsu;	/* fits into supplied decNumber */
   else {				/* rounding needed */
    Intint32_t needbytes=D2U(d)((d)<=49?d2utable[d]:((d)+3 -1)/3)*sizeof(Unituint16_t);/* bytes needed */
    res=resbuff;			/* assume use local buffer */
    if (needbytes>(Intint32_t)sizeof(resbuff)) { /* too big for local */
allocres=(Unituint16_t *)malloc(needbytes);
if (allocres==NULL((void*)0)) {status|=DEC_Insufficient_storage0x00000010; break;}
res=allocres;
}
    }
  /* res now -> number lsu, buffer, or allocated storage for Unit array */

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

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

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

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

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

719/* ================================================================== */
720/* Operators							      */
721/* ================================================================== */

723/* ------------------------------------------------------------------ */
724/* decNumberAbs -- absolute value operator			      */
725/*								      */
726/*   This computes C = abs(A)					      */
727/*								      */
728/*   res is C, the result.  C may be A				      */
729/*   rhs is A							      */
730/*   set is the context 					      */
731/*								      */
732/* See also decNumberCopyAbs for a quiet bitwise version of this.     */
733/* C must have space for set->digits digits.			      */
734/* ------------------------------------------------------------------ */
735/* This has the same effect as decNumberPlus unless A is negative,    */
736/* in which case it has the same effect as decNumberMinus.	      */
737/* ------------------------------------------------------------------ */
738decNumber * decNumberAbs(decNumber *res, const decNumber *rhs,
	 decContext *set) {
decNumber dzero;			/* for 0 */
uIntuint32_t status=0;			/* accumulator */

#if DECCHECK0
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, (uByteuint8_t)(rhs->bits & DECNEG0x80), &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberAbs */

757/* ------------------------------------------------------------------ */
758/* decNumberAdd -- add two Numbers				      */
759/*								      */
760/*   This computes C = A + B					      */
761/*								      */
762/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)	      */
763/*   lhs is A							      */
764/*   rhs is B							      */
765/*   set is the context 					      */
766/*								      */
767/* C must have space for set->digits digits.			      */
768/* ------------------------------------------------------------------ */
769/* This just calls the routine shared with Subtract		      */
770decNumber * decNumberAdd(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decAddOp(res, lhs, rhs, set, 0, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberAdd */

781/* ------------------------------------------------------------------ */
782/* decNumberAnd -- AND two Numbers, digitwise			      */
783/*								      */
784/*   This computes C = A & B					      */
785/*								      */
786/*   res is C, the result.  C may be A and/or B (e.g., X=X&X)	      */
787/*   lhs is A							      */
788/*   rhs is B							      */
789/*   set is the context (used for result length and error report)     */
790/*								      */
791/* C must have space for set->digits digits.			      */
792/*								      */
793/* Logical function restrictions apply (see above); a NaN is	      */
794/* returned with Invalid_operation if a restriction is violated.      */
795/* ------------------------------------------------------------------ */
796decNumber * decNumberAnd(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
const Unituint16_t *ua, *ub;			/* -> operands */
const Unituint16_t *msua, *msub;		/* -> operand msus */
Unituint16_t *uc,  *msuc;			/* -> result and its msu */
Intint32_t	msudigs;			/* digits in res msu */
#if DECCHECK0
if (decCheckOperands(res, lhs, rhs, set)) return res;
#endif

if (lhs->exponent!=0 || decNumberIsSpecial(lhs)(((lhs)->bits&(0x40|0x20|0x10))!=0) || decNumberIsNegative(lhs)(((lhs)->bits&0x80)!=0)
 || rhs->exponent!=0 || decNumberIsSpecial(rhs)(((rhs)->bits&(0x40|0x20|0x10))!=0) || decNumberIsNegative(rhs)(((rhs)->bits&0x80)!=0)) {
  decStatus(res, DEC_Invalid_operation0x00000080, set);
  return res;
  }

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

852/* ------------------------------------------------------------------ */
853/* decNumberCompare -- compare two Numbers			      */
854/*								      */
855/*   This computes C = A ? B					      */
856/*								      */
857/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
858/*   lhs is A							      */
859/*   rhs is B							      */
860/*   set is the context 					      */
861/*								      */
862/* C must have space for one digit (or NaN).			      */
863/* ------------------------------------------------------------------ */
864decNumber * decNumberCompare(decNumber *res, const decNumber *lhs,
	     const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPARE0x01, &status);
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberCompare */

872/* ------------------------------------------------------------------ */
873/* decNumberCompareSignal -- compare, signalling on all NaNs	      */
874/*								      */
875/*   This computes C = A ? B					      */
876/*								      */
877/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
878/*   lhs is A							      */
879/*   rhs is B							      */
880/*   set is the context 					      */
881/*								      */
882/* C must have space for one digit (or NaN).			      */
883/* ------------------------------------------------------------------ */
884decNumber * decNumberCompareSignal(decNumber *res, const decNumber *lhs,
		   const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPSIG0x06, &status);
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberCompareSignal */

892/* ------------------------------------------------------------------ */
893/* decNumberCompareTotal -- compare two Numbers, using total ordering */
894/*								      */
895/*   This computes C = A ? B, under total ordering		      */
896/*								      */
897/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
898/*   lhs is A							      */
899/*   rhs is B							      */
900/*   set is the context 					      */
901/*								      */
902/* C must have space for one digit; the result will always be one of  */
903/* -1, 0, or 1. 						      */
904/* ------------------------------------------------------------------ */
905decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs,
		  const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPTOTAL0x04, &status);
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberCompareTotal */

913/* ------------------------------------------------------------------ */
914/* decNumberCompareTotalMag -- compare, total ordering of magnitudes  */
915/*								      */
916/*   This computes C = |A| ? |B|, under total ordering		      */
917/*								      */
918/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
919/*   lhs is A							      */
920/*   rhs is B							      */
921/*   set is the context 					      */
922/*								      */
923/* C must have space for one digit; the result will always be one of  */
924/* -1, 0, or 1. 						      */
925/* ------------------------------------------------------------------ */
926decNumber * decNumberCompareTotalMag(decNumber *res, const decNumber *lhs,
		     const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;		   /* accumulator */
uIntuint32_t needbytes;		   /* for space calculations */
decNumber bufa[D2N(DECBUFFER+1)(((((((36 +1)+3 -1)/3)-1)*sizeof(uint16_t))+sizeof(decNumber)
*2-1)/sizeof(decNumber))];/* +1 in case DECBUFFER=0 */
decNumber *allocbufa=NULL((void*)0);	   /* -> allocated bufa, iff allocated */
decNumber bufb[D2N(DECBUFFER+1)(((((((36 +1)+3 -1)/3)-1)*sizeof(uint16_t))+sizeof(decNumber)
*2-1)/sizeof(decNumber))];
decNumber *allocbufb=NULL((void*)0);	   /* -> allocated bufb, iff allocated */
decNumber *a, *b;		   /* temporary pointers */

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

do {					/* protect allocated storage */
  /* if either is negative, take a copy and absolute */
  if (decNumberIsNegative(lhs)(((lhs)->bits&0x80)!=0)) {	/* lhs<0 */
    a=bufa;
    needbytes=sizeof(decNumber)+(D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3)-1)*sizeof(Unituint16_t);
    if (needbytes>sizeof(bufa)) {	/* need malloc space */
allocbufa=(decNumber *)malloc(needbytes);
if (allocbufa==NULL((void*)0)) {		/* hopeless -- abandon */
 status|=DEC_Insufficient_storage0x00000010;
 break;}
a=allocbufa;			/* use the allocated space */
}
    decNumberCopy(a, lhs);		/* copy content */
    a->bits&=~DECNEG0x80; 		/* .. and clear the sign */
    lhs=a;				/* use copy from here on */
    }
  if (decNumberIsNegative(rhs)(((rhs)->bits&0x80)!=0)) {	/* rhs<0 */
    b=bufb;
    needbytes=sizeof(decNumber)+(D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3)-1)*sizeof(Unituint16_t);
    if (needbytes>sizeof(bufb)) {	/* need malloc space */
allocbufb=(decNumber *)malloc(needbytes);
if (allocbufb==NULL((void*)0)) {		/* hopeless -- abandon */
 status|=DEC_Insufficient_storage0x00000010;
 break;}
b=allocbufb;			/* use the allocated space */
}
    decNumberCopy(b, rhs);		/* copy content */
    b->bits&=~DECNEG0x80; 		/* .. and clear the sign */
    rhs=b;				/* use copy from here on */
    }
  decCompareOp(res, lhs, rhs, set, COMPTOTAL0x04, &status);
  } while(0); 			/* end protected */

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

979/* ------------------------------------------------------------------ */
980/* decNumberDivide -- divide one number by another		      */
981/*								      */
982/*   This computes C = A / B					      */
983/*								      */
984/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)	      */
985/*   lhs is A							      */
986/*   rhs is B							      */
987/*   set is the context 					      */
988/*								      */
989/* C must have space for set->digits digits.			      */
990/* ------------------------------------------------------------------ */
991decNumber * decNumberDivide(decNumber *res, const decNumber *lhs,
	    const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decDivideOp(res, lhs, rhs, set, DIVIDE0x80, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberDivide */

1002/* ------------------------------------------------------------------ */
1003/* decNumberDivideInteger -- divide and return integer quotient       */
1004/*								      */
1005/*   This computes C = A # B, where # is the integer divide operator  */
1006/*								      */
1007/*   res is C, the result.  C may be A and/or B (e.g., X=X#X)	      */
1008/*   lhs is A							      */
1009/*   rhs is B							      */
1010/*   set is the context 					      */
1011/*								      */
1012/* C must have space for set->digits digits.			      */
1013/* ------------------------------------------------------------------ */
1014decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs,
		   const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decDivideOp(res, lhs, rhs, set, DIVIDEINT0x20, &status);
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberDivideInteger */

1022/* ------------------------------------------------------------------ */
1023/* decNumberExp -- exponentiation				      */
1024/*								      */
1025/*   This computes C = exp(A)					      */
1026/*								      */
1027/*   res is C, the result.  C may be A				      */
1028/*   rhs is A							      */
1029/*   set is the context; note that rounding mode has no effect	      */
1030/*								      */
1031/* C must have space for set->digits digits.			      */
1032/*								      */
1033/* Mathematical function restrictions apply (see above); a NaN is     */
1034/* returned with Invalid_operation if a restriction is violated.      */
1035/*								      */
1036/* Finite results will always be full precision and Inexact, except   */
1037/* when A is a zero or -Infinity (giving 1 or 0 respectively).	      */
1038/*								      */
1039/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1040/* almost always be correctly rounded, but may be up to 1 ulp in      */
1041/* error in rare cases. 					      */
1042/* ------------------------------------------------------------------ */
1043/* This is a wrapper for decExpOp which can handle the slightly wider */
1044/* (double) range needed by Ln (which has to be able to calculate     */
1045/* exp(-a) where a can be the tiniest number (Ntiny).		      */
1046/* ------------------------------------------------------------------ */
1047decNumber * decNumberExp(decNumber *res, const decNumber *rhs,
	 decContext *set) {
uIntuint32_t status=0;			/* accumulator */
#if DECSUBSET0
decNumber *allocrhs=NULL((void*)0);	   /* non-NULL if rounded rhs allocated */
#endif

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

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

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

1087/* ------------------------------------------------------------------ */
1088/* decNumberFMA -- fused multiply add				      */
1089/*								      */
1090/*   This computes D = (A * B) + C with only one rounding	      */
1091/*								      */
1092/*   res is D, the result.  D may be A or B or C (e.g., X=FMA(X,X,X)) */
1093/*   lhs is A							      */
1094/*   rhs is B							      */
1095/*   fhs is C [far hand side]					      */
1096/*   set is the context 					      */
1097/*								      */
1098/* Mathematical function restrictions apply (see above); a NaN is     */
1099/* returned with Invalid_operation if a restriction is violated.      */
1100/*								      */
1101/* C must have space for set->digits digits.			      */
1102/* ------------------------------------------------------------------ */
1103decNumber * decNumberFMA(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, const decNumber *fhs,
	 decContext *set) {
uIntuint32_t status=0;		   /* accumulator */
decContext dcmul;		   /* context for the multiplication */
uIntuint32_t needbytes;		   /* for space calculations */
decNumber bufa[D2N(DECBUFFER*2+1)(((((((36*2+1)+3 -1)/3)-1)*sizeof(uint16_t))+sizeof(decNumber
)*2-1)/sizeof(decNumber))];
decNumber *allocbufa=NULL((void*)0);	   /* -> allocated bufa, iff allocated */
decNumber *acc;		   /* accumulator pointer */
decNumber dzero;		   /* work */

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

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

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

1180/* ------------------------------------------------------------------ */
1181/* decNumberInvert -- invert a Number, digitwise		      */
1182/*								      */
1183/*   This computes C = ~A					      */
1184/*								      */
1185/*   res is C, the result.  C may be A (e.g., X=~X)		      */
1186/*   rhs is A							      */
1187/*   set is the context (used for result length and error report)     */
1188/*								      */
1189/* C must have space for set->digits digits.			      */
1190/*								      */
1191/* Logical function restrictions apply (see above); a NaN is	      */
1192/* returned with Invalid_operation if a restriction is violated.      */
1193/* ------------------------------------------------------------------ */
1194decNumber * decNumberInvert(decNumber *res, const decNumber *rhs,
	    decContext *set) {
const Unituint16_t *ua, *msua;		/* -> operand and its msu */
Unituint16_t	*uc, *msuc;			/* -> result and its msu */
Intint32_t	msudigs;			/* digits in res msu */
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
#endif

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

1239/* ------------------------------------------------------------------ */
1240/* decNumberLn -- natural logarithm				      */
1241/*								      */
1242/*   This computes C = ln(A)					      */
1243/*								      */
1244/*   res is C, the result.  C may be A				      */
1245/*   rhs is A							      */
1246/*   set is the context; note that rounding mode has no effect	      */
1247/*								      */
1248/* C must have space for set->digits digits.			      */
1249/*								      */
1250/* Notable cases:						      */
1251/*   A<0 -> Invalid						      */
1252/*   A=0 -> -Infinity (Exact)					      */
1253/*   A=+Infinity -> +Infinity (Exact)				      */
1254/*   A=1 exactly -> 0 (Exact)					      */
1255/*								      */
1256/* Mathematical function restrictions apply (see above); a NaN is     */
1257/* returned with Invalid_operation if a restriction is violated.      */
1258/*								      */
1259/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1260/* almost always be correctly rounded, but may be up to 1 ulp in      */
1261/* error in rare cases. 					      */
1262/* ------------------------------------------------------------------ */
1263/* This is a wrapper for decLnOp which can handle the slightly wider  */
1264/* (+11) range needed by Ln, Log10, etc. (which may have to be able   */
1265/* to calculate at p+e+2).					      */
1266/* ------------------------------------------------------------------ */
1267decNumber * decNumberLn(decNumber *res, const decNumber *rhs,
	decContext *set) {
uIntuint32_t status=0;		   /* accumulator */
#if DECSUBSET0
decNumber *allocrhs=NULL((void*)0);	   /* non-NULL if rounded rhs allocated */
#endif

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

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

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

1309/* ------------------------------------------------------------------ */
1310/* decNumberLogB - get adjusted exponent, by 754 rules		      */
1311/*								      */
1312/*   This computes C = adjustedexponent(A)			      */
1313/*								      */
1314/*   res is C, the result.  C may be A				      */
1315/*   rhs is A							      */
1316/*   set is the context, used only for digits and status	      */
1317/*								      */
1318/* C must have space for 10 digits (A might have 10**9 digits and     */
1319/* an exponent of +999999999, or one digit and an exponent of	      */
1320/* -1999999999).						      */
1321/*								      */
1322/* This returns the adjusted exponent of A after (in theory) padding  */
1323/* with zeros on the right to set->digits digits while keeping the    */
1324/* same value.	The exponent is not limited by emin/emax.	      */
1325/*								      */
1326/* Notable cases:						      */
1327/*   A<0 -> Use |A|						      */
1328/*   A=0 -> -Infinity (Division by zero)			      */
1329/*   A=Infinite -> +Infinity (Exact)				      */
1330/*   A=1 exactly -> 0 (Exact)					      */
1331/*   NaNs are propagated as usual				      */
1332/* ------------------------------------------------------------------ */
1333decNumber * decNumberLogB(decNumber *res, const decNumber *rhs,
	  decContext *set) {
uIntuint32_t status=0;		   /* accumulator */

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

/* NaNs as usual; Infinities return +Infinity; 0->oops */
if (decNumberIsNaN(rhs)(((rhs)->bits&(0x20|0x10))!=0)) decNaNs(res, rhs, NULL((void*)0), set, &status);
 else if (decNumberIsInfinite(rhs)(((rhs)->bits&0x40)!=0)) decNumberCopyAbs(res, rhs);
 else if (decNumberIsZero(rhs)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0))) {
  decNumberZero(res); 		/* prepare for Infinity */
  res->bits=DECNEG0x80|DECINF0x40;		/* -Infinity */
  status|=DEC_Division_by_zero0x00000002;	/* as per 754 */
  }
 else { /* finite non-zero */
  Intint32_t ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */
  decNumberFromInt32(res, ae);	/* lay it out */
  }

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

1358/* ------------------------------------------------------------------ */
1359/* decNumberLog10 -- logarithm in base 10			      */
1360/*								      */
1361/*   This computes C = log10(A) 				      */
1362/*								      */
1363/*   res is C, the result.  C may be A				      */
1364/*   rhs is A							      */
1365/*   set is the context; note that rounding mode has no effect	      */
1366/*								      */
1367/* C must have space for set->digits digits.			      */
1368/*								      */
1369/* Notable cases:						      */
1370/*   A<0 -> Invalid						      */
1371/*   A=0 -> -Infinity (Exact)					      */
1372/*   A=+Infinity -> +Infinity (Exact)				      */
1373/*   A=10**n (if n is an integer) -> n (Exact)			      */
1374/*								      */
1375/* Mathematical function restrictions apply (see above); a NaN is     */
1376/* returned with Invalid_operation if a restriction is violated.      */
1377/*								      */
1378/* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will    */
1379/* almost always be correctly rounded, but may be up to 1 ulp in      */
1380/* error in rare cases. 					      */
1381/* ------------------------------------------------------------------ */
1382/* This calculates ln(A)/ln(10) using appropriate precision.  For     */
1383/* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the      */
1384/* requested digits and t is the number of digits in the exponent     */
1385/* (maximum 6).  For ln(10) it is p + 3; this is often handled by the */
1386/* fastpath in decLnOp.  The final division is done to the requested  */
1387/* precision.							      */
1388/* ------------------------------------------------------------------ */
1389decNumber * decNumberLog10(decNumber *res, const decNumber *rhs,
	  decContext *set) {
uIntuint32_t status=0, ignore=0;	   /* status accumulators */
uIntuint32_t needbytes;		   /* for space calculations */
Intint32_t p;			   /* working precision */
Intint32_t t;			   /* digits in exponent of A */

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

decContext aset;		   /* working context */

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

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

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

  /* handle exact powers of 10; only check if +ve finite */
  if (!(rhs->bits&(DECNEG0x80|DECSPECIAL(0x40|0x20|0x10))) && !ISZERO(rhs)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0))) {
    Intint32_t residue=0;		   /* (no residue) */
    uIntuint32_t copystat=0;		   /* clean status */

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

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

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

  /* skip the division if the result so far is infinite, NaN, or */
  /* zero, or there was an error; note NaN from sNaN needs copy */
  if (status&DEC_NaNs(0x00000001 | 0x00000004 | 0x00000008 | 0x00000010 | 0x00000040
 | 0x00000080) && !(status&DEC_sNaN0x40000000)) break;
  if (a->bits&DECSPECIAL(0x40|0x20|0x10) || ISZERO(a)(*(a)->lsu==0 && (a)->digits==1 && (((a
)->bits&(0x40|0x20|0x10))==0))) {
    decNumberCopy(res, a);		/* [will fit] */
    break;}

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

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

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

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

1525/* ------------------------------------------------------------------ */
1526/* decNumberMax -- compare two Numbers and return the maximum	      */
1527/*								      */
1528/*   This computes C = A ? B, returning the maximum by 754 rules      */
1529/*								      */
1530/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
1531/*   lhs is A							      */
1532/*   rhs is B							      */
1533/*   set is the context 					      */
1534/*								      */
1535/* C must have space for set->digits digits.			      */
1536/* ------------------------------------------------------------------ */
1537decNumber * decNumberMax(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPMAX0x02, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMax */

1548/* ------------------------------------------------------------------ */
1549/* decNumberMaxMag -- compare and return the maximum by magnitude     */
1550/*								      */
1551/*   This computes C = A ? B, returning the maximum by 754 rules      */
1552/*								      */
1553/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
1554/*   lhs is A							      */
1555/*   rhs is B							      */
1556/*   set is the context 					      */
1557/*								      */
1558/* C must have space for set->digits digits.			      */
1559/* ------------------------------------------------------------------ */
1560decNumber * decNumberMaxMag(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPMAXMAG0x07, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMaxMag */

1571/* ------------------------------------------------------------------ */
1572/* decNumberMin -- compare two Numbers and return the minimum	      */
1573/*								      */
1574/*   This computes C = A ? B, returning the minimum by 754 rules      */
1575/*								      */
1576/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
1577/*   lhs is A							      */
1578/*   rhs is B							      */
1579/*   set is the context 					      */
1580/*								      */
1581/* C must have space for set->digits digits.			      */
1582/* ------------------------------------------------------------------ */
1583decNumber * decNumberMin(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPMIN0x03, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMin */

1594/* ------------------------------------------------------------------ */
1595/* decNumberMinMag -- compare and return the minimum by magnitude     */
1596/*								      */
1597/*   This computes C = A ? B, returning the minimum by 754 rules      */
1598/*								      */
1599/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)	      */
1600/*   lhs is A							      */
1601/*   rhs is B							      */
1602/*   set is the context 					      */
1603/*								      */
1604/* C must have space for set->digits digits.			      */
1605/* ------------------------------------------------------------------ */
1606decNumber * decNumberMinMag(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decCompareOp(res, lhs, rhs, set, COMPMINMAG0x08, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMinMag */

1617/* ------------------------------------------------------------------ */
1618/* decNumberMinus -- prefix minus operator			      */
1619/*								      */
1620/*   This computes C = 0 - A					      */
1621/*								      */
1622/*   res is C, the result.  C may be A				      */
1623/*   rhs is A							      */
1624/*   set is the context 					      */
1625/*								      */
1626/* See also decNumberCopyNegate for a quiet bitwise version of this.  */
1627/* C must have space for set->digits digits.			      */
1628/* ------------------------------------------------------------------ */
1629/* Simply use AddOp for the subtract, which will do the necessary.    */
1630/* ------------------------------------------------------------------ */
1631decNumber * decNumberMinus(decNumber *res, const decNumber *rhs,
	   decContext *set) {
decNumber dzero;
uIntuint32_t status=0;			/* accumulator */

#if DECCHECK0
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, DECNEG0x80, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMinus */

1650/* ------------------------------------------------------------------ */
1651/* decNumberNextMinus -- next towards -Infinity 		      */
1652/*								      */
1653/*   This computes C = A - infinitesimal, rounded towards -Infinity   */
1654/*								      */
1655/*   res is C, the result.  C may be A				      */
1656/*   rhs is A							      */
1657/*   set is the context 					      */
1658/*								      */
1659/* This is a generalization of 754 NextDown.			      */
1660/* ------------------------------------------------------------------ */
1661decNumber * decNumberNextMinus(decNumber *res, const decNumber *rhs,
	       decContext *set) {
decNumber dtiny;			     /* constant */
decContext workset=*set;		     /* work */
uIntuint32_t status=0;			     /* accumulator */
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
#endif

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

1686/* ------------------------------------------------------------------ */
1687/* decNumberNextPlus -- next towards +Infinity			      */
1688/*								      */
1689/*   This computes C = A + infinitesimal, rounded towards +Infinity   */
1690/*								      */
1691/*   res is C, the result.  C may be A				      */
1692/*   rhs is A							      */
1693/*   set is the context 					      */
1694/*								      */
1695/* This is a generalization of 754 NextUp.			      */
1696/* ------------------------------------------------------------------ */
1697decNumber * decNumberNextPlus(decNumber *res, const decNumber *rhs,
	      decContext *set) {
decNumber dtiny;			     /* constant */
decContext workset=*set;		     /* work */
uIntuint32_t status=0;			     /* accumulator */
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
#endif

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

1723/* ------------------------------------------------------------------ */
1724/* decNumberNextToward -- next towards rhs			      */
1725/*								      */
1726/*   This computes C = A +/- infinitesimal, rounded towards	      */
1727/*   +/-Infinity in the direction of B, as per 754-1985 nextafter     */
1728/*   modified during revision but dropped from 754-2008.	      */
1729/*								      */
1730/*   res is C, the result.  C may be A or B.			      */
1731/*   lhs is A							      */
1732/*   rhs is B							      */
1733/*   set is the context 					      */
1734/*								      */
1735/* This is a generalization of 754-1985 NextAfter.		      */
1736/* ------------------------------------------------------------------ */
1737decNumber * decNumberNextToward(decNumber *res, const decNumber *lhs,
		const decNumber *rhs, decContext *set) {
decNumber dtiny;			     /* constant */
decContext workset=*set;		     /* work */
Intint32_t result;				     /* .. */
uIntuint32_t status=0;			     /* accumulator */
#if DECCHECK0
if (decCheckOperands(res, lhs, rhs, set)) return res;
#endif

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

1790/* ------------------------------------------------------------------ */
1791/* decNumberOr -- OR two Numbers, digitwise			      */
1792/*								      */
1793/*   This computes C = A | B					      */
1794/*								      */
1795/*   res is C, the result.  C may be A and/or B (e.g., X=X|X)	      */
1796/*   lhs is A							      */
1797/*   rhs is B							      */
1798/*   set is the context (used for result length and error report)     */
1799/*								      */
1800/* C must have space for set->digits digits.			      */
1801/*								      */
1802/* Logical function restrictions apply (see above); a NaN is	      */
1803/* returned with Invalid_operation if a restriction is violated.      */
1804/* ------------------------------------------------------------------ */
1805decNumber * decNumberOr(decNumber *res, const decNumber *lhs,
	const decNumber *rhs, decContext *set) {
const Unituint16_t *ua, *ub;			/* -> operands */
const Unituint16_t *msua, *msub;		/* -> operand msus */
Unituint16_t	*uc, *msuc;			/* -> result and its msu */
Intint32_t	msudigs;			/* digits in res msu */
#if DECCHECK0
if (decCheckOperands(res, lhs, rhs, set)) return res;
#endif

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

1859/* ------------------------------------------------------------------ */
1860/* decNumberPlus -- prefix plus operator			      */
1861/*								      */
1862/*   This computes C = 0 + A					      */
1863/*								      */
1864/*   res is C, the result.  C may be A				      */
1865/*   rhs is A							      */
1866/*   set is the context 					      */
1867/*								      */
1868/* See also decNumberCopy for a quiet bitwise version of this.	      */
1869/* C must have space for set->digits digits.			      */
1870/* ------------------------------------------------------------------ */
1871/* This simply uses AddOp; Add will take fast path after preparing A. */
1872/* Performance is a concern here, as this routine is often used to    */
1873/* check operands and apply rounding and overflow/underflow testing.  */
1874/* ------------------------------------------------------------------ */
1875decNumber * decNumberPlus(decNumber *res, const decNumber *rhs,
	  decContext *set) {
decNumber dzero;
uIntuint32_t status=0;			/* accumulator */
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
#endif

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

1893/* ------------------------------------------------------------------ */
1894/* decNumberMultiply -- multiply two Numbers			      */
1895/*								      */
1896/*   This computes C = A x B					      */
1897/*								      */
1898/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)	      */
1899/*   lhs is A							      */
1900/*   rhs is B							      */
1901/*   set is the context 					      */
1902/*								      */
1903/* C must have space for set->digits digits.			      */
1904/* ------------------------------------------------------------------ */
1905decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs,
	      const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;		   /* accumulator */
decMultiplyOp(res, lhs, rhs, set, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberMultiply */

1916/* ------------------------------------------------------------------ */
1917/* decNumberPower -- raise a number to a power			      */
1918/*								      */
1919/*   This computes C = A ** B					      */
1920/*								      */
1921/*   res is C, the result.  C may be A and/or B (e.g., X=X**X)	      */
1922/*   lhs is A							      */
1923/*   rhs is B							      */
1924/*   set is the context 					      */
1925/*								      */
1926/* C must have space for set->digits digits.			      */
1927/*								      */
1928/* Mathematical function restrictions apply (see above); a NaN is     */
1929/* returned with Invalid_operation if a restriction is violated.      */
1930/*								      */
1931/* However, if 1999999997<=B<=999999999 and B is an integer then the  */
1932/* restrictions on A and the context are relaxed to the usual bounds, */
1933/* for compatibility with the earlier (integer power only) version    */
1934/* of this function.						      */
1935/*								      */
1936/* When B is an integer, the result may be exact, even if rounded.    */
1937/*								      */
1938/* The final result is rounded according to the context; it will      */
1939/* almost always be correctly rounded, but may be up to 1 ulp in      */
1940/* error in rare cases. 					      */
1941/* ------------------------------------------------------------------ */
1942decNumber * decNumberPower(decNumber *res, const decNumber *lhs,
	   const decNumber *rhs, decContext *set) {
#if DECSUBSET0
decNumber *alloclhs=NULL((void*)0);	   /* non-NULL if rounded lhs allocated */
decNumber *allocrhs=NULL((void*)0);	   /* .., rhs */
#endif
decNumber *allocdac=NULL((void*)0);	   /* -> allocated acc buffer, iff used */
decNumber *allocinv=NULL((void*)0);	   /* -> allocated 1/x buffer, iff used */
Intint32_t	reqdigits=set->digits;	   /* requested DIGITS */
Intint32_t	n;			   /* rhs in binary */
Flaguint8_t	rhsint=0;		   /* 1 if rhs is an integer */
Flaguint8_t	useint=0;		   /* 1 if can use integer calculation */
Flaguint8_t	isoddint=0;		   /* 1 if rhs is an integer and odd */
Intint32_t	i;			   /* work */
#if DECSUBSET0
Intint32_t	dropped;		   /* .. */
#endif
uIntuint32_t	needbytes;		   /* buffer size needed */
Flaguint8_t	seenbit;		   /* seen a bit while powering */
Intint32_t	residue=0;		   /* rounding residue */
uIntuint32_t	status=0;		   /* accumulators */
uByteuint8_t bits=0; 		   /* result sign if errors */
decContext aset;		   /* working context */
decNumber dnOne;		   /* work value 1... */
/* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
decNumber dacbuff[D2N(DECBUFFER+9)(((((((36 +9)+3 -1)/3)-1)*sizeof(uint16_t))+sizeof(decNumber)
*2-1)/sizeof(decNumber))];
decNumber *dac=dacbuff;	   /* -> result accumulator */
/* same again for possible 1/lhs calculation */
decNumber invbuff[D2N(DECBUFFER+9)(((((((36 +9)+3 -1)/3)-1)*sizeof(uint16_t))+sizeof(decNumber)
*2-1)/sizeof(decNumber))];

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

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

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

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

  if (decNumberIsNegative(lhs)(((lhs)->bits&0x80)!=0)	/* -x .. */
    && isoddint) bits=DECNEG0x80; 	/* .. to an odd power */

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

  /* similarly handle LHS zero */
  if (decNumberIsZero(lhs)(*(lhs)->lsu==0 && (lhs)->digits==1 && (
((lhs)->bits&(0x40|0x20|0x10))==0))) {
    if (n==0) {			     /* 0**0 => Error */
#if DECSUBSET0
if (!set->extended) {		     /* [unless subset] */
 decNumberZero(res);
 *res->lsu=1;			     /* return 1 */
 break;}
#endif
status|=DEC_Invalid_operation0x00000080;
}
     else {				     /* 0**x */
uByteuint8_t rbits=rhs->bits;		     /* save */
if (rbits & DECNEG0x80) {		     /* was a 0**(-n) */
 #if DECSUBSET0
 if (!set->extended) { 	     /* [bad if subset] */
   status|=DEC_Invalid_operation0x00000080;
   break;}
 #endif
 bits|=DECINF0x40;
 }
decNumberZero(res);		     /* prepare */
/* [otherwise will be 0 or -0] */
res->bits=bits;
}
    break;}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2295/* ------------------------------------------------------------------ */
2296/* decNumberReduce -- remove trailing zeros			      */
2297/*								      */
2298/*   This computes C = 0 + A, and normalizes the result 	      */
2299/*								      */
2300/*   res is C, the result.  C may be A				      */
2301/*   rhs is A							      */
2302/*   set is the context 					      */
2303/*								      */
2304/* C must have space for set->digits digits.			      */
2305/* ------------------------------------------------------------------ */
2306/* Previously known as Normalize */
2307decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs,
	       decContext *set) {
return decNumberReduce(res, rhs, set);
} /* decNumberNormalize */

2312decNumber * decNumberReduce(decNumber *res, const decNumber *rhs,
	    decContext *set) {
#if DECSUBSET0
decNumber *allocrhs=NULL((void*)0);	   /* non-NULL if rounded rhs allocated */
#endif
uIntuint32_t status=0;		   /* as usual */
Intint32_t  residue=0;		   /* as usual */
Intint32_t  dropped; 		   /* work */

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

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

  /* Infinities copy through; NaNs need usual treatment */
  if (decNumberIsNaN(rhs)(((rhs)->bits&(0x20|0x10))!=0)) {
    decNaNs(res, rhs, NULL((void*)0), set, &status);
    break;
    }

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

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

2358/* ------------------------------------------------------------------ */
2359/* decNumberRescale -- force exponent to requested value	      */
2360/*								      */
2361/*   This computes C = op(A, B), where op adjusts the coefficient     */
2362/*   of C (by rounding or shifting) such that the exponent (-scale)   */
2363/*   of C has the value B.  The numerical value of C will equal A,    */
2364/*   except for the effects of any rounding that occurred.	      */
2365/*								      */
2366/*   res is C, the result.  C may be A or B			      */
2367/*   lhs is A, the number to adjust				      */
2368/*   rhs is B, the requested exponent				      */
2369/*   set is the context 					      */
2370/*								      */
2371/* C must have space for set->digits digits.			      */
2372/*								      */
2373/* Unless there is an error or the result is infinite, the exponent   */
2374/* after the operation is guaranteed to be equal to B.		      */
2375/* ------------------------------------------------------------------ */
2376decNumber * decNumberRescale(decNumber *res, const decNumber *lhs,
	     const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decQuantizeOp(res, lhs, rhs, set, 0, &status);
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberRescale */

2384/* ------------------------------------------------------------------ */
2385/* decNumberRemainder -- divide and return remainder		      */
2386/*								      */
2387/*   This computes C = A % B					      */
2388/*								      */
2389/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)	      */
2390/*   lhs is A							      */
2391/*   rhs is B							      */
2392/*   set is the context 					      */
2393/*								      */
2394/* C must have space for set->digits digits.			      */
2395/* ------------------------------------------------------------------ */
2396decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs,
	       const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decDivideOp(res, lhs, rhs, set, REMAINDER0x40, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberRemainder */

2407/* ------------------------------------------------------------------ */
2408/* decNumberRemainderNear -- divide and return remainder from nearest */
2409/*								      */
2410/*   This computes C = A % B, where % is the IEEE remainder operator  */
2411/*								      */
2412/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)	      */
2413/*   lhs is A							      */
2414/*   rhs is B							      */
2415/*   set is the context 					      */
2416/*								      */
2417/* C must have space for set->digits digits.			      */
2418/* ------------------------------------------------------------------ */
2419decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs,
		   const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */
decDivideOp(res, lhs, rhs, set, REMNEAR0x10, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberRemainderNear */

2430/* ------------------------------------------------------------------ */
2431/* decNumberRotate -- rotate the coefficient of a Number left/right   */
2432/*								      */
2433/*   This computes C = A rot B	(in base ten and rotating set->digits */
2434/*   digits).							      */
2435/*								      */
2436/*   res is C, the result.  C may be A and/or B (e.g., X=XrotX)       */
2437/*   lhs is A							      */
2438/*   rhs is B, the number of digits to rotate (-ve to right)	      */
2439/*   set is the context 					      */
2440/*								      */
2441/* The digits of the coefficient of A are rotated to the left (if B   */
2442/* is positive) or to the right (if B is negative) without adjusting  */
2443/* the exponent or the sign of A.  If lhs->digits is less than	      */
2444/* set->digits the coefficient is padded with zeros on the left       */
2445/* before the rotate.  Any leading zeros in the result are removed    */
2446/* as usual.							      */
2447/*								      */
2448/* B must be an integer (q=0) and in the range -set->digits through   */
2449/* +set->digits.						      */
2450/* C must have space for set->digits digits.			      */
2451/* NaNs are propagated as usual.  Infinities are unaffected (but      */
2452/* B must be valid).  No status is set unless B is invalid or an      */
2453/* operand is an sNaN.						      */
2454/* ------------------------------------------------------------------ */
2455decNumber * decNumberRotate(decNumber *res, const decNumber *lhs,
	   const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;	      /* accumulator */
Intint32_t  rotate;		      /* rhs as an Int */

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

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

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

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

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

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

2567/* ------------------------------------------------------------------ */
2568/* decNumberSameQuantum -- test for equal exponents		      */
2569/*								      */
2570/*   res is the result number, which will contain either 0 or 1       */
2571/*   lhs is a number to test					      */
2572/*   rhs is the second (usually a pattern)			      */
2573/*								      */
2574/* No errors are possible and no context is needed.		      */
2575/* ------------------------------------------------------------------ */
2576decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs,
		 const decNumber *rhs) {
Unituint16_t ret=0;			   /* return value */

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

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

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

2596/* ------------------------------------------------------------------ */
2597/* decNumberScaleB -- multiply by a power of 10 		      */
2598/*								      */
2599/* This computes C = A x 10**B where B is an integer (q=0) with       */
2600/* maximum magnitude 2*(emax+digits)				      */
2601/*								      */
2602/*   res is C, the result.  C may be A or B			      */
2603/*   lhs is A, the number to adjust				      */
2604/*   rhs is B, the requested power of ten to use		      */
2605/*   set is the context 					      */
2606/*								      */
2607/* C must have space for set->digits digits.			      */
2608/*								      */
2609/* The result may underflow or overflow.			      */
2610/* ------------------------------------------------------------------ */
2611decNumber * decNumberScaleB(decNumber *res, const decNumber *lhs,
	    const decNumber *rhs, decContext *set) {
Intint32_t  reqexp;		      /* requested exponent change [B] */
uIntuint32_t status=0;	      /* accumulator */
Intint32_t  residue; 	      /* work */

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

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

2647/* ------------------------------------------------------------------ */
2648/* decNumberShift -- shift the coefficient of a Number left or right  */
2649/*								      */
2650/*   This computes C = A << B or C = A >> -B  (in base ten).	      */
2651/*								      */
2652/*   res is C, the result.  C may be A and/or B (e.g., X=X<<X)	      */
2653/*   lhs is A							      */
2654/*   rhs is B, the number of digits to shift (-ve to right)	      */
2655/*   set is the context 					      */
2656/*								      */
2657/* The digits of the coefficient of A are shifted to the left (if B   */
2658/* is positive) or to the right (if B is negative) without adjusting  */
2659/* the exponent or the sign of A.				      */
2660/*								      */
2661/* B must be an integer (q=0) and in the range -set->digits through   */
2662/* +set->digits.						      */
2663/* C must have space for set->digits digits.			      */
2664/* NaNs are propagated as usual.  Infinities are unaffected (but      */
2665/* B must be valid).  No status is set unless B is invalid or an      */
2666/* operand is an sNaN.						      */
2667/* ------------------------------------------------------------------ */
2668decNumber * decNumberShift(decNumber *res, const decNumber *lhs,
	   const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;	      /* accumulator */
Intint32_t  shift;		      /* rhs as an Int */

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

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

2725/* ------------------------------------------------------------------ */
2726/* decNumberSquareRoot -- square root operator			      */
2727/*								      */
2728/*   This computes C = squareroot(A)				      */
2729/*								      */
2730/*   res is C, the result.  C may be A				      */
2731/*   rhs is A							      */
2732/*   set is the context; note that rounding mode has no effect	      */
2733/*								      */
2734/* C must have space for set->digits digits.			      */
2735/* ------------------------------------------------------------------ */
2736/* This uses the following varying-precision algorithm in:	      */
2737/*								      */
2738/*   Properly Rounded Variable Precision Square Root, T. E. Hull and  */
2739/*   A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
2740/*   pp229-237, ACM, September 1985.				      */
2741/*								      */
2742/* The square-root is calculated using Newton's method, after which   */
2743/* a check is made to ensure the result is correctly rounded.	      */
2744/*								      */
2745/* % [Reformatted original Numerical Turing source code follows.]     */
2746/* function sqrt(x : real) : real				      */
2747/* % sqrt(x) returns the properly rounded approximation to the square */
2748/* % root of x, in the precision of the calling environment, or it    */
2749/* % fails if x < 0.						      */
2750/* % t e hull and a abrham, august, 1984			      */
2751/* if x <= 0 then						      */
2752/*   if x < 0 then						      */
2753/*     assert false						      */
2754/*   else							      */
2755/*     result 0 						      */
2756/*   end if							      */
2757/* end if							      */
2758/* var f := setexp(x, 0)  % fraction part of x	 [0.1 <= x < 1]       */
2759/* var e := getexp(x)	  % exponent part of x			      */
2760/* var approx : real						      */
2761/* if e mod 2 = 0  then 					      */
2762/*   approx := .259 + .819 * f	 % approx to root of f		      */
2763/* else 							      */
2764/*   f := f/l0			 % adjustments			      */
2765/*   e := e + 1 		 %   for odd			      */
2766/*   approx := .0819 + 2.59 * f  %   exponent			      */
2767/* end if							      */
2768/*								      */
2769/* var p:= 3							      */
2770/* const maxp := currentprecision + 2				      */
2771/* loop 							      */
2772/*   p := min(2*p - 2, maxp)	 % p = 4,6,10, . . . , maxp	      */
2773/*   precision p						      */
2774/*   approx := .5 * (approx + f/approx) 			      */
2775/*   exit when p = maxp 					      */
2776/* end loop							      */
2777/*								      */
2778/* % approx is now within 1 ulp of the properly rounded square root   */
2779/* % of f; to ensure proper rounding, compare squares of (approx -    */
2780/* % l/2 ulp) and (approx + l/2 ulp) with f.			      */
2781/* p := currentprecision					      */
2782/* begin							      */
2783/*   precision p + 2						      */
2784/*   const approxsubhalf := approx - setexp(.5, -p)		      */
2785/*   if mulru(approxsubhalf, approxsubhalf) > f then		      */
2786/*     approx := approx - setexp(.l, -p + 1)			      */
2787/*   else							      */
2788/*     const approxaddhalf := approx + setexp(.5, -p)		      */
2789/*     if mulrd(approxaddhalf, approxaddhalf) < f then		      */
2790/*	 approx := approx + setexp(.l, -p + 1)			      */
2791/*     end if							      */
2792/*   end if							      */
2793/* end								      */
2794/* result setexp(approx, e div 2)  % fix exponent		      */
2795/* end sqrt							      */
2796/* ------------------------------------------------------------------ */
2797decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs,
		decContext *set) {
decContext workset, approxset;   /* work contexts */
decNumber dzero;		   /* used for constant zero */
Intint32_t  maxp;			   /* largest working precision */
Intint32_t  workp;			   /* working precision */
Intint32_t  residue=0;		   /* rounding residue */
uIntuint32_t status=0, ignore=0;	   /* status accumulators */
uIntuint32_t rstatus; 		   /* .. */
Intint32_t  exp;			   /* working exponent */
Intint32_t  ideal;			   /* ideal (preferred) exponent */
Intint32_t  needbytes;		   /* work */
Intint32_t  dropped; 		   /* .. */

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

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

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

  /* handle infinities and NaNs */
  if (SPECIALARG(rhs->bits & (0x40|0x20|0x10))) {
    if (decNumberIsInfinite(rhs)(((rhs)->bits&0x40)!=0)) {	      /* an infinity */
if (decNumberIsNegative(rhs)(((rhs)->bits&0x80)!=0)) status|=DEC_Invalid_operation0x00000080;
else decNumberCopy(res, rhs);	      /* +Infinity */
}
     else decNaNs(res, rhs, NULL((void*)0), set, &status); /* a NaN */
    break;
    }

  /* calculate the ideal (preferred) exponent [floor(exp/2)] */
  /* [It would be nicer 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)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0))) {
    decNumberCopy(res, rhs);		/* could be 0 or -0 */
    res->exponent=ideal;		/* use the ideal [safe] */
    /* use decFinish to clamp any out-of-range exponent, etc. */
    decFinish(res, set, &residue, &status)decFinalize(res,set,&residue,&status);
    break;
    }

  /* any other -x is an oops */
  if (decNumberIsNegative(rhs)(((rhs)->bits&0x80)!=0)) {
    status|=DEC_Invalid_operation0x00000080;
    break;
    }

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

  needbytes=sizeof(decNumber)+(D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3)-1)*sizeof(Unituint16_t);
  if (needbytes>(Intint32_t)sizeof(buff)) {
    allocbuff=(decNumber *)malloc(needbytes);
    if (allocbuff==NULL((void*)0)) {  /* hopeless -- abandon */
status|=DEC_Insufficient_storage0x00000010;
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)((maxp)<=49?d2utable[maxp]:((maxp)+3 -1)/3)-1)*sizeof(Unituint16_t);
  if (needbytes>(Intint32_t)sizeof(bufa)) {		  /* [same applies to b] */
    allocbufa=(decNumber *)malloc(needbytes);
    allocbufb=(decNumber *)malloc(needbytes);
    if (allocbufa==NULL((void*)0) || allocbufb==NULL((void*)0)) {   /* hopeless */
status|=DEC_Insufficient_storage0x00000010;
break;}
    a=allocbufa;	      /* use the allocated spaces */
    b=allocbufb;	      /* .. */
    }

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

  /* set up working context */
  decContextDefault(&workset, DEC_INIT_DECIMAL6464);
  workset.emax=DEC_MAX_EMAX999999999;
  workset.emin=DEC_MIN_EMIN-999999999;

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

  /* Calculate initial approximation, and allow for odd exponent */
  workset.digits=workp;		     /* p for initial calculation */
  t->bits=0; t->digits=3;
  a->bits=0; a->digits=3;
  if ((exp & 1)==0) { 		     /* even exponent */
    /* Set t=0.259, a=0.819 */
    t->exponent=-3;
    a->exponent=-3;
    #if DECDPUN3>=3
t->lsu[0]=259;
a->lsu[0]=819;
    #elif DECDPUN3==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 DECDPUN3>=3
t->lsu[0]=819;
a->lsu[0]=259;
    #elif DECDPUN3==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 (; workset.digits<maxp;) {
    /* set p to min(2*p - 2, maxp)  [hence 3; or: 4, 6, 10, ... , maxp] */
    workset.digits=MINI(workset.digits*2-2, maxp)((workset.digits*2-2)>(maxp)?(maxp):(workset.digits*2-2));
    /* a = 0.5 * (a + f/a) */
    /* [calculated at p then rounded to currentprecision] */
    decDivideOp(b, f, a, &workset, DIVIDE0x80, &ignore); /* b=f/a */
    decAddOp(b, b, a, &workset, 0, &ignore);	       /* b=b+a */
    decMultiplyOp(a, b, t, &workset, &ignore);       /* a=b*0.5 */
    } /* loop */

  /* Here, 0.1 <= a < 1 [Hull], and a has maxp digits */
  /* now reduce to length, etc.; this needs to be done with a */
  /* having the correct exponent so as to handle subnormals */
  /* correctly */
  approxset=*set;			     /* get emin, emax, etc. */
  approxset.round=DEC_ROUND_HALF_EVEN;
  a->exponent+=exp/2; 		     /* set correct exponent */
  rstatus=0;				     /* clear status */
  residue=0;				     /* .. and accumulator */
  decCopyFit(a, a, &approxset, &residue, &rstatus);  /* reduce (if needed) */
  decFinish(a, &approxset, &residue, &rstatus)decFinalize(a,&approxset,&residue,&rstatus);      /* clean and finalize */

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

  /* Preserve status except Inexact/Rounded */
  status|=(rstatus & ~(DEC_Rounded0x00000800|DEC_Inexact0x00000020));

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

  /* a is now at final precision and within 1 ulp of the properly */
  /* rounded square root of f; to ensure proper rounding, compare */
  /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
  /* Here workset.digits=maxp and t=0.5, and a->digits determines */
  /* the ulp */
  workset.digits--;				  /* maxp-1 is OK now */
  t->exponent=-a->digits-1;			  /* make 0.5 ulp */
  decAddOp(b, a, t, &workset, DECNEG0x80, &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, COMPARE0x01, &ignore); /* b ? f, reversed */
  if (decNumberIsNegative(b)(((b)->bits&0x80)!=0)) {		  /* 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, DECNEG0x80, &ignore); /* a = a - 1 ulp */
    /* assign to approx [round to length] */
    approxset.emin-=exp/2;			  /* adjust to match a */
    approxset.emax-=exp/2;
    decAddOp(a, &dzero, a, &approxset, 0, &ignore);
    }
   else {
    decAddOp(b, a, 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, COMPARE0x01, &ignore);   /* b ? f */
    if (decNumberIsNegative(b)(((b)->bits&0x80)!=0)) {		  /* 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] */
approxset.emin-=exp/2;			  /* adjust to match a */
approxset.emax-=exp/2;
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  (still), so adjust back */
  a->exponent+=exp/2; 		     /* set correct exponent */

  /* count droppable zeros [after any subnormal rounding] by */
  /* trimming a copy */
  decNumberCopy(b, a);
  decTrim(b, set, 1, 1, &dropped);	     /* [drops trailing zeros] */

  /* Set Inexact and Rounded.  The answer can only be exact if */
  /* it is short enough so that squaring it could fit in workp */
  /* digits, so this is the only (relatively rare) condition that */
  /* a careful check is needed */
  if (b->digits*2-1 > workp) {	     /* cannot fit */
    status|=DEC_Inexact0x00000020|DEC_Rounded0x00000800;
    }
   else {				     /* could be exact/unrounded */
    uIntuint32_t mstatus=0;			     /* local status */
    decMultiplyOp(b, b, b, &workset, &mstatus); /* try the multiply */
    if (mstatus&DEC_Overflow0x00000200) {	     /* result just won't fit */
status|=DEC_Inexact0x00000020|DEC_Rounded0x00000800;
}
     else {				     /* plausible */
decCompareOp(t, b, rhs, &workset, COMPARE0x01, &mstatus); /* b ? rhs */
if (!ISZERO(t)(*(t)->lsu==0 && (t)->digits==1 && (((t
)->bits&(0x40|0x20|0x10))==0))) status|=DEC_Inexact0x00000020|DEC_Rounded0x00000800; /* not equal */
else { 			     /* is Exact */
 /* here, dropped is the count of trailing zeros in 'a' */
 /* use closest exponent to ideal... */
 Intint32_t todrop=ideal-a->exponent;      /* most that can be dropped */
 if (todrop<0) status|=DEC_Rounded0x00000800; /* ideally would add 0s */
  else {			     /* unrounded */
   /* there are some to drop, but emax may not allow all */
   Intint32_t maxexp=set->emax-set->digits+1;
   Intint32_t maxdrop=maxexp-a->exponent;
   if (todrop>maxdrop && set->clamp) { /* apply clamping */
     todrop=maxdrop;
     status|=DEC_Clamped0x00000400;
     }
   if (dropped<todrop) {	     /* clamp to those available */
     todrop=dropped;
     status|=DEC_Clamped0x00000400;
     }
   if (todrop>0) {		     /* have some to drop */
     decShiftToLeast(a->lsu, D2U(a->digits)((a->digits)<=49?d2utable[a->digits]:((a->digits)
+3 -1)/3), todrop);
     a->exponent+=todrop;	     /* maintain numerical value */
     a->digits-=todrop;	     /* new length */
     }
   }
 }
}
    }

  /* double-check Underflow, as perhaps the result could not have */
  /* been subnormal (initial argument too big), or it is now Exact */
  if (status&DEC_Underflow0x00002000) {
    Intint32_t ae=rhs->exponent+rhs->digits-1;    /* adjusted exponent */
    /* check if truly subnormal */
    #if DECEXTFLAG1			     /* DEC_Subnormal too */
if (ae>=set->emin*2) status&=~(DEC_Subnormal0x00001000|DEC_Underflow0x00002000);
    #else
if (ae>=set->emin*2) status&=~DEC_Underflow0x00002000;
    #endif
    /* check if truly inexact */
    if (!(status&DEC_Inexact0x00000020)) status&=~DEC_Underflow0x00002000;
    }

  decNumberCopy(res, a);		     /* a is now the result */
  } while(0); 			     /* end protected */

free(allocbuff);      /* drop any storage used */
free(allocbufa);      /* .. */
free(allocbufb);      /* .. */
#if DECSUBSET0
free(allocrhs);	     /* .. */
#endif
if (status!=0) decStatus(res, status, set);/* then report status */
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberSquareRoot */

3127/* ------------------------------------------------------------------ */
3128/* decNumberSubtract -- subtract two Numbers			      */
3129/*								      */
3130/*   This computes C = A - B					      */
3131/*								      */
3132/*   res is C, the result.  C may be A and/or B (e.g., X=X-X)	      */
3133/*   lhs is A							      */
3134/*   rhs is B							      */
3135/*   set is the context 					      */
3136/*								      */
3137/* C must have space for set->digits digits.			      */
3138/* ------------------------------------------------------------------ */
3139decNumber * decNumberSubtract(decNumber *res, const decNumber *lhs,
	      const decNumber *rhs, decContext *set) {
uIntuint32_t status=0;			/* accumulator */

decAddOp(res, lhs, rhs, set, DECNEG0x80, &status);
if (status!=0) decStatus(res, status, set);
#if DECCHECK0
decCheckInexact(res, set);
#endif
return res;
} /* decNumberSubtract */

3151/* ------------------------------------------------------------------ */
3152/* decNumberToIntegralExact -- round-to-integral-value with InExact   */
3153/* decNumberToIntegralValue -- round-to-integral-value		      */
3154/*								      */
3155/*   res is the result						      */
3156/*   rhs is input number					      */
3157/*   set is the context 					      */
3158/*								      */
3159/* res must have space for any value of rhs.			      */
3160/*								      */
3161/* This implements the IEEE special operators and therefore treats    */
3162/* special values as valid.  For finite numbers it returns	      */
3163/* rescale(rhs, 0) if rhs->exponent is <0.			      */
3164/* Otherwise the result is rhs (so no error is possible, except for   */
3165/* sNaN).							      */
3166/*								      */
3167/* The context is used for rounding mode and status after sNaN, but   */
3168/* the digits setting is ignored.  The Exact version will signal      */
3169/* Inexact if the result differs numerically from rhs; the other      */
3170/* never signals Inexact.					      */
3171/* ------------------------------------------------------------------ */
3172decNumber * decNumberToIntegralExact(decNumber *res, const decNumber *rhs,
		     decContext *set) {
decNumber dn;
decContext workset;		   /* working context */
uIntuint32_t status=0;		   /* accumulator */

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

/* handle infinities and NaNs */
if (SPECIALARG(rhs->bits & (0x40|0x20|0x10))) {
  if (decNumberIsInfinite(rhs)(((rhs)->bits&0x40)!=0)) decNumberCopy(res, rhs); /* an Infinity */
   else decNaNs(res, rhs, NULL((void*)0), set, &status); /* a NaN */
  }
 else { /* finite */
  /* have a finite number; no error possible (res must be big enough) */
  if (rhs->exponent>=0) return decNumberCopy(res, rhs);
  /* that was easy, but if negative exponent there is 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 */
  decNumberQuantize(res, rhs, &dn, &workset);
  status|=workset.status;
  }
if (status!=0) decStatus(res, status, set);
return res;
} /* decNumberToIntegralExact */

3202decNumber * decNumberToIntegralValue(decNumber *res, const decNumber *rhs,
		     decContext *set) {
decContext workset=*set;	   /* working context */
workset.traps=0;		   /* no traps */
decNumberToIntegralExact(res, rhs, &workset);
/* this never affects set, except for sNaNs; NaN will have been set */
/* or propagated already, so no need to call decStatus */
set->status|=workset.status&DEC_Invalid_operation0x00000080;
return res;
} /* decNumberToIntegralValue */

3213/* ------------------------------------------------------------------ */
3214/* decNumberXor -- XOR two Numbers, digitwise			      */
3215/*								      */
3216/*   This computes C = A ^ B					      */
3217/*								      */
3218/*   res is C, the result.  C may be A and/or B (e.g., X=X^X)	      */
3219/*   lhs is A							      */
3220/*   rhs is B							      */
3221/*   set is the context (used for result length and error report)     */
3222/*								      */
3223/* C must have space for set->digits digits.			      */
3224/*								      */
3225/* Logical function restrictions apply (see above); a NaN is	      */
3226/* returned with Invalid_operation if a restriction is violated.      */
3227/* ------------------------------------------------------------------ */
3228decNumber * decNumberXor(decNumber *res, const decNumber *lhs,
	 const decNumber *rhs, decContext *set) {
const Unituint16_t *ua, *ub;			/* -> operands */
const Unituint16_t *msua, *msub;		/* -> operand msus */
Unituint16_t	*uc, *msuc;			/* -> result and its msu */
Intint32_t	msudigs;			/* digits in res msu */
#if DECCHECK0
if (decCheckOperands(res, lhs, rhs, set)) return res;
#endif

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


3283/* ================================================================== */
3284/* Utility routines						      */
3285/* ================================================================== */

3287/* ------------------------------------------------------------------ */
3288/* decNumberClass -- return the decClass of a decNumber 	      */
3289/*   dn -- the decNumber to test				      */
3290/*   set -- the context to use for Emin 			      */
3291/*   returns the decClass enum					      */
3292/* ------------------------------------------------------------------ */
3293enum decClass decNumberClass(const decNumber *dn, decContext *set) {
if (decNumberIsSpecial(dn)(((dn)->bits&(0x40|0x20|0x10))!=0)) {
  if (decNumberIsQNaN(dn)(((dn)->bits&(0x20))!=0)) return DEC_CLASS_QNAN;
  if (decNumberIsSNaN(dn)(((dn)->bits&(0x10))!=0)) return DEC_CLASS_SNAN;
  /* must be an infinity */
  if (decNumberIsNegative(dn)(((dn)->bits&0x80)!=0)) return DEC_CLASS_NEG_INF;
  return DEC_CLASS_POS_INF;
  }
/* is finite */
if (decNumberIsNormal(dn, set)) { /* most common */
  if (decNumberIsNegative(dn)(((dn)->bits&0x80)!=0)) return DEC_CLASS_NEG_NORMAL;
  return DEC_CLASS_POS_NORMAL;
  }
/* is subnormal or zero */
if (decNumberIsZero(dn)(*(dn)->lsu==0 && (dn)->digits==1 && ((
(dn)->bits&(0x40|0x20|0x10))==0))) {	/* most common */
  if (decNumberIsNegative(dn)(((dn)->bits&0x80)!=0)) return DEC_CLASS_NEG_ZERO;
  return DEC_CLASS_POS_ZERO;
  }
if (decNumberIsNegative(dn)(((dn)->bits&0x80)!=0)) return DEC_CLASS_NEG_SUBNORMAL;
return DEC_CLASS_POS_SUBNORMAL;
} /* decNumberClass */

3315/* ------------------------------------------------------------------ */
3316/* decNumberClassToString -- convert decClass to a string	      */
3317/*								      */
3318/*  eclass is a valid decClass					      */
3319/*  returns a constant string describing the class (max 13+1 chars)   */
3320/* ------------------------------------------------------------------ */
3321const char *decNumberClassToString(enum decClass eclass) {
if (eclass==DEC_CLASS_POS_NORMAL)    return DEC_ClassString_PN"+Normal";
if (eclass==DEC_CLASS_NEG_NORMAL)    return DEC_ClassString_NN"-Normal";
if (eclass==DEC_CLASS_POS_ZERO)      return DEC_ClassString_PZ"+Zero";
if (eclass==DEC_CLASS_NEG_ZERO)      return DEC_ClassString_NZ"-Zero";
if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS"+Subnormal";
if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS"-Subnormal";
if (eclass==DEC_CLASS_POS_INF)       return DEC_ClassString_PI"+Infinity";
if (eclass==DEC_CLASS_NEG_INF)       return DEC_ClassString_NI"-Infinity";
if (eclass==DEC_CLASS_QNAN)	       return DEC_ClassString_QN"NaN";
if (eclass==DEC_CLASS_SNAN)	       return DEC_ClassString_SN"sNaN";
return DEC_ClassString_UN"Invalid";	       /* Unknown */
} /* decNumberClassToString */

3335/* ------------------------------------------------------------------ */
3336/* decNumberCopy -- copy a number				      */
3337/*								      */
3338/*   dest is the target decNumber				      */
3339/*   src  is the source decNumber				      */
3340/*   returns dest						      */
3341/*								      */
3342/* (dest==src is allowed and is a no-op)			      */
3343/* All fields are updated as required.	This is a utility operation,  */
3344/* so special values are unchanged and no error is possible.	      */
3345/* ------------------------------------------------------------------ */
3346decNumber * decNumberCopy(decNumber *dest, const decNumber *src) {

#if DECCHECK0
if (src==NULL((void*)0)) return decNumberZero(dest);
#endif

if (dest==src) return dest;		     /* no copy required */

/* Use explicit assignments here as structure assignment could copy */
/* more than just the lsu (for small DECDPUN).  This would not affect */
/* the value of the results, but could 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>DECDPUN3) {		     /* more Units to come */
  const Unituint16_t *smsup, *s;		     /* work */
  Unituint16_t  *d;				     /* .. */
  /* 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)((src->digits)<=49?d2utable[src->digits]:((src->digits
)+3 -1)/3);	     /* -> source msu+1 */
  for (s=src->lsu+1; s<smsup; s++, d++) *d=*s;
  }
return dest;
} /* decNumberCopy */

3374/* ------------------------------------------------------------------ */
3375/* decNumberCopyAbs -- quiet absolute value operator		      */
3376/*								      */
3377/*   This sets C = abs(A)					      */
3378/*								      */
3379/*   res is C, the result.  C may be A				      */
3380/*   rhs is A							      */
3381/*								      */
3382/* C must have space for set->digits digits.			      */
3383/* No exception or error can occur; this is a quiet bitwise operation.*/
3384/* See also decNumberAbs for a checking version of this.	      */
3385/* ------------------------------------------------------------------ */
3386decNumber * decNumberCopyAbs(decNumber *res, const decNumber *rhs) {
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
#endif
decNumberCopy(res, rhs);
res->bits&=~DECNEG0x80;			/* turn off sign */
return res;
} /* decNumberCopyAbs */

3395/* ------------------------------------------------------------------ */
3396/* decNumberCopyNegate -- quiet negate value operator		      */
3397/*								      */
3398/*   This sets C = negate(A)					      */
3399/*								      */
3400/*   res is C, the result.  C may be A				      */
3401/*   rhs is A							      */
3402/*								      */
3403/* C must have space for set->digits digits.			      */
3404/* No exception or error can occur; this is a quiet bitwise operation.*/
3405/* See also decNumberMinus for a checking version of this.	      */
3406/* ------------------------------------------------------------------ */
3407decNumber * decNumberCopyNegate(decNumber *res, const decNumber *rhs) {
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
#endif
decNumberCopy(res, rhs);
res->bits^=DECNEG0x80;			/* invert the sign */
return res;
} /* decNumberCopyNegate */

3416/* ------------------------------------------------------------------ */
3417/* decNumberCopySign -- quiet copy and set sign operator	      */
3418/*								      */
3419/*   This sets C = A with the sign of B 			      */
3420/*								      */
3421/*   res is C, the result.  C may be A				      */
3422/*   lhs is A							      */
3423/*   rhs is B							      */
3424/*								      */
3425/* C must have space for set->digits digits.			      */
3426/* No exception or error can occur; this is a quiet bitwise operation.*/
3427/* ------------------------------------------------------------------ */
3428decNumber * decNumberCopySign(decNumber *res, const decNumber *lhs,
	      const decNumber *rhs) {
uByteuint8_t sign;				/* rhs sign */
#if DECCHECK0
if (decCheckOperands(res, DECUNUSED, rhs, DECUNCONT)) return res;
#endif
sign=rhs->bits & DECNEG0x80;		/* save sign bit */
decNumberCopy(res, lhs);
res->bits&=~DECNEG0x80;			/* clear the sign */
res->bits|=sign;			/* set from rhs */
return res;
} /* decNumberCopySign */

3441/* ------------------------------------------------------------------ */
3442/* decNumberGetBCD -- get the coefficient in BCD8		      */
3443/*   dn is the source decNumber 				      */
3444/*   bcd is the uInt array that will receive dn->digits BCD bytes,    */
3445/*     most-significant at offset 0				      */
3446/*   returns bcd						      */
3447/*								      */
3448/* bcd must have at least dn->digits bytes.  No error is possible; if */
3449/* dn is a NaN or Infinite, digits must be 1 and the coefficient 0.   */
3450/* ------------------------------------------------------------------ */
3451uByteuint8_t * decNumberGetBCD(const decNumber *dn, uByteuint8_t *bcd) {
uByteuint8_t *ub=bcd+dn->digits-1;	   /* -> lsd */
const Unituint16_t *up=dn->lsu;	   /* Unit pointer, -> lsu */

#if DECDPUN3==1		   /* trivial simple copy */
  for (; ub>=bcd; ub--, up++) *ub=*up;
#else 			   /* chopping needed */
  uIntuint32_t u=*up; 		   /* work */
  uIntuint32_t cut=DECDPUN3;		   /* downcounter through unit */
  for (; ub>=bcd; ub--) {
    *ub=(uByteuint8_t)(u%10);	   /* [*6554 trick inhibits, here] */
    u=u/10;
    cut--;
    if (cut>0) continue;	   /* more in this unit */
    up++;
    u=*up;
    cut=DECDPUN3;
    }
#endif
return bcd;
} /* decNumberGetBCD */

3473/* ------------------------------------------------------------------ */
3474/* decNumberSetBCD -- set (replace) the coefficient from BCD8	      */
3475/*   dn is the target decNumber 				      */
3476/*   bcd is the uInt array that will source n BCD bytes, most-	      */
3477/*     significant at offset 0					      */
3478/*   n is the number of digits in the source BCD array (bcd)	      */
3479/*   returns dn 						      */
3480/*								      */
3481/* dn must have space for at least n digits.  No error is possible;   */
3482/* if dn is a NaN, or Infinite, or is to become a zero, n must be 1   */
3483/* and bcd[0] zero.						      */
3484/* ------------------------------------------------------------------ */
3485decNumber * decNumberSetBCD(decNumber *dn, const uByteuint8_t *bcd, uIntuint32_t n) {
Unituint16_t *up=dn->lsu+D2U(dn->digits)((dn->digits)<=49?d2utable[dn->digits]:((dn->digits
)+3 -1)/3)-1;	/* -> msu [target pointer] */
const uByteuint8_t *ub=bcd;			/* -> source msd */

#if DECDPUN3==1			/* trivial simple copy */
  for (; ub<bcd+n; ub++, up--) *up=*ub;
#else 				/* some assembly needed */
  /* calculate how many digits in msu, and hence first cut */
  Intint32_t cut=MSUDIGITS(n)((n)-(((n)<=49?d2utable[n]:((n)+3 -1)/3)-1)*3);		/* [faster than remainder] */
  for (;up>=dn->lsu; up--) {		/* each Unit from msu */
    *up=0;				/* will take <=DECDPUN digits */
    for (; cut>0; ub++, cut--) *up=X10(*up)(((*up)<<1)+((*up)<<3))+*ub;
    cut=DECDPUN3;			/* next Unit has all digits */
    }
#endif
dn->digits=n; 			/* set digit count */
return dn;
} /* decNumberSetBCD */

3504/* ------------------------------------------------------------------ */
3505/* decNumberIsNormal -- test normality of a decNumber		      */
3506/*   dn is the decNumber to test				      */
3507/*   set is the context to use for Emin 			      */
3508/*   returns 1 if |dn| is finite and >=Nmin, 0 otherwise	      */
3509/* ------------------------------------------------------------------ */
3510Intint32_t decNumberIsNormal(const decNumber *dn, decContext *set) {
Intint32_t ae;				/* adjusted exponent */
#if DECCHECK0
if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
#endif

if (decNumberIsSpecial(dn)(((dn)->bits&(0x40|0x20|0x10))!=0)) return 0; /* not finite */
if (decNumberIsZero(dn)(*(dn)->lsu==0 && (dn)->digits==1 && ((
(dn)->bits&(0x40|0x20|0x10))==0))) return 0;	/* not non-zero */

ae=dn->exponent+dn->digits-1; 	/* adjusted exponent */
if (ae<set->emin) return 0;		/* is subnormal */
return 1;
} /* decNumberIsNormal */

3524/* ------------------------------------------------------------------ */
3525/* decNumberIsSubnormal -- test subnormality of a decNumber	      */
3526/*   dn is the decNumber to test				      */
3527/*   set is the context to use for Emin 			      */
3528/*   returns 1 if |dn| is finite, non-zero, and <Nmin, 0 otherwise    */
3529/* ------------------------------------------------------------------ */
3530Intint32_t decNumberIsSubnormal(const decNumber *dn, decContext *set) {
Intint32_t ae;				/* adjusted exponent */
#if DECCHECK0
if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
#endif

if (decNumberIsSpecial(dn)(((dn)->bits&(0x40|0x20|0x10))!=0)) return 0; /* not finite */
if (decNumberIsZero(dn)(*(dn)->lsu==0 && (dn)->digits==1 && ((
(dn)->bits&(0x40|0x20|0x10))==0))) return 0;	/* not non-zero */

ae=dn->exponent+dn->digits-1; 	/* adjusted exponent */
if (ae<set->emin) return 1;		/* is subnormal */
return 0;
} /* decNumberIsSubnormal */

3544/* ------------------------------------------------------------------ */
3545/* decNumberTrim -- remove insignificant zeros			      */
3546/*								      */
3547/*   dn is the number to trim					      */
3548/*   returns dn 						      */
3549/*								      */
3550/* All fields are updated as required.	This is a utility operation,  */
3551/* so special values are unchanged and no error is possible.  The     */
3552/* zeros are removed unconditionally.				      */
3553/* ------------------------------------------------------------------ */
3554decNumber * decNumberTrim(decNumber *dn) {
Intint32_t  dropped; 		   /* work */
decContext set;		   /* .. */
#if DECCHECK0
if (decCheckOperands(DECUNRESU, DECUNUSED, dn, DECUNCONT)) return dn;
#endif
decContextDefault(&set, DEC_INIT_BASE0);    /* clamp=0 */
return decTrim(dn, &set, 0, 1, &dropped);
} /* decNumberTrim */

3564/* ------------------------------------------------------------------ */
3565/* decNumberVersion -- return the name and version of this module     */
3566/*								      */
3567/* No error is possible.					      */
3568/* ------------------------------------------------------------------ */
3569const char * decNumberVersion(void) {
return DECVERSION"decNumber 3.61";
} /* decNumberVersion */

3573/* ------------------------------------------------------------------ */
3574/* decNumberZero -- set a number to 0				      */
3575/*								      */
3576/*   dn is the number to set, with space for one digit		      */
3577/*   returns dn 						      */
3578/*								      */
3579/* No error is possible.					      */
3580/* ------------------------------------------------------------------ */
3581/* Memset is not used as it is much slower in some environments. */
3582decNumber * decNumberZero(decNumber *dn) {

#if DECCHECK0
if (decCheckOperands(dn, DECUNUSED, DECUNUSED, DECUNCONT)) return dn;
#endif

dn->bits=0;
dn->exponent=0;
dn->digits=1;
dn->lsu[0]=0;
return dn;
} /* decNumberZero */

3595/* ================================================================== */
3596/* Local routines						      */
3597/* ================================================================== */

3599/* ------------------------------------------------------------------ */
3600/* decToString -- lay out a number into a string		      */
3601/*								      */
3602/*   dn     is the number to lay out				      */
3603/*   string is where to lay out the number			      */
3604/*   eng    is 1 if Engineering, 0 if Scientific		      */
3605/*								      */
3606/* string must be at least dn->digits+14 characters long	      */
3607/* No error is possible.					      */
3608/*								      */
3609/* Note that this routine can generate a -0 or 0.000.  These are      */
3610/* never generated in subset to-number or arithmetic, but can occur   */
3611/* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234).	      */
3612/* ------------------------------------------------------------------ */
3613/* If DECCHECK is enabled the string "?" is returned if a number is */
3614/* invalid. */
3615static void decToString(const decNumber *dn, char *string, Flaguint8_t eng) {
Intint32_t exp=dn->exponent;       /* local copy */
Intint32_t e;		      /* E-part value */
Intint32_t pre;		      /* digits before the '.' */
Intint32_t cut;		      /* for counting digits in a Unit */
char *c=string;	      /* work [output pointer] */
const Unituint16_t *up=dn->lsu+D2U(dn->digits)((dn->digits)<=49?d2utable[dn->digits]:((dn->digits
)+3 -1)/3)-1; /* -> msu [input pointer] */
uIntuint32_t u, pow;		      /* work */

#if DECCHECK0
if (decCheckOperands(DECUNRESU, dn, DECUNUSED, DECUNCONT)) {
  strcpy(string, "?");
  return;}
#endif

if (decNumberIsNegative(dn)(((dn)->bits&0x80)!=0)) {   /* Negatives get a minus */
  *c='-';
  c++;
  }
if (dn->bits&DECSPECIAL(0x40|0x20|0x10)) {	   /* Is a special value */
  if (decNumberIsInfinite(dn)(((dn)->bits&0x40)!=0)) {
    strcpy(c,   "Inf");
    strcpy(c+3, "inity");
    return;}
  /* a NaN */
  if (dn->bits&DECSNAN0x10) {	   /* signalling NaN */
    *c='s';
    c++;
    }
  strcpy(c, "NaN");
  c+=3;			   /* step past */
  /* if not a clean non-zero coefficient, that's all there is 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=MSUDIGITS(dn->digits)((dn->digits)-(((dn->digits)<=49?d2utable[dn->digits
]:((dn->digits)+3 -1)/3)-1)*3);	   /* [faster than remainder] */
cut--;			   /* power of ten for digit */

if (exp==0) { 		   /* simple integer [common fastpath] */
  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, pow){ *(c)='0'; pow=DECPOWERS[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;} };
    cut=DECDPUN3-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)) {	   /* engineering: may need to adjust */
    Intint32_t adj;			   /* adjustment */
    /* The C remainder operator is undefined for negative numbers, so */
    /* a positive remainder calculation must be used here */
    if (e<0) {
adj=(-e)%3;
if (adj!=0) adj=3-adj;
}
     else { /* e>0 */
adj=e%3;
}
    e=e-adj;
    /* if dealing with zero still produce an 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)(*(dn)->lsu==0 && (dn)->digits==1 && ((
(dn)->bits&(0x40|0x20|0x10))==0))) pre+=adj;
     else {  /* is zero */
if (adj!=0) {		   /* 0.00Esnn needed */
 e=e+3;
 pre=-(2-adj);
 }
} /* zero */
    } /* eng */
  } /* need exponent */

/* lay out the digits of the coefficient, adding 0s and . as needed */
u=*up;
if (pre>0) {			   /* xxx.xxx or xx00 (engineering) form */
  Intint32_t n=pre;
  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=DECDPUN3-1;
u=*up;
}
    TODIGIT(u, cut, c, pow){ *(c)='0'; pow=DECPOWERS[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;} };
    }
  if (n<dn->digits) { 	   /* 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=DECDPUN3-1;
 u=*up;
 }
TODIGIT(u, cut, c, pow){ *(c)='0'; pow=DECPOWERS[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;} };
}
    }
   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=DECDPUN3-1;
u=*up;
}
    TODIGIT(u, cut, c, pow){ *(c)='0'; pow=DECPOWERS[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;} };
    }
  }

/* 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 anormal numbers */
if (e!=0) {
  Flaguint8_t 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 */
    }
  /* lay out the exponent [_itoa or equivalent is not ANSI C] */
  for (cut=9; cut>=0; cut--) {
    TODIGIT(u, cut, c, pow){ *(c)='0'; pow=DECPOWERS[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;} };
    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;
} /* decToString */

3762/* ------------------------------------------------------------------ */
3763/* decAddOp -- add/subtract operation				      */
3764/*								      */
3765/*   This computes C = A + B					      */
3766/*								      */
3767/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)	      */
3768/*   lhs is A							      */
3769/*   rhs is B							      */
3770/*   set is the context 					      */
3771/*   negate is DECNEG if rhs should be negated, or 0 otherwise	      */
3772/*   status accumulates status for the caller			      */
3773/*								      */
3774/* C must have space for set->digits digits.			      */
3775/* Inexact in status must be 0 for correct Exact zero sign in result  */
3776/* ------------------------------------------------------------------ */
3777/* If possible, the coefficient is calculated directly into C.	      */
3778/* However, if: 						      */
3779/*   -- a digits+1 calculation is needed because the numbers are      */
3780/*	unaligned and span more than set->digits digits 	      */
3781/*   -- a carry to digits+1 digits looks possible		      */
3782/*   -- C is the same as A or B, and the result would destructively   */
3783/*	overlap the A or B coefficient				      */
3784/* then the result must be calculated into a temporary buffer.	In    */
3785/* this case a local (stack) buffer is used if possible, and only if  */
3786/* too long for that does malloc become the final resort.	      */
3787/*								      */
3788/* Misalignment is handled as follows:				      */
3789/*   Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp.    */
3790/*   BPad: Apply the padding by a combination of shifting (whole      */
3791/*	   units) and multiplication (part units).		      */
3792/*								      */
3793/* Addition, especially x=x+1, is speed-critical.		      */
3794/* The static buffer is larger than might be expected to allow for    */
3795/* calls from higher-level funtions (notable exp).		      */
3796/* ------------------------------------------------------------------ */
3797static decNumber * decAddOp(decNumber *res, const decNumber *lhs,
	    const decNumber *rhs, decContext *set,
	    uByteuint8_t negate, uIntuint32_t *status) {
#if DECSUBSET0
decNumber *alloclhs=NULL((void*)0);	   /* non-NULL if rounded lhs allocated */
decNumber *allocrhs=NULL((void*)0);	   /* .., rhs */
#endif
Intint32_t	rhsshift;		   /* working shift (in Units) */
Intint32_t	maxdigits;		   /* longest logical length */
Intint32_t	mult;			   /* multiplier */
Intint32_t	residue;		   /* rounding accumulator */
uByteuint8_t bits;			   /* result bits */
Flaguint8_t	diffsign;		   /* non-0 if arguments have different sign */
Unituint16_t	*acc;			   /* accumulator for result */
Unituint16_t	accbuff[SD2U(DECBUFFER*2+20)(((36*2+20)+3 -1)/3)]; /* local buffer [*2+20 reduces many */
		   /* allocations when called from */
		   /* other operations, notable exp] */
Unituint16_t	*allocacc=NULL((void*)0); 	   /* -> allocated acc buffer, iff allocated */
Intint32_t	reqdigits=set->digits;	   /* local copy; requested DIGITS */
Intint32_t	padding;		   /* work */

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

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

  /* note whether signs differ [used all paths] */
  diffsign=(Flaguint8_t)((lhs->bits^rhs->bits^negate)&DECNEG0x80);

  /* handle infinities and NaNs */
  if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10))) {			/* a special bit set */
    if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10)) & (DECSNAN0x10 | DECNAN0x20))  /* a NaN */
decNaNs(res, lhs, rhs, set, status);
     else { /* one or two infinities */
if (decNumberIsInfinite(lhs)(((lhs)->bits&0x40)!=0)) { /* LHS is infinity */
 /* two infinities with different signs is invalid */
 if (decNumberIsInfinite(rhs)(((rhs)->bits&0x40)!=0) && diffsign) {
   *status|=DEC_Invalid_operation0x00000080;
   break;
   }
 bits=lhs->bits & DECNEG0x80;	/* get sign from LHS */
 }
else bits=(rhs->bits^negate) & DECNEG0x80;/* RHS must be Infinity */
bits|=DECINF0x40;
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)(*(lhs)->lsu==0 && (lhs)->digits==1 && (
((lhs)->bits&(0x40|0x20|0x10))==0))) {
    Intint32_t adjust;			/* work */
    Intint32_t 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 DECSUBSET0
    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)(*(res)->lsu==0 && (res)->digits==1 && (
((res)->bits&(0x40|0x20|0x10))==0))) {		/* both 0: special IEEE 754 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=DECNEG0x80;	/* 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_Rounded0x00000800;		  /* [but exact] */
     }
   res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
   res->exponent+=adjust;		  /* set the exponent. */
   }
 } /* non-0 res */
    #if DECSUBSET0
} /* extended */
    #endif
    decFinish(res, set, &residue, status)decFinalize(res,set,&residue,status);	  /* clean and finalize */
    break;}

  if (ISZERO(rhs)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0))) {			/* [lhs is non-zero] */
    Intint32_t adjust;			/* work */
    Intint32_t 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 DECSUBSET0
    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_Rounded0x00000800;		/* [but exact] */
   }
 res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
 res->exponent+=adjust;		/* set the exponent. */
 }
    #if DECSUBSET0
} /* extended */
    #endif
    decFinish(res, set, &residue, status)decFinalize(res,set,&residue,status);	  /* clean and finalize */
    break;}

  /* [NB: 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 (padding==0
&& rhs->digits<=DECDPUN3
&& rhs->exponent>=set->emin	/* [some normals drop through] */
&& rhs->exponent<=set->emax-set->digits+1 /* [could clamp] */
&& rhs->digits<=reqdigits
&& lhs->digits<=reqdigits) {
    Intint32_t partial=*lhs->lsu;
    if (!diffsign) {			/* adding */
partial+=*rhs->lsu;
if ((partial<=DECDPUNMAX999)	/* result fits in unit */
&& (lhs->digits>=DECDPUN3 ||	/* .. and no digits-count change */
    partial<(Intint32_t)powersDECPOWERS[lhs->digits])) { /* .. */
 if (res!=lhs) decNumberCopy(res, lhs);  /* not in place */
 *res->lsu=(Unituint16_t)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=(Unituint16_t)partial;
 /* this could have reduced digits [but result>0] */
 res->digits=decGetDigits(res->lsu, D2U(res->digits)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3));
 break;
 }
/* else drop out for careful subtract */
}
    }

  /* Now align (pad) the lhs or rhs so they can be added or */
  /* subtracted, 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 of one and the highest of the */
  /* other) padding with up to DIGITS-1 trailing zeros may be */
  /* needed; 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 is needed] */
  if (padding!=0) {
    /* some padding needed; always pad the RHS, as any required */
    /* padding can then be effected by a simple combination of */
    /* shifts and a multiply */
    Flaguint8_t swapped=0;
    if (padding<0) {			/* LHS needs the padding */
const decNumber *t;
padding=-padding;		/* will be +ve */
bits=(uByteuint8_t)(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 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 use the relative sign indication... */
Intint32_t 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)decFinalize(res,set,&residue,status);	  /* done */
break;}

    /* LHS digits may affect result */
    rhsshift=D2U(padding+1)((padding+1)<=49?d2utable[padding+1]:((padding+1)+3 -1)/3)-1;	/* this much by Unit shift .. */
    mult=powersDECPOWERS[padding-(rhsshift*DECDPUN3)]; /* .. 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 add direct to result */
  /* If destructive overlap, or the number is too long, or a carry or */
  /* borrow to DIGITS+1 might be possible, a buffer must be used. */
  /* [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; units for maxdigits digits will be */
    /* needed, +1 Unit for carry or borrow */
    Intint32_t need=D2U(maxdigits)((maxdigits)<=49?d2utable[maxdigits]:((maxdigits)+3 -1)/3)+1;
    acc=accbuff;			/* assume use local buffer */
    if (need*sizeof(Unituint16_t)>sizeof(accbuff)) {
/* printf("malloc add %ld %ld\n", need, sizeof(accbuff)); */
allocacc=(Unituint16_t *)malloc(need*sizeof(Unituint16_t));
if (allocacc==NULL((void*)0)) {		/* hopeless -- abandon */
 *status|=DEC_Insufficient_storage0x00000010;
 break;}
acc=allocacc;
}
    }

  res->bits=(uByteuint8_t)(bits&DECNEG0x80);	/* it's now safe to overwrite.. */
  res->exponent=lhs->exponent;	/* .. operands (even if aliased) */

  #if DECTRACE0
    decDumpAr('A', lhs->lsu, D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3));
    decDumpAr('B', rhs->lsu, D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3));
    printf("	:h: %ld %ld\n", rhsshift, mult);
  #endif

  /* add [A+B*m] or subtract [A+B*(-m)] */
  res->digits=decUnitAddSub(lhs->lsu, D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3),
	      rhs->lsu, D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3),
	      rhsshift, acc, mult)
      *DECDPUN3;	   /* [units -> digits] */
  if (res->digits<0) {	   /* borrowed... */
    res->digits=-res->digits;
    res->bits^=DECNEG0x80;	   /* flip the sign */
    }
  #if DECTRACE0
    decDumpAr('+', acc, D2U(res->digits)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3));
  #endif

  /* If a buffer was used the result must be copied back, possibly */
  /* shortening.  (If no buffer was used then the result must have */
  /* fit, so can't need rounding and residue must be 0.) */
  residue=0;			   /* clear accumulator */
  if (acc!=res->lsu) {
    #if DECSUBSET0
    if (set->extended) {	   /* round from first significant digit */
    #endif
/* remove leading zeros that were added due to rounding up to */
/* integral Units -- before the test for rounding. */
if (res->digits>reqdigits)
 res->digits=decGetDigits(acc, D2U(res->digits)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3));
decSetCoeff(res, set, acc, res->digits, &residue, status);
    #if DECSUBSET0
}
     else { /* subset arithmetic rounds from original significant digit */
/* May have an underestimate.  This only occurs when both */
/* numbers fit in DECDPUN digits and are padding with a */
/* negative multiple (-10, -100...) and the top digit(s) become */
/* 0.  (This only matters when using X3.274 rules where the */
/* leading zero could be included in the rounding.) */
if (res->digits<maxdigits) {
 *(acc+D2U(res->digits)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3))=0; /* ensure leading 0 is there */
 res->digits=maxdigits;
 }
else {
 /* remove leading zeros that 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)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3));
   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;		     /* did what needed to be done */
 }
} /* subset */
    #endif
    } /* used buffer */

  /* strip leading zeros [these were left on in case of subset subtract] */
  res->digits=decGetDigits(res->lsu, D2U(res->digits)((res->digits)<=49?d2utable[res->digits]:((res->digits
)+3 -1)/3));

  /* apply checks and rounding */
  decFinish(res, set, &residue, status)decFinalize(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)(*(res)->lsu==0 && (res)->digits==1 && (
((res)->bits&(0x40|0x20|0x10))==0)) && diffsign
   #if DECSUBSET0
   && set->extended
   #endif
   && (*status&DEC_Inexact0x00000020)==0) {
    if (set->round==DEC_ROUND_FLOOR) res->bits|=DECNEG0x80;   /* sign - */
		  else res->bits&=~DECNEG0x80;  /* sign + */
    }
  } while(0); 			     /* end protected */

free(allocacc);	     /* drop any storage used */
#if DECSUBSET0
free(allocrhs);	     /* .. */
free(alloclhs);	     /* .. */
#endif
return res;
} /* decAddOp */

4141/* ------------------------------------------------------------------ */
4142/* decDivideOp -- division operation				      */
4143/*								      */
4144/*  This routine performs the calculations for all four division      */
4145/*  operators (divide, divideInteger, remainder, remainderNear).      */
4146/*								      */
4147/*  C=A op B							      */
4148/*								      */
4149/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)	      */
4150/*   lhs is A							      */
4151/*   rhs is B							      */
4152/*   set is the context 					      */
4153/*   op  is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively.    */
4154/*   status is the usual accumulator				      */
4155/*								      */
4156/* C must have space for set->digits digits.			      */
4157/*								      */
4158/* ------------------------------------------------------------------ */
4159/*   The underlying algorithm of this routine is the same as in the   */
4160/*   1981 S/370 implementation, that is, non-restoring long division  */
4161/*   with bi-unit (rather than bi-digit) estimation for each unit     */
4162/*   multiplier.  In this pseudocode overview, complications for the  */
4163/*   Remainder operators and division residues for exact rounding are */
4164/*   omitted for clarity.					      */
4165/*								      */
4166/*     Prepare operands and handle special values		      */
4167/*     Test for x/0 and then 0/x				      */
4168/*     Exp =Exp1 - Exp2 					      */
4169/*     Exp =Exp +len(var1) -len(var2)				      */
4170/*     Sign=Sign1 * Sign2					      */
4171/*     Pad accumulator (Var1) to double-length with 0's (pad1)	      */
4172/*     Pad Var2 to same length as Var1				      */
4173/*     msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round  */
4174/*     have=0							      */
4175/*     Do until (have=digits+1 OR residue=0)			      */
4176/*	 if exp<0 then if integer divide/residue then leave	      */
4177/*	 this_unit=0						      */
4178/*	 Do forever						      */
4179/*	    compare numbers					      */
4180/*	    if <0 then leave inner_loop 			      */
4181/*	    if =0 then (* quick exit without subtract *) do	      */
4182/*	       this_unit=this_unit+1; output this_unit		      */
4183/*	       leave outer_loop; end				      */
4184/*	    Compare lengths of numbers (mantissae):		      */
4185/*	    If same then tops2=msu2pair -- {units 1&2 of var2}	      */
4186/*		    else tops2=msu2plus -- {0, unit 1 of var2}	      */
4187/*	    tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
4188/*	    mult=tops1/tops2  -- Good and safe guess at divisor       */
4189/*	    if mult=0 then mult=1				      */
4190/*	    this_unit=this_unit+mult				      */
4191/*	    subtract						      */
4192/*	    end inner_loop					      */
4193/*	  if have\=0 | this_unit\=0 then do			      */
4194/*	    output this_unit					      */
4195/*	    have=have+1; end					      */
4196/*	  var2=var2/10						      */
4197/*	  exp=exp-1						      */
4198/*	  end outer_loop					      */
4199/*     exp=exp+1   -- set the proper exponent			      */
4200/*     if have=0 then generate answer=0 			      */
4201/*     Return (Result is defined by Var1)			      */
4202/*								      */
4203/* ------------------------------------------------------------------ */
4204/* Two working buffers are needed during the division; one (digits+   */
4205/* 1) to accumulate the result, and the other (up to 2*digits+1) for  */
4206/* long subtractions.  These are acc and var1 respectively.	      */
4207/* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
4208/* The static buffers may be larger than might be expected to allow   */
4209/* for calls from higher-level funtions (notable exp).		      */
4210/* ------------------------------------------------------------------ */
4211static decNumber * decDivideOp(decNumber *res,
	       const decNumber *lhs, const decNumber *rhs,
	       decContext *set, Flaguint8_t op, uIntuint32_t *status) {
#if DECSUBSET0
decNumber *alloclhs=NULL((void*)0);	   /* non-NULL if rounded lhs allocated */
decNumber *allocrhs=NULL((void*)0);	   /* .., rhs */
#endif
Unituint16_t	accbuff[SD2U(DECBUFFER+DECDPUN+10)(((36 +3 +10)+3 -1)/3)]; /* local buffer */
Unituint16_t	*acc=accbuff;		   /* -> accumulator array for result */
Unituint16_t	*allocacc=NULL((void*)0); 	   /* -> allocated buffer, iff allocated */
Unituint16_t	*accnext;		   /* -> where next digit will go */
Intint32_t	acclength;		   /* length of acc needed [Units] */
Intint32_t	accunits;		   /* count of units accumulated */
Intint32_t	accdigits;		   /* count of digits accumulated */

Unituint16_t	varbuff[SD2U(DECBUFFER*2+DECDPUN)(((36*2+3)+3 -1)/3)];  /* buffer for var1 */
Unituint16_t	*var1=varbuff;		   /* -> var1 array for long subtraction */
Unituint16_t	*varalloc=NULL((void*)0); 	   /* -> allocated buffer, iff used */
Unituint16_t	*msu1;			   /* -> msu of var1 */

const Unituint16_t *var2;		   /* -> var2 array */
const Unituint16_t *msu2;		   /* -> msu of var2 */
Intint32_t	msu2plus;		   /* msu2 plus one [does not vary] */
eIntint32_t	msu2pair;		   /* msu2 pair plus one [does not vary] */

Intint32_t	var1units, var2units;	   /* actual lengths */
Intint32_t	var2ulen;		   /* logical length (units) */
Intint32_t	var1initpad=0;		   /* var1 initial padding (digits) */
Intint32_t	maxdigits;		   /* longest LHS or required acc length */
Intint32_t	mult;			   /* multiplier for subtraction */
Unituint16_t	thisunit;		   /* current unit being accumulated */
Intint32_t	residue;		   /* for rounding */
Intint32_t	reqdigits=set->digits;	   /* requested DIGITS */
Intint32_t	exponent;		   /* working exponent */
Intint32_t	maxexponent=0;		   /* DIVIDE maximum exponent if unrounded */
uByteuint8_t bits;			   /* working sign */
Unituint16_t	*target;		   /* work */
const Unituint16_t *source;		   /* .. */
uIntuint32_t	const *pow;		   /* .. */
Intint32_t	shift, cut;		   /* .. */
#if DECSUBSET0
Intint32_t	dropped;		   /* work */
#endif

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

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

  bits=(lhs->bits^rhs->bits)&DECNEG0x80;	/* assumed sign for divisions */

  /* handle infinities and NaNs */
  if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10))) {			/* a special bit set */
1
Assuming the condition is false→
    if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10)) & (DECSNAN0x10 | DECNAN0x20)) { /* one or two NaNs */
decNaNs(res, lhs, rhs, set, status);
break;
}
    /* one or two infinities */
    if (decNumberIsInfinite(lhs)(((lhs)->bits&0x40)!=0)) {	/* LHS (dividend) is infinite */
if (decNumberIsInfinite(rhs)(((rhs)->bits&0x40)!=0) || /* two infinities are invalid .. */
   op & (REMAINDER0x40 | REMNEAR0x10)) { /* as is remainder of infinity */
 *status|=DEC_Invalid_operation0x00000080;
 break;
 }
/* [Note that infinity/0 raises no exceptions] */
decNumberZero(res);
res->bits=bits|DECINF0x40;		/* set +/- infinity */
break;
}
     else {				/* RHS (divisor) is infinite */
residue=0;
if (op&(REMAINDER0x40|REMNEAR0x10)) {
 /* 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&DIVIDE0x80) {
   res->exponent=set->emin-set->digits+1;
   *status|=DEC_Clamped0x00000400;
   }
 }
decFinish(res, set, &residue, status)decFinalize(res,set,&residue,status);
break;
}
    }

  /* handle 0 rhs (x/0) */
  if (ISZERO(rhs)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0))) {			/* x/0 is always exceptional */
2
←
Taking false branch→
3
←
Assuming the condition is false→
    if (ISZERO(lhs)(*(lhs)->lsu==0 && (lhs)->digits==1 && (
((lhs)->bits&(0x40|0x20|0x10))==0))) {
decNumberZero(res);		/* [after lhs test] */
*status|=DEC_Division_undefined0x00000008;/* 0/0 will become NaN */
}
     else {
decNumberZero(res);
if (op&(REMAINDER0x40|REMNEAR0x10)) *status|=DEC_Invalid_operation0x00000080;
else {
 *status|=DEC_Division_by_zero0x00000002; /* x/0 */
 res->bits=bits|DECINF0x40;	 /* .. is +/- Infinity */
 }
}
    break;}

  /* handle 0 lhs (0/x) */
  if (ISZERO(lhs)(*(lhs)->lsu==0 && (lhs)->digits==1 && (
((lhs)->bits&(0x40|0x20|0x10))==0))) {			/* 0/x [x!=0] */
4
←
Assuming the condition is false→
    #if DECSUBSET0
    if (!set->extended) decNumberZero(res);
     else {
    #endif
if (op&DIVIDE0x80) {
 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&DIVIDEINT0x20) {
 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 DECSUBSET0
}
    #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==DIVIDE0x80)) {
5
←
Assuming 'exponent' is >= 0→
    if (op&DIVIDEINT0x20) {
decNumberZero(res);		     /* integer part is 0 */
#if DECSUBSET0
if (set->extended)
#endif
 res->bits=bits;		     /* set +/- zero */
break;}
    /* fastpath remainders so long as the lhs has the smaller */
    /* (or equal) exponent */
    if (lhs->exponent<=rhs->exponent) {
if (op&REMAINDER0x40 || 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)decFinalize(res,set,&residue,status);
 break;
 }
/* [unsafe REMNEAR drops through] */
}
    } /* fastpaths */

  /* Long (slow) division is needed; 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)((reqdigits+3)<=49?d2utable[reqdigits+3]:((reqdigits+3)+3 -
1)/3);	/* in Units */
6
←
Assuming the condition is false→
7
←
'?' condition is false→
  if (acclength*sizeof(Unituint16_t)>sizeof(accbuff)) {
8
←
Assuming the condition is false→
9
←
Taking false branch→
    /* printf("malloc dvacc %ld units\n", acclength); */
    allocacc=(Unituint16_t *)malloc(acclength*sizeof(Unituint16_t));
    if (allocacc==NULL((void*)0)) {		/* hopeless -- abandon */
*status|=DEC_Insufficient_storage0x00000010;
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 needed 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;
10
←
Assuming 'maxdigits' is >= field 'digits'→
  var1units=D2U(maxdigits)((maxdigits)<=49?d2utable[maxdigits]:((maxdigits)+3 -1)/3)+2;
11
←
Taking false branch→
12
←
Assuming 'maxdigits' is > 49→
13
←
'?' condition is false→
  /* allocate a guard unit above msu1 for REMAINDERNEAR */
  if (!(op&DIVIDE0x80)) var1units++;
14
←
Assuming the condition is false→
15
←
Taking false branch→
  if ((var1units+1)*sizeof(Unituint16_t)>sizeof(varbuff)) {
16
←
Assuming the condition is false→
17
←
Taking false branch→
    /* printf("malloc dvvar %ld units\n", var1units+1); */
    varalloc=(Unituint16_t *)malloc((var1units+1)*sizeof(Unituint16_t));
    if (varalloc==NULL((void*)0)) {		/* hopeless -- abandon */
*status|=DEC_Insufficient_storage0x00000010;
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 a */
  /* subtract in place is always possible.  The rhs (var2) has */
  /* virtual padding (implemented by decUnitAddSub). */
  /* One guard unit was allocated above msu1 for rem=rem+rem in */
  /* REMAINDERNEAR. */
  msu1=var1+var1units-1;		/* msu of var1 */
  source=lhs->lsu+D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3)-1; /* msu of input array */
18
←
Assuming field 'digits' is > 49→
19
←
'?' condition is false→
  for (target=msu1; source>=lhs->lsu; source--, target--) *target=*source;
20
←
Assuming 'source' is < field 'lsu'→
21
←
Loop condition is false. Execution continues on line 4442→
  for (; target>=var1; target--) *target=0;
22
←
Assuming 'target' is < 'var1'→
23
←
Loop condition is false. Execution continues on line 4445→

  /* rhs (var2) is left-aligned with var1 at the start */
  var2ulen=var1units; 		/* rhs logical length (units) */
  var2units=D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3); 	/* rhs actual length (units) */
24
←
Assuming field 'digits' is > 49→
25
←
'?' condition is false→
  var2=rhs->lsu;			/* -> rhs array */
  msu2=var2+var2units-1;		/* -> msu of var2 [never changes] */
  /* now set up the variables which will be used for estimating the */
  /* multiplication factor.  If these variables are not exact, add */
  /* 1 to make sure that the multiplier is never overestimated. */
  msu2plus=*msu2;			/* it's value .. */
  if (var2units>1) msu2plus++;	/* .. +1 if any more */
26
←
Assuming 'var2units' is <= 1→
27
←
Taking false branch→
  msu2pair=(eIntint32_t)*msu2*(DECDPUNMAX999+1);/* top two pair .. */
  if (var2units27.1
'var2units' is <= 1
>1) {			/* .. [else treat 2nd as 0] */
28
←
Taking false branch→
    msu2pair+=*(msu2-1);		/* .. */
    if (var2units>2) msu2pair++;	/* .. +1 if any more */
    }

  /* The calculation is working in units, which may have leading zeros, */
  /* but the exponent was calculated 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. */
  /* [This is actually done by counting the digits and negating, as */
  /* lead1=DECDPUN-digits1, and similarly for lead2.] */
  for (pow=&powersDECPOWERS[1]; *msu1>=*pow; pow++) exponent--;
29
←
Loop condition is false. Execution continues on line 4467→
  for (pow=&powersDECPOWERS[1]; *msu2>=*pow; pow++) exponent++;
30
←
Loop condition is false. Execution continues on line 4474→

  /* Now, if doing an integer divide or remainder, ensure that */
  /* the result will be Unit-aligned.  To do this, shift the var1 */
  /* accumulator towards least if need be.  (It's much easier to */
  /* do this now than to reassemble the residue afterwards, if */
  /* doing a remainder.)  Also ensure the exponent is not negative. */
  if (!(op&DIVIDE0x80)) {
31
←
Taking false branch→
    Unituint16_t *u;				/* work */
    /* save the initial 'false' padding of var1, in digits */
    var1initpad=(var1units-D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3))*DECDPUN3;
    /* Determine the shift to do. */
    if (exponent<0) cut=-exponent;
     else cut=DECDPUN3-exponent%DECDPUN3;
    decShiftToLeast(var1, var1units, cut);
    exponent+=cut;			/* maintain numerical value */
    var1initpad-=cut; 		/* .. and reduce padding */
    /* clean any most-significant units which were just emptied */
    for (u=msu1; cut>=DECDPUN3; cut-=DECDPUN3, 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) {
32
←
Assuming the condition is true→
33
←
Taking true branch→
var2ulen--;			/* shift down */
exponent-=DECDPUN3;		/* 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 */
34
←
Loop condition is true.  Entering loop body→
    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--;
35
←
Assuming the condition is false→

if (var1units<var2ulen) break;	     /* var1 too low for subtract */
36
←
Assuming 'var1units' is < 'var2ulen'→
37
←
Taking true branch→
if (var1units==var2ulen) {	     /* unit-by-unit compare needed */
 /* compare the two numbers, from msu */
 const Unituint16_t *pv1, *pv2;
 Unituint16_t 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 */
   /* the calculation is done; leave the loop with residue=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=(Intint32_t)(((eIntint32_t)*msu1*(DECDPUNMAX999+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 only one unit for var2 can be used. */
 mult=(Intint32_t)(((eIntint32_t)*msu1*(DECDPUNMAX999+1)+*(msu1-1))/msu2plus);
 }
if (mult==0) mult=1;		     /* must always be at least 1 */
/* subtraction needed; var1 is > var2 */
thisunit=(Unituint16_t)(thisunit+mult);      /* accumulate */
/* subtract var1-var2, into var1; only the overlap needs */
/* processing, as this is an in-place calculation */
shift=var2ulen-var2units;
#if DECTRACE0
 decDumpAr('1', &var1[shift], var1units-shift);
 decDumpAr('2', var2, var2units);
 printf("m=%ld\n", -mult);
#endif
decUnitAddSub(&var1[shift], var1units-shift,
      var2, var2units, 0,
      &var1[shift], -mult);
#if DECTRACE0
 decDumpAr('#', &var1[shift], var1units-shift);
#endif
/* var1 now probably has leading zeros; these are removed at the */
/* top of the inner loop. */
} /* inner loop */

    /* The next unit has been calculated in full; unless it's a */
    /* leading zero, add to acc */
    if (accunits37.1
'accunits' is equal to 0
!=0 || thisunit37.2
'thisunit' is equal to 0
!=0) {      /* is first or non-zero */
*accnext=thisunit;		     /* store in accumulator */
/* account exactly for the new digits */
if (accunits==0) {
 accdigits++;			     /* at least one */
 for (pow=&powersDECPOWERS[1]; thisunit>=*pow; pow++) accdigits++;
 }
else accdigits+=DECDPUN3;
accunits++;			     /* update count */
accnext--;			     /* ready for next */
if (accdigits>reqdigits) break;      /* have enough digits */
}

    /* if the residue is zero, the operation is done (unless divide */
    /* or divideInteger and still not enough digits yet) */
    if (*var1==0 && var1units==1) {	     /* residue is 0 */
38
←
The left operand of '==' is a garbage value
if (op&(REMAINDER0x40|REMNEAR0x10)) break;
if ((op&DIVIDE0x80) && (exponent<=maxexponent)) break;
/* [drop through if divideInteger] */
}
    /* also done enough if calculating remainder or integer */
    /* divide and just did the last ('units') unit */
    if (exponent==0 && !(op&DIVIDE0x80)) 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-=DECDPUN3;		     /* 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 collected in acc */
  if (accunits==0) {		   /* acc is 0 */
    accunits=1;		   /* show have a unit .. */
    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&DIVIDE0x80) {
    /* record the presence of any residue, for rounding */
    if (*var1!=0 || var1units>1) residue=1;
     else { /* no residue */
/* 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 DECDPUN3>1
Unituint16_t lsu=*accnext;
if (!(lsu&0x01) && (lsu!=0)) {
 /* count the trailing zeros */
 Intint32_t drop=0;
 for (;; drop++) {    /* [will terminate because lsu!=0] */
   if (exponent>=maxexponent) break;	  /* don't chop real 0s */
   #if DECDPUN3<=4
     if ((lsu-QUOT10(lsu, drop+1)((((uint32_t)(lsu)>>(drop+1))*multies[drop+1])>>17
)
  *powersDECPOWERS[drop+1])!=0) break;	  /* found non-0 digit */
   #else
     if (lsu%powersDECPOWERS[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)((accdigits)<=49?d2utable[accdigits]:((accdigits)+3 -1)/3);
   /* [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_impossible0x00000004;
break;
}
    if (op & (REMAINDER0x40|REMNEAR0x10)) {
/* [Here, the exponent will be 0, because var1 was adjusted */
/* appropriately.] */
Intint32_t postshift;			     /* work */
Flaguint8_t wasodd=0;			     /* integer was odd */
Unituint16_t *quotlsu;			     /* for save */
Intint32_t  quotdigits;		     /* .. */

bits=lhs->bits; 		     /* remainder sign is always as lhs */

/* Fastpath when residue is truly 0 is worthwhile [and */
/* simplifies the code below] */
if (*var1==0 && var1units==1) {      /* residue is 0 */
 Intint32_t exp=lhs->exponent;	     /* save min(exponents) */
 if (rhs->exponent<exp) exp=rhs->exponent;
 decNumberZero(res);		     /* 0 coefficient */
 #if DECSUBSET0
 if (set->extended)
 #endif
 res->exponent=exp;		     /* .. with proper exponent */
 res->bits=(uByteuint8_t)(bits&DECNEG0x80);	   /* [cleaned] */
 decFinish(res, set, &residue, status)decFinalize(res,set,&residue,status);   /* might clamp */
 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)((accdigits)<=49?d2utable[accdigits]:((accdigits)+3 -1)/3);

exponent=lhs->exponent; 	/* exponent is smaller of lhs & rhs */
if (rhs->exponent<exponent) exponent=rhs->exponent;

/* Now correct the result if 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&REMNEAR0x10) {
 Intint32_t compare, tarunits;	/* work */
 Unituint16_t *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 must now be compared to the */
 /* RHS.  The remainder's exponent may be smaller than the RHS's. */
 compare=decUnitCompare(accnext, tarunits, rhs->lsu, D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3),
		 rhs->exponent-exponent);
 if (compare==BADINT(int32_t)0x80000000) {	     /* deep trouble */
   *status|=DEC_Insufficient_storage0x00000010;
   break;}

 /* now restore the remainder by dividing by two; the lsu */
 /* is known to be even. */
 for (up=accnext; up<accnext+tarunits; up++) {
   Intint32_t half;		   /* half to add to lower unit */
   half=*up & 0x01;
   *up/=2;		   /* [shift] */
   if (!half) continue;
   *(up-1)+=(DECDPUNMAX999+1)/2;
   }
 /* [accunits still describes the original remainder length] */

 if (compare>0 || (compare==0 && wasodd)) { /* adjustment needed */
   Intint32_t 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 */
   Flaguint8_t allnines=0;		     /* 1 if quotient all nines */
   if (quotdigits==reqdigits) {     /* could be borderline */
     for (up=quotlsu; ; up++) {
if (quotdigits>DECDPUN3) {
  if (*up!=DECDPUNMAX999) break;/* non-nines */
  }
 else { 		     /* this is the last Unit */
  if (*up==powersDECPOWERS[quotdigits]-1) allnines=1;
  break;
  }
quotdigits-=DECDPUN3;	     /* checked those digits */
} /* up */
     } /* borderline check */
   if (allnines) {
     *status|=DEC_Division_impossible0x00000004;
     break;}

   /* rem-rhs is needed; the sign will invert.  Again, var1 */
   /* can safely be used for the working Units array. */
   exp=rhs->exponent-exponent;      /* RHS padding needed */
   /* Calculate units and remainder from exponent. */
   expunits=exp/DECDPUN3;
   exprem=exp%DECDPUN3;
   /* subtract [A+B*(-m)]; the result will always be negative */
   accunits=-decUnitAddSub(accnext, accunits,
		    rhs->lsu, D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3),
		    expunits, accnext, -(Intint32_t)powersDECPOWERS[exprem]);
   accdigits=decGetDigits(accnext, accunits); /* count digits exactly */
   accunits=D2U(accdigits)((accdigits)<=49?d2utable[accdigits]:((accdigits)+3 -1)/3);	/* and recalculate the units for copy */
   /* [exponent is as for original remainder] */
   bits^=DECNEG0x80;		/* flip the sign */
   }
 } /* REMNEAR */
} /* REMAINDER or REMNEAR */
    } /* not DIVIDE */

  /* Set exponent and bits */
  res->exponent=exponent;
  res->bits=(uByteuint8_t)(bits&DECNEG0x80);	     /* [cleaned] */

  /* Now the coefficient. */
  decSetCoeff(res, set, accnext, accdigits, &residue, status);

  decFinish(res, set, &residue, status)decFinalize(res,set,&residue,status);   /* final cleanup */

  #if DECSUBSET0
  /* If a divide then strip trailing zeros if subset [after round] */
  if (!set->extended && (op==DIVIDE0x80)) decTrim(res, set, 0, 1, &dropped);
  #endif
  } while(0); 			     /* end protected */

free(varalloc);	/* drop any storage used */
free(allocacc);	/* .. */
#if DECSUBSET0
free(allocrhs);	/* .. */
free(alloclhs);	/* .. */
#endif
return res;
} /* decDivideOp */

4794/* ------------------------------------------------------------------ */
4795/* decMultiplyOp -- multiplication operation			      */
4796/*								      */
4797/*  This routine performs the multiplication C=A x B.		      */
4798/*								      */
4799/*   res is C, the result.  C may be A and/or B (e.g., X=X*X)	      */
4800/*   lhs is A							      */
4801/*   rhs is B							      */
4802/*   set is the context 					      */
4803/*   status is the usual accumulator				      */
4804/*								      */
4805/* C must have space for set->digits digits.			      */
4806/*								      */
4807/* ------------------------------------------------------------------ */
4808/* 'Classic' multiplication is used rather than Karatsuba, as the     */
4809/* latter would give only a minor improvement for the short numbers   */
4810/* expected to be handled most (and uses much more memory).	      */
4811/*								      */
4812/* There are two major paths here: the general-purpose ('old code')   */
4813/* path which handles all DECDPUN values, and a fastpath version      */
4814/* which is used if 64-bit ints are available, DECDPUN<=4, and more   */
4815/* than two calls to decUnitAddSub would be made.		      */
4816/*								      */
4817/* The fastpath version lumps units together into 8-digit or 9-digit  */
4818/* chunks, and also uses a lazy carry strategy to minimise expensive  */
4819/* 64-bit divisions.  The chunks are then broken apart again into     */
4820/* units for continuing processing.  Despite this overhead, the       */
4821/* fastpath can speed up some 16-digit operations by 10x (and much    */
4822/* more for higher-precision calculations).			      */
4823/*								      */
4824/* A buffer always has to be used for the accumulator; in the	      */
4825/* fastpath, buffers are also always needed for the chunked copies of */
4826/* of the operand coefficients. 				      */
4827/* Static buffers are larger than needed just for multiply, to allow  */
4828/* for calls from other operations (notably exp).		      */
4829/* ------------------------------------------------------------------ */
4830#define FASTMUL(1 && 3<5) (DECUSE641 && DECDPUN3<5)
4831static decNumber * decMultiplyOp(decNumber *res, const decNumber *lhs,
		 const decNumber *rhs, decContext *set,
		 uIntuint32_t *status) {
Intint32_t	 accunits;		   /* Units of accumulator in use */
Intint32_t	 exponent;		   /* work */
Intint32_t	 residue=0;		   /* rounding residue */
uByteuint8_t  bits;			   /* result sign */
Unituint16_t	*acc;			   /* -> accumulator Unit array */
Intint32_t	 needbytes;		   /* size calculator */
void	*allocacc=NULL((void*)0); 	   /* -> allocated accumulator, iff allocated */
Unituint16_t	accbuff[SD2U(DECBUFFER*4+1)(((36*4+1)+3 -1)/3)]; /* buffer (+1 for DECBUFFER==0, */
		   /* *4 for calls from other operations) */
const Unituint16_t *mer, *mermsup;	   /* work */
Intint32_t	madlength;		   /* Units in multiplicand */
Intint32_t	shift;			   /* Units to shift multiplicand by */

#if FASTMUL(1 && 3<5)
  /* if DECDPUN is 1 or 3 work in base 10**9, otherwise */
  /* (DECDPUN is 2 or 4) then work in base 10**8 */
  #if DECDPUN3 & 1		   /* odd */
    #define FASTBASE1000000000 1000000000  /* base */
    #define FASTDIGS9		9  /* digits in base */
    #define FASTLAZY18	       18  /* carry resolution point [1->18] */
  #else
    #define FASTBASE1000000000	100000000
    #define FASTDIGS9		8
    #define FASTLAZY18	     1844  /* carry resolution point [1->1844] */
  #endif
  /* three buffers are used, two for chunked copies of the operands */
  /* (base 10**8 or base 10**9) and one base 2**64 accumulator with */
  /* lazy carry evaluation */
  uIntuint32_t   zlhibuff[(DECBUFFER36*2+1)/8+1]; /* buffer (+1 for DECBUFFER==0) */
  uIntuint32_t  *zlhi=zlhibuff;		  /* -> lhs array */
  uIntuint32_t  *alloclhi=NULL((void*)0);		  /* -> allocated buffer, iff allocated */
  uIntuint32_t   zrhibuff[(DECBUFFER36*2+1)/8+1]; /* buffer (+1 for DECBUFFER==0) */
  uIntuint32_t  *zrhi=zrhibuff;		  /* -> rhs array */
  uIntuint32_t  *allocrhi=NULL((void*)0);		  /* -> allocated buffer, iff allocated */
  uLonguint64_t  zaccbuff[(DECBUFFER36*2+1)/4+2]; /* buffer (+1 for DECBUFFER==0) */
  /* [allocacc is shared for both paths, as only one will run] */
  uLonguint64_t *zacc=zaccbuff;	   /* -> accumulator array for exact result */
  #if DECDPUN3==1
  Intint32_t    zoff;		   /* accumulator offset */
  #endif
  uIntuint32_t  *lip, *rip;		   /* item pointers */
  uIntuint32_t  *lmsi, *rmsi; 	   /* most significant items */
  Intint32_t    ilhs, irhs, iacc;	   /* item counts in the arrays */
  Intint32_t    lazy;		   /* lazy carry counter */
  uLonguint64_t  lcarry;		   /* uLong carry */
  uIntuint32_t   carry;		   /* carry (NB not uLong) */
  Intint32_t    count;		   /* work */
  const  Unituint16_t *cup;		   /* .. */
  Unituint16_t  *up;			   /* .. */
  uLonguint64_t *lp;			   /* .. */
  Intint32_t    p;			   /* .. */
#endif

#if DECSUBSET0
  decNumber *alloclhs=NULL((void*)0);	   /* -> allocated buffer, iff allocated */
  decNumber *allocrhs=NULL((void*)0);	   /* -> allocated buffer, iff allocated */
#endif

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

/* precalculate result sign */
bits=(uByteuint8_t)((lhs->bits^rhs->bits)&DECNEG0x80);

/* handle infinities and NaNs */
if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10))) {		   /* a special bit set */
  if (SPECIALARGS((lhs->bits | rhs->bits) & (0x40|0x20|0x10)) & (DECSNAN0x10 | DECNAN0x20)) { /* one or two NaNs */
    decNaNs(res, lhs, rhs, set, status);
    return res;}
  /* one or two infinities; Infinity * 0 is invalid */
  if (((lhs->bits & DECINF0x40)==0 && ISZERO(lhs)(*(lhs)->lsu==0 && (lhs)->digits==1 && (
((lhs)->bits&(0x40|0x20|0x10))==0)))
    ||((rhs->bits & DECINF0x40)==0 && ISZERO(rhs)(*(rhs)->lsu==0 && (rhs)->digits==1 && (
((rhs)->bits&(0x40|0x20|0x10))==0)))) {
    *status|=DEC_Invalid_operation0x00000080;
    return res;}
  decNumberZero(res);
  res->bits=bits|DECINF0x40;	   /* infinity */
  return res;}

/* For best speed, as in DMSRCN [the original Rexx numerics */
/* module], use the shorter number as the multiplier (rhs) and */
/* the longer as the multiplicand (lhs) to minimise the number of */
/* adds (partial products) */
if (lhs->digits<rhs->digits) {   /* swap... */
  const decNumber *hold=lhs;
  lhs=rhs;
  rhs=hold;
  }

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

  #if FASTMUL(1 && 3<5) 		   /* fastpath can be used */
  /* use the fast path if there are enough digits in the shorter */
  /* operand to make the setup and takedown worthwhile */
  #define NEEDTWO(3*2) (DECDPUN3*2)    /* within two decUnitAddSub calls */
  if (rhs->digits>NEEDTWO(3*2)) {	   /* use fastpath... */
    /* calculate the number of elements in each array */
    ilhs=(lhs->digits+FASTDIGS9-1)/FASTDIGS9; /* [ceiling] */
    irhs=(rhs->digits+FASTDIGS9-1)/FASTDIGS9; /* .. */
    iacc=ilhs+irhs;

    /* allocate buffers if required, as usual */
    needbytes=ilhs*sizeof(uIntuint32_t);
    if (needbytes>(Intint32_t)sizeof(zlhibuff)) {
alloclhi=(uIntuint32_t *)malloc(needbytes);
zlhi=alloclhi;}
    needbytes=irhs*sizeof(uIntuint32_t);
    if (needbytes>(Intint32_t)sizeof(zrhibuff)) {
allocrhi=(uIntuint32_t *)malloc(needbytes);
zrhi=allocrhi;}

    /* Allocating the accumulator space needs a special case when */
    /* DECDPUN=1 because when converting the accumulator to Units */
    /* after the multiplication each 8-byte item becomes 9 1-byte */
    /* units.  Therefore iacc extra bytes are needed at the front */
    /* (rounded up to a multiple of 8 bytes), and the uLong */
    /* accumulator starts offset the appropriate number of units */
    /* to the right to avoid overwrite during the unchunking. */
    needbytes=iacc*sizeof(uLonguint64_t);
    #if DECDPUN3==1
    zoff=(iacc+7)/8;	      /* items to offset by */
    needbytes+=zoff*8;
    #endif
    if (needbytes>(Intint32_t)sizeof(zaccbuff)) {
allocacc=(uLonguint64_t *)malloc(needbytes);
zacc=(uLonguint64_t *)allocacc;}
    if (zlhi==NULL((void*)0)||zrhi==NULL((void*)0)||zacc==NULL((void*)0)) {
*status|=DEC_Insufficient_storage0x00000010;
break;}

    acc=(Unituint16_t *)zacc;       /* -> target Unit array */
    #if DECDPUN3==1
    zacc+=zoff;	      /* start uLong accumulator to right */
    #endif

    /* assemble the chunked copies of the left and right sides */
    for (count=lhs->digits, cup=lhs->lsu, lip=zlhi; count>0; lip++)
for (p=0, *lip=0; p<FASTDIGS9 && count>0;
    p+=DECDPUN3, cup++, count-=DECDPUN3)
 *lip+=*cup*powersDECPOWERS[p];
    lmsi=lip-1;     /* save -> msi */
    for (count=rhs->digits, cup=rhs->lsu, rip=zrhi; count>0; rip++)
for (p=0, *rip=0; p<FASTDIGS9 && count>0;
    p+=DECDPUN3, cup++, count-=DECDPUN3)
 *rip+=*cup*powersDECPOWERS[p];
    rmsi=rip-1;     /* save -> msi */

    /* zero the accumulator */
    for (lp=zacc; lp<zacc+iacc; lp++) *lp=0;

    /* Start the multiplication */
    /* Resolving carries can dominate the cost of accumulating the */
    /* partial products, so this is only done when necessary. */
    /* Each uLong item in the accumulator can hold values up to */
    /* 2**64-1, and each partial product can be as large as */
    /* (10**FASTDIGS-1)**2.  When FASTDIGS=9, this can be added to */
    /* itself 18.4 times in a uLong without overflowing, so during */
    /* the main calculation resolution is carried out every 18th */
    /* add -- every 162 digits.  Similarly, when FASTDIGS=8, the */
    /* partial products can be added to themselves 1844.6 times in */
    /* a uLong without overflowing, so intermediate carry */
    /* resolution occurs only every 14752 digits.  Hence for common */
    /* short numbers usually only the one final carry resolution */
    /* occurs. */
    /* (The count is set via FASTLAZY to simplify experiments to */
    /* measure the value of this approach: a 35% improvement on a */
    /* [34x34] multiply.) */
    lazy=FASTLAZY18;			     /* carry delay count */
    for (rip=zrhi; rip<=rmsi; rip++) {     /* over each item in rhs */
lp=zacc+(rip-zrhi);		     /* where to add the lhs */
for (lip=zlhi; lip<=lmsi; lip++, lp++) { /* over each item in lhs */
 *lp+=(uLonguint64_t)(*lip)*(*rip);	     /* [this should in-line] */
 } /* lip loop */
lazy--;
if (lazy>0 && rip!=rmsi) continue;
lazy=FASTLAZY18;			     /* reset delay count */
/* spin up the accumulator resolving overflows */
for (lp=zacc; lp<zacc+iacc; lp++) {
 if (*lp<FASTBASE1000000000) continue;	     /* it fits */
 lcarry=*lp/FASTBASE1000000000;		     /* top part [slow divide] */
 /* lcarry can exceed 2**32-1, so check again; this check */
 /* and occasional extra divide (slow) is well worth it, as */
 /* it allows FASTLAZY to be increased to 18 rather than 4 */
 /* in the FASTDIGS=9 case */
 if (lcarry<FASTBASE1000000000) carry=(uIntuint32_t)lcarry;  /* [usual] */
  else { /* two-place carry [fairly rare] */
   uIntuint32_t carry2=(uIntuint32_t)(lcarry/FASTBASE1000000000);    /* top top part */
   *(lp+2)+=carry2;			    /* add to item+2 */
   *lp-=((uLonguint64_t)FASTBASE1000000000*FASTBASE1000000000*carry2); /* [slow] */
   carry=(uIntuint32_t)(lcarry-((uLonguint64_t)FASTBASE1000000000*carry2)); /* [inline] */
   }
 *(lp+1)+=carry;		     /* add to item above [inline] */
 *lp-=((uLonguint64_t)FASTBASE1000000000*carry);      /* [inline] */
 } /* carry resolution */
} /* rip loop */

    /* The multiplication is complete; time to convert back into */
    /* units.  This can be done in-place in the accumulator and in */
    /* 32-bit operations, because carries were resolved after the */
    /* final add.  This needs N-1 divides and multiplies for */
    /* each item in the accumulator (which will become up to N */
    /* units, where 2<=N<=9). */
    for (lp=zacc, up=acc; lp<zacc+iacc; lp++) {
uIntuint32_t item=(uIntuint32_t)*lp;		     /* decapitate to uInt */
for (p=0; p<FASTDIGS9-DECDPUN3; p+=DECDPUN3, up++) {
 uIntuint32_t part=item/(DECDPUNMAX999+1);
 *up=(Unituint16_t)(item-(part*(DECDPUNMAX999+1)));
 item=part;
 } /* p */
*up=(Unituint16_t)item; up++;		     /* [final needs no division] */
} /* lp */
    accunits=up-acc;			     /* count of units */
    }
   else { /* here to use units directly, without chunking ['old code'] */
  #endif

    /* if accumulator will be too long for local storage, then allocate */
    acc=accbuff;		   /* -> assume buffer for accumulator */
    needbytes=(D2U(lhs->digits)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3)+D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3))*sizeof(Unituint16_t);
    if (needbytes>(Intint32_t)sizeof(accbuff)) {
allocacc=(Unituint16_t *)malloc(needbytes);
if (allocacc==NULL((void*)0)) {*status|=DEC_Insufficient_storage0x00000010; break;}
acc=(Unituint16_t *)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, 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)((lhs->digits)<=49?d2utable[lhs->digits]:((lhs->digits
)+3 -1)/3);  /* this won't change */
    mermsup=rhs->lsu+D2U(rhs->digits)((rhs->digits)<=49?d2utable[rhs->digits]:((rhs->digits
)+3 -1)/3); /* -> 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; it will be used shortly */
 *(acc+accunits)=0;	   /* [this avoids length of <=0 later] */
 accunits++;
 }
/* multiply multiplicand by 10**DECDPUN for next Unit to left */
shift++;		   /* add this for 'logical length' */
} /* n */
  #if FASTMUL(1 && 3<5)
    } /* unchunked units */
  #endif
  /* common end-path */
  #if DECTRACE0
    decDumpAr('*', acc, accunits);	     /* Show exact result */
  #endif

  /* acc now contains the exact result of the multiplication, */
  /* possibly with a leading zero unit; build the decNumber from */
  /* it, noting if any residue */
  res->bits=bits;			     /* set sign */
  res->digits=decGetDigits(acc, accunits); /* count digits exactly */

  /* There can be 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 there was a wrap, set a */
  /* safe very negative exponent, from which decFinalize() will */
  /* raise a hard underflow shortly. */
  exponent=lhs->exponent+rhs->exponent;    /* calculate exponent */
  if (lhs->exponent<0 && rhs->exponent<0 && exponent>0)
    exponent=-2*DECNUMMAXE999999999;		     /* force underflow */
  res->exponent=exponent;		     /* OK to overwrite now */
