Bug Summary

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

Annotated Source Code

Press '?' to see keyboard shortcuts

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