decNumber.c · samitouri/QOSamiQemu

Code Commits Issues Pulls Sessions CI/CD

master

c 8,322 lines 328 KB

Raw

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

Showing first 5,000 of 8,322 lines. View raw