1 /* ------------------------------------------------------------------ */
2 /* Decimal Number arithmetic module */
3 /* ------------------------------------------------------------------ */
4 /* Copyright (c) IBM Corporation, 2000, 2006. All rights reserved. */
5 /* */
6 /* This software is made available under the terms of the */
7 /* ICU License -- ICU 1.8.1 and later. */
8 /* */
9 /* The description and User's Guide ("The decNumber C Library") for */
10 /* this software is called decNumber.pdf. This document is */
11 /* available, together with arithmetic and format specifications, */
12 /* testcases, and Web links, at: http://www2.hursley.ibm.com/decimal */
13 /* */
14 /* Please send comments, suggestions, and corrections to the author: */
15 /* mfc@uk.ibm.com */
16 /* Mike Cowlishaw, IBM Fellow */
17 /* IBM UK, PO Box 31, Birmingham Road, Warwick CV34 5JL, UK */
18 /* ------------------------------------------------------------------ */
19 /* This module comprises the routines for General Decimal Arithmetic */
20 /* as defined in the specification which may be found on the */
21 /* http://www2.hursley.ibm.com/decimal web pages. It implements both */
22 /* the full ('extended') arithmetic and the simpler ('subset') */
23 /* arithmetic. */
24 /* */
25 /* Usage notes: */
26 /* */
27 /* 1. This code is ANSI C89 except: */
28 /* */
29 /* a) C99 line comments (double forward slash) are used. (Most C */
30 /* compilers accept these. If yours does not, a simple script */
31 /* can be used to convert them to ANSI C comments.) */
32 /* */
33 /* b) Types from C99 stdint.h are used. If you do not have this */
34 /* header file, see the User's Guide section of the decNumber */
35 /* documentation; this lists the necessary definitions. */
36 /* */
37 /* c) If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and */
38 /* uint64_t types may be used. To avoid these, set DECUSE64=0 */
39 /* and DECDPUN<=4 (see documentation). */
40 /* */
41 /* 2. The decNumber format which this library uses is optimized for */
42 /* efficient processing of relatively short numbers; in particular */
43 /* it allows the use of fixed sized structures and minimizes copy */
44 /* and move operations. It does, however, support arbitrary */
45 /* precision (up to 999,999,999 digits) and arbitrary exponent */
46 /* range (Emax in the range 0 through 999,999,999 and Emin in the */
47 /* range -999,999,999 through 0). Mathematical functions (for */
48 /* example decNumberExp) as identified below are restricted more */
49 /* tightly: digits, emax, and -emin in the context must be <= */
50 /* DEC_MAX_MATH (999999), and their operand(s) must be within */
51 /* these bounds. */
52 /* */
53 /* 3. Operands to operator functions are never modified unless they */
54 /* are also specified to be the result number (which is always */
55 /* permitted). Other than that case, operands must not overlap. */
56 /* */
57 /* 4. Error handling: the type of the error is ORed into the status */
58 /* flags in the current context (decContext structure). The */
59 /* SIGFPE signal is then raised if the corresponding trap-enabler */
60 /* flag in the decContext is set (is 1). */
61 /* */
62 /* It is the responsibility of the caller to clear the status */
63 /* flags as required. */
64 /* */
65 /* The result of any routine which returns a number will always */
66 /* be a valid number (which may be a special value, such as an */
67 /* Infinity or NaN). */
68 /* */
69 /* 5. The decNumber format is not an exchangeable concrete */
70 /* representation as it comprises fields which may be machine- */
71 /* dependent (packed or unpacked, or special length, for example). */
72 /* Canonical conversions to and from strings are provided; other */
73 /* conversions are available in separate modules. */
74 /* */
75 /* 6. Normally, input operands are assumed to be valid. Set DECCHECK */
76 /* to 1 for extended operand checking (including NULL operands). */
77 /* Results are undefined if a badly-formed structure (or a NULL */
78 /* pointer to a structure) is provided, though with DECCHECK */
79 /* enabled the operator routines are protected against exceptions. */
80 /* (Except if the result pointer is NULL, which is unrecoverable.) */
81 /* */
82 /* However, the routines will never cause exceptions if they are */
83 /* given well-formed operands, even if the value of the operands */
84 /* is inappropriate for the operation and DECCHECK is not set. */
85 /* (Except for SIGFPE, as and where documented.) */
86 /* */
87 /* 7. Subset arithmetic is available only if DECSUBSET is set to 1. */
88 /* ------------------------------------------------------------------ */
89 /* Implementation notes for maintenance of this module: */
90 /* */
91 /* 1. Storage leak protection: Routines which use malloc are not */
92 /* permitted to use return for fastpath or error exits (i.e., */
93 /* they follow strict structured programming conventions). */
94 /* Instead they have a do{}while(0); construct surrounding the */
95 /* code which is protected -- break may be used to exit this. */
96 /* Other routines can safely use the return statement inline. */
97 /* */
98 /* Storage leak accounting can be enabled using DECALLOC. */
99 /* */
100 /* 2. All loops use the for(;;) construct. Any do construct does */
101 /* not loop; it is for allocation protection as just described. */
102 /* */
103 /* 3. Setting status in the context must always be the very last */
104 /* action in a routine, as non-0 status may raise a trap and hence */
105 /* the call to set status may not return (if the handler uses long */
106 /* jump). Therefore all cleanup must be done first. In general, */
107 /* to achieve this status is accumulated and is only applied just */
108 /* before return by calling decContextSetStatus (via decStatus). */
109 /* */
110 /* Routines which allocate storage cannot, in general, use the */
111 /* 'top level' routines which could cause a non-returning */
112 /* transfer of control. The decXxxxOp routines are safe (do not */
113 /* call decStatus even if traps are set in the context) and should */
114 /* be used instead (they are also a little faster). */
115 /* */
116 /* 4. Exponent checking is minimized by allowing the exponent to */
117 /* grow outside its limits during calculations, provided that */
118 /* the decFinalize function is called later. Multiplication and */
119 /* division, and intermediate calculations in exponentiation, */
120 /* require more careful checks because of the risk of 31-bit */
121 /* overflow (the most negative valid exponent is -1999999997, for */
122 /* a 999999999-digit number with adjusted exponent of -999999999). */
123 /* */
124 /* 5. Rounding is deferred until finalization of results, with any */
125 /* 'off to the right' data being represented as a single digit */
126 /* residue (in the range -1 through 9). This avoids any double- */
127 /* rounding when more than one shortening takes place (for */
128 /* example, when a result is subnormal). */
129 /* */
130 /* 6. The digits count is allowed to rise to a multiple of DECDPUN */
131 /* during many operations, so whole Units are handled and exact */
132 /* accounting of digits is not needed. The correct digits value */
133 /* is found by decGetDigits, which accounts for leading zeros. */
134 /* This must be called before any rounding if the number of digits */
135 /* is not known exactly. */
136 /* */
137 /* 7. The multiply-by-reciprocal 'trick' is used for partitioning */
138 /* numbers up to four digits, using appropriate constants. This */
139 /* is not useful for longer numbers because overflow of 32 bits */
140 /* would lead to 4 multiplies, which is almost as expensive as */
141 /* a divide (unless a floating-point or 64-bit multiply is */
142 /* assumed to be available). */
143 /* */
144 /* 8. Unusual abbreviations that may be used in the commentary: */
145 /* lhs -- left hand side (operand, of an operation) */
146 /* lsd -- least significant digit (of coefficient) */
147 /* lsu -- least significant Unit (of coefficient) */
148 /* msd -- most significant digit (of coefficient) */
149 /* msi -- most significant item (in an array) */
150 /* msu -- most significant Unit (of coefficient) */
151 /* rhs -- right hand side (operand, of an operation) */
152 /* +ve -- positive */
153 /* -ve -- negative */
154 /* ** -- raise to the power */
155 /* ------------------------------------------------------------------ */
156
157 #include <stdlib.h> // for malloc, free, etc.
158 #include <stdio.h> // for printf [if needed]
159 #include <string.h> // for strcpy
160 #include <ctype.h> // for lower
161 #if defined(_MSVC_)
162 #pragma warning(disable:4244) // [for win64]
163 #endif /*defined(_MSVC_)*/
164 #include "decNumber.h" // base number library
165 #include "decNumberLocal.h" // decNumber local types, etc.
166
167 /* Constants */
168 // Public constant array: powers of ten (powers[n]==10**n, 0<=n<=9)
169 const uInt powers[10]={1, 10, 100, 1000, 10000, 100000, 1000000,
170 10000000, 100000000, 1000000000};
171 // Public lookup table used by the D2U macro
172 const uByte d2utable[DECMAXD2U+1]=D2UTABLE;
173
174 // Local constants
175 #define DIVIDE 0x80 // Divide operators
176 #define REMAINDER 0x40 // ..
177 #define DIVIDEINT 0x20 // ..
178 #define REMNEAR 0x10 // ..
179 #define COMPARE 0x01 // Compare operators
180 #define COMPMAX 0x02 // ..
181 #define COMPMIN 0x03 // ..
182 #define COMPTOTAL 0x04 // ..
183 #define COMPNAN 0x05 // .. [NaN processing]
184
185 #define DEC_sNaN 0x40000000 // local status: sNaN signal
186 #define BADINT (Int)0x80000000 // most-negative Int; error indicator
187 // Next two indicate an integer >= 10**6, and its parity (bottom bit)
188 #define BIGEVEN (Int)0x80000002
189 #define BIGODD (Int)0x80000003
190
191 static Unit uarrone[1]={1}; // Unit array of 1, used for incrementing
192
193 /* Granularity-dependent code */
194 #if DECDPUN<=4
195 #define eInt Int // extended integer
196 #define ueInt uInt // unsigned extended integer
197 // Constant multipliers for divide-by-power-of five using reciprocal
198 // multiply, after removing powers of 2 by shifting, and final shift
199 // of 17 [we only need up to **4]
200 static const uInt multies[]={131073, 26215, 5243, 1049, 210};
201 // QUOT10 -- macro to return the quotient of unit u divided by 10**n
202 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
203 #else
204 // For DECDPUN>4 non-ANSI-89 64-bit types are needed.
205 #if !DECUSE64
206 #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4
207 #endif
208 #define eInt Long // extended integer
209 #define ueInt uLong // unsigned extended integer
210 #endif
211
212 /* Local routines */
213 static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *,
214 decContext *, uByte, uInt *);
215 static Flag decBiStr(const char *, const char *, const char *);
216 static uInt decCheckMath(const decNumber *, decContext *, uInt *);
217 static void decApplyRound(decNumber *, decContext *, Int, uInt *);
218 static Int decCompare(const decNumber *lhs, const decNumber *rhs, Flag);
219 static decNumber * decCompareOp(decNumber *, const decNumber *,
220 const decNumber *, decContext *,
221 Flag, uInt *);
222 static void decCopyFit(decNumber *, const decNumber *, decContext *,
223 Int *, uInt *);
224 static decNumber * decDivideOp(decNumber *, const decNumber *,
225 const decNumber *, decContext *, Flag, uInt *);
226 static decNumber * decExpOp(decNumber *, const decNumber *,
227 decContext *, uInt *);
228 static void decFinalize(decNumber *, decContext *, Int *, uInt *);
229 static Int decGetDigits(Unit *, Int);
230 static Int decGetInt(const decNumber *);
231 static decNumber * decLnOp(decNumber *, const decNumber *,
232 decContext *, uInt *);
233 static decNumber * decMultiplyOp(decNumber *, const decNumber *,
234 const decNumber *, decContext *,
235 uInt *);
236 static decNumber * decNaNs(decNumber *, const decNumber *,
237 const decNumber *, uInt *);
238 static decNumber * decPutInt(decNumber *, Int);
239 static decNumber * decQuantizeOp(decNumber *, const decNumber *,
240 const decNumber *, decContext *, Flag,
241 uInt *);
242 static void decSetCoeff(decNumber *, decContext *, const Unit *,
243 Int, Int *, uInt *);
244 static void decSetOverflow(decNumber *, decContext *, uInt *);
245 static void decSetSubnormal(decNumber *, decContext *, Int *, uInt *);
246 static Int decShiftToLeast(Unit *, Int, Int);
247 static Int decShiftToMost(Unit *, Int, Int);
248 static void decStatus(decNumber *, uInt, decContext *);
249 static void decToString(const decNumber *, char[], Flag);
250 static decNumber * decTrim(decNumber *, Flag, Int *);
251 static Int decUnitAddSub(const Unit *, Int, const Unit *, Int, Int,
252 Unit *, Int);
253 static Int decUnitCompare(const Unit *, Int, const Unit *, Int, Int);
254
255 #if !DECSUBSET
256 /* decFinish == decFinalize when no subset arithmetic needed */
257 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
258 #else
259 static void decFinish(decNumber *, decContext *, Int *, uInt *);
260 static decNumber * decRoundOperand(const decNumber *, decContext *, uInt *);
261 #endif
262
263 /* Local macros */
264 // masked special-values bits
265 #define SPECIALARG (rhs->bits & DECSPECIAL)
266 #define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)
267
268 /* Diagnostic macros, etc. */
269 #if DECALLOC
270 // Handle malloc/free accounting. If enabled, our accountable routines
271 // are used; otherwise the code just goes straight to the system malloc
272 // and free routines.
273 #define malloc(a) decMalloc(a)
274 #define free(a) decFree(a)
275 #define DECFENCE 0x5a // corruption detector
276 // 'Our' malloc and free:
277 static void *decMalloc(size_t);
278 static void decFree(void *);
279 uInt decAllocBytes=0; // count of bytes allocated
280 // Note that DECALLOC code only checks for storage buffer overflow.
281 // To check for memory leaks, the decAllocBytes variable must be
282 // checked to be 0 at appropriate times (e.g., after the test
283 // harness completes a set of tests). This checking may be unreliable
284 // if the testing is done in a multi-thread environment.
285 #endif
286
287 #if DECCHECK
288 // Optional checking routines. Enabling these means that decNumber
289 // and decContext operands to operator routines are checked for
290 // correctness. This roughly doubles the execution time of the
291 // fastest routines (and adds 600+ bytes), so should not normally be
292 // used in 'production'.
293 // decCheckInexact is used to check that inexact results have a full
294 // complement of digits (where appropriate -- this is not the case
295 // for Quantize, for example)
296 #define DECUNUSED (void *)(0xffffffff)
297 static Flag decCheckOperands(decNumber *, const decNumber *,
298 const decNumber *, decContext *);
299 static Flag decCheckNumber(const decNumber *, decContext *);
300 static void decCheckInexact(const decNumber *, decContext *);
301 #endif
302
303 #if DECTRACE || DECCHECK
304 // Optional trace/debugging routines (may or may not be used)
305 void decNumberShow(const decNumber *); // displays the components of a number
306 static void decDumpAr(char, const Unit *, Int);
307 #endif
308
309 /* ================================================================== */
310 /* Conversions */
311 /* ================================================================== */
312
313 /* ------------------------------------------------------------------ */
314 /* to-scientific-string -- conversion to numeric string */
315 /* to-engineering-string -- conversion to numeric string */
316 /* */
317 /* decNumberToString(dn, string); */
318 /* decNumberToEngString(dn, string); */
319 /* */
320 /* dn is the decNumber to convert */
321 /* string is the string where the result will be laid out */
322 /* */
323 /* string must be at least dn->digits+14 characters long */
324 /* */
325 /* No error is possible, and no status can be set. */
326 /* ------------------------------------------------------------------ */
decNumberToString(const decNumber * dn,char * string)327 char * decNumberToString(const decNumber *dn, char *string){
328 decToString(dn, string, 0);
329 return string;
330 } // DecNumberToString
331
decNumberToEngString(const decNumber * dn,char * string)332 char * decNumberToEngString(const decNumber *dn, char *string){
333 decToString(dn, string, 1);
334 return string;
335 } // DecNumberToEngString
336
337 /* ------------------------------------------------------------------ */
338 /* to-number -- conversion from numeric string */
339 /* */
340 /* decNumberFromString -- convert string to decNumber */
341 /* dn -- the number structure to fill */
342 /* chars[] -- the string to convert ('\0' terminated) */
343 /* set -- the context used for processing any error, */
344 /* determining the maximum precision available */
345 /* (set.digits), determining the maximum and minimum */
346 /* exponent (set.emax and set.emin), determining if */
347 /* extended values are allowed, and checking the */
348 /* rounding mode if overflow occurs or rounding is */
349 /* needed. */
350 /* */
351 /* The length of the coefficient and the size of the exponent are */
352 /* checked by this routine, so the correct error (Underflow or */
353 /* Overflow) can be reported or rounding applied, as necessary. */
354 /* */
355 /* If bad syntax is detected, the result will be a quiet NaN. */
356 /* ------------------------------------------------------------------ */
decNumberFromString(decNumber * dn,const char chars[],decContext * set)357 decNumber * decNumberFromString(decNumber *dn, const char chars[],
358 decContext *set) {
359 Int exponent=0; // working exponent [assume 0]
360 uByte bits=0; // working flags [assume +ve]
361 Unit *res; // where result will be built
362 Unit resbuff[SD2U(DECBUFFER+1)];// local buffer in case need temporary
363 Unit *allocres=NULL; // -> allocated result, iff allocated
364 Int d=0; // count of digits found in decimal part
365 const char *dotchar=NULL; // where dot was found
366 const char *cfirst=chars; // -> first character of decimal part
367 const char *last=NULL; // -> last digit of decimal part
368 const char *c; // work
369 Unit *up; // ..
370 #if DECDPUN>1
371 Int cut, out; // ..
372 #endif
373 Int residue; // rounding residue
374 uInt status=0; // error code
375
376 #if DECCHECK
377 if (decCheckOperands(DECUNUSED, DECUNUSED, DECUNUSED, set))
378 return decNumberZero(dn);
379 #endif
380
381 do { // status & malloc protection
382 for (c=chars;; c++) { // -> input character
383 if (*c>='0' && *c<='9') { // test for Arabic digit
384 last=c;
385 d++; // count of real digits
386 continue; // still in decimal part
387 }
388 if (*c=='.' && dotchar==NULL) { // first '.'
389 dotchar=c; // record offset into decimal part
390 if (c==cfirst) cfirst++; // first digit must follow
391 continue;}
392 if (c==chars) { // first in string...
393 if (*c=='-') { // valid - sign
394 cfirst++;
395 bits=DECNEG;
396 continue;}
397 if (*c=='+') { // valid + sign
398 cfirst++;
399 continue;}
400 }
401 // *c is not a digit, or a valid +, -, or '.'
402 break;
403 } // c
404
405 if (last==NULL) { // no digits yet
406 status=DEC_Conversion_syntax;// assume the worst
407 if (*c=='\0') break; // and no more to come...
408 #if DECSUBSET
409 // if subset then infinities and NaNs are not allowed
410 if (!set->extended) break; // hopeless
411 #endif
412 // Infinities and NaNs are possible, here
413 if (dotchar!=NULL) break; // .. unless had a dot
414 decNumberZero(dn); // be optimistic
415 if (decBiStr(c, "infinity", "INFINITY")
416 || decBiStr(c, "inf", "INF")) {
417 dn->bits=bits | DECINF;
418 status=0; // is OK
419 break; // all done
420 }
421 // a NaN expected
422 // 2003.09.10 NaNs are now permitted to have a sign
423 dn->bits=bits | DECNAN; // assume simple NaN
424 if (*c=='s' || *c=='S') { // looks like an sNaN
425 c++;
426 dn->bits=bits | DECSNAN;
427 }
428 if (*c!='n' && *c!='N') break; // check caseless "NaN"
429 c++;
430 if (*c!='a' && *c!='A') break; // ..
431 c++;
432 if (*c!='n' && *c!='N') break; // ..
433 c++;
434 // now either nothing, or nnnn payload, expected
435 // -> start of integer and skip leading 0s [including plain 0]
436 for (cfirst=c; *cfirst=='0';) cfirst++;
437 if (*cfirst=='\0') { // "NaN" or "sNaN", maybe with all 0s
438 status=0; // it's good
439 break; // ..
440 }
441 // something other than 0s; setup last and d as usual [no dots]
442 for (c=cfirst;; c++, d++) {
443 if (*c<'0' || *c>'9') break; // test for Arabic digit
444 last=c;
445 }
446 if (*c!='\0') break; // not all digits
447 if (d>set->digits-1) {
448 // [NB: payload in a decNumber can be full length unless
449 // clamped, in which case can only be digits-1]
450 if (set->clamp) break;
451 if (d>set->digits) break;
452 } // too many digits?
453 // good; drop through to convert the integer to coefficient
454 status=0; // syntax is OK
455 bits=dn->bits; // for copy-back
456 } // last==NULL
457
458 else if (*c!='\0') { // more to process...
459 // had some digits; exponent is only valid sequence now
460 Flag nege; // 1=negative exponent
461 const char *firstexp; // -> first significant exponent digit
462 status=DEC_Conversion_syntax;// assume the worst
463 if (*c!='e' && *c!='E') break;
464 /* Found 'e' or 'E' -- now process explicit exponent */
465 // 1998.07.11: sign no longer required
466 nege=0;
467 c++; // to (possible) sign
468 if (*c=='-') {nege=1; c++;}
469 else if (*c=='+') c++;
470 if (*c=='\0') break;
471
472 for (; *c=='0' && *(c+1)!='\0';) c++; // strip insignificant zeros
473 firstexp=c; // save exponent digit place
474 for (; ;c++) {
475 if (*c<'0' || *c>'9') break; // not a digit
476 exponent=X10(exponent)+(Int)*c-(Int)'0';
477 } // c
478 // if not now on a '\0', *c must not be a digit
479 if (*c!='\0') break;
480
481 // (this next test must be after the syntax checks)
482 // if it was too long the exponent may have wrapped, so check
483 // carefully and set it to a certain overflow if wrap possible
484 if (c>=firstexp+9+1) {
485 if (c>firstexp+9+1 || *firstexp>'1') exponent=DECNUMMAXE*2;
486 // [up to 1999999999 is OK, for example 1E-1000000998]
487 }
488 if (nege) exponent=-exponent; // was negative
489 status=0; // is OK
490 } // stuff after digits
491
492 // Here when whole string has been inspected; syntax is good
493 // cfirst->first digit (never dot), last->last digit (ditto)
494
495 // strip leading zeros/dot [leave final 0 if all 0's]
496 if (*cfirst=='0') { // [cfirst has stepped over .]
497 for (c=cfirst; c<last; c++, cfirst++) {
498 if (*c=='.') continue; // ignore dots
499 if (*c!='0') break; // non-zero found
500 d--; // 0 stripped
501 } // c
502 #if DECSUBSET
503 // make a rapid exit for easy zeros if !extended
504 if (*cfirst=='0' && !set->extended) {
505 decNumberZero(dn); // clean result
506 break; // [could be return]
507 }
508 #endif
509 } // at least one leading 0
510
511 // Handle decimal point...
512 if (dotchar!=NULL && dotchar<last) // non-trailing '.' found?
513 exponent-=(last-dotchar); // adjust exponent
514 // [we can now ignore the .]
515
516 // OK, the digits string is good. Assemble in the decNumber, or in
517 // a temporary units array if rounding is needed
518 if (d<=set->digits) res=dn->lsu; // fits into supplied decNumber
519 else { // rounding needed
520 Int needbytes=D2U(d)*sizeof(Unit);// bytes needed
521 res=resbuff; // assume use local buffer
522 if (needbytes>(Int)sizeof(resbuff)) { // too big for local
523 allocres=(Unit *)malloc(needbytes);
524 if (allocres==NULL) {status|=DEC_Insufficient_storage; break;}
525 res=allocres;
526 }
527 }
528 // res now -> number lsu, buffer, or allocated storage for Unit array
529
530 // Place the coefficient into the selected Unit array
531 // [this is often 70% of the cost of this function when DECDPUN>1]
532 #if DECDPUN>1
533 out=0; // accumulator
534 up=res+D2U(d)-1; // -> msu
535 cut=d-(up-res)*DECDPUN; // digits in top unit
536 for (c=cfirst;; c++) { // along the digits
537 if (*c=='.') continue; // ignore '.' [don't decrement cut]
538 out=X10(out)+(Int)*c-(Int)'0';
539 if (c==last) break; // done [never get to trailing '.']
540 cut--;
541 if (cut>0) continue; // more for this unit
542 *up=(Unit)out; // write unit
543 up--; // prepare for unit below..
544 cut=DECDPUN; // ..
545 out=0; // ..
546 } // c
547 *up=(Unit)out; // write lsu
548
549 #else
550 // DECDPUN==1
551 up=res; // -> lsu
552 for (c=last; c>=cfirst; c--) { // over each character, from least
553 if (*c=='.') continue; // ignore . [don't step up]
554 *up=(Unit)((Int)*c-(Int)'0');
555 up++;
556 } // c
557 #endif
558
559 dn->bits=bits;
560 dn->exponent=exponent;
561 dn->digits=d;
562
563 // if not in number (too long) shorten into the number
564 if (d>set->digits) {
565 residue=0;
566 decSetCoeff(dn, set, res, d, &residue, &status);
567 // always check for overflow or subnormal and round as needed
568 decFinalize(dn, set, &residue, &status);
569 }
570 else { // no rounding, but may still have overflow or subnormal
571 // [these tests are just for performance; finalize repeats them]
572 if ((dn->exponent-1<set->emin-dn->digits)
573 || (dn->exponent-1>set->emax-set->digits)) {
574 residue=0;
575 decFinalize(dn, set, &residue, &status);
576 }
577 }
578 // decNumberShow(dn);
579 } while(0); // [for break]
580
581 if (allocres!=NULL) free(allocres); // drop any storage used
582 if (status!=0) decStatus(dn, status, set);
583 return dn;
584 } /* decNumberFromString */
585
586 /* ================================================================== */
587 /* Operators */
588 /* ================================================================== */
589
590 /* ------------------------------------------------------------------ */
591 /* decNumberAbs -- absolute value operator */
592 /* */
593 /* This computes C = abs(A) */
594 /* */
595 /* res is C, the result. C may be A */
596 /* rhs is A */
597 /* set is the context */
598 /* */
599 /* C must have space for set->digits digits. */
600 /* ------------------------------------------------------------------ */
601 /* This has the same effect as decNumberPlus unless A is negative, */
602 /* in which case it has the same effect as decNumberMinus. */
603 /* ------------------------------------------------------------------ */
decNumberAbs(decNumber * res,const decNumber * rhs,decContext * set)604 decNumber * decNumberAbs(decNumber *res, const decNumber *rhs,
605 decContext *set) {
606 decNumber dzero; // for 0
607 uInt status=0; // accumulator
608
609 #if DECCHECK
610 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
611 #endif
612
613 decNumberZero(&dzero); // set 0
614 dzero.exponent=rhs->exponent; // [no coefficient expansion]
615 decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status);
616 if (status!=0) decStatus(res, status, set);
617 #if DECCHECK
618 decCheckInexact(res, set);
619 #endif
620 return res;
621 } // decNumberAbs
622
623 /* ------------------------------------------------------------------ */
624 /* decNumberAdd -- add two Numbers */
625 /* */
626 /* This computes C = A + B */
627 /* */
628 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
629 /* lhs is A */
630 /* rhs is B */
631 /* set is the context */
632 /* */
633 /* C must have space for set->digits digits. */
634 /* ------------------------------------------------------------------ */
635 /* This just calls the routine shared with Subtract */
decNumberAdd(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)636 decNumber * decNumberAdd(decNumber *res, const decNumber *lhs,
637 const decNumber *rhs, decContext *set) {
638 uInt status=0; // accumulator
639 decAddOp(res, lhs, rhs, set, 0, &status);
640 if (status!=0) decStatus(res, status, set);
641 #if DECCHECK
642 decCheckInexact(res, set);
643 #endif
644 return res;
645 } // decNumberAdd
646
647 /* ------------------------------------------------------------------ */
648 /* decNumberCompare -- compare two Numbers */
649 /* */
650 /* This computes C = A ? B */
651 /* */
652 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
653 /* lhs is A */
654 /* rhs is B */
655 /* set is the context */
656 /* */
657 /* C must have space for one digit (or NaN). */
658 /* ------------------------------------------------------------------ */
decNumberCompare(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)659 decNumber * decNumberCompare(decNumber *res, const decNumber *lhs,
660 const decNumber *rhs, decContext *set) {
661 uInt status=0; // accumulator
662 decCompareOp(res, lhs, rhs, set, COMPARE, &status);
663 if (status!=0) decStatus(res, status, set);
664 return res;
665 } // decNumberCompare
666
667 /* ------------------------------------------------------------------ */
668 /* decNumberCompareTotal -- compare two Numbers, using total ordering */
669 /* */
670 /* This computes C = A ? B, under total ordering */
671 /* */
672 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
673 /* lhs is A */
674 /* rhs is B */
675 /* set is the context */
676 /* */
677 /* C must have space for one digit; the result will always be one of */
678 /* -1, 0, or 1. */
679 /* ------------------------------------------------------------------ */
decNumberCompareTotal(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)680 decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs,
681 const decNumber *rhs, decContext *set) {
682 uInt status=0; // accumulator
683 decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
684 if (status!=0) decStatus(res, status, set);
685 return res;
686 } // decNumberCompareTotal
687
688 /* ------------------------------------------------------------------ */
689 /* decNumberDivide -- divide one number by another */
690 /* */
691 /* This computes C = A / B */
692 /* */
693 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
694 /* lhs is A */
695 /* rhs is B */
696 /* set is the context */
697 /* */
698 /* C must have space for set->digits digits. */
699 /* ------------------------------------------------------------------ */
decNumberDivide(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)700 decNumber * decNumberDivide(decNumber *res, const decNumber *lhs,
701 const decNumber *rhs, decContext *set) {
702 uInt status=0; // accumulator
703 decDivideOp(res, lhs, rhs, set, DIVIDE, &status);
704 if (status!=0) decStatus(res, status, set);
705 #if DECCHECK
706 decCheckInexact(res, set);
707 #endif
708 return res;
709 } // decNumberDivide
710
711 /* ------------------------------------------------------------------ */
712 /* decNumberDivideInteger -- divide and return integer quotient */
713 /* */
714 /* This computes C = A # B, where # is the integer divide operator */
715 /* */
716 /* res is C, the result. C may be A and/or B (e.g., X=X#X) */
717 /* lhs is A */
718 /* rhs is B */
719 /* set is the context */
720 /* */
721 /* C must have space for set->digits digits. */
722 /* ------------------------------------------------------------------ */
decNumberDivideInteger(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)723 decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs,
724 const decNumber *rhs, decContext *set) {
725 uInt status=0; // accumulator
726 decDivideOp(res, lhs, rhs, set, DIVIDEINT, &status);
727 if (status!=0) decStatus(res, status, set);
728 return res;
729 } // decNumberDivideInteger
730
731 /* ------------------------------------------------------------------ */
732 /* decNumberExp -- exponentiation */
733 /* */
734 /* This computes C = exp(A) */
735 /* */
736 /* res is C, the result. C may be A */
737 /* rhs is A */
738 /* set is the context; note that rounding mode has no effect */
739 /* */
740 /* C must have space for set->digits digits. */
741 /* */
742 /* Mathematical function restrictions apply (see above); a NaN is */
743 /* returned with Invalid_operation if a restriction is violated. */
744 /* */
745 /* Finite results will always be full precision and Inexact, except */
746 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
747 /* */
748 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
749 /* almost always be correctly rounded, but may be up to 1 ulp in */
750 /* error in rare cases. */
751 /* ------------------------------------------------------------------ */
752 /* This is a wrapper for decExpOp which can handle the slightly wider */
753 /* (double) range needed by Ln (which has to be able to calculate */
754 /* exp(-a) where a can be the tiniest number (Ntiny). */
755 /* ------------------------------------------------------------------ */
decNumberExp(decNumber * res,const decNumber * rhs,decContext * set)756 decNumber * decNumberExp(decNumber *res, const decNumber *rhs,
757 decContext *set) {
758 uInt status=0; // accumulator
759 #if DECSUBSET
760 decNumber *allocrhs=NULL; // non-NULL if rounded rhs allocated
761 #endif
762
763 #if DECCHECK
764 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
765 #endif
766
767 // Check restrictions; these restrictions ensure that if h=8 (see
768 // decExpOp) then the result will either overflow or underflow to 0.
769 // Other math functions restrict the input range, too, for inverses.
770 // If not violated then carry out the operation.
771 if (!decCheckMath(rhs, set, &status)) do { // protect allocation
772 #if DECSUBSET
773 if (!set->extended) {
774 // reduce operand and set lostDigits status, as needed
775 if (rhs->digits>set->digits) {
776 allocrhs=decRoundOperand(rhs, set, &status);
777 if (allocrhs==NULL) break;
778 rhs=allocrhs;
779 }
780 }
781 #endif
782 decExpOp(res, rhs, set, &status);
783 } while(0); // end protected
784
785 #if DECSUBSET
786 if (allocrhs !=NULL) free(allocrhs); // drop any storage used
787 #endif
788 // apply significant status
789 if (status!=0) decStatus(res, status, set);
790 #if DECCHECK
791 decCheckInexact(res, set);
792 #endif
793 return res;
794 } // decNumberExp
795
796 /* ------------------------------------------------------------------ */
797 /* decNumberLn -- natural logarithm */
798 /* */
799 /* This computes C = ln(A) */
800 /* */
801 /* res is C, the result. C may be A */
802 /* rhs is A */
803 /* set is the context; note that rounding mode has no effect */
804 /* */
805 /* C must have space for set->digits digits. */
806 /* */
807 /* Notable cases: */
808 /* A<0 -> Invalid */
809 /* A=0 -> -Infinity (Exact) */
810 /* A=+Infinity -> +Infinity (Exact) */
811 /* A=1 exactly -> 0 (Exact) */
812 /* */
813 /* Mathematical function restrictions apply (see above); a NaN is */
814 /* returned with Invalid_operation if a restriction is violated. */
815 /* */
816 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
817 /* almost always be correctly rounded, but may be up to 1 ulp in */
818 /* error in rare cases. */
819 /* ------------------------------------------------------------------ */
820 /* This is a wrapper for decLnOp which can handle the slightly wider */
821 /* (+11) range needed by Ln, Log10, etc. (which may have to be able */
822 /* to calculate at p+e+2). */
823 /* ------------------------------------------------------------------ */
decNumberLn(decNumber * res,const decNumber * rhs,decContext * set)824 decNumber * decNumberLn(decNumber *res, const decNumber *rhs,
825 decContext *set) {
826 uInt status=0; // accumulator
827 #if DECSUBSET
828 decNumber *allocrhs=NULL; // non-NULL if rounded rhs allocated
829 #endif
830
831 #if DECCHECK
832 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
833 #endif
834
835 // Check restrictions; this is a math function; if not violated
836 // then carry out the operation.
837 if (!decCheckMath(rhs, set, &status)) do { // protect allocation
838 #if DECSUBSET
839 if (!set->extended) {
840 // reduce operand and set lostDigits status, as needed
841 if (rhs->digits>set->digits) {
842 allocrhs=decRoundOperand(rhs, set, &status);
843 if (allocrhs==NULL) break;
844 rhs=allocrhs;
845 }
846 // special check in subset for rhs=0
847 if (ISZERO(rhs)) { // +/- zeros -> error
848 status|=DEC_Invalid_operation;
849 break;}
850 } // extended=0
851 #endif
852 decLnOp(res, rhs, set, &status);
853 } while(0); // end protected
854
855 #if DECSUBSET
856 if (allocrhs !=NULL) free(allocrhs); // drop any storage used
857 #endif
858 // apply significant status
859 if (status!=0) decStatus(res, status, set);
860 #if DECCHECK
861 decCheckInexact(res, set);
862 #endif
863 return res;
864 } // decNumberLn
865
866 /* ------------------------------------------------------------------ */
867 /* decNumberLog10 -- logarithm in base 10 */
868 /* */
869 /* This computes C = log10(A) */
870 /* */
871 /* res is C, the result. C may be A */
872 /* rhs is A */
873 /* set is the context; note that rounding mode has no effect */
874 /* */
875 /* C must have space for set->digits digits. */
876 /* */
877 /* Notable cases: */
878 /* A<0 -> Invalid */
879 /* A=0 -> -Infinity (Exact) */
880 /* A=+Infinity -> +Infinity (Exact) */
881 /* A=10**n (if n is an integer) -> n (Exact) */
882 /* */
883 /* Mathematical function restrictions apply (see above); a NaN is */
884 /* returned with Invalid_operation if a restriction is violated. */
885 /* */
886 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
887 /* almost always be correctly rounded, but may be up to 1 ulp in */
888 /* error in rare cases. */
889 /* ------------------------------------------------------------------ */
890 /* This calculates ln(A)/ln(10) using appropriate precision. For */
891 /* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the */
892 /* requested digits and t is the number of digits in the exponent */
893 /* (maximum 6). For ln(10) it is p + 3; this is often handled by the */
894 /* fastpath in decLnOp. The final division is done to the requested */
895 /* precision. */
896 /* ------------------------------------------------------------------ */
decNumberLog10(decNumber * res,const decNumber * rhs,decContext * set)897 decNumber * decNumberLog10(decNumber *res, const decNumber *rhs,
898 decContext *set) {
899 uInt status=0, ignore=0; // status accumulators
900 uInt needbytes; // for space calculations
901 Int p; // working precision
902 Int t; // digits in exponent of A
903
904 // buffers for a and b working decimals
905 // (adjustment calculator, same size)
906 decNumber bufa[D2N(DECBUFFER+2)];
907 decNumber *allocbufa=NULL; // -> allocated bufa, iff allocated
908 decNumber *a=bufa; // temporary a
909 decNumber bufb[D2N(DECBUFFER+2)];
910 decNumber *allocbufb=NULL; // -> allocated bufa, iff allocated
911 decNumber *b=bufb; // temporary b
912 decNumber bufw[D2N(10)]; // working 2-10 digit number
913 decNumber *w=bufw; // ..
914 #if DECSUBSET
915 decNumber *allocrhs=NULL; // non-NULL if rounded rhs allocated
916 #endif
917
918 decContext aset; // working context
919
920 #if DECCHECK
921 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
922 #endif
923
924 // Check restrictions; this is a math function; if not violated
925 // then carry out the operation.
926 if (!decCheckMath(rhs, set, &status)) do { // protect malloc
927 #if DECSUBSET
928 if (!set->extended) {
929 // reduce operand and set lostDigits status, as needed
930 if (rhs->digits>set->digits) {
931 allocrhs=decRoundOperand(rhs, set, &status);
932 if (allocrhs==NULL) break;
933 rhs=allocrhs;
934 }
935 // special check in subset for rhs=0
936 if (ISZERO(rhs)) { // +/- zeros -> error
937 status|=DEC_Invalid_operation;
938 break;}
939 } // extended=0
940 #endif
941
942 decContextDefault(&aset, DEC_INIT_DECIMAL64); // clean context
943
944 // handle exact powers of 10; only check if +ve finite
945 if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) {
946 Int residue=0; // (no residue)
947 uInt copystat=0; // clean status
948
949 // round to a single digit...
950 aset.digits=1;
951 decCopyFit(w, rhs, &aset, &residue, ©stat); // copy & shorten
952 // if exact and the digit is 1, rhs is a power of 10
953 if (!(copystat&DEC_Inexact) && w->lsu[0]==1) {
954 // the exponent, conveniently, is the power of 10; making
955 // this the result needs a little care as it might not fit,
956 // so first convert it into the working number, and then move
957 // to res
958 decPutInt(w, w->exponent);
959 residue=0;
960 decCopyFit(res, w, set, &residue, &status); // copy & round
961 decFinish(res, set, &residue, &status); // cleanup/set flags
962 break;
963 } // not a power of 10
964 } // not a candidate for exact
965
966 // simplify the information-content calculation to use 'total
967 // number of digits in a, including exponent' as compared to the
968 // requested digits, as increasing this will only rarely cost an
969 // iteration in ln(a) anyway
970 t=6; // it can never be >6
971
972 // allocate space when needed...
973 p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3;
974 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
975 if (needbytes>sizeof(bufa)) { // need malloc space
976 allocbufa=(decNumber *)malloc(needbytes);
977 if (allocbufa==NULL) { // hopeless -- abandon
978 status|=DEC_Insufficient_storage;
979 break;}
980 a=allocbufa; // use the allocated space
981 }
982 aset.digits=p; // as calculated
983 aset.emax=DEC_MAX_MATH; // usual bounds
984 aset.emin=-DEC_MAX_MATH; // ..
985 aset.clamp=0; // and no concrete format
986 decLnOp(a, rhs, &aset, &status); // a=ln(rhs)
987
988 // skip the division if the result so far is infinite, NaN, or
989 // zero, or there was an error; note NaN from sNaN needs copy
990 if (status&DEC_NaNs && !(status&DEC_sNaN)) break;
991 if (a->bits&DECSPECIAL || ISZERO(a)) {
992 decNumberCopy(res, a); // [will fit]
993 break;}
994
995 // for ln(10) an extra 3 digits of precision are needed
996 p=set->digits+3;
997 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
998 if (needbytes>sizeof(bufb)) { // need malloc space
999 allocbufb=(decNumber *)malloc(needbytes);
1000 if (allocbufb==NULL) { // hopeless -- abandon
1001 status|=DEC_Insufficient_storage;
1002 break;}
1003 b=allocbufb; // use the allocated space
1004 }
1005 decNumberZero(w); // set up 10...
1006 #if DECDPUN==1
1007 w->lsu[1]=1; w->lsu[0]=0; // ..
1008 #else
1009 w->lsu[0]=10; // ..
1010 #endif
1011 w->digits=2; // ..
1012
1013 aset.digits=p;
1014 decLnOp(b, w, &aset, &ignore); // b=ln(10)
1015
1016 aset.digits=set->digits; // for final divide
1017 decDivideOp(res, a, b, &aset, DIVIDE, &status); // into result
1018 } while(0); // [for break]
1019
1020 if (allocbufa!=NULL) free(allocbufa); // drop any storage used
1021 if (allocbufb!=NULL) free(allocbufb); // ..
1022 #if DECSUBSET
1023 if (allocrhs !=NULL) free(allocrhs); // ..
1024 #endif
1025 // apply significant status
1026 if (status!=0) decStatus(res, status, set);
1027 #if DECCHECK
1028 decCheckInexact(res, set);
1029 #endif
1030 return res;
1031 } // decNumberLog10
1032
1033 /* ------------------------------------------------------------------ */
1034 /* decNumberMax -- compare two Numbers and return the maximum */
1035 /* */
1036 /* This computes C = A ? B, returning the maximum or A if equal */
1037 /* */
1038 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1039 /* lhs is A */
1040 /* rhs is B */
1041 /* set is the context */
1042 /* */
1043 /* C must have space for set->digits digits. */
1044 /* ------------------------------------------------------------------ */
decNumberMax(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1045 decNumber * decNumberMax(decNumber *res, const decNumber *lhs,
1046 const decNumber *rhs, decContext *set) {
1047 uInt status=0; // accumulator
1048 decCompareOp(res, lhs, rhs, set, COMPMAX, &status);
1049 if (status!=0) decStatus(res, status, set);
1050 #if DECCHECK
1051 decCheckInexact(res, set);
1052 #endif
1053 return res;
1054 } // decNumberMax
1055
1056 /* ------------------------------------------------------------------ */
1057 /* decNumberMin -- compare two Numbers and return the minimum */
1058 /* */
1059 /* This computes C = A ? B, returning the minimum or A if equal */
1060 /* */
1061 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1062 /* lhs is A */
1063 /* rhs is B */
1064 /* set is the context */
1065 /* */
1066 /* C must have space for set->digits digits. */
1067 /* ------------------------------------------------------------------ */
decNumberMin(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1068 decNumber * decNumberMin(decNumber *res, const decNumber *lhs,
1069 const decNumber *rhs, decContext *set) {
1070 uInt status=0; // accumulator
1071 decCompareOp(res, lhs, rhs, set, COMPMIN, &status);
1072 if (status!=0) decStatus(res, status, set);
1073 #if DECCHECK
1074 decCheckInexact(res, set);
1075 #endif
1076 return res;
1077 } // decNumberMin
1078
1079 /* ------------------------------------------------------------------ */
1080 /* decNumberMinus -- prefix minus operator */
1081 /* */
1082 /* This computes C = 0 - A */
1083 /* */
1084 /* res is C, the result. C may be A */
1085 /* rhs is A */
1086 /* set is the context */
1087 /* */
1088 /* C must have space for set->digits digits. */
1089 /* ------------------------------------------------------------------ */
1090 /* Simply use AddOp for the subtract, which will do the necessary. */
1091 /* ------------------------------------------------------------------ */
decNumberMinus(decNumber * res,const decNumber * rhs,decContext * set)1092 decNumber * decNumberMinus(decNumber *res, const decNumber *rhs,
1093 decContext *set) {
1094 decNumber dzero;
1095 uInt status=0; // accumulator
1096
1097 #if DECCHECK
1098 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1099 #endif
1100
1101 decNumberZero(&dzero); // make 0
1102 dzero.exponent=rhs->exponent; // [no coefficient expansion]
1103 decAddOp(res, &dzero, rhs, set, DECNEG, &status);
1104 if (status!=0) decStatus(res, status, set);
1105 #if DECCHECK
1106 decCheckInexact(res, set);
1107 #endif
1108 return res;
1109 } // decNumberMinus
1110
1111 /* ------------------------------------------------------------------ */
1112 /* decNumberPlus -- prefix plus operator */
1113 /* */
1114 /* This computes C = 0 + A */
1115 /* */
1116 /* res is C, the result. C may be A */
1117 /* rhs is A */
1118 /* set is the context */
1119 /* */
1120 /* C must have space for set->digits digits. */
1121 /* ------------------------------------------------------------------ */
1122 /* This simply uses AddOp; Add will take fast path after preparing A. */
1123 /* Performance is a concern here, as this routine is often used to */
1124 /* check operands and apply rounding and overflow/underflow testing. */
1125 /* ------------------------------------------------------------------ */
decNumberPlus(decNumber * res,const decNumber * rhs,decContext * set)1126 decNumber * decNumberPlus(decNumber *res, const decNumber *rhs,
1127 decContext *set) {
1128 decNumber dzero;
1129 uInt status=0; // accumulator
1130
1131 #if DECCHECK
1132 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1133 #endif
1134
1135 decNumberZero(&dzero); // make 0
1136 dzero.exponent=rhs->exponent; // [no coefficient expansion]
1137 decAddOp(res, &dzero, rhs, set, 0, &status);
1138 if (status!=0) decStatus(res, status, set);
1139 #if DECCHECK
1140 decCheckInexact(res, set);
1141 #endif
1142 return res;
1143 } // decNumberPlus
1144
1145 /* ------------------------------------------------------------------ */
1146 /* decNumberMultiply -- multiply two Numbers */
1147 /* */
1148 /* This computes C = A x B */
1149 /* */
1150 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
1151 /* lhs is A */
1152 /* rhs is B */
1153 /* set is the context */
1154 /* */
1155 /* C must have space for set->digits digits. */
1156 /* ------------------------------------------------------------------ */
decNumberMultiply(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1157 decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs,
1158 const decNumber *rhs, decContext *set) {
1159 uInt status=0; // accumulator
1160 decMultiplyOp(res, lhs, rhs, set, &status);
1161 if (status!=0) decStatus(res, status, set);
1162 #if DECCHECK
1163 decCheckInexact(res, set);
1164 #endif
1165 return res;
1166 } // decNumberMultiply
1167
1168 /* ------------------------------------------------------------------ */
1169 /* decNumberNormalize -- remove trailing zeros */
1170 /* */
1171 /* This computes C = 0 + A, and normalizes the result */
1172 /* */
1173 /* res is C, the result. C may be A */
1174 /* rhs is A */
1175 /* set is the context */
1176 /* */
1177 /* C must have space for set->digits digits. */
1178 /* ------------------------------------------------------------------ */
decNumberNormalize(decNumber * res,const decNumber * rhs,decContext * set)1179 decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs,
1180 decContext *set) {
1181 #if DECSUBSET
1182 decNumber *allocrhs=NULL; // non-NULL if rounded rhs allocated
1183 #endif
1184 uInt status=0; // as usual
1185 Int residue=0; // as usual
1186 Int dropped; // work
1187
1188 #if DECCHECK
1189 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1190 #endif
1191
1192 do { // protect allocated storage
1193 #if DECSUBSET
1194 if (!set->extended) {
1195 // reduce operand and set lostDigits status, as needed
1196 if (rhs->digits>set->digits) {
1197 allocrhs=decRoundOperand(rhs, set, &status);
1198 if (allocrhs==NULL) break;
1199 rhs=allocrhs;
1200 }
1201 }
1202 #endif
1203 // [following code does not require input rounding]
1204
1205 // specials copy through, except NaNs need care
1206 if (decNumberIsNaN(rhs)) {
1207 decNaNs(res, rhs, NULL, &status);
1208 break;
1209 }
1210
1211 // reduce result to the requested length and copy to result
1212 decCopyFit(res, rhs, set, &residue, &status); // copy & round
1213 decFinish(res, set, &residue, &status); // cleanup/set flags
1214 decTrim(res, 1, &dropped); // normalize in place
1215 } while(0); // end protected
1216
1217 #if DECSUBSET
1218 if (allocrhs !=NULL) free(allocrhs); // ..
1219 #endif
1220 if (status!=0) decStatus(res, status, set);// then report status
1221 return res;
1222 } // decNumberNormalize
1223
1224 /* ------------------------------------------------------------------ */
1225 /* decNumberPower -- raise a number to a power */
1226 /* */
1227 /* This computes C = A ** B */
1228 /* */
1229 /* res is C, the result. C may be A and/or B (e.g., X=X**X) */
1230 /* lhs is A */
1231 /* rhs is B */
1232 /* set is the context */
1233 /* */
1234 /* C must have space for set->digits digits. */
1235 /* */
1236 /* Mathematical function restrictions apply (see above); a NaN is */
1237 /* returned with Invalid_operation if a restriction is violated. */
1238 /* */
1239 /* However, if 1999999997<=B<=999999999 and B is an integer then the */
1240 /* restrictions on A and the context are relaxed to the usual bounds, */
1241 /* for compatibility with the earlier (integer power only) version */
1242 /* of this function. */
1243 /* */
1244 /* When B is an integer, the result may be exact, even if rounded. */
1245 /* */
1246 /* The final result is rounded according to the context; it will */
1247 /* almost always be correctly rounded, but may be up to 1 ulp in */
1248 /* error in rare cases. */
1249 /* ------------------------------------------------------------------ */
decNumberPower(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1250 decNumber * decNumberPower(decNumber *res, const decNumber *lhs,
1251 const decNumber *rhs, decContext *set) {
1252 #if DECSUBSET
1253 decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated
1254 decNumber *allocrhs=NULL; // .., rhs
1255 #endif
1256 decNumber *allocdac=NULL; // -> allocated acc buffer, iff used
1257 decNumber *allocinv=NULL; // -> allocated 1/x buffer, iff used
1258 Int reqdigits=set->digits; // requested DIGITS
1259 Int n; // rhs in binary
1260 Flag rhsint=0; // 1 if rhs is an integer
1261 Flag useint=0; // 1 if can use integer calculation
1262 Flag isoddint=0; // 1 if rhs is an integer and odd
1263 Int i; // work
1264 #if DECSUBSET
1265 Int dropped; // ..
1266 #endif
1267 uInt needbytes; // buffer size needed
1268 Flag seenbit; // seen a bit while powering
1269 Int residue=0; // rounding residue
1270 uInt status=0; // accumulators
1271 uByte bits=0; // result sign if errors
1272 decContext aset; // working context
1273 decNumber dnOne; // work value 1...
1274 // local accumulator buffer [a decNumber, with digits+elength+1 digits]
1275 decNumber dacbuff[D2N(DECBUFFER+9)];
1276 decNumber *dac=dacbuff; // -> result accumulator
1277 // same again for possible 1/lhs calculation
1278 decNumber invbuff[D2N(DECBUFFER+9)];
1279
1280 #if DECCHECK
1281 if (decCheckOperands(res, lhs, rhs, set)) return res;
1282 #endif
1283
1284 do { // protect allocated storage
1285 #if DECSUBSET
1286 if (!set->extended) { // reduce operands and set status, as needed
1287 if (lhs->digits>reqdigits) {
1288 alloclhs=decRoundOperand(lhs, set, &status);
1289 if (alloclhs==NULL) break;
1290 lhs=alloclhs;
1291 }
1292 if (rhs->digits>reqdigits) {
1293 allocrhs=decRoundOperand(rhs, set, &status);
1294 if (allocrhs==NULL) break;
1295 rhs=allocrhs;
1296 }
1297 }
1298 #endif
1299 // [following code does not require input rounding]
1300
1301 // handle NaNs and rhs Infinity (lhs infinity is harder)
1302 if (SPECIALARGS) {
1303 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { // NaNs
1304 decNaNs(res, lhs, rhs, &status);
1305 break;}
1306 if (decNumberIsInfinite(rhs)) { // rhs Infinity
1307 Flag rhsneg=rhs->bits&DECNEG; // save rhs sign
1308 if (decNumberIsNegative(lhs) // lhs<0
1309 && !decNumberIsZero(lhs)) // ..
1310 status|=DEC_Invalid_operation;
1311 else { // lhs >=0
1312 decNumberZero(&dnOne); // set up 1
1313 dnOne.lsu[0]=1;
1314 decNumberCompare(dac, lhs, &dnOne, set); // lhs ? 1
1315 decNumberZero(res); // prepare for 0/1/Infinity
1316 if (decNumberIsNegative(dac)) { // lhs<1
1317 if (rhsneg) res->bits|=DECINF; // +Infinity [else is +0]
1318 }
1319 else if (dac->lsu[0]==0) { // lhs=1
1320 // 1**Infinity is inexact, so return fully-padded 1.0000
1321 Int shift=set->digits-1;
1322 *res->lsu=1; // was 0, make int 1
1323 res->digits=decShiftToMost(res->lsu, 1, shift);
1324 res->exponent=-shift; // make 1.0000...
1325 status|=DEC_Inexact|DEC_Rounded; // deemed inexact
1326 }
1327 else { // lhs>1
1328 if (!rhsneg) res->bits|=DECINF; // +Infinity [else is +0]
1329 }
1330 } // lhs>=0
1331 break;}
1332 // [lhs infinity drops through]
1333 } // specials
1334
1335 // Original rhs may be an integer that fits and is in range
1336 n=decGetInt(rhs);
1337 if (n!=BADINT) { // it is an integer
1338 rhsint=1; // record the fact for 1**n
1339 isoddint=(Flag)n&1; // [works even if big]
1340 if (n!=BIGEVEN && n!=BIGODD) // can use integer path?
1341 useint=1; // looks good
1342 }
1343
1344 if (decNumberIsNegative(lhs) // -x ..
1345 && isoddint) bits=DECNEG; // .. to an odd power
1346
1347 // handle LHS infinity
1348 if (decNumberIsInfinite(lhs)) { // [NaNs already handled]
1349 uByte rbits=rhs->bits; // save
1350 decNumberZero(res); // prepare
1351 if (n==0) *res->lsu=1; // [-]Inf**0 => 1
1352 else {
1353 // -Inf**nonint -> error
1354 if (!rhsint && decNumberIsNegative(lhs)) {
1355 status|=DEC_Invalid_operation; // -Inf**nonint is error
1356 break;}
1357 if (!(rbits & DECNEG)) bits|=DECINF; // was not a **-n
1358 // [otherwise will be 0 or -0]
1359 res->bits=bits;
1360 }
1361 break;}
1362
1363 // similarly handle LHS zero
1364 if (decNumberIsZero(lhs)) {
1365 if (n==0) { // 0**0 => Error
1366 #if DECSUBSET
1367 if (!set->extended) { // [unless subset]
1368 decNumberZero(res);
1369 *res->lsu=1; // return 1
1370 break;}
1371 #endif
1372 status|=DEC_Invalid_operation;
1373 }
1374 else { // 0**x
1375 uByte rbits=rhs->bits; // save
1376 if (rbits & DECNEG) { // was a 0**(-n)
1377 #if DECSUBSET
1378 if (!set->extended) { // [bad if subset]
1379 status|=DEC_Invalid_operation;
1380 break;}
1381 #endif
1382 bits|=DECINF;
1383 }
1384 decNumberZero(res); // prepare
1385 // [otherwise will be 0 or -0]
1386 res->bits=bits;
1387 }
1388 break;}
1389
1390 // here both lhs and rhs are finite; rhs==0 is handled in the
1391 // integer path. Next handle the non-integer cases
1392 if (!useint) { // non-integral rhs
1393 // any -ve lhs is bad, as is either operand or context out of
1394 // bounds
1395 if (decNumberIsNegative(lhs)) {
1396 status|=DEC_Invalid_operation;
1397 break;}
1398 if (decCheckMath(lhs, set, &status)
1399 || decCheckMath(rhs, set, &status)) break; // variable status
1400
1401 decContextDefault(&aset, DEC_INIT_DECIMAL64); // clean context
1402 aset.emax=DEC_MAX_MATH; // usual bounds
1403 aset.emin=-DEC_MAX_MATH; // ..
1404 aset.clamp=0; // and no concrete format
1405
1406 // calculate the result using exp(ln(lhs)*rhs), which can
1407 // all be done into the accumulator, dac. The precision needed
1408 // is enough to contain the full information in the lhs (which
1409 // is the total digits, including exponent), or the requested
1410 // precision, if larger, + 4; 6 is used for the exponent
1411 // maximum length, and this is also used when it is shorter
1412 // than the requested digits as it greatly reduces the >0.5 ulp
1413 // cases at little cost (because Ln doubles digits each
1414 // iteration so a few extra digits rarely causes an extra
1415 // iteration)
1416 aset.digits=MAX(lhs->digits, set->digits)+6+4;
1417 } // non-integer rhs
1418
1419 else { // rhs is in-range integer
1420 if (n==0) { // x**0 = 1
1421 // (0**0 was handled above)
1422 decNumberZero(res); // result=1
1423 *res->lsu=1; // ..
1424 break;}
1425 // rhs is a non-zero integer
1426 if (n<0) n=-n; // use abs(n)
1427
1428 aset=*set; // clone the context
1429 aset.round=DEC_ROUND_HALF_EVEN; // internally use balanced
1430 // calculate the working DIGITS
1431 aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2;
1432 #if DECSUBSET
1433 if (!set->extended) aset.digits--; // use classic precision
1434 #endif
1435 // it's an error if this is more than can be handled
1436 if (aset.digits>DECNUMMAXP) {status|=DEC_Invalid_operation; break;}
1437 } // integer path
1438
1439 // aset.digits is the count of digits for the accumulator needed
1440 // if accumulator is too long for local storage, then allocate
1441 needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit);
1442 // [needbytes also used below if 1/lhs needed]
1443 if (needbytes>sizeof(dacbuff)) {
1444 allocdac=(decNumber *)malloc(needbytes);
1445 if (allocdac==NULL) { // hopeless -- abandon
1446 status|=DEC_Insufficient_storage;
1447 break;}
1448 dac=allocdac; // use the allocated space
1449 }
1450 // here, aset is set up and accumulator is ready for use
1451
1452 if (!useint) { // non-integral rhs
1453 // x ** y; special-case x=1 here as it will otherwise always
1454 // reduce to integer 1; decLnOp has a fastpath which detects
1455 // the case of x=1
1456 decLnOp(dac, lhs, &aset, &status); // dac=ln(lhs)
1457 // [no error possible, as lhs 0 already handled]
1458 if (ISZERO(dac)) { // x==1, 1.0, etc.
1459 // need to return fully-padded 1.0000 etc., but rhsint->1
1460 *dac->lsu=1; // was 0, make int 1
1461 if (!rhsint) { // add padding
1462 Int shift=set->digits-1;
1463 dac->digits=decShiftToMost(dac->lsu, 1, shift);
1464 dac->exponent=-shift; // make 1.0000...
1465 status|=DEC_Inexact|DEC_Rounded; // deemed inexact
1466 }
1467 }
1468 else {
1469 decMultiplyOp(dac, dac, rhs, &aset, &status); // dac=dac*rhs
1470 decExpOp(dac, dac, &aset, &status); // dac=exp(dac)
1471 }
1472 // and drop through for final rounding
1473 } // non-integer rhs
1474
1475 else { // carry on with integer
1476 decNumberZero(dac); // acc=1
1477 *dac->lsu=1; // ..
1478
1479 // if a negative power the constant 1 is needed, and if not subset
1480 // invert the lhs now rather than inverting the result later
1481 if (decNumberIsNegative(rhs)) { // was a **-n [hence digits>0]
1482 decNumber *inv=invbuff; // asssume use fixed buffer
1483 decNumberCopy(&dnOne, dac); // dnOne=1; [needed now or later]
1484 #if DECSUBSET
1485 if (set->extended) { // need to calculate 1/lhs
1486 #endif
1487 // divide lhs into 1, putting result in dac [dac=1/dac]
1488 decDivideOp(dac, &dnOne, lhs, &aset, DIVIDE, &status);
1489 // now locate or allocate space for the inverted lhs
1490 if (needbytes>sizeof(invbuff)) {
1491 allocinv=(decNumber *)malloc(needbytes);
1492 if (allocinv==NULL) { // hopeless -- abandon
1493 status|=DEC_Insufficient_storage;
1494 break;}
1495 inv=allocinv; // use the allocated space
1496 }
1497 // [inv now points to big-enough buffer or allocated storage]
1498 decNumberCopy(inv, dac); // copy the 1/lhs
1499 decNumberCopy(dac, &dnOne); // restore acc=1
1500 lhs=inv; // .. and go forward with new lhs
1501 #if DECSUBSET
1502 }
1503 #endif
1504 }
1505
1506 // Raise-to-the-power loop...
1507 seenbit=0; // set once a 1-bit is encountered
1508 for (i=1;;i++){ // for each bit [top bit ignored]
1509 // abandon if had overflow or terminal underflow
1510 if (status & (DEC_Overflow|DEC_Underflow)) { // interesting?
1511 if (status&DEC_Overflow || ISZERO(dac)) break;
1512 }
1513 // [the following two lines revealed an optimizer bug in a C++
1514 // compiler, with symptom: 5**3 -> 25, when n=n+n was used]
1515 n=n<<1; // move next bit to testable position
1516 if (n<0) { // top bit is set
1517 seenbit=1; // OK, significant bit seen
1518 decMultiplyOp(dac, dac, lhs, &aset, &status); // dac=dac*x
1519 }
1520 if (i==31) break; // that was the last bit
1521 if (!seenbit) continue; // no need to square 1
1522 decMultiplyOp(dac, dac, dac, &aset, &status); // dac=dac*dac [square]
1523 } /*i*/ // 32 bits
1524
1525 // complete internal overflow or underflow processing
1526 if (status & (DEC_Overflow|DEC_Subnormal)) {
1527 #if DECSUBSET
1528 // If subset, and power was negative, reverse the kind of -erflow
1529 // [1/x not yet done]
1530 if (!set->extended && decNumberIsNegative(rhs)) {
1531 if (status & DEC_Overflow)
1532 status^=DEC_Overflow | DEC_Underflow | DEC_Subnormal;
1533 else { // trickier -- Underflow may or may not be set
1534 status&=~(DEC_Underflow | DEC_Subnormal); // [one or both]
1535 status|=DEC_Overflow;
1536 }
1537 }
1538 #endif
1539 dac->bits=(dac->bits & ~DECNEG) | bits; // force correct sign
1540 // round subnormals [to set.digits rather than aset.digits]
1541 // or set overflow result similarly as required
1542 decFinalize(dac, set, &residue, &status);
1543 decNumberCopy(res, dac); // copy to result (is now OK length)
1544 break;
1545 }
1546
1547 #if DECSUBSET
1548 if (!set->extended && // subset math
1549 decNumberIsNegative(rhs)) { // was a **-n [hence digits>0]
1550 // so divide result into 1 [dac=1/dac]
1551 decDivideOp(dac, &dnOne, dac, &aset, DIVIDE, &status);
1552 }
1553 #endif
1554 } // rhs integer path
1555
1556 // reduce result to the requested length and copy to result
1557 decCopyFit(res, dac, set, &residue, &status);
1558 decFinish(res, set, &residue, &status); // final cleanup
1559 #if DECSUBSET
1560 if (!set->extended) decTrim(res, 0, &dropped); // trailing zeros
1561 #endif
1562 } while(0); // end protected
1563
1564 if (allocdac!=NULL) free(allocdac); // drop any storage used
1565 if (allocinv!=NULL) free(allocinv); // ..
1566 #if DECSUBSET
1567 if (alloclhs!=NULL) free(alloclhs); // ..
1568 if (allocrhs!=NULL) free(allocrhs); // ..
1569 #endif
1570 if (status!=0) decStatus(res, status, set);
1571 #if DECCHECK
1572 decCheckInexact(res, set);
1573 #endif
1574 return res;
1575 } // decNumberPower
1576
1577 /* ------------------------------------------------------------------ */
1578 /* decNumberQuantize -- force exponent to requested value */
1579 /* */
1580 /* This computes C = op(A, B), where op adjusts the coefficient */
1581 /* of C (by rounding or shifting) such that the exponent (-scale) */
1582 /* of C has exponent of B. The numerical value of C will equal A, */
1583 /* except for the effects of any rounding that occurred. */
1584 /* */
1585 /* res is C, the result. C may be A or B */
1586 /* lhs is A, the number to adjust */
1587 /* rhs is B, the number with exponent to match */
1588 /* set is the context */
1589 /* */
1590 /* C must have space for set->digits digits. */
1591 /* */
1592 /* Unless there is an error or the result is infinite, the exponent */
1593 /* after the operation is guaranteed to be equal to that of B. */
1594 /* ------------------------------------------------------------------ */
decNumberQuantize(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1595 decNumber * decNumberQuantize(decNumber *res, const decNumber *lhs,
1596 const decNumber *rhs, decContext *set) {
1597 uInt status=0; // accumulator
1598 decQuantizeOp(res, lhs, rhs, set, 1, &status);
1599 if (status!=0) decStatus(res, status, set);
1600 return res;
1601 } // decNumberQuantize
1602
1603 /* ------------------------------------------------------------------ */
1604 /* decNumberRescale -- force exponent to requested value */
1605 /* */
1606 /* This computes C = op(A, B), where op adjusts the coefficient */
1607 /* of C (by rounding or shifting) such that the exponent (-scale) */
1608 /* of C has the value B. The numerical value of C will equal A, */
1609 /* except for the effects of any rounding that occurred. */
1610 /* */
1611 /* res is C, the result. C may be A or B */
1612 /* lhs is A, the number to adjust */
1613 /* rhs is B, the requested exponent */
1614 /* set is the context */
1615 /* */
1616 /* C must have space for set->digits digits. */
1617 /* */
1618 /* Unless there is an error or the result is infinite, the exponent */
1619 /* after the operation is guaranteed to be equal to B. */
1620 /* ------------------------------------------------------------------ */
decNumberRescale(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1621 decNumber * decNumberRescale(decNumber *res, const decNumber *lhs,
1622 const decNumber *rhs, decContext *set) {
1623 uInt status=0; // accumulator
1624 decQuantizeOp(res, lhs, rhs, set, 0, &status);
1625 if (status!=0) decStatus(res, status, set);
1626 return res;
1627 } // decNumberRescale
1628
1629 /* ------------------------------------------------------------------ */
1630 /* decNumberRemainder -- divide and return remainder */
1631 /* */
1632 /* This computes C = A % B */
1633 /* */
1634 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1635 /* lhs is A */
1636 /* rhs is B */
1637 /* set is the context */
1638 /* */
1639 /* C must have space for set->digits digits. */
1640 /* ------------------------------------------------------------------ */
decNumberRemainder(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1641 decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs,
1642 const decNumber *rhs, decContext *set) {
1643 uInt status=0; // accumulator
1644 decDivideOp(res, lhs, rhs, set, REMAINDER, &status);
1645 if (status!=0) decStatus(res, status, set);
1646 #if DECCHECK
1647 decCheckInexact(res, set);
1648 #endif
1649 return res;
1650 } // decNumberRemainder
1651
1652 /* ------------------------------------------------------------------ */
1653 /* decNumberRemainderNear -- divide and return remainder from nearest */
1654 /* */
1655 /* This computes C = A % B, where % is the IEEE remainder operator */
1656 /* */
1657 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
1658 /* lhs is A */
1659 /* rhs is B */
1660 /* set is the context */
1661 /* */
1662 /* C must have space for set->digits digits. */
1663 /* ------------------------------------------------------------------ */
decNumberRemainderNear(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)1664 decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs,
1665 const decNumber *rhs, decContext *set) {
1666 uInt status=0; // accumulator
1667 decDivideOp(res, lhs, rhs, set, REMNEAR, &status);
1668 if (status!=0) decStatus(res, status, set);
1669 #if DECCHECK
1670 decCheckInexact(res, set);
1671 #endif
1672 return res;
1673 } // decNumberRemainderNear
1674
1675 /* ------------------------------------------------------------------ */
1676 /* decNumberSameQuantum -- test for equal exponents */
1677 /* */
1678 /* res is the result number, which will contain either 0 or 1 */
1679 /* lhs is a number to test */
1680 /* rhs is the second (usually a pattern) */
1681 /* */
1682 /* No errors are possible and no context is needed. */
1683 /* ------------------------------------------------------------------ */
decNumberSameQuantum(decNumber * res,const decNumber * lhs,const decNumber * rhs)1684 decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs,
1685 const decNumber *rhs) {
1686 Unit ret=0; // return value
1687
1688 #if DECCHECK
1689 if (decCheckOperands(res, lhs, rhs, DECUNUSED)) return res;
1690 #endif
1691
1692 if (SPECIALARGS) {
1693 if (decNumberIsNaN(lhs) && decNumberIsNaN(rhs)) ret=1;
1694 else if (decNumberIsInfinite(lhs) && decNumberIsInfinite(rhs)) ret=1;
1695 // [anything else with a special gives 0]
1696 }
1697 else if (lhs->exponent==rhs->exponent) ret=1;
1698
1699 decNumberZero(res); // OK to overwrite an operand now
1700 *res->lsu=ret;
1701 return res;
1702 } // decNumberSameQuantum
1703
1704 /* ------------------------------------------------------------------ */
1705 /* decNumberSquareRoot -- square root operator */
1706 /* */
1707 /* This computes C = squareroot(A) */
1708 /* */
1709 /* res is C, the result. C may be A */
1710 /* rhs is A */
1711 /* set is the context; note that rounding mode has no effect */
1712 /* */
1713 /* C must have space for set->digits digits. */
1714 /* ------------------------------------------------------------------ */
1715 /* This uses the following varying-precision algorithm in: */
1716 /* */
1717 /* Properly Rounded Variable Precision Square Root, T. E. Hull and */
1718 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
1719 /* pp229-237, ACM, September 1985. */
1720 /* */
1721 /* The square-root is calculated using Newton's method, after which */
1722 /* a check is made to ensure the result is correctly rounded. */
1723 /* */
1724 /* % [Reformatted original Numerical Turing source code follows.] */
1725 /* function sqrt(x : real) : real */
1726 /* % sqrt(x) returns the properly rounded approximation to the square */
1727 /* % root of x, in the precision of the calling environment, or it */
1728 /* % fails if x < 0. */
1729 /* % t e hull and a abrham, august, 1984 */
1730 /* if x <= 0 then */
1731 /* if x < 0 then */
1732 /* assert false */
1733 /* else */
1734 /* result 0 */
1735 /* end if */
1736 /* end if */
1737 /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
1738 /* var e := getexp(x) % exponent part of x */
1739 /* var approx : real */
1740 /* if e mod 2 = 0 then */
1741 /* approx := .259 + .819 * f % approx to root of f */
1742 /* else */
1743 /* f := f/l0 % adjustments */
1744 /* e := e + 1 % for odd */
1745 /* approx := .0819 + 2.59 * f % exponent */
1746 /* end if */
1747 /* */
1748 /* var p:= 3 */
1749 /* const maxp := currentprecision + 2 */
1750 /* loop */
1751 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
1752 /* precision p */
1753 /* approx := .5 * (approx + f/approx) */
1754 /* exit when p = maxp */
1755 /* end loop */
1756 /* */
1757 /* % approx is now within 1 ulp of the properly rounded square root */
1758 /* % of f; to ensure proper rounding, compare squares of (approx - */
1759 /* % l/2 ulp) and (approx + l/2 ulp) with f. */
1760 /* p := currentprecision */
1761 /* begin */
1762 /* precision p + 2 */
1763 /* const approxsubhalf := approx - setexp(.5, -p) */
1764 /* if mulru(approxsubhalf, approxsubhalf) > f then */
1765 /* approx := approx - setexp(.l, -p + 1) */
1766 /* else */
1767 /* const approxaddhalf := approx + setexp(.5, -p) */
1768 /* if mulrd(approxaddhalf, approxaddhalf) < f then */
1769 /* approx := approx + setexp(.l, -p + 1) */
1770 /* end if */
1771 /* end if */
1772 /* end */
1773 /* result setexp(approx, e div 2) % fix exponent */
1774 /* end sqrt */
1775 /* ------------------------------------------------------------------ */
decNumberSquareRoot(decNumber * res,const decNumber * rhs,decContext * set)1776 decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs,
1777 decContext *set) {
1778 decContext workset, approxset; // work contexts
1779 decNumber dzero; // used for constant zero
1780 Int maxp=set->digits+2; // largest working precision
1781 Int residue=0; // rounding residue
1782 uInt status=0, ignore=0; // status accumulators
1783 Int exp; // working exponent
1784 Int ideal; // ideal (preferred) exponent
1785 Int needbytes; // work
1786 Int dropped; // ..
1787
1788 #if DECSUBSET
1789 decNumber *allocrhs=NULL; // non-NULL if rounded rhs allocated
1790 #endif
1791 // buffer for f [needs +1 in case DECBUFFER 0]
1792 decNumber buff[D2N(DECBUFFER+1)];
1793 // buffer for a [needs +2 to match maxp]
1794 decNumber bufa[D2N(DECBUFFER+2)];
1795 // buffer for temporary, b [must be same size as a]
1796 decNumber bufb[D2N(DECBUFFER+2)];
1797 decNumber *allocbuff=NULL; // -> allocated buff, iff allocated
1798 decNumber *allocbufa=NULL; // -> allocated bufa, iff allocated
1799 decNumber *allocbufb=NULL; // -> allocated bufb, iff allocated
1800 decNumber *f=buff; // reduced fraction
1801 decNumber *a=bufa; // approximation to result
1802 decNumber *b=bufb; // intermediate result
1803 // buffer for temporary variable, up to 3 digits
1804 decNumber buft[D2N(3)];
1805 decNumber *t=buft; // up-to-3-digit constant or work
1806
1807 #if DECCHECK
1808 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1809 #endif
1810
1811 do { // protect allocated storage
1812 #if DECSUBSET
1813 if (!set->extended) {
1814 // reduce operand and set lostDigits status, as needed
1815 if (rhs->digits>set->digits) {
1816 allocrhs=decRoundOperand(rhs, set, &status);
1817 if (allocrhs==NULL) break;
1818 // [Note: 'f' allocation below could reuse this buffer if
1819 // used, but as this is rare they are kept separate for clarity.]
1820 rhs=allocrhs;
1821 }
1822 }
1823 #endif
1824 // [following code does not require input rounding]
1825
1826 // handle infinities and NaNs
1827 if (SPECIALARG) {
1828 if (decNumberIsInfinite(rhs)) { // an infinity
1829 if (decNumberIsNegative(rhs)) status|=DEC_Invalid_operation;
1830 else decNumberCopy(res, rhs); // +Infinity
1831 }
1832 else decNaNs(res, rhs, NULL, &status); // a NaN
1833 break;
1834 }
1835
1836 // calculate the ideal (preferred) exponent [floor(exp/2)]
1837 // [We would like to write: ideal=rhs->exponent>>1, but this
1838 // generates a compiler warning. Generated code is the same.]
1839 ideal=(rhs->exponent&~1)/2; // target
1840
1841 // handle zeros
1842 if (ISZERO(rhs)) {
1843 decNumberCopy(res, rhs); // could be 0 or -0
1844 res->exponent=ideal; // use the ideal [safe]
1845 break;
1846 }
1847
1848 // any other -x is an oops
1849 if (decNumberIsNegative(rhs)) {
1850 status|=DEC_Invalid_operation;
1851 break;
1852 }
1853
1854 // space is needed for three working variables
1855 // f -- the same precision as the RHS, reduced to 0.01->0.99...
1856 // a -- Hull's approximation -- precision, when assigned, is
1857 // currentprecision (we allow +2 for use as temporary)
1858 // b -- intermediate temporary result
1859 // if any is too long for local storage, then allocate
1860 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
1861 if (needbytes>(Int)sizeof(buff)) {
1862 allocbuff=(decNumber *)malloc(needbytes);
1863 if (allocbuff==NULL) { // hopeless -- abandon
1864 status|=DEC_Insufficient_storage;
1865 break;}
1866 f=allocbuff; // use the allocated space
1867 }
1868 // a and b both need to be able to hold a maxp-length number
1869 needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit);
1870 if (needbytes>(Int)sizeof(bufa)) { // [same applies to b]
1871 allocbufa=(decNumber *)malloc(needbytes);
1872 allocbufb=(decNumber *)malloc(needbytes);
1873 if (allocbufa==NULL || allocbufb==NULL) { // hopeless
1874 status|=DEC_Insufficient_storage;
1875 break;}
1876 a=allocbufa; // use the allocated space
1877 b=allocbufb; // ..
1878 }
1879
1880 // copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1
1881 decNumberCopy(f, rhs);
1882 exp=f->exponent+f->digits; // adjusted to Hull rules
1883 f->exponent=-(f->digits); // to range
1884
1885 // set up working contexts (the second is used for Numerical
1886 // Turing assignment)
1887 decContextDefault(&workset, DEC_INIT_DECIMAL64);
1888 decContextDefault(&approxset, DEC_INIT_DECIMAL64);
1889 approxset.digits=set->digits; // approx's length
1890
1891 // [Until further notice, no error is possible and status bits
1892 // (Rounded, etc.) should be ignored, not accumulated.]
1893
1894 // Calculate initial approximation, and allow for odd exponent
1895 workset.digits=set->digits; // p for initial calculation
1896 t->bits=0; t->digits=3;
1897 a->bits=0; a->digits=3;
1898 if ((exp & 1)==0) { // even exponent
1899 // Set t=0.259, a=0.819
1900 t->exponent=-3;
1901 a->exponent=-3;
1902 #if DECDPUN>=3
1903 t->lsu[0]=259;
1904 a->lsu[0]=819;
1905 #elif DECDPUN==2
1906 t->lsu[0]=59; t->lsu[1]=2;
1907 a->lsu[0]=19; a->lsu[1]=8;
1908 #else
1909 t->lsu[0]=9; t->lsu[1]=5; t->lsu[2]=2;
1910 a->lsu[0]=9; a->lsu[1]=1; a->lsu[2]=8;
1911 #endif
1912 }
1913 else { // odd exponent
1914 // Set t=0.0819, a=2.59
1915 f->exponent--; // f=f/10
1916 exp++; // e=e+1
1917 t->exponent=-4;
1918 a->exponent=-2;
1919 #if DECDPUN>=3
1920 t->lsu[0]=819;
1921 a->lsu[0]=259;
1922 #elif DECDPUN==2
1923 t->lsu[0]=19; t->lsu[1]=8;
1924 a->lsu[0]=59; a->lsu[1]=2;
1925 #else
1926 t->lsu[0]=9; t->lsu[1]=1; t->lsu[2]=8;
1927 a->lsu[0]=9; a->lsu[1]=5; a->lsu[2]=2;
1928 #endif
1929 }
1930 decMultiplyOp(a, a, f, &workset, &ignore); // a=a*f
1931 decAddOp(a, a, t, &workset, 0, &ignore); // ..+t
1932 // [a is now the initial approximation for sqrt(f), calculated with
1933 // currentprecision, which is also a's precision.]
1934
1935 // the main calculation loop
1936 decNumberZero(&dzero); // make 0
1937 decNumberZero(t); // set t = 0.5
1938 t->lsu[0]=5; // ..
1939 t->exponent=-1; // ..
1940 workset.digits=3; // initial p
1941 for (;;) {
1942 // set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp]
1943 workset.digits=workset.digits*2-2;
1944 if (workset.digits>maxp) workset.digits=maxp;
1945 // a = 0.5 * (a + f/a)
1946 // [calculated at p then rounded to currentprecision]
1947 decDivideOp(b, f, a, &workset, DIVIDE, &ignore); // b=f/a
1948 decAddOp(b, b, a, &workset, 0, &ignore); // b=b+a
1949 decMultiplyOp(a, b, t, &workset, &ignore); // a=b*0.5
1950 // assign to approx [round to length]
1951 decAddOp(a, &dzero, a, &approxset, 0, &ignore);
1952 if (workset.digits==maxp) break; // just did final
1953 } // loop
1954
1955 // a is now at currentprecision and within 1 ulp of the properly
1956 // rounded square root of f; to ensure proper rounding, compare
1957 // squares of (a - l/2 ulp) and (a + l/2 ulp) with f.
1958 // Here workset.digits=maxp and t=0.5
1959 workset.digits--; // maxp-1 is OK now
1960 t->exponent=-set->digits-1; // make 0.5 ulp
1961 decNumberCopy(b, a);
1962 decAddOp(b, b, t, &workset, DECNEG, &ignore); // b = a - 0.5 ulp
1963 workset.round=DEC_ROUND_UP;
1964 decMultiplyOp(b, b, b, &workset, &ignore); // b = mulru(b, b)
1965 decCompareOp(b, f, b, &workset, COMPARE, &ignore); // b ? f, reversed
1966 if (decNumberIsNegative(b)) { // f < b [i.e., b > f]
1967 // this is the more common adjustment, though both are rare
1968 t->exponent++; // make 1.0 ulp
1969 t->lsu[0]=1; // ..
1970 decAddOp(a, a, t, &workset, DECNEG, &ignore); // a = a - 1 ulp
1971 // assign to approx [round to length]
1972 decAddOp(a, &dzero, a, &approxset, 0, &ignore);
1973 }
1974 else {
1975 decNumberCopy(b, a);
1976 decAddOp(b, b, t, &workset, 0, &ignore); // b = a + 0.5 ulp
1977 workset.round=DEC_ROUND_DOWN;
1978 decMultiplyOp(b, b, b, &workset, &ignore); // b = mulrd(b, b)
1979 decCompareOp(b, b, f, &workset, COMPARE, &ignore); // b ? f
1980 if (decNumberIsNegative(b)) { // b < f
1981 t->exponent++; // make 1.0 ulp
1982 t->lsu[0]=1; // ..
1983 decAddOp(a, a, t, &workset, 0, &ignore); // a = a + 1 ulp
1984 // assign to approx [round to length]
1985 decAddOp(a, &dzero, a, &approxset, 0, &ignore);
1986 }
1987 }
1988 // [no errors are possible in the above, and rounding/inexact during
1989 // estimation are irrelevant, so status was not accumulated]
1990
1991 // Here, 0.1 <= a < 1 [Hull]
1992 a->exponent+=exp/2; // set correct exponent
1993
1994 // Process Subnormals
1995 decFinalize(a, set, &residue, &status);
1996
1997 // count droppable zeros [after any subnormal rounding] by
1998 // trimming a copy
1999 decNumberCopy(b, a);
2000 decTrim(b, 1, &dropped); // [drops trailing zeros]
2001
2002 // Finally set Inexact and Rounded. The answer can only be exact if
2003 // it is short enough so that squaring it could fit in set->digits,
2004 // so this is the only (relatively rare) time a careful check is
2005 // needed
2006 if (b->digits*2-1 > set->digits) { // cannot fit
2007 status|=DEC_Inexact|DEC_Rounded;
2008 }
2009 else { // could be exact/unrounded
2010 uInt mstatus=0; // local status
2011 decMultiplyOp(b, b, b, &workset, &mstatus); // try the multiply
2012 if (mstatus!=0) { // result won't fit
2013 status|=DEC_Inexact|DEC_Rounded;
2014 }
2015 else { // plausible
2016 decCompareOp(t, b, rhs, &workset, COMPARE, &mstatus); // b ? rhs
2017 if (!ISZERO(t)) {
2018 status|=DEC_Inexact|DEC_Rounded;
2019 }
2020 else { // is Exact
2021 Int todrop; // work
2022 // here, dropped is the count of trailing zeros in 'a'
2023 // use closest exponent to ideal...
2024 todrop=ideal-a->exponent; // most that can be dropped
2025 if (todrop<0) { // ideally would add 0s
2026 status|=DEC_Rounded;
2027 }
2028 else { // unrounded
2029 if (dropped<todrop) todrop=dropped; // clamp to those available
2030 if (todrop>0) { // OK, some to drop
2031 decShiftToLeast(a->lsu, D2U(a->digits), todrop);
2032 a->exponent+=todrop; // maintain numerical value
2033 a->digits-=todrop; // new length
2034 }
2035 }
2036 }
2037 }
2038 }
2039
2040 // make sure there is a full complement of digits for normal
2041 // inexact results
2042 if ((status & (DEC_Inexact|DEC_Subnormal))==DEC_Inexact) {
2043 Int shift=set->digits-a->digits;
2044 if (shift>0) {
2045 a->digits=decShiftToMost(a->lsu, a->digits, shift);
2046 a->exponent-=shift; // adjust the exponent.
2047 }
2048 }
2049 decNumberCopy(res, a); // a is now the result
2050 } while(0); // end protected
2051
2052 if (allocbuff!=NULL) free(allocbuff); // drop any storage used
2053 if (allocbufa!=NULL) free(allocbufa); // ..
2054 if (allocbufb!=NULL) free(allocbufb); // ..
2055 #if DECSUBSET
2056 if (allocrhs !=NULL) free(allocrhs); // ..
2057 #endif
2058 if (status!=0) decStatus(res, status, set);// then report status
2059 #if DECCHECK
2060 decCheckInexact(res, set);
2061 #endif
2062 return res;
2063 } // decNumberSquareRoot
2064
2065 /* ------------------------------------------------------------------ */
2066 /* decNumberSubtract -- subtract two Numbers */
2067 /* */
2068 /* This computes C = A - B */
2069 /* */
2070 /* res is C, the result. C may be A and/or B (e.g., X=X-X) */
2071 /* lhs is A */
2072 /* rhs is B */
2073 /* set is the context */
2074 /* */
2075 /* C must have space for set->digits digits. */
2076 /* ------------------------------------------------------------------ */
decNumberSubtract(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set)2077 decNumber * decNumberSubtract(decNumber *res, const decNumber *lhs,
2078 const decNumber *rhs, decContext *set) {
2079 uInt status=0; // accumulator
2080
2081 decAddOp(res, lhs, rhs, set, DECNEG, &status);
2082 if (status!=0) decStatus(res, status, set);
2083 #if DECCHECK
2084 decCheckInexact(res, set);
2085 #endif
2086 return res;
2087 } // decNumberSubtract
2088
2089 /* ------------------------------------------------------------------ */
2090 /* decNumberToIntegralValue -- round-to-integral-value */
2091 /* */
2092 /* res is the result */
2093 /* rhs is input number */
2094 /* set is the context */
2095 /* */
2096 /* res must have space for any value of rhs. */
2097 /* */
2098 /* This implements the IEEE special operator and therefore treats */
2099 /* special values as valid, and also never sets Inexact. For finite */
2100 /* numbers it returns rescale(rhs, 0) if rhs->exponent is <0. */
2101 /* Otherwise the result is rhs (so no error is possible). */
2102 /* */
2103 /* The context is used for rounding mode and status after sNaN, but */
2104 /* the digits setting is ignored. */
2105 /* ------------------------------------------------------------------ */
decNumberToIntegralValue(decNumber * res,const decNumber * rhs,decContext * set)2106 decNumber * decNumberToIntegralValue(decNumber *res, const decNumber *rhs,
2107 decContext *set) {
2108 decNumber dn;
2109 decContext workset; // working context
2110
2111 #if DECCHECK
2112 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
2113 #endif
2114
2115 // handle infinities and NaNs
2116 if (SPECIALARG) {
2117 uInt status=0;
2118 if (decNumberIsInfinite(rhs)) decNumberCopy(res, rhs); // an Infinity
2119 else decNaNs(res, rhs, NULL, &status); // a NaN
2120 if (status!=0) decStatus(res, status, set);
2121 return res;
2122 }
2123
2124 // have a finite number; no error possible (res must be big enough)
2125 if (rhs->exponent>=0) return decNumberCopy(res, rhs);
2126 // that was easy, but if negative exponent there is work to do...
2127 workset=*set; // clone rounding, etc.
2128 workset.digits=rhs->digits; // no length rounding
2129 workset.traps=0; // no traps
2130 decNumberZero(&dn); // make a number with exponent 0
2131 return decNumberQuantize(res, rhs, &dn, &workset);
2132 } // decNumberToIntegralValue
2133
2134 /* ================================================================== */
2135 /* Utility routines */
2136 /* ================================================================== */
2137
2138 /* ------------------------------------------------------------------ */
2139 /* decNumberCopy -- copy a number */
2140 /* */
2141 /* dest is the target decNumber */
2142 /* src is the source decNumber */
2143 /* returns dest */
2144 /* */
2145 /* (dest==src is allowed and is a no-op) */
2146 /* All fields are updated as required. This is a utility operation, */
2147 /* so special values are unchanged and no error is possible. */
2148 /* ------------------------------------------------------------------ */
decNumberCopy(decNumber * dest,const decNumber * src)2149 decNumber * decNumberCopy(decNumber *dest, const decNumber *src) {
2150
2151 #if DECCHECK
2152 if (src==NULL) return decNumberZero(dest);
2153 #endif
2154
2155 if (dest==src) return dest; // no copy required
2156
2157 // Use explicit assignments here as structure assignment could copy
2158 // more than just the lsu (for small DECDPUN). This would not affect
2159 // the value of the results, but could disturb test harness spill
2160 // checking.
2161 dest->bits=src->bits;
2162 dest->exponent=src->exponent;
2163 dest->digits=src->digits;
2164 dest->lsu[0]=src->lsu[0];
2165 if (src->digits>DECDPUN) { // more Units to come
2166 const Unit *smsup, *s; // work
2167 Unit *d; // ..
2168 // memcpy for the remaining Units would be safe as they cannot
2169 // overlap. However, this explicit loop is faster in short cases.
2170 d=dest->lsu+1; // -> first destination
2171 smsup=src->lsu+D2U(src->digits); // -> source msu+1
2172 for (s=src->lsu+1; s<smsup; s++, d++) *d=*s;
2173 }
2174 return dest;
2175 } // decNumberCopy
2176
2177 /* ------------------------------------------------------------------ */
2178 /* decNumberTrim -- remove insignificant zeros */
2179 /* */
2180 /* dn is the number to trim */
2181 /* returns dn */
2182 /* */
2183 /* All fields are updated as required. This is a utility operation, */
2184 /* so special values are unchanged and no error is possible. */
2185 /* ------------------------------------------------------------------ */
decNumberTrim(decNumber * dn)2186 decNumber * decNumberTrim(decNumber *dn) {
2187 Int dropped; // work
2188 return decTrim(dn, 0, &dropped);
2189 } // decNumberTrim
2190
2191 /* ------------------------------------------------------------------ */
2192 /* decNumberVersion -- return the name and version of this module */
2193 /* */
2194 /* No error is possible. */
2195 /* ------------------------------------------------------------------ */
decNumberVersion(void)2196 const char * decNumberVersion(void) {
2197 return DECVERSION;
2198 } // decNumberVersion
2199
2200 /* ------------------------------------------------------------------ */
2201 /* decNumberZero -- set a number to 0 */
2202 /* */
2203 /* dn is the number to set, with space for one digit */
2204 /* returns dn */
2205 /* */
2206 /* No error is possible. */
2207 /* ------------------------------------------------------------------ */
2208 // Memset is not used as it is much slower in some environments.
decNumberZero(decNumber * dn)2209 decNumber * decNumberZero(decNumber *dn) {
2210
2211 #if DECCHECK
2212 if (decCheckOperands(dn, DECUNUSED, DECUNUSED, DECUNUSED)) return dn;
2213 #endif
2214
2215 dn->bits=0;
2216 dn->exponent=0;
2217 dn->digits=1;
2218 dn->lsu[0]=0;
2219 return dn;
2220 } // decNumberZero
2221
2222 /* ================================================================== */
2223 /* Local routines */
2224 /* ================================================================== */
2225
2226 /* ------------------------------------------------------------------ */
2227 /* decToString -- lay out a number into a string */
2228 /* */
2229 /* dn is the number to lay out */
2230 /* string is where to lay out the number */
2231 /* eng is 1 if Engineering, 0 if Scientific */
2232 /* */
2233 /* string must be at least dn->digits+14 characters long */
2234 /* No error is possible. */
2235 /* */
2236 /* Note that this routine can generate a -0 or 0.000. These are */
2237 /* never generated in subset to-number or arithmetic, but can occur */
2238 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234). */
2239 /* ------------------------------------------------------------------ */
2240 // If DECCHECK is enabled the string "?" is returned if a number is
2241 // invalid.
decToString(const decNumber * dn,char * string,Flag eng)2242 static void decToString(const decNumber *dn, char *string, Flag eng) {
2243 Int exp=dn->exponent; // local copy
2244 Int e; // E-part value
2245 Int pre; // digits before the '.'
2246 Int cut; // for counting digits in a Unit
2247 char *c=string; // work [output pointer]
2248 const Unit *up=dn->lsu+D2U(dn->digits)-1; // -> msu [input pointer]
2249 uInt u, pow; // work
2250
2251 #if DECCHECK
2252 if (decCheckOperands(DECUNUSED, dn, DECUNUSED, DECUNUSED)) {
2253 strcpy(string, "?");
2254 return;}
2255 #endif
2256
2257 if (decNumberIsNegative(dn)) { // Negatives get a minus
2258 *c='-';
2259 c++;
2260 }
2261 if (dn->bits&DECSPECIAL) { // Is a special value
2262 if (decNumberIsInfinite(dn)) {
2263 strcpy(c, "Infinity");
2264 return;}
2265 // a NaN
2266 if (dn->bits&DECSNAN) { // signalling NaN
2267 *c='s';
2268 c++;
2269 }
2270 strcpy(c, "NaN");
2271 c+=3; // step past
2272 // if not a clean non-zero coefficient, that's all there is in a
2273 // NaN string
2274 if (exp!=0 || (*dn->lsu==0 && dn->digits==1)) return;
2275 // [drop through to add integer]
2276 }
2277
2278 // calculate how many digits in msu, and hence first cut
2279 cut=MSUDIGITS(dn->digits); // [faster than remainder]
2280 cut--; // power of ten for digit
2281
2282 if (exp==0) { // simple integer [common fastpath]
2283 for (;up>=dn->lsu; up--) { // each Unit from msu
2284 u=*up; // contains DECDPUN digits to lay out
2285 for (; cut>=0; c++, cut--) TODIGIT(u, cut, c, pow);
2286 cut=DECDPUN-1; // next Unit has all digits
2287 }
2288 *c='\0'; // terminate the string
2289 return;}
2290
2291 /* non-0 exponent -- assume plain form */
2292 pre=dn->digits+exp; // digits before '.'
2293 e=0; // no E
2294 if ((exp>0) || (pre<-5)) { // need exponential form
2295 e=exp+dn->digits-1; // calculate E value
2296 pre=1; // assume one digit before '.'
2297 if (eng && (e!=0)) { // engineering: may need to adjust
2298 Int adj; // adjustment
2299 // The C remainder operator is undefined for negative numbers, so
2300 // a positive remainder calculation must be used here
2301 if (e<0) {
2302 adj=(-e)%3;
2303 if (adj!=0) adj=3-adj;
2304 }
2305 else { // e>0
2306 adj=e%3;
2307 }
2308 e=e-adj;
2309 // if dealing with zero still produce an exponent which is a
2310 // multiple of three, as expected, but there will only be the
2311 // one zero before the E, still. Otherwise note the padding.
2312 if (!ISZERO(dn)) pre+=adj;
2313 else { // is zero
2314 if (adj!=0) { // 0.00Esnn needed
2315 e=e+3;
2316 pre=-(2-adj);
2317 }
2318 } // zero
2319 } // eng
2320 } // need exponent
2321
2322 /* lay out the digits of the coefficient, adding 0s and . as needed */
2323 u=*up;
2324 if (pre>0) { // xxx.xxx or xx00 (engineering) form
2325 Int n=pre;
2326 for (; pre>0; pre--, c++, cut--) {
2327 if (cut<0) { // need new Unit
2328 if (up==dn->lsu) break; // out of input digits (pre>digits)
2329 up--;
2330 cut=DECDPUN-1;
2331 u=*up;
2332 }
2333 TODIGIT(u, cut, c, pow);
2334 }
2335 if (n<dn->digits) { // more to come, after '.'
2336 *c='.'; c++;
2337 for (;; c++, cut--) {
2338 if (cut<0) { // need new Unit
2339 if (up==dn->lsu) break; // out of input digits
2340 up--;
2341 cut=DECDPUN-1;
2342 u=*up;
2343 }
2344 TODIGIT(u, cut, c, pow);
2345 }
2346 }
2347 else for (; pre>0; pre--, c++) *c='0'; // 0 padding (for engineering) needed
2348 }
2349 else { // 0.xxx or 0.000xxx form
2350 *c='0'; c++;
2351 *c='.'; c++;
2352 for (; pre<0; pre++, c++) *c='0'; // add any 0's after '.'
2353 for (; ; c++, cut--) {
2354 if (cut<0) { // need new Unit
2355 if (up==dn->lsu) break; // out of input digits
2356 up--;
2357 cut=DECDPUN-1;
2358 u=*up;
2359 }
2360 TODIGIT(u, cut, c, pow);
2361 }
2362 }
2363
2364 /* Finally add the E-part, if needed. It will never be 0, has a
2365 base maximum and minimum of +999999999 through -999999999, but
2366 could range down to -1999999998 for anormal numbers */
2367 if (e!=0) {
2368 Flag had=0; // 1=had non-zero
2369 *c='E'; c++;
2370 *c='+'; c++; // assume positive
2371 u=e; // ..
2372 if (e<0) {
2373 *(c-1)='-'; // oops, need -
2374 u=-e; // uInt, please
2375 }
2376 // layout the exponent [_itoa or equivalent is not ANSI C]
2377 for (cut=9; cut>=0; cut--) {
2378 TODIGIT(u, cut, c, pow);
2379 if (*c=='0' && !had) continue; // skip leading zeros
2380 had=1; // had non-0
2381 c++; // step for next
2382 } // cut
2383 }
2384 *c='\0'; // terminate the string (all paths)
2385 return;
2386 } // decToString
2387
2388 /* ------------------------------------------------------------------ */
2389 /* decAddOp -- add/subtract operation */
2390 /* */
2391 /* This computes C = A + B */
2392 /* */
2393 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
2394 /* lhs is A */
2395 /* rhs is B */
2396 /* set is the context */
2397 /* negate is DECNEG if rhs should be negated, or 0 otherwise */
2398 /* status accumulates status for the caller */
2399 /* */
2400 /* C must have space for set->digits digits. */
2401 /* ------------------------------------------------------------------ */
2402 /* If possible, the coefficient is calculated directly into C. */
2403 /* However, if: */
2404 /* -- a digits+1 calculation is needed because the numbers are */
2405 /* unaligned and span more than set->digits digits */
2406 /* -- a carry to digits+1 digits looks possible */
2407 /* -- C is the same as A or B, and the result would destructively */
2408 /* overlap the A or B coefficient */
2409 /* then the result must be calculated into a temporary buffer. In */
2410 /* this case a local (stack) buffer is used if possible, and only if */
2411 /* too long for that does malloc become the last resort. */
2412 /* */
2413 /* Misalignment is handled as follows: */
2414 /* Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp. */
2415 /* BPad: Apply the padding by a combination of shifting (whole */
2416 /* units) and multiplication (part units). */
2417 /* */
2418 /* Addition, especially x=x+1, is speed-critical. */
2419 /* ------------------------------------------------------------------ */
decAddOp(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set,uByte negate,uInt * status)2420 static decNumber * decAddOp(decNumber *res, const decNumber *lhs,
2421 const decNumber *rhs, decContext *set,
2422 uByte negate, uInt *status) {
2423 #if DECSUBSET
2424 decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated
2425 decNumber *allocrhs=NULL; // .., rhs
2426 #endif
2427 Int rhsshift; // working shift (in Units)
2428 Int maxdigits; // longest logical length
2429 Int mult; // multiplier
2430 Int residue; // rounding accumulator
2431 uByte bits; // result bits
2432 Flag diffsign; // non-0 if arguments have different sign
2433 Unit *acc; // accumulator for result
2434 Unit accbuff[SD2U(DECBUFFER+20)]; // local buffer [+20 reduces many
2435 // allocations when called from
2436 // other operations]
2437 Unit *allocacc=NULL; // -> allocated acc buffer, iff allocated
2438 Int reqdigits=set->digits; // local copy; requested DIGITS
2439 Int padding; // work
2440
2441 #if DECCHECK
2442 if (decCheckOperands(res, lhs, rhs, set)) return res;
2443 #endif
2444
2445 do { // protect allocated storage
2446 #if DECSUBSET
2447 if (!set->extended) {
2448 // reduce operands and set lostDigits status, as needed
2449 if (lhs->digits>reqdigits) {
2450 alloclhs=decRoundOperand(lhs, set, status);
2451 if (alloclhs==NULL) break;
2452 lhs=alloclhs;
2453 }
2454 if (rhs->digits>reqdigits) {
2455 allocrhs=decRoundOperand(rhs, set, status);
2456 if (allocrhs==NULL) break;
2457 rhs=allocrhs;
2458 }
2459 }
2460 #endif
2461 // [following code does not require input rounding]
2462
2463 // note whether signs differ [used all paths]
2464 diffsign=(Flag)((lhs->bits^rhs->bits^negate)&DECNEG);
2465
2466 // handle infinities and NaNs
2467 if (SPECIALARGS) { // a special bit set
2468 if (SPECIALARGS & (DECSNAN | DECNAN)) // a NaN
2469 decNaNs(res, lhs, rhs, status);
2470 else { // one or two infinities
2471 if (decNumberIsInfinite(lhs)) { // LHS is infinity
2472 // two infinities with different signs is invalid
2473 if (decNumberIsInfinite(rhs) && diffsign) {
2474 *status|=DEC_Invalid_operation;
2475 break;
2476 }
2477 bits=lhs->bits & DECNEG; // get sign from LHS
2478 }
2479 else bits=(rhs->bits^negate) & DECNEG;// RHS must be Infinity
2480 bits|=DECINF;
2481 decNumberZero(res);
2482 res->bits=bits; // set +/- infinity
2483 } // an infinity
2484 break;
2485 }
2486
2487 // Quick exit for add 0s; return the non-0, modified as need be
2488 if (ISZERO(lhs)) {
2489 Int adjust; // work
2490 Int lexp=lhs->exponent; // save in case LHS==RES
2491 bits=lhs->bits; // ..
2492 residue=0; // clear accumulator
2493 decCopyFit(res, rhs, set, &residue, status); // copy (as needed)
2494 res->bits^=negate; // flip if rhs was negated
2495 #if DECSUBSET
2496 if (set->extended) { // exponents on zeros count
2497 #endif
2498 // exponent will be the lower of the two
2499 adjust=lexp-res->exponent; // adjustment needed [if -ve]
2500 if (ISZERO(res)) { // both 0: special IEEE 854 rules
2501 if (adjust<0) res->exponent=lexp; // set exponent
2502 // 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0
2503 if (diffsign) {
2504 if (set->round!=DEC_ROUND_FLOOR) res->bits=0;
2505 else res->bits=DECNEG; // preserve 0 sign
2506 }
2507 }
2508 else { // non-0 res
2509 if (adjust<0) { // 0-padding needed
2510 if ((res->digits-adjust)>set->digits) {
2511 adjust=res->digits-set->digits; // to fit exactly
2512 *status|=DEC_Rounded; // [but exact]
2513 }
2514 res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
2515 res->exponent+=adjust; // set the exponent.
2516 }
2517 } // non-0 res
2518 #if DECSUBSET
2519 } // extended
2520 #endif
2521 decFinish(res, set, &residue, status); // clean and finalize
2522 break;}
2523
2524 if (ISZERO(rhs)) { // [lhs is non-zero]
2525 Int adjust; // work
2526 Int rexp=rhs->exponent; // save in case RHS==RES
2527 bits=rhs->bits; // be clean
2528 residue=0; // clear accumulator
2529 decCopyFit(res, lhs, set, &residue, status); // copy (as needed)
2530 #if DECSUBSET
2531 if (set->extended) { // exponents on zeros count
2532 #endif
2533 // exponent will be the lower of the two
2534 // [0-0 case handled above]
2535 adjust=rexp-res->exponent; // adjustment needed [if -ve]
2536 if (adjust<0) { // 0-padding needed
2537 if ((res->digits-adjust)>set->digits) {
2538 adjust=res->digits-set->digits; // to fit exactly
2539 *status|=DEC_Rounded; // [but exact]
2540 }
2541 res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
2542 res->exponent+=adjust; // set the exponent.
2543 }
2544 #if DECSUBSET
2545 } // extended
2546 #endif
2547 decFinish(res, set, &residue, status); // clean and finalize
2548 break;}
2549 // [NB: both fastpath and mainpath code below assume these cases
2550 // (notably 0-0) have already been handled]
2551
2552 // calculate the padding needed to align the operands
2553 padding=rhs->exponent-lhs->exponent;
2554
2555 // Fastpath cases where the numbers are aligned and normal, the RHS
2556 // is all in one unit, no operand rounding is needed, and no carry,
2557 // lengthening, or borrow is needed
2558 if (padding==0
2559 && rhs->digits<=DECDPUN
2560 && rhs->exponent>=set->emin // [some normals drop through]
2561 && rhs->digits<=reqdigits
2562 && lhs->digits<=reqdigits) {
2563 Int partial=*lhs->lsu;
2564 if (!diffsign) { // adding
2565 partial+=*rhs->lsu;
2566 if ((partial<=DECDPUNMAX) // result fits in unit
2567 && (lhs->digits>=DECDPUN || // .. and no digits-count change
2568 partial<(Int)powers[lhs->digits])) { // ..
2569 if (res!=lhs) decNumberCopy(res, lhs); // not in place
2570 *res->lsu=(Unit)partial; // [copy could have overwritten RHS]
2571 break;
2572 }
2573 // else drop out for careful add
2574 }
2575 else { // signs differ
2576 partial-=*rhs->lsu;
2577 if (partial>0) { // no borrow needed, and non-0 result
2578 if (res!=lhs) decNumberCopy(res, lhs); // not in place
2579 *res->lsu=(Unit)partial;
2580 // this could have reduced digits [but result>0]
2581 res->digits=decGetDigits(res->lsu, D2U(res->digits));
2582 break;
2583 }
2584 // else drop out for careful subtract
2585 }
2586 }
2587
2588 // Now align (pad) the lhs or rhs so they can be added or
2589 // subtracted, as necessary. If one number is much larger than
2590 // the other (that is, if in plain form there is a least one
2591 // digit between the lowest digit of one and the highest of the
2592 // other) padding with up to DIGITS-1 trailing zeros may be
2593 // needed; then apply rounding (as exotic rounding modes may be
2594 // affected by the residue).
2595 rhsshift=0; // rhs shift to left (padding) in Units
2596 bits=lhs->bits; // assume sign is that of LHS
2597 mult=1; // likely multiplier
2598
2599 // if padding==0 the operands are aligned; no padding needed
2600 if (padding!=0) {
2601 // some padding needed; always pad the RHS, as any required
2602 // padding can then be effected by a simple combination of
2603 // shifts and a multiply
2604 Flag swapped=0;
2605 if (padding<0) { // LHS needs the padding
2606 const decNumber *t;
2607 padding=-padding; // will be +ve
2608 bits=(uByte)(rhs->bits^negate); // assumed sign is now that of RHS
2609 t=lhs; lhs=rhs; rhs=t;
2610 swapped=1;
2611 }
2612
2613 // If, after pad, rhs would be longer than lhs by digits+1 or
2614 // more then lhs cannot affect the answer, except as a residue,
2615 // so only need to pad up to a length of DIGITS+1.
2616 if (rhs->digits+padding > lhs->digits+reqdigits+1) {
2617 // The RHS is sufficient
2618 // for residue use the relative sign indication...
2619 Int shift=reqdigits-rhs->digits; // left shift needed
2620 residue=1; // residue for rounding
2621 if (diffsign) residue=-residue; // signs differ
2622 // copy, shortening if necessary
2623 decCopyFit(res, rhs, set, &residue, status);
2624 // if it was already shorter, then need to pad with zeros
2625 if (shift>0) {
2626 res->digits=decShiftToMost(res->lsu, res->digits, shift);
2627 res->exponent-=shift; // adjust the exponent.
2628 }
2629 // flip the result sign if unswapped and rhs was negated
2630 if (!swapped) res->bits^=negate;
2631 decFinish(res, set, &residue, status); // done
2632 break;}
2633
2634 // LHS digits may affect result
2635 rhsshift=D2U(padding+1)-1; // this much by Unit shift ..
2636 mult=powers[padding-(rhsshift*DECDPUN)]; // .. this by multiplication
2637 } // padding needed
2638
2639 if (diffsign) mult=-mult; // signs differ
2640
2641 // determine the longer operand
2642 maxdigits=rhs->digits+padding; // virtual length of RHS
2643 if (lhs->digits>maxdigits) maxdigits=lhs->digits;
2644
2645 // Decide on the result buffer to use; if possible place directly
2646 // into result.
2647 acc=res->lsu; // assume add direct to result
2648 // If destructive overlap, or the number is too long, or a carry or
2649 // borrow to DIGITS+1 might be possible, a buffer must be used.
2650 // [Might be worth more sophisticated tests when maxdigits==reqdigits]
2651 if ((maxdigits>=reqdigits) // is, or could be, too large
2652 || (res==rhs && rhsshift>0)) { // destructive overlap
2653 // buffer needed, choose it; units for maxdigits digits will be
2654 // needed, +1 Unit for carry or borrow
2655 Int need=D2U(maxdigits)+1;
2656 acc=accbuff; // assume use local buffer
2657 if (need*sizeof(Unit)>sizeof(accbuff)) {
2658 allocacc=(Unit *)malloc(need*sizeof(Unit));
2659 if (allocacc==NULL) { // hopeless -- abandon
2660 *status|=DEC_Insufficient_storage;
2661 break;}
2662 acc=allocacc;
2663 }
2664 }
2665
2666 res->bits=(uByte)(bits&DECNEG); // it's now safe to overwrite..
2667 res->exponent=lhs->exponent; // .. operands (even if aliased)
2668
2669 #if DECTRACE
2670 decDumpAr('A', lhs->lsu, D2U(lhs->digits));
2671 decDumpAr('B', rhs->lsu, D2U(rhs->digits));
2672 printf(" :h: %d %d\n", rhsshift, mult);
2673 #endif
2674
2675 // add [A+B*m] or subtract [A+B*(-m)]
2676 res->digits=decUnitAddSub(lhs->lsu, D2U(lhs->digits),
2677 rhs->lsu, D2U(rhs->digits),
2678 rhsshift, acc, mult)
2679 *DECDPUN; // [units -> digits]
2680 if (res->digits<0) { // borrowed...
2681 res->digits=-res->digits;
2682 res->bits^=DECNEG; // flip the sign
2683 }
2684 #if DECTRACE
2685 decDumpAr('+', acc, D2U(res->digits));
2686 #endif
2687
2688 // If a buffer was used the result must be copied back, possibly
2689 // shortening. (If no buffer was used then the result must have
2690 // fit, so can't need rounding and residue must be 0.)
2691 residue=0; // clear accumulator
2692 if (acc!=res->lsu) {
2693 #if DECSUBSET
2694 if (set->extended) { // round from first significant digit
2695 #endif
2696 // remove leading zeros that were added due to rounding up to
2697 // integral Units -- before the test for rounding.
2698 if (res->digits>reqdigits)
2699 res->digits=decGetDigits(acc, D2U(res->digits));
2700 decSetCoeff(res, set, acc, res->digits, &residue, status);
2701 #if DECSUBSET
2702 }
2703 else { // subset arithmetic rounds from original significant digit
2704 // May have an underestimate. This only occurs when both
2705 // numbers fit in DECDPUN digits and are padding with a
2706 // negative multiple (-10, -100...) and the top digit(s) become
2707 // 0. (This only matters when using X3.274 rules where the
2708 // leading zero could be included in the rounding.)
2709 if (res->digits<maxdigits) {
2710 *(acc+D2U(res->digits))=0; // ensure leading 0 is there
2711 res->digits=maxdigits;
2712 }
2713 else {
2714 // remove leading zeros that added due to rounding up to
2715 // integral Units (but only those in excess of the original
2716 // maxdigits length, unless extended) before test for rounding.
2717 if (res->digits>reqdigits) {
2718 res->digits=decGetDigits(acc, D2U(res->digits));
2719 if (res->digits<maxdigits) res->digits=maxdigits;
2720 }
2721 }
2722 decSetCoeff(res, set, acc, res->digits, &residue, status);
2723 // Now apply rounding if needed before removing leading zeros.
2724 // This is safe because subnormals are not a possibility
2725 if (residue!=0) {
2726 decApplyRound(res, set, residue, status);
2727 residue=0; // did what needed to be done
2728 }
2729 } // subset
2730 #endif
2731 } // used buffer
2732
2733 // strip leading zeros [these were left on in case of subset subtract]
2734 res->digits=decGetDigits(res->lsu, D2U(res->digits));
2735
2736 // apply checks and rounding
2737 decFinish(res, set, &residue, status);
2738
2739 // "When the sum of two operands with opposite signs is exactly
2740 // zero, the sign of that sum shall be '+' in all rounding modes
2741 // except round toward -Infinity, in which mode that sign shall be
2742 // '-'." [Subset zeros also never have '-', set by decFinish.]
2743 if (ISZERO(res)
2744 && diffsign
2745 #if DECSUBSET
2746 && set->extended
2747 #endif
2748 && (*status&DEC_Inexact)==0) {
2749 if (set->round==DEC_ROUND_FLOOR) res->bits|=DECNEG; // sign -
2750 else res->bits&=~DECNEG; // sign +
2751 }
2752 } while(0); // end protected
2753
2754 if (allocacc!=NULL) free(allocacc); // drop any storage used
2755 #if DECSUBSET
2756 if (allocrhs!=NULL) free(allocrhs); // ..
2757 if (alloclhs!=NULL) free(alloclhs); // ..
2758 #endif
2759 return res;
2760 } // decAddOp
2761
2762 /* ------------------------------------------------------------------ */
2763 /* decDivideOp -- division operation */
2764 /* */
2765 /* This routine performs the calculations for all four division */
2766 /* operators (divide, divideInteger, remainder, remainderNear). */
2767 /* */
2768 /* C=A op B */
2769 /* */
2770 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
2771 /* lhs is A */
2772 /* rhs is B */
2773 /* set is the context */
2774 /* op is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively. */
2775 /* status is the usual accumulator */
2776 /* */
2777 /* C must have space for set->digits digits. */
2778 /* */
2779 /* ------------------------------------------------------------------ */
2780 /* The underlying algorithm of this routine is the same as in the */
2781 /* 1981 S/370 implementation, that is, non-restoring long division */
2782 /* with bi-unit (rather than bi-digit) estimation for each unit */
2783 /* multiplier. In this pseudocode overview, complications for the */
2784 /* Remainder operators and division residues for exact rounding are */
2785 /* omitted for clarity. */
2786 /* */
2787 /* Prepare operands and handle special values */
2788 /* Test for x/0 and then 0/x */
2789 /* Exp =Exp1 - Exp2 */
2790 /* Exp =Exp +len(var1) -len(var2) */
2791 /* Sign=Sign1 * Sign2 */
2792 /* Pad accumulator (Var1) to double-length with 0's (pad1) */
2793 /* Pad Var2 to same length as Var1 */
2794 /* msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round */
2795 /* have=0 */
2796 /* Do until (have=digits+1 OR residue=0) */
2797 /* if exp<0 then if integer divide/residue then leave */
2798 /* this_unit=0 */
2799 /* Do forever */
2800 /* compare numbers */
2801 /* if <0 then leave inner_loop */
2802 /* if =0 then (* quick exit without subtract *) do */
2803 /* this_unit=this_unit+1; output this_unit */
2804 /* leave outer_loop; end */
2805 /* Compare lengths of numbers (mantissae): */
2806 /* If same then tops2=msu2pair -- {units 1&2 of var2} */
2807 /* else tops2=msu2plus -- {0, unit 1 of var2} */
2808 /* tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
2809 /* mult=tops1/tops2 -- Good and safe guess at divisor */
2810 /* if mult=0 then mult=1 */
2811 /* this_unit=this_unit+mult */
2812 /* subtract */
2813 /* end inner_loop */
2814 /* if have\=0 | this_unit\=0 then do */
2815 /* output this_unit */
2816 /* have=have+1; end */
2817 /* var2=var2/10 */
2818 /* exp=exp-1 */
2819 /* end outer_loop */
2820 /* exp=exp+1 -- set the proper exponent */
2821 /* if have=0 then generate answer=0 */
2822 /* Return (Result is defined by Var1) */
2823 /* */
2824 /* ------------------------------------------------------------------ */
2825 /* Two working buffers are needed during the division; one (digits+ */
2826 /* 1) to accumulate the result, and the other (up to 2*digits+1) for */
2827 /* long subtractions. These are acc and var1 respectively. */
2828 /* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
2829 /* ------------------------------------------------------------------ */
decDivideOp(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set,Flag op,uInt * status)2830 static decNumber * decDivideOp(decNumber *res,
2831 const decNumber *lhs, const decNumber *rhs,
2832 decContext *set, Flag op, uInt *status) {
2833 #if DECSUBSET
2834 decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated
2835 decNumber *allocrhs=NULL; // .., rhs
2836 #endif
2837 Unit accbuff[SD2U(DECBUFFER+DECDPUN)]; // local buffer
2838 Unit *acc=accbuff; // -> accumulator array for result
2839 Unit *allocacc=NULL; // -> allocated buffer, iff allocated
2840 Unit *accnext; // -> where next digit will go
2841 Int acclength; // length of acc needed [Units]
2842 Int accunits; // count of units accumulated
2843 Int accdigits; // count of digits accumulated
2844
2845 Unit varbuff[SD2U(DECBUFFER*2+DECDPUN)*sizeof(Unit)]; // buffer for var1
2846 Unit *var1=varbuff; // -> var1 array for long subtraction
2847 Unit *varalloc=NULL; // -> allocated buffer, iff used
2848 Unit *msu1; // -> msu of var1
2849
2850 const Unit *var2; // -> var2 array
2851 const Unit *msu2; // -> msu of var2
2852 Int msu2plus; // msu2 plus one [does not vary]
2853 eInt msu2pair; // msu2 pair plus one [does not vary]
2854
2855 Int var1units, var2units; // actual lengths
2856 Int var2ulen; // logical length (units)
2857 Int var1initpad=0; // var1 initial padding (digits)
2858 Int maxdigits; // longest LHS or required acc length
2859 Int mult; // multiplier for subtraction
2860 Unit thisunit; // current unit being accumulated
2861 Int residue; // for rounding
2862 Int reqdigits=set->digits; // requested DIGITS
2863 Int exponent; // working exponent
2864 Int maxexponent=0; // DIVIDE maximum exponent if unrounded
2865 uByte bits; // working sign
2866 Unit *target; // work
2867 const Unit *source; // ..
2868 uInt const *pow; // ..
2869 Int shift, cut; // ..
2870 #if DECSUBSET
2871 Int dropped; // work
2872 #endif
2873
2874 #if DECCHECK
2875 if (decCheckOperands(res, lhs, rhs, set)) return res;
2876 #endif
2877
2878 do { // protect allocated storage
2879 #if DECSUBSET
2880 if (!set->extended) {
2881 // reduce operands and set lostDigits status, as needed
2882 if (lhs->digits>reqdigits) {
2883 alloclhs=decRoundOperand(lhs, set, status);
2884 if (alloclhs==NULL) break;
2885 lhs=alloclhs;
2886 }
2887 if (rhs->digits>reqdigits) {
2888 allocrhs=decRoundOperand(rhs, set, status);
2889 if (allocrhs==NULL) break;
2890 rhs=allocrhs;
2891 }
2892 }
2893 #endif
2894 // [following code does not require input rounding]
2895
2896 bits=(lhs->bits^rhs->bits)&DECNEG; // assumed sign for divisions
2897
2898 // handle infinities and NaNs
2899 if (SPECIALARGS) { // a special bit set
2900 if (SPECIALARGS & (DECSNAN | DECNAN)) { // one or two NaNs
2901 decNaNs(res, lhs, rhs, status);
2902 break;
2903 }
2904 // one or two infinities
2905 if (decNumberIsInfinite(lhs)) { // LHS (dividend) is infinite
2906 if (decNumberIsInfinite(rhs) || // two infinities are invalid ..
2907 op & (REMAINDER | REMNEAR)) { // as is remainder of infinity
2908 *status|=DEC_Invalid_operation;
2909 break;
2910 }
2911 // [Note that infinity/0 raises no exceptions]
2912 decNumberZero(res);
2913 res->bits=bits|DECINF; // set +/- infinity
2914 break;
2915 }
2916 else { // RHS (divisor) is infinite
2917 residue=0;
2918 if (op&(REMAINDER|REMNEAR)) {
2919 // result is [finished clone of] lhs
2920 decCopyFit(res, lhs, set, &residue, status);
2921 }
2922 else { // a division
2923 decNumberZero(res);
2924 res->bits=bits; // set +/- zero
2925 // for DIVIDEINT the exponent is always 0. For DIVIDE, result
2926 // is a 0 with infinitely negative exponent, clamped to minimum
2927 if (op&DIVIDE) {
2928 res->exponent=set->emin-set->digits+1;
2929 *status|=DEC_Clamped;
2930 }
2931 }
2932 decFinish(res, set, &residue, status);
2933 break;
2934 }
2935 }
2936
2937 // handle 0 rhs (x/0)
2938 if (ISZERO(rhs)) { // x/0 is always exceptional
2939 if (ISZERO(lhs)) {
2940 decNumberZero(res); // [after lhs test]
2941 *status|=DEC_Division_undefined;// 0/0 will become NaN
2942 }
2943 else {
2944 decNumberZero(res);
2945 if (op&(REMAINDER|REMNEAR)) *status|=DEC_Invalid_operation;
2946 else {
2947 *status|=DEC_Division_by_zero; // x/0
2948 res->bits=bits|DECINF; // .. is +/- Infinity
2949 }
2950 }
2951 break;}
2952
2953 // handle 0 lhs (0/x)
2954 if (ISZERO(lhs)) { // 0/x [x!=0]
2955 #if DECSUBSET
2956 if (!set->extended) decNumberZero(res);
2957 else {
2958 #endif
2959 if (op&DIVIDE) {
2960 residue=0;
2961 exponent=lhs->exponent-rhs->exponent; // ideal exponent
2962 decNumberCopy(res, lhs); // [zeros always fit]
2963 res->bits=bits; // sign as computed
2964 res->exponent=exponent; // exponent, too
2965 decFinalize(res, set, &residue, status); // check exponent
2966 }
2967 else if (op&DIVIDEINT) {
2968 decNumberZero(res); // integer 0
2969 res->bits=bits; // sign as computed
2970 }
2971 else { // a remainder
2972 exponent=rhs->exponent; // [save in case overwrite]
2973 decNumberCopy(res, lhs); // [zeros always fit]
2974 if (exponent<res->exponent) res->exponent=exponent; // use lower
2975 }
2976 #if DECSUBSET
2977 }
2978 #endif
2979 break;}
2980
2981 // Precalculate exponent. This starts off adjusted (and hence fits
2982 // in 31 bits) and becomes the usual unadjusted exponent as the
2983 // division proceeds. The order of evaluation is important, here,
2984 // to avoid wrap.
2985 exponent=(lhs->exponent+lhs->digits)-(rhs->exponent+rhs->digits);
2986
2987 // If the working exponent is -ve, then some quick exits are
2988 // possible because the quotient is known to be <1
2989 // [for REMNEAR, it needs to be < -1, as -0.5 could need work]
2990 if (exponent<0 && !(op==DIVIDE)) {
2991 if (op&DIVIDEINT) {
2992 decNumberZero(res); // integer part is 0
2993 #if DECSUBSET
2994 if (set->extended)
2995 #endif
2996 res->bits=bits; // set +/- zero
2997 break;}
2998 // fastpath remainders so long as the lhs has the smaller
2999 // (or equal) exponent
3000 if (lhs->exponent<=rhs->exponent) {
3001 if (op&REMAINDER || exponent<-1) {
3002 // It is REMAINDER or safe REMNEAR; result is [finished
3003 // clone of] lhs (r = x - 0*y)
3004 residue=0;
3005 decCopyFit(res, lhs, set, &residue, status);
3006 decFinish(res, set, &residue, status);
3007 break;
3008 }
3009 // [unsafe REMNEAR drops through]
3010 }
3011 } // fastpaths
3012
3013 /* Long (slow) division is needed; roll up the sleeves... */
3014
3015 // The accumulator will hold the quotient of the division.
3016 // If it needs to be too long for stack storage, then allocate.
3017 acclength=D2U(reqdigits+DECDPUN); // in Units
3018 if (acclength*sizeof(Unit)>sizeof(accbuff)) {
3019 allocacc=(Unit *)malloc(acclength*sizeof(Unit));
3020 if (allocacc==NULL) { // hopeless -- abandon
3021 *status|=DEC_Insufficient_storage;
3022 break;}
3023 acc=allocacc; // use the allocated space
3024 }
3025
3026 // var1 is the padded LHS ready for subtractions.
3027 // If it needs to be too long for stack storage, then allocate.
3028 // The maximum units needed for var1 (long subtraction) is:
3029 // Enough for
3030 // (rhs->digits+reqdigits-1) -- to allow full slide to right
3031 // or (lhs->digits) -- to allow for long lhs
3032 // whichever is larger
3033 // +1 -- for rounding of slide to right
3034 // +1 -- for leading 0s
3035 // +1 -- for pre-adjust if a remainder or DIVIDEINT
3036 // [Note: unused units do not participate in decUnitAddSub data]
3037 maxdigits=rhs->digits+reqdigits-1;
3038 if (lhs->digits>maxdigits) maxdigits=lhs->digits;
3039 var1units=D2U(maxdigits)+2;
3040 // allocate a guard unit above msu1 for REMAINDERNEAR
3041 if (!(op&DIVIDE)) var1units++;
3042 if ((var1units+1)*sizeof(Unit)>sizeof(varbuff)) {
3043 varalloc=(Unit *)malloc((var1units+1)*sizeof(Unit));
3044 if (varalloc==NULL) { // hopeless -- abandon
3045 *status|=DEC_Insufficient_storage;
3046 break;}
3047 var1=varalloc; // use the allocated space
3048 }
3049
3050 // Extend the lhs and rhs to full long subtraction length. The lhs
3051 // is truly extended into the var1 buffer, with 0 padding, so a
3052 // subtract in place is always possible. The rhs (var2) has
3053 // virtual padding (implemented by decUnitAddSub).
3054 // One guard unit was allocated above msu1 for rem=rem+rem in
3055 // REMAINDERNEAR.
3056 msu1=var1+var1units-1; // msu of var1
3057 source=lhs->lsu+D2U(lhs->digits)-1; // msu of input array
3058 for (target=msu1; source>=lhs->lsu; source--, target--) *target=*source;
3059 for (; target>=var1; target--) *target=0;
3060
3061 // rhs (var2) is left-aligned with var1 at the start
3062 var2ulen=var1units; // rhs logical length (units)
3063 var2units=D2U(rhs->digits); // rhs actual length (units)
3064 var2=rhs->lsu; // -> rhs array
3065 msu2=var2+var2units-1; // -> msu of var2 [never changes]
3066 // now set up the variables which will be used for estimating the
3067 // multiplication factor. If these variables are not exact, add
3068 // 1 to make sure that the multiplier is never overestimated.
3069 msu2plus=*msu2; // it's value ..
3070 if (var2units>1) msu2plus++; // .. +1 if any more
3071 msu2pair=(eInt)*msu2*(DECDPUNMAX+1);// top two pair ..
3072 if (var2units>1) { // .. [else treat 2nd as 0]
3073 msu2pair+=*(msu2-1); // ..
3074 if (var2units>2) msu2pair++; // .. +1 if any more
3075 }
3076
3077 // The calculation is working in units, which may have leading zeros,
3078 // but the exponent was calculated on the assumption that they are
3079 // both left-aligned. Adjust the exponent to compensate: add the
3080 // number of leading zeros in var1 msu and subtract those in var2 msu.
3081 // [This is actually done by counting the digits and negating, as
3082 // lead1=DECDPUN-digits1, and similarly for lead2.]
3083 for (pow=&powers[1]; *msu1>=*pow; pow++) exponent--;
3084 for (pow=&powers[1]; *msu2>=*pow; pow++) exponent++;
3085
3086 // Now, if doing an integer divide or remainder, ensure that
3087 // the result will be Unit-aligned. To do this, shift the var1
3088 // accumulator towards least if need be. (It's much easier to
3089 // do this now than to reassemble the residue afterwards, if
3090 // doing a remainder.) Also ensure the exponent is not negative.
3091 if (!(op&DIVIDE)) {
3092 Unit *u; // work
3093 // save the initial 'false' padding of var1, in digits
3094 var1initpad=(var1units-D2U(lhs->digits))*DECDPUN;
3095 // Determine the shift to do.
3096 if (exponent<0) cut=-exponent;
3097 else cut=DECDPUN-exponent%DECDPUN;
3098 decShiftToLeast(var1, var1units, cut);
3099 exponent+=cut; // maintain numerical value
3100 var1initpad-=cut; // .. and reduce padding
3101 // clean any most-significant units which were just emptied
3102 for (u=msu1; cut>=DECDPUN; cut-=DECDPUN, u--) *u=0;
3103 } // align
3104 else { // is DIVIDE
3105 maxexponent=lhs->exponent-rhs->exponent; // save
3106 // optimization: if the first iteration will just produce 0,
3107 // preadjust to skip it [valid for DIVIDE only]
3108 if (*msu1<*msu2) {
3109 var2ulen--; // shift down
3110 exponent-=DECDPUN; // update the exponent
3111 }
3112 }
3113
3114 // ---- start the long-division loops ------------------------------
3115 accunits=0; // no units accumulated yet
3116 accdigits=0; // .. or digits
3117 accnext=acc+acclength-1; // -> msu of acc [NB: allows digits+1]
3118 for (;;) { // outer forever loop
3119 thisunit=0; // current unit assumed 0
3120 // find the next unit
3121 for (;;) { // inner forever loop
3122 // strip leading zero units [from either pre-adjust or from
3123 // subtract last time around]. Leave at least one unit.
3124 for (; *msu1==0 && msu1>var1; msu1--) var1units--;
3125
3126 if (var1units<var2ulen) break; // var1 too low for subtract
3127 if (var1units==var2ulen) { // unit-by-unit compare needed
3128 // compare the two numbers, from msu
3129 const Unit *pv1, *pv2;
3130 Unit v2; // units to compare
3131 pv2=msu2; // -> msu
3132 for (pv1=msu1; ; pv1--, pv2--) {
3133 // v1=*pv1 -- always OK
3134 v2=0; // assume in padding
3135 if (pv2>=var2) v2=*pv2; // in range
3136 if (*pv1!=v2) break; // no longer the same
3137 if (pv1==var1) break; // done; leave pv1 as is
3138 }
3139 // here when all inspected or a difference seen
3140 if (*pv1<v2) break; // var1 too low to subtract
3141 if (*pv1==v2) { // var1 == var2
3142 // reach here if var1 and var2 are identical; subtraction
3143 // would increase digit by one, and the residue will be 0 so
3144 // the calculation is done; leave the loop with residue=0.
3145 thisunit++; // as though subtracted
3146 *var1=0; // set var1 to 0
3147 var1units=1; // ..
3148 break; // from inner
3149 } // var1 == var2
3150 // *pv1>v2. Prepare for real subtraction; the lengths are equal
3151 // Estimate the multiplier (there's always a msu1-1)...
3152 // Bring in two units of var2 to provide a good estimate.
3153 mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2pair);
3154 } // lengths the same
3155 else { // var1units > var2ulen, so subtraction is safe
3156 // The var2 msu is one unit towards the lsu of the var1 msu,
3157 // so only one unit for var2 can be used.
3158 mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2plus);
3159 }
3160 if (mult==0) mult=1; // must always be at least 1
3161 // subtraction needed; var1 is > var2
3162 thisunit=(Unit)(thisunit+mult); // accumulate
3163 // subtract var1-var2, into var1; only the overlap needs
3164 // processing, as this is an in-place calculation
3165 shift=var2ulen-var2units;
3166 #if DECTRACE
3167 decDumpAr('1', &var1[shift], var1units-shift);
3168 decDumpAr('2', var2, var2units);
3169 printf("m=%d\n", -mult);
3170 #endif
3171 decUnitAddSub(&var1[shift], var1units-shift,
3172 var2, var2units, 0,
3173 &var1[shift], -mult);
3174 #if DECTRACE
3175 decDumpAr('#', &var1[shift], var1units-shift);
3176 #endif
3177 // var1 now probably has leading zeros; these are removed at the
3178 // top of the inner loop.
3179 } // inner loop
3180
3181 // The next unit has been calculated in full; unless it's a
3182 // leading zero, add to acc
3183 if (accunits!=0 || thisunit!=0) { // is first or non-zero
3184 *accnext=thisunit; // store in accumulator
3185 // account exactly for the new digits
3186 if (accunits==0) {
3187 accdigits++; // at least one
3188 for (pow=&powers[1]; thisunit>=*pow; pow++) accdigits++;
3189 }
3190 else accdigits+=DECDPUN;
3191 accunits++; // update count
3192 accnext--; // ready for next
3193 if (accdigits>reqdigits) break; // have enough digits
3194 }
3195
3196 // if the residue is zero, the operation is done (unless divide
3197 // or divideInteger and still not enough digits yet)
3198 if (*var1==0 && var1units==1) { // residue is 0
3199 if (op&(REMAINDER|REMNEAR)) break;
3200 if ((op&DIVIDE) && (exponent<=maxexponent)) break;
3201 // [drop through if divideInteger]
3202 }
3203 // also done enough if calculating remainder or integer
3204 // divide and just did the last ('units') unit
3205 if (exponent==0 && !(op&DIVIDE)) break;
3206
3207 // to get here, var1 is less than var2, so divide var2 by the per-
3208 // Unit power of ten and go for the next digit
3209 var2ulen--; // shift down
3210 exponent-=DECDPUN; // update the exponent
3211 } // outer loop
3212
3213 // ---- division is complete ---------------------------------------
3214 // here: acc has at least reqdigits+1 of good results (or fewer
3215 // if early stop), starting at accnext+1 (its lsu)
3216 // var1 has any residue at the stopping point
3217 // accunits is the number of digits collected in acc
3218 if (accunits==0) { // acc is 0
3219 accunits=1; // show have a unit ..
3220 accdigits=1; // ..
3221 *accnext=0; // .. whose value is 0
3222 }
3223 else accnext++; // back to last placed
3224 // accnext now -> lowest unit of result
3225
3226 residue=0; // assume no residue
3227 if (op&DIVIDE) {
3228 // record the presence of any residue, for rounding
3229 if (*var1!=0 || var1units>1) residue=1;
3230 else { // no residue
3231 // Had an exact division; clean up spurious trailing 0s.
3232 // There will be at most DECDPUN-1, from the final multiply,
3233 // and then only if the result is non-0 (and even) and the
3234 // exponent is 'loose'.
3235 #if DECDPUN>1
3236 Unit lsu=*accnext;
3237 if (!(lsu&0x01) && (lsu!=0)) {
3238 // count the trailing zeros
3239 Int drop=0;
3240 for (;; drop++) { // [will terminate because lsu!=0]
3241 if (exponent>=maxexponent) break; // don't chop real 0s
3242 #if DECDPUN<=4
3243 if ((lsu-QUOT10(lsu, drop+1)
3244 *powers[drop+1])!=0) break; // found non-0 digit
3245 #else
3246 if (lsu%powers[drop+1]!=0) break; // found non-0 digit
3247 #endif
3248 exponent++;
3249 }
3250 if (drop>0) {
3251 accunits=decShiftToLeast(accnext, accunits, drop);
3252 accdigits=decGetDigits(accnext, accunits);
3253 accunits=D2U(accdigits);
3254 // [exponent was adjusted in the loop]
3255 }
3256 } // neither odd nor 0
3257 #endif
3258 } // exact divide
3259 } // divide
3260 else /* op!=DIVIDE */ {
3261 // check for coefficient overflow
3262 if (accdigits+exponent>reqdigits) {
3263 *status|=DEC_Division_impossible;
3264 break;
3265 }
3266 if (op & (REMAINDER|REMNEAR)) {
3267 // [Here, the exponent will be 0, because var1 was adjusted
3268 // appropriately.]
3269 Int postshift; // work
3270 Flag wasodd=0; // integer was odd
3271 Unit *quotlsu; // for save
3272 Int quotdigits; // ..
3273
3274 // Fastpath when residue is truly 0 is worthwhile [and
3275 // simplifies the code below]
3276 if (*var1==0 && var1units==1) { // residue is 0
3277 Int exp=lhs->exponent; // save min(exponents)
3278 if (rhs->exponent<exp) exp=rhs->exponent;
3279 decNumberZero(res); // 0 coefficient
3280 #if DECSUBSET
3281 if (set->extended)
3282 #endif
3283 res->exponent=exp; // .. with proper exponent
3284 break;
3285 }
3286 // note if the quotient was odd
3287 if (*accnext & 0x01) wasodd=1; // acc is odd
3288 quotlsu=accnext; // save in case need to reinspect
3289 quotdigits=accdigits; // ..
3290
3291 // treat the residue, in var1, as the value to return, via acc
3292 // calculate the unused zero digits. This is the smaller of:
3293 // var1 initial padding (saved above)
3294 // var2 residual padding, which happens to be given by:
3295 postshift=var1initpad+exponent-lhs->exponent+rhs->exponent;
3296 // [the 'exponent' term accounts for the shifts during divide]
3297 if (var1initpad<postshift) postshift=var1initpad;
3298
3299 // shift var1 the requested amount, and adjust its digits
3300 var1units=decShiftToLeast(var1, var1units, postshift);
3301 accnext=var1;
3302 accdigits=decGetDigits(var1, var1units);
3303 accunits=D2U(accdigits);
3304
3305 exponent=lhs->exponent; // exponent is smaller of lhs & rhs
3306 if (rhs->exponent<exponent) exponent=rhs->exponent;
3307 bits=lhs->bits; // remainder sign is always as lhs
3308
3309 // Now correct the result if doing remainderNear; if it
3310 // (looking just at coefficients) is > rhs/2, or == rhs/2 and
3311 // the integer was odd then the result should be rem-rhs.
3312 if (op&REMNEAR) {
3313 Int compare, tarunits; // work
3314 Unit *up; // ..
3315
3316
3317 // calculate remainder*2 into the var1 buffer (which has
3318 // 'headroom' of an extra unit and hence enough space)
3319 // [a dedicated 'double' loop would be faster, here]
3320 tarunits=decUnitAddSub(accnext, accunits, accnext, accunits,
3321 0, accnext, 1);
3322 // decDumpAr('r', accnext, tarunits);
3323
3324 // Here, accnext (var1) holds tarunits Units with twice the
3325 // remainder's coefficient, which must now be compared to the
3326 // RHS. The remainder's exponent may be smaller than the RHS's.
3327 compare=decUnitCompare(accnext, tarunits, rhs->lsu, D2U(rhs->digits),
3328 rhs->exponent-exponent);
3329 if (compare==BADINT) { // deep trouble
3330 *status|=DEC_Insufficient_storage;
3331 break;}
3332
3333 // now restore the remainder by dividing by two; the lsu
3334 // is known to be even.
3335 for (up=accnext; up<accnext+tarunits; up++) {
3336 Int half; // half to add to lower unit
3337 half=*up & 0x01;
3338 *up/=2; // [shift]
3339 if (!half) continue;
3340 *(up-1)+=(DECDPUNMAX+1)/2;
3341 }
3342 // [accunits still describes the original remainder length]
3343
3344 if (compare>0 || (compare==0 && wasodd)) { // adjustment needed
3345 Int exp, expunits, exprem; // work
3346 // This is effectively causing round-up of the quotient,
3347 // so if it was the rare case where it was full and all
3348 // nines, it would overflow and hence division-impossible
3349 // should be raised
3350 Flag allnines=0; // 1 if quotient all nines
3351 if (quotdigits==reqdigits) { // could be borderline
3352 for (up=quotlsu; ; up++) {
3353 if (quotdigits>DECDPUN) {
3354 if (*up!=DECDPUNMAX) break;// non-nines
3355 }
3356 else { // this is the last Unit
3357 if (*up==powers[quotdigits]-1) allnines=1;
3358 break;
3359 }
3360 quotdigits-=DECDPUN; // checked those digits
3361 } // up
3362 } // borderline check
3363 if (allnines) {
3364 *status|=DEC_Division_impossible;
3365 break;}
3366
3367 // rem-rhs is needed; the sign will invert. Again, var1
3368 // can safely be used for the working Units array.
3369 exp=rhs->exponent-exponent; // RHS padding needed
3370 // Calculate units and remainder from exponent.
3371 expunits=exp/DECDPUN;
3372 exprem=exp%DECDPUN;
3373 // subtract [A+B*(-m)]; the result will always be negative
3374 accunits=-decUnitAddSub(accnext, accunits,
3375 rhs->lsu, D2U(rhs->digits),
3376 expunits, accnext, -(Int)powers[exprem]);
3377 accdigits=decGetDigits(accnext, accunits); // count digits exactly
3378 accunits=D2U(accdigits); // and recalculate the units for copy
3379 // [exponent is as for original remainder]
3380 bits^=DECNEG; // flip the sign
3381 }
3382 } // REMNEAR
3383 } // REMAINDER or REMNEAR
3384 } // not DIVIDE
3385
3386 // Set exponent and bits
3387 res->exponent=exponent;
3388 res->bits=(uByte)(bits&DECNEG); // [cleaned]
3389
3390 // Now the coefficient.
3391 decSetCoeff(res, set, accnext, accdigits, &residue, status);
3392
3393 decFinish(res, set, &residue, status); // final cleanup
3394
3395 #if DECSUBSET
3396 // If a divide then strip trailing zeros if subset [after round]
3397 if (!set->extended && (op==DIVIDE)) decTrim(res, 0, &dropped);
3398 #endif
3399 } while(0); // end protected
3400
3401 if (varalloc!=NULL) free(varalloc); // drop any storage used
3402 if (allocacc!=NULL) free(allocacc); // ..
3403 #if DECSUBSET
3404 if (allocrhs!=NULL) free(allocrhs); // ..
3405 if (alloclhs!=NULL) free(alloclhs); // ..
3406 #endif
3407 return res;
3408 } // decDivideOp
3409
3410 /* ------------------------------------------------------------------ */
3411 /* decMultiplyOp -- multiplication operation */
3412 /* */
3413 /* This routine performs the multiplication C=A x B. */
3414 /* */
3415 /* res is C, the result. C may be A and/or B (e.g., X=X*X) */
3416 /* lhs is A */
3417 /* rhs is B */
3418 /* set is the context */
3419 /* status is the usual accumulator */
3420 /* */
3421 /* C must have space for set->digits digits. */
3422 /* */
3423 /* ------------------------------------------------------------------ */
3424 /* 'Classic' multiplication is used rather than Karatsuba, as the */
3425 /* latter would give only a minor improvement for the short numbers */
3426 /* expected to be handled most (and uses much more memory). */
3427 /* */
3428 /* There are two major paths here: the general-purpose ('old code') */
3429 /* path which handles all DECDPUN values, and a fastpath version */
3430 /* which is used if 64-bit ints are available, DECDPUN<=4, and more */
3431 /* than two calls to decUnitAddSub would be made. */
3432 /* */
3433 /* The fastpath version lumps units together into 8-digit or 9-digit */
3434 /* chunks, and also uses a lazy carry strategy to minimise expensive */
3435 /* 64-bit divisions. The chunks are then broken apart again into */
3436 /* units for continuing processing. Despite this overhead, the */
3437 /* fastpath can speed up some 16-digit operations by 10x (and much */
3438 /* more for higher-precision calculations). */
3439 /* */
3440 /* A buffer always has to be used for the accumulator; in the */
3441 /* fastpath, buffers are also always needed for the chunked copies of */
3442 /* of the operand coefficients. */
3443 /* ------------------------------------------------------------------ */
3444 #define FASTMUL (DECUSE64 && DECDPUN<5)
decMultiplyOp(decNumber * res,const decNumber * lhs,const decNumber * rhs,decContext * set,uInt * status)3445 static decNumber * decMultiplyOp(decNumber *res, const decNumber *lhs,
3446 const decNumber *rhs, decContext *set,
3447 uInt *status) {
3448 Int accunits; // Units of accumulator in use
3449 Int exponent; // work
3450 Int residue=0; // rounding residue
3451 uByte bits; // result sign
3452 Unit *acc; // -> accumulator Unit array
3453 Int needbytes; // size calculator
3454 void *allocacc=NULL; // -> allocated accumulator, iff allocated
3455 Unit accbuff[SD2U(DECBUFFER*2+5)]; // buffer (+1 for DECBUFFER==0,
3456 // + 4 for calls from other operations)
3457 const Unit *mer, *mermsup; // work
3458 Int madlength; // Units in multiplicand
3459 Int shift; // Units to shift multiplicand by
3460
3461 #if FASTMUL
3462 // if DECDPUN is 1 or 3 work in base 10**9, otherwise
3463 // (DECDPUN is 2 or 4) then work in base 10**8
3464 #if DECDPUN & 1 // odd
3465 #define FASTBASE 1000000000 // base
3466 #define FASTDIGS 9 // digits in base
3467 #define FASTLAZY 18 // carry resolution point [1->18]
3468 #else
3469 #define FASTBASE 100000000
3470 #define FASTDIGS 8
3471 #define FASTLAZY 1844 // carry resolution point [1->1844]
3472 #endif
3473 // three buffers are used, two for chunked copies of the operands
3474 // (base 10**8 or base 10**9) and one base 2**64 accumulator with
3475 // lazy carry evaluation
3476 uInt zlhibuff[(DECBUFFER+7)/8+1]; // buffer (+1 for DECBUFFER==0)
3477 uInt *zlhi=zlhibuff; // -> lhs array
3478 uInt *alloclhi=NULL; // -> allocated buffer, iff allocated
3479 uInt zrhibuff[(DECBUFFER+7)/8+1]; // buffer (+1 for DECBUFFER==0)
3480 uInt *zrhi=zrhibuff; // -> rhs array
3481 uInt *allocrhi=NULL; // -> allocated buffer, iff allocated
3482 uLong zaccbuff[(DECBUFFER+3)/4+2]; // buffer (+1 for DECBUFFER==0)
3483 // + 1 for calls from other operations)
3484 // [allocacc is shared for both paths, as only one will run]
3485 uLong *zacc=zaccbuff; // -> accumulator array for exact result
3486 #if DECDPUN==1
3487 Int zoff; // accumulator offset
3488 #endif
3489 uInt *lip, *rip; // item pointers
3490 uInt *lmsi, *rmsi; // most significant items
3491 Int ilhs, irhs, iacc; // item counts in the arrays
3492 Int lazy; // lazy carry counter
3493 uLong lcarry; // uLong carry
3494 uInt carry; // carry (NB not uLong)
3495 Int count; // work
3496 const Unit *cup; // ..
3497 Unit *up; // ..
3498 uLong *lp; // ..
3499 Int p; // ..
3500 #endif
3501
3502 #if DECSUBSET
3503 decNumber *alloclhs=NULL; // -> allocated buffer, iff allocated
3504 decNumber *allocrhs=NULL; // -> allocated buffer, iff allocated
3505 #endif
3506
3507 #if DECCHECK
3508 if (decCheckOperands(res, lhs, rhs, set)) return res;
3509 #endif
3510
3511 // precalculate result sign
3512 bits=(uByte)((lhs->bits^rhs->bits)&DECNEG);
3513
3514 // handle infinities and NaNs
3515 if (SPECIALARGS) { // a special bit set
3516 if (SPECIALARGS & (DECSNAN | DECNAN)) { // one or two NaNs
3517 decNaNs(res, lhs, rhs, status);
3518 return res;}
3519 // one or two infinities; Infinity * 0 is invalid
3520 if (((lhs->bits & DECINF)==0 && ISZERO(lhs))
3521 ||((rhs->bits & DECINF)==0 && ISZERO(rhs))) {
3522 *status|=DEC_Invalid_operation;
3523 return res;}
3524 decNumberZero(res);
3525 res->bits=bits|DECINF; // infinity
3526 return res;}
3527
3528 // For best speed, as in DMSRCN [the original Rexx numerics
3529 // module], use the shorter number as the multiplier (rhs) and
3530 // the longer as the multiplicand (lhs) to minimise the number of
3531 // adds (partial products)
3532 if (lhs->digits<rhs->digits) { // swap...
3533 const decNumber *hold=lhs;
3534 lhs=rhs;
3535 rhs=hold;
3536 }
3537
3538 do { // protect allocated storage
3539 #if DECSUBSET
3540 if (!set->extended) {
3541 // reduce operands and set lostDigits status, as needed
3542 if (lhs->digits>set->digits) {
3543 alloclhs=decRoundOperand(lhs, set, status);
3544 if (alloclhs==NULL) break;
3545 lhs=alloclhs;
3546 }
3547 if (rhs->digits>set->digits) {
3548 allocrhs=decRoundOperand(rhs, set, status);
3549 if (allocrhs==NULL) break;
3550 rhs=allocrhs;
3551 }
3552 }
3553 #endif
3554 // [following code does not require input rounding]
3555
3556 #if FASTMUL // fastpath can be used
3557 // use the fast path if there are enough digits in the shorter
3558 // operand to make the setup and takedown worthwhile
3559 #define NEEDTWO (DECDPUN*2) // within two decUnitAddSub calls
3560 if (rhs->digits>NEEDTWO) { // use fastpath...
3561 // calculate the number of elements in each array
3562 ilhs=(lhs->digits+FASTDIGS-1)/FASTDIGS; // [ceiling]
3563 irhs=(rhs->digits+FASTDIGS-1)/FASTDIGS; // ..
3564 iacc=ilhs+irhs;
3565
3566 // allocate buffers if required, as usual
3567 needbytes=ilhs*sizeof(uInt);
3568 if (needbytes>(Int)sizeof(zlhibuff)) {
3569 alloclhi=(uInt *)malloc(needbytes);
3570 zlhi=alloclhi;}
3571 needbytes=irhs*sizeof(uInt);
3572 if (needbytes>(Int)sizeof(zrhibuff)) {
3573 allocrhi=(uInt *)malloc(needbytes);
3574 zrhi=allocrhi;}
3575
3576 // Allocating the accumulator space needs a special case when
3577 // DECDPUN=1 because when converting the accumulator to Units
3578 // after the multiplication each 8-byte item becomes 9 1-byte
3579 // units. Therefore iacc extra bytes are needed at the front
3580 // (rounded up to a multiple of 8 bytes), and the uLong
3581 // accumulator starts offset the appropriate number of units
3582 // to the right to avoid overwrite during the unchunking.
3583 needbytes=iacc*sizeof(uLong);
3584 #if DECDPUN==1
3585 zoff=(iacc+7)/8; // items to offset by
3586 needbytes+=zoff*8;
3587 #endif
3588 if (needbytes>(Int)sizeof(zaccbuff)) {
3589 allocacc=(uLong *)malloc(needbytes);
3590 zacc=(uLong *)allocacc;}
3591 if (zlhi==NULL||zrhi==NULL||zacc==NULL) {
3592 *status|=DEC_Insufficient_storage;
3593 break;}
3594
3595 acc=(Unit *)zacc; // -> target Unit array
3596 #if DECDPUN==1
3597 zacc+=zoff; // start uLong accumulator to right
3598 #endif
3599
3600 // assemble the chunked copies of the left and right sides
3601 for (count=lhs->digits, cup=lhs->lsu, lip=zlhi; count>0; lip++)
3602 for (p=0, *lip=0; p<FASTDIGS && count>0;
3603 p+=DECDPUN, cup++, count-=DECDPUN)
3604 *lip+=*cup*powers[p];
3605 lmsi=lip-1; // save -> msi
3606 for (count=rhs->digits, cup=rhs->lsu, rip=zrhi; count>0; rip++)
3607 for (p=0, *rip=0; p<FASTDIGS && count>0;
3608 p+=DECDPUN, cup++, count-=DECDPUN)
3609 *rip+=*cup*powers[p];
3610 rmsi=rip-1; // save -> msi
3611
3612 // zero the accumulator
3613 for (lp=zacc; lp<zacc+iacc; lp++) *lp=0;
3614
3615 /* Start the multiplication */
3616 // Resolving carries can dominate the cost of accumulating the
3617 // partial products, so this is only done when necessary.
3618 // Each uLong item in the accumulator can hold values up to
3619 // 2**64-1, and each partial product can be as large as
3620 // (10**FASTDIGS-1)**2. When FASTDIGS=9, this can be added to
3621 // itself 18.4 times in a uLong without overflowing, so during
3622 // the main calculation resolution is carried out every 18th
3623 // add -- every 162 digits. Similarly, when FASTDIGS=8, the
3624 // partial products can be added to themselves 1844.6 times in
3625 // a uLong without overflowing, so intermediate carry
3626 // resolution occurs only every 14752 digits. Hence for common
3627 // short numbers usually only the one final carry resolution
3628 // occurs.
3629 // (The count is set via FASTLAZY to simplify experiments to
3630 // measure the value of this approach: a 35% improvement on a
3631 // [34x34] multiply.)
3632 lazy=FASTLAZY; // carry delay count
3633 for (rip=zrhi; rip<=rmsi; rip++) { // over each item in rhs
3634 lp=zacc+(rip-zrhi); // where to add the lhs
3635 for (lip=zlhi; lip<=lmsi; lip++, lp++) { // over each item in lhs
3636 *lp+=(uLong)(*lip)*(*rip); // [this should in-line]
3637 } // lip loop
3638 lazy--;
3639 if (lazy>0 && rip!=rmsi) continue;
3640 lazy=FASTLAZY; // reset delay count
3641 // spin up the accumulator resolving overflows
3642 for (lp=zacc; lp<zacc+iacc; lp++) {
3643 if (*lp<FASTBASE) continue; // it fits
3644 lcarry=*lp/FASTBASE; // top part [slow divide]
3645 // lcarry can exceed 2**32-1, so check again; this check
3646 // and occasional extra divide (slow) is well worth it, as
3647 // it allows FASTLAZY to be increased to 18 rather than 4
3648 // in the FASTDIGS=9 case
3649 if (lcarry<FASTBASE) carry=(uInt)lcarry; // [usual]
3650 else { // two-place carry [fairly rare]
3651 uInt carry2=(uInt)(lcarry/FASTBASE); // top top part
3652 *(lp+2)+=carry2; // add to item+2
3653 *lp-=((uLong)FASTBASE*FASTBASE*carry2); // [slow]
3654 carry=(uInt)(lcarry-((uLong)FASTBASE*carry2)); // [inline]
3655 }
3656 *(lp+1)+=carry; // add to item above [inline]
3657 *lp-=((uLong)FASTBASE*carry); // [inline]
3658 } // carry resolution
3659 } // rip loop
3660
3661 // The multiplication is complete; time to convert back into
3662 // units. This can be done in-place in the accumulator and in
3663 // 32-bit operations, because carries were resolved after the
3664 // final add. This needs N-1 divides and multiplies for
3665 // each item in the accumulator (which will become up to N
3666 // units, where 2<=N<=9).
3667 for (lp=zacc, up=acc; lp<zacc+iacc; lp++) {
3668 uInt item=(uInt)*lp; // decapitate to uInt
3669 for (p=0; p<FASTDIGS-DECDPUN; p+=DECDPUN, up++) {
3670 uInt part=item/(DECDPUNMAX+1);
3671 *up=(Unit)(item-(part*(DECDPUNMAX+1)));
3672 item=part;
3673 } // p
3674 *up=(Unit)item; up++; // [final needs no division]
3675 } // lp
3676 accunits=up-acc; // count of units
3677 }
3678 else { // here to use units directly, without chunking ['old code']
3679 #endif
3680
3681 // if accumulator will be too long for local storage, then allocate
3682 acc=accbuff; // -> assume buffer for accumulator
3683 needbytes=(D2U(lhs->digits)+D2U(rhs->digits))*sizeof(Unit);
3684 if (needbytes>(Int)sizeof(accbuff)) {
3685 allocacc=(Unit *)malloc(needbytes);
3686 if (allocacc==NULL) {*status|=DEC_Insufficient_storage; break;}
3687 acc=(Unit *)allocacc; // use the allocated space
3688 }
3689
3690 /* Now the main long multiplication loop */
3691 // Unlike the equivalent in the IBM Java implementation, there
3692 // is no advantage in calculating from msu to lsu. So, do it
3693 // by the book, as it were.
3694 // Each iteration calculates ACC=ACC+MULTAND*MULT
3695 accunits=1; // accumulator starts at '0'
3696 *acc=0; // .. (lsu=0)
3697 shift=0; // no multiplicand shift at first
3698 madlength=D2U(lhs->digits); // this won't change
3699 mermsup=rhs->lsu+D2U(rhs->digits); // -> msu+1 of multiplier
3700
3701 for (mer=rhs->lsu; mer<mermsup; mer++) {
3702 // Here, *mer is the next Unit in the multiplier to use
3703 // If non-zero [optimization] add it...
3704 if (*mer!=0) accunits=decUnitAddSub(&acc[shift], accunits-shift,
3705 lhs->lsu, madlength, 0,
3706 &acc[shift], *mer)
3707 + shift;
3708 else { // extend acc with a 0; it will be used shortly
3709 *(acc+accunits)=0; // [this avoids length of <=0 later]
3710 accunits++;
3711 }
3712 // multiply multiplicand by 10**DECDPUN for next Unit to left
3713 shift++; // add this for 'logical length'
3714 } // n
3715 #if FASTMUL
3716 } // unchunked units
3717 #endif
3718 // common end-path
3719 #if DECTRACE
3720 decDumpAr('*', acc, accunits); // Show exact result
3721 #endif
3722
3723 // acc now contains the exact result of the multiplication,
3724 // possibly with a leading zero unit; build the decNumber from
3725 // it, noting if any residue
3726 res->bits=bits; // set sign
3727 res->digits=decGetDigits(acc, accunits); // count digits exactly
3728
3729 // There can be a 31-bit wrap in calculating the exponent.
3730 // This can only happen if both input exponents are negative and
3731 // both their magnitudes are large. If there was a wrap, set a
3732 // safe very negative exponent, from which decFinalize() will
3733 // raise a hard underflow shortly.
3734 exponent=lhs->exponent+rhs->exponent; // calculate exponent
3735 if (lhs->exponent<0 && rhs->exponent<0 && exponent>0)
3736 exponent=-2*DECNUMMAXE; // force underflow
3737 res->exponent=exponent; // OK to overwrite now
3738
3739 // Set the coefficient. If any rounding, residue records
3740 decSetCoeff(res, set, acc, res->digits, &residue, status);
3741 decFinish(res, set, &residue, status); // final cleanup
3742 } while(0); // end protected
3743
3744 if (allocacc!=NULL) free(allocacc); // drop any storage used
3745 #if DECSUBSET
3746 if (allocrhs!=NULL) free(allocrhs); // ..
3747 if (alloclhs!=NULL) free(alloclhs); // ..
3748 #endif
3749 #if FASTMUL
3750 if (allocrhi!=NULL) free(allocrhi); // ..
3751 if (alloclhi!=NULL) free(alloclhi); // ..
3752 #endif
3753 return res;
3754 } // decMultiplyOp
3755
3756 /* ------------------------------------------------------------------ */
3757 /* decExpOp -- effect exponentiation */
3758 /* */
3759 /* This computes C = exp(A) */
3760 /* */
3761 /* res is C, the result. C may be A */
3762 /* rhs is A */
3763 /* set is the context; note that rounding mode has no effect */
3764 /* */
3765 /* C must have space for set->digits digits. status is updated but */
3766 /* not set. */
3767 /* */
3768 /* Restrictions: */
3769 /* */
3770 /* digits, emax, and -emin in the context must be less than */
3771 /* 2*DEC_MAX_MATH (1999998), and the rhs must be within these */
3772 /* bounds or a zero. This is an internal routine, so these */
3773 /* restrictions are contractual and not enforced. */
3774 /* */
3775 /* A finite result is rounded using DEC_ROUND_HALF_EVEN; it will */
3776 /* almost always be correctly rounded, but may be up to 1 ulp in */
3777 /* error in rare cases. */
3778 /* */
3779 /* Finite results will always be full precision and Inexact, except */
3780 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
3781 /* ------------------------------------------------------------------ */
3782 /* This approach used here is similar to the algorithm described in */
3783 /* */
3784 /* Variable Precision Exponential Function, T. E. Hull and */
3785 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 12 #2, */
3786 /* pp79-91, ACM, June 1986. */
3787 /* */
3788 /* with the main difference being that the iterations in the series */
3789 /* evaluation are terminated dynamically (which does not require the */
3790 /* extra variable-precision variables which are expensive in this */
3791 /* context). */
3792 /* */
3793 /* The error analysis in Hull & Abrham's paper applies except for the */
3794 /* round-off error accumulation during the series evaluation. This */
3795 /* code does not precalculate the number of iterations and so cannot */
3796 /* use Horner's scheme. Instead, the accumulation is done at double- */
3797 /* precision, which ensures that the additions of the terms are exact */
3798 /* and do not accumulate round-off (and any round-off errors in the */
3799 /* terms themselves move 'to the right' faster than they can */
3800 /* accumulate). This code also extends the calculation by allowing, */
3801 /* in the spirit of other decNumber operators, the input to be more */
3802 /* precise than the result (the precision used is based on the more */
3803 /* precise of the input or requested result). */
3804 /* */
3805 /* Implementation notes: */
3806 /* */
3807 /* 1. This is separated out as decExpOp so it can be called from */
3808 /* other Mathematical functions (notably Ln) with a wider range */
3809 /* than normal. In particular, it can handle the slightly wider */
3810 /* (double) range needed by Ln (which has to be able to calculate */
3811 /* exp(-x) where x can be the tiniest number (Ntiny). */
3812 /* */
3813 /* 2. Normalizing x to be <=0.1 (instead of <=1) reduces loop */
3814 /* iterations by appoximately a third with additional (although */
3815 /* diminishing) returns as the range is reduced to even smaller */
3816 /* fractions. However, h (the power of 10 used to correct the */
3817 /* result at the end, see below) must be kept <=8 as otherwise */
3818 /* the final result cannot be computed. Hence the leverage is a */
3819 /* sliding value (8-h), where potentially the range is reduced */
3820 /* more for smaller values. */
3821 /* */
3822 /* The leverage that can be applied in this way is severely */
3823 /* limited by the cost of the raise-to-the power at the end, */
3824 /* which dominates when the number of iterations is small (less */
3825 /* than ten) or when rhs is short. As an example, the adjustment */
3826 /* x**10,000,000 needs 31 multiplications, all but one full-width. */
3827 /* */
3828 /* 3. The restrictions (especially precision) could be raised with */
3829 /* care, but the full decNumber range seems very hard within the */
3830 /* 32-bit limits. */
3831 /* ------------------------------------------------------------------ */
decExpOp(decNumber * res,const decNumber * rhs,decContext * set,uInt * status)3832 decNumber * decExpOp(decNumber *res, const decNumber *rhs,
3833 decContext *set, uInt *status) {
3834 uInt ignore=0; // working status
3835 Int h; // adjusted exponent for 0.xxxx
3836 Int p; // working precision
3837 Int residue; // rounding residue
3838 uInt needbytes; // for space calculations
3839 const decNumber *x=rhs; // (may point to safe copy later)
3840 decContext aset, tset, dset; // working contexts
3841
3842 // the argument is often copied to normalize it, so (unusually) it
3843 // is treated like other buffers, using DECBUFFER, +1 in case
3844 // DECBUFFER is 0
3845 decNumber bufr[D2N(DECBUFFER+1)];
3846 decNumber *allocrhs=NULL; // non-NULL if rhs buffer allocated
3847
3848 // the working precision will be no more than set->digits+8+1
3849 // so for on-stack buffers DECBUFFER+9 is used, +1 in case DECBUFFER
3850 // is 0 (and twice that for the accumulator)
3851
3852 // buffer for t, term (working precision plus)
3853 decNumber buft[D2N(DECBUFFER+9+1)];
3854 decNumber *allocbuft=NULL; // -> allocated buft, iff allocated
3855 decNumber *t=buft; // term
3856 // buffer for a, accumulator (working precision * 2), at least 9
3857 decNumber bufa[D2N(DECBUFFER*2+18+1)];
3858 decNumber *allocbufa=NULL; // -> allocated bufa, iff allocated
3859 decNumber *a=bufa; // accumulator
3860 // decNumber for the divisor term; this needs at most 9 digits
3861 // and so can be fixed size [16 so can use standard context]
3862 decNumber bufd[D2N(16)];
3863 decNumber *d=bufd; // divisor
3864 decNumber numone; // constant 1
3865
3866 #if DECCHECK
3867 Int iterations=0; // for later sanity check
3868 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
3869 #endif
3870
3871 do { // protect allocated storage
3872 if (SPECIALARG) { // handle infinities and NaNs
3873 if (decNumberIsInfinite(rhs)) { // an infinity
3874 if (decNumberIsNegative(rhs)) // -Infinity -> +0
3875 decNumberZero(res);
3876 else decNumberCopy(res, rhs); // +Infinity -> self
3877 }
3878 else decNaNs(res, rhs, NULL, status); // a NaN
3879 break;}
3880
3881 if (ISZERO(rhs)) { // zeros -> exact 1
3882 decNumberZero(res); // make clean 1
3883 *res->lsu=1; // ..
3884 break;} // [no status to set]
3885
3886 // e**x when 0 < x < 0.66 is < 1+3x/2, hence can fast-path
3887 // positive and negative tiny cases which will result in inexact
3888 // 1. This also allows the later add-accumulate to always be
3889 // exact (because its length will never be more than twice the
3890 // working precision).
3891 // The comparator (tiny) needs just one digit, so use the
3892 // decNumber d for it (reused as the divisor, etc., below); its
3893 // exponent is such that if x is positive it will have
3894 // set->digits-1 zeros between the decimal point and the digit,
3895 // which is 4, and if x is negative one more zero there as the
3896 // more precise result will be of the form 0.9999999 rather than
3897 // 1.0000001. Hence, tiny will be 0.0000004 if digits=7 and x>0
3898 // or 0.00000004 if digits=7 and x<0. If RHS not larger than
3899 // this then the result will be 1.000000
3900 decNumberZero(d); // clean
3901 *d->lsu=4; // set 4 ..
3902 d->exponent=-set->digits; // * 10**(-d)
3903 if (decNumberIsNegative(rhs)) d->exponent--; // negative case
3904 if (decCompare(d, rhs, 1)>=0) { // signless compare
3905 Int shift=set->digits-1;
3906 decNumberZero(res); // set 1
3907 *res->lsu=1; // ..
3908 res->digits=decShiftToMost(res->lsu, 1, shift);
3909 res->exponent=-shift; // make 1.0000...
3910 *status|=DEC_Inexact | DEC_Rounded; // .. inexactly
3911 break;} // tiny
3912
3913 // set up the context to be used for calculating a, as this is
3914 // used on both paths below
3915 decContextDefault(&aset, DEC_INIT_DECIMAL64);
3916 // accumulator bounds are as requested (could underflow)
3917 aset.emax=set->emax; // usual bounds
3918 aset.emin=set->emin; // ..
3919 aset.clamp=0; // and no concrete format
3920
3921 // calculate the adjusted (Hull & Abrham) exponent (where the
3922 // decimal point is just to the left of the coefficient msd)
3923 h=rhs->exponent+rhs->digits;
3924 // if h>8 then 10**h cannot be calculated safely; however, when
3925 // h=8 then exp(|rhs|) will be at least exp(1E+7) which is at
3926 // least 6.59E+4342944, so (due to the restriction on Emax/Emin)
3927 // overflow (or underflow to 0) is guaranteed -- so this case can
3928 // be handled by simply forcing the appropriate excess
3929 if (h>8) { // overflow/underflow
3930 // set up here so Power call below will over or underflow to
3931 // zero; set accumulator to either 2 or 0.02
3932 // [stack buffer for a is always big enough for this]
3933 decNumberZero(a);
3934 *a->lsu=2; // not 1 but < exp(1)
3935 if (decNumberIsNegative(rhs)) a->exponent=-2; // make 0.02
3936 h=8; // clamp so 10**h computable
3937 p=9; // set a working precision
3938 }
3939 else { // h<=8
3940 Int maxlever=(rhs->digits>8?1:0);
3941 // [could/should increase this for precisions >40 or so, too]
3942
3943 // if h is 8, cannot normalize to a lower upper limit because
3944 // the final result will not be computable (see notes above),
3945 // but leverage can be applied whenever h is less than 8.
3946 // Apply as much as possible, up to a MAXLEVER digits, which
3947 // sets the tradeoff against the cost of the later a**(10**h).
3948 // As h is increased, the working precision below also
3949 // increases to compensate for the "constant digits at the
3950 // front" effect.
3951 Int lever=MIN(8-h, maxlever); // leverage attainable
3952 Int use=-rhs->digits-lever; // exponent to use for RHS
3953 h+=lever; // apply leverage selected
3954 if (h<0) { // clamp
3955 use+=h; // [may end up subnormal]
3956 h=0;
3957 }
3958 // Take a copy of RHS if it needs normalization (true whenever x>=1)
3959 if (rhs->exponent!=use) {
3960 decNumber *newrhs=bufr; // assume will fit on stack
3961 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
3962 if (needbytes>sizeof(bufr)) { // need malloc space
3963 allocrhs=(decNumber *)malloc(needbytes);
3964 if (allocrhs==NULL) { // hopeless -- abandon
3965 *status|=DEC_Insufficient_storage;
3966 break;}
3967 newrhs=allocrhs; // use the allocated space
3968 }
3969 decNumberCopy(newrhs, rhs); // copy to safe space
3970 newrhs->exponent=use; // normalize; now <1
3971 x=newrhs; // ready for use
3972 // decNumberShow(x);
3973 }
3974
3975 // Now use the usual power series to evaluate exp(x). The
3976 // series starts as 1 + x + x^2/2 ... so prime ready for the
3977 // third term by setting the term variable t=x, the accumulator
3978 // a=1, and the divisor d=2.
3979
3980 // First determine the working precision. From Hull & Abrham
3981 // this is set->digits+h+2. However, if x is 'over-precise' we
3982 // need to allow for all its digits to potentially participate
3983 // (consider an x where all the excess digits are 9s) so in
3984 // this case use x->digits+h+2
3985 p=MAX(x->digits, set->digits)+h+2;
3986
3987 // a and t are variable precision, and depend on p, so space
3988 // must be allocated for them if necessary
3989
3990 // the accumulator needs to be able to hold 2p digits so that
3991 // the additions on the second and subsequent iterations are
3992 // sufficiently exact.
3993 needbytes=sizeof(decNumber)+(D2U(p*2)-1)*sizeof(Unit);
3994 if (needbytes>sizeof(bufa)) { // need malloc space
3995 allocbufa=(decNumber *)malloc(needbytes);
3996 if (allocbufa==NULL) { // hopeless -- abandon
3997 *status|=DEC_Insufficient_storage;
3998 break;}
3999 a=allocbufa; // use the allocated space
4000 }
4001 // the term needs to be able to hold p digits (which is
4002 // guaranteed to be larger than x->digits, so the initial copy
4003 // is safe); it may also be used for the raise-to-power
4004 // calculation below, which needs an extra two digits
4005 needbytes=sizeof(decNumber)+(D2U(p+2)-1)*sizeof(Unit);
4006 if (needbytes>sizeof(buft)) { // need malloc space
4007 allocbuft=(decNumber *)malloc(needbytes);
4008 if (allocbuft==NULL) { // hopeless -- abandon
4009 *status|=DEC_Insufficient_storage;
4010 break;}
4011 t=allocbuft; // use the allocated space
4012 }
4013
4014 decNumberCopy(t, x); // term=x
4015 decNumberZero(a); *a->lsu=1; // accumulator=1
4016 decNumberZero(d); *d->lsu=2; // divisor=2
4017 decNumberZero(&numone); *numone.lsu=1; // constant 1 for increment
4018
4019 // set up the contexts for calculating a, t, and d
4020 decContextDefault(&tset, DEC_INIT_DECIMAL64);
4021 dset=tset;
4022 // accumulator bounds are set above, set precision now
4023 aset.digits=p*2; // double
4024 // term bounds avoid any underflow or overflow
4025 tset.digits=p;
4026 tset.emin=DEC_MIN_EMIN; // [emax is plenty]
4027 // [dset.digits=16, etc., are sufficient]
4028
4029 // finally ready to roll
4030 for (;;) {
4031 #if DECCHECK
4032 iterations++;
4033 #endif
4034 // only the status from the accumulation is interesting
4035 // [but it should remain unchanged after first add]
4036 decAddOp(a, a, t, &aset, 0, status); // a=a+t
4037 decMultiplyOp(t, t, x, &tset, &ignore); // t=t*x
4038 decDivideOp(t, t, d, &tset, DIVIDE, &ignore); // t=t/d
4039 // the iteration ends when the term cannot affect the result,
4040 // if rounded to p digits, which is when its value is smaller
4041 // than the accumulator by p+1 digits. There must also be
4042 // full precision in a.
4043 if (((a->digits+a->exponent)>=(t->digits+t->exponent+p+1))
4044 && (a->digits>=p)) break;
4045 decAddOp(d, d, &numone, &dset, 0, &ignore); // d=d+1
4046 } // iterate
4047
4048 #if DECCHECK
4049 // just a sanity check; comment out test to show always
4050 if (iterations>p+3)
4051 printf("Exp iterations=%d, status=%08x, p=%d, d=%d\n",
4052 iterations, *status, p, x->digits);
4053 #endif
4054 } // h<=8
4055
4056 // apply postconditioning: a=a**(10**h) -- this is calculated
4057 // at a slightly higher precision than Hull & Abrham suggest
4058 if (h>0) {
4059 Int seenbit=0; // set once a 1-bit is seen
4060 Int i; // counter
4061 Int n=powers[h]; // always positive
4062 aset.digits=p+2; // sufficient precision
4063 // avoid the overhead and many extra digits of decNumberPower
4064 // as all that is needed is the short 'multipliers' loop; here
4065 // accumulate the answer into t
4066 decNumberZero(t); *t->lsu=1; // acc=1
4067 for (i=1;;i++){ // for each bit [top bit ignored]
4068 // abandon if have had overflow or terminal underflow
4069 if (*status & (DEC_Overflow|DEC_Underflow)) { // interesting?
4070 if (*status&DEC_Overflow || ISZERO(t)) break;}
4071 n=n<<1; // move next bit to testable position
4072 if (n<0) { // top bit is set
4073 seenbit=1; // OK, have a significant bit
4074 decMultiplyOp(t, t, a, &aset, status); // acc=acc*x
4075 }
4076 if (i==31) break; // that was the last bit
4077 if (!seenbit) continue; // no need to square 1
4078 decMultiplyOp(t, t, t, &aset, status); // acc=acc*acc [square]
4079 } /*i*/ // 32 bits
4080 // decNumberShow(t);
4081 a=t; // and carry on using t instead of a
4082 }
4083
4084 // Copy and round the result to res
4085 residue=1; // indicate dirt to right ..
4086 if (ISZERO(a)) residue=0; // .. unless underflowed to 0
4087 aset.digits=set->digits; // [use default rounding]
4088 decCopyFit(res, a, &aset, &residue, status); // copy & shorten
4089 decFinish(res, set, &residue, status); // cleanup/set flags
4090 } while(0); // end protected
4091
4092 if (allocrhs !=NULL) free(allocrhs); // drop any storage used
4093 if (allocbufa!=NULL) free(allocbufa); // ..
4094 if (allocbuft!=NULL) free(allocbuft); // ..
4095 // [status is handled by caller]
4096 return res;
4097 } // decExpOp
4098
4099 /* ------------------------------------------------------------------ */
4100 /* Initial-estimate natural logarithm table */
4101 /* */
4102 /* LNnn -- 90-entry 16-bit table for values from .10 through .99. */
4103 /* The result is a 4-digit encode of the coefficient (c=the */
4104 /* top 14 bits encoding 0-9999) and a 2-digit encode of the */
4105 /* exponent (e=the bottom 2 bits encoding 0-3) */
4106 /* */
4107 /* The resulting value is given by: */
4108 /* */
4109 /* v = -c * 10**(-e-3) */
4110 /* */
4111 /* where e and c are extracted from entry k = LNnn[x-10] */
4112 /* where x is truncated (NB) into the range 10 through 99, */
4113 /* and then c = k>>2 and e = k&3. */
4114 /* ------------------------------------------------------------------ */
4115 const uShort LNnn[90]={9016, 8652, 8316, 8008, 7724, 7456, 7208,
4116 6972, 6748, 6540, 6340, 6148, 5968, 5792, 5628, 5464, 5312,
4117 5164, 5020, 4884, 4748, 4620, 4496, 4376, 4256, 4144, 4032,
4118 39233, 38181, 37157, 36157, 35181, 34229, 33297, 32389, 31501, 30629,
4119 29777, 28945, 28129, 27329, 26545, 25777, 25021, 24281, 23553, 22837,
4120 22137, 21445, 20769, 20101, 19445, 18801, 18165, 17541, 16925, 16321,
4121 15721, 15133, 14553, 13985, 13421, 12865, 12317, 11777, 11241, 10717,
4122 10197, 9685, 9177, 8677, 8185, 7697, 7213, 6737, 6269, 5801,
4123 5341, 4889, 4437, 39930, 35534, 31186, 26886, 22630, 18418, 14254,
4124 10130, 6046, 20055};
4125
4126 /* ------------------------------------------------------------------ */
4127 /* decLnOp -- effect natural logarithm */
4128 /* */
4129 /* This computes C = ln(A) */
4130 /* */
4131 /* res is C, the result. C may be A */
4132 /* rhs is A */
4133 /* set is the context; note that rounding mode has no effect */
4134 /* */
4135 /* C must have space for set->digits digits. */
4136 /* */
4137 /* Notable cases: */
4138 /* A<0 -> Invalid */
4139 /* A=0 -> -Infinity (Exact) */
4140 /* A=+Infinity -> +Infinity (Exact) */
4141 /* A=1 exactly -> 0 (Exact) */
4142 /* */
4143 /* Restrictions (as for Exp): */
4144 /* */
4145 /* digits, emax, and -emin in the context must be less than */
4146 /* DEC_MAX_MATH+11 (1000010), and the rhs must be within these */
4147 /* bounds or a zero. This is an internal routine, so these */
4148 /* restrictions are contractual and not enforced. */
4149 /* */
4150 /* A finite result is rounded using DEC_ROUND_HALF_EVEN; it will */
4151 /* almost always be correctly rounded, but may be up to 1 ulp in */
4152 /* error in rare cases. */
4153 /* ------------------------------------------------------------------ */
4154 /* The result is calculated using Newton's method, with each */
4155 /* iteration calculating a' = a + x * exp(-a) - 1. See, for example, */
4156 /* Epperson 1989. */
4157 /* */
4158 /* The iteration ends when the adjustment x*exp(-a)-1 is tiny enough. */
4159 /* This has to be calculated at the sum of the precision of x and the */
4160 /* working precision. */
4161 /* */
4162 /* Implementation notes: */
4163 /* */
4164 /* 1. This is separated out as decLnOp so it can be called from */
4165 /* other Mathematical functions (e.g., Log 10) with a wider range */
4166 /* than normal. In particular, it can handle the slightly wider */
4167 /* (+9+2) range needed by a power function. */
4168 /* */
4169 /* 2. The speed of this function is about 10x slower than exp, as */
4170 /* it typically needs 4-6 iterations for short numbers, and the */
4171 /* extra precision needed adds a squaring effect, twice. */
4172 /* */
4173 /* 3. Fastpaths are included for ln(10) and ln(2), up to length 40, */
4174 /* as these are common requests. ln(10) is used by log10(x). */
4175 /* */
4176 /* 4. An iteration might be saved by widening the LNnn table, and */
4177 /* would certainly save at least one if it were made ten times */
4178 /* bigger, too (for truncated fractions 0.100 through 0.999). */
4179 /* However, for most practical evaluations, at least four or five */
4180 /* iterations will be neede -- so this would only speed up by */
4181 /* 20-25% and that probably does not justify increasing the table */
4182 /* size. */
4183 /* ------------------------------------------------------------------ */
decLnOp(decNumber * res,const decNumber * rhs,decContext * set,uInt * status)4184 decNumber * decLnOp(decNumber *res, const decNumber *rhs,
4185 decContext *set, uInt *status) {
4186 uInt ignore=0; // working status accumulator
4187 uInt needbytes; // for space calculations
4188 Int residue; // rounding residue
4189 Int r; // rhs=f*10**r [see below]
4190 Int p; // working precision
4191 Int pp; // precision for iteration
4192 Int t; // work
4193
4194 // buffers for a (accumulator, typically precision+2) and b
4195 // (adjustment calculator, same size)
4196 decNumber bufa[D2N(DECBUFFER+2)];
4197 decNumber *allocbufa=NULL; // -> allocated bufa, iff allocated
4198 decNumber *a=bufa; // accumulator/work
4199 decNumber bufb[D2N(DECBUFFER+2)];
4200 decNumber *allocbufb=NULL; // -> allocated bufa, iff allocated
4201 decNumber *b=bufb; // adjustment/work
4202
4203 decNumber numone; // constant 1
4204 decNumber cmp; // work
4205 decContext aset, bset; // working contexts
4206
4207 #if DECCHECK
4208 Int iterations=0; // for later sanity check
4209 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
4210 #endif
4211
4212 do { // protect allocated storage
4213 if (SPECIALARG) { // handle infinities and NaNs
4214 if (decNumberIsInfinite(rhs)) { // an infinity
4215 if (decNumberIsNegative(rhs)) // -Infinity -> error
4216 *status|=DEC_Invalid_operation;
4217 else decNumberCopy(res, rhs); // +Infinity -> self
4218 }
4219 else decNaNs(res, rhs, NULL, status); // a NaN
4220 break;}
4221
4222 if (ISZERO(rhs)) { // +/- zeros -> -Infinity
4223 decNumberZero(res); // make clean
4224 res->bits=DECINF|DECNEG; // set - infinity
4225 break;} // [no status to set]
4226
4227 // Non-zero negatives are bad...
4228 if (decNumberIsNegative(rhs)) { // -x -> error
4229 *status|=DEC_Invalid_operation;
4230 break;}
4231
4232 // Here, rhs is positive, finite, and in range
4233
4234 // lookaside fastpath code for ln(2) and ln(10) at common lengths
4235 if (rhs->exponent==0 && set->digits<=40) {
4236 #if DECDPUN==1
4237 if (rhs->lsu[0]==0 && rhs->lsu[1]==1 && rhs->digits==2) { // ln(10)
4238 #else
4239 if (rhs->lsu[0]==10 && rhs->digits==2) { // ln(10)
4240 #endif
4241 aset=*set; aset.round=DEC_ROUND_HALF_EVEN;
4242 #define LN10 "2.302585092994045684017991454684364207601"
4243 decNumberFromString(res, LN10, &aset);
4244 *status|=(DEC_Inexact | DEC_Rounded); // is inexact
4245 break;}
4246 if (rhs->lsu[0]==2 && rhs->digits==1) { // ln(2)
4247 aset=*set; aset.round=DEC_ROUND_HALF_EVEN;
4248 #define LN2 "0.6931471805599453094172321214581765680755"
4249 decNumberFromString(res, LN2, &aset);
4250 *status|=(DEC_Inexact | DEC_Rounded);
4251 break;}
4252 } // integer and short
4253
4254 // Determine the working precision. This is normally the
4255 // requested precision + 2, with a minimum of 9. However, if
4256 // the rhs is 'over-precise' then allow for all its digits to
4257 // potentially participate (consider an rhs where all the excess
4258 // digits are 9s) so in this case use rhs->digits+2.
4259 p=MAX(rhs->digits, MAX(set->digits, 7))+2;
4260
4261 // Allocate space for the accumulator and the high-precision
4262 // adjustment calculator, if necessary. The accumulator must
4263 // be able to hold p digits, and the adjustment up to
4264 // rhs->digits+p digits. They are also made big enough for 16
4265 // digits so that they can be used for calculating the initial
4266 // estimate.
4267 needbytes=sizeof(decNumber)+(D2U(MAX(p,16))-1)*sizeof(Unit);
4268 if (needbytes>sizeof(bufa)) { // need malloc space
4269 allocbufa=(decNumber *)malloc(needbytes);
4270 if (allocbufa==NULL) { // hopeless -- abandon
4271 *status|=DEC_Insufficient_storage;
4272 break;}
4273 a=allocbufa; // use the allocated space
4274 }
4275 pp=p+rhs->digits;
4276 needbytes=sizeof(decNumber)+(D2U(MAX(pp,16))-1)*sizeof(Unit);
4277 if (needbytes>sizeof(bufb)) { // need malloc space
4278 allocbufb=(decNumber *)malloc(needbytes);
4279 if (allocbufb==NULL) { // hopeless -- abandon
4280 *status|=DEC_Insufficient_storage;
4281 break;}
4282 b=allocbufb; // use the allocated space
4283 }
4284
4285 // Prepare an initial estimate in acc. Calculate this by
4286 // considering the coefficient of x to be a normalized fraction,
4287 // f, with the decimal point at far left and multiplied by
4288 // 10**r. Then, rhs=f*10**r and 0.1<=f<1, and
4289 // ln(x) = ln(f) + ln(10)*r
4290 // Get the initial estimate for ln(f) from a small lookup
4291 // table (see above) indexed by the first two digits of f,
4292 // truncated.
4293
4294 decContextDefault(&aset, DEC_INIT_DECIMAL64); // 16-digit extended
4295 r=rhs->exponent+rhs->digits; // 'normalised' exponent
4296 decPutInt(a, r); // a=r
4297 decPutInt(b, 2302585); // b=ln(10) (2.302585)
4298 b->exponent=-6; // ..
4299 decMultiplyOp(a, a, b, &aset, &ignore); // a=a*b
4300 // now get top two digits of rhs into b by simple truncate and
4301 // force to integer
4302 residue=0; // (no residue)
4303 aset.digits=2; aset.round=DEC_ROUND_DOWN;
4304 decCopyFit(b, rhs, &aset, &residue, &ignore); // copy & shorten
4305 b->exponent=0; // make integer
4306 t=decGetInt(b); // [cannot fail]
4307 if (t<10) t=X10(t); // adjust single-digit b
4308 t=LNnn[t-10]; // look up ln(b)
4309 decPutInt(b, t>>2); // b=ln(b) coefficient
4310 b->exponent=-(t&3)-3; // set exponent
4311 b->bits=DECNEG; // ln(0.10)->ln(0.99) always -ve
4312 aset.digits=16; aset.round=DEC_ROUND_HALF_EVEN; // restore
4313 decAddOp(a, a, b, &aset, 0, &ignore); // acc=a+b
4314 // the initial estimate is now in a, with up to 4 digits correct.
4315 // When rhs is at or near Nmax the estimate will be low, so we
4316 // will approach it from below, avoiding overflow when calling exp.
4317
4318 decNumberZero(&numone); *numone.lsu=1; // constant 1 for adjustment
4319
4320 // accumulator bounds are as requested (could underflow, but
4321 // cannot overflow)
4322 aset.emax=set->emax;
4323 aset.emin=set->emin;
4324 aset.clamp=0; // no concrete format
4325 // set up a context to be used for the multiply and subtract
4326 bset=aset;
4327 bset.emax=DEC_MAX_MATH*2; // use double bounds for the
4328 bset.emin=-DEC_MAX_MATH*2; // adjustment calculation
4329 // [see decExpOp call below]
4330 // for each iteration double the number of digits to calculate,
4331 // up to a maximum of p
4332 pp=9; // initial precision
4333 // [initially 9 as then the sequence starts 7+2, 16+2, and
4334 // 34+2, which is ideal for standard-sized numbers]
4335 aset.digits=pp; // working context
4336 bset.digits=pp+rhs->digits; // wider context
4337 for (;;) { // iterate
4338 #if DECCHECK
4339 iterations++;
4340 if (iterations>24) break; // consider 9 * 2**24
4341 #endif
4342 // calculate the adjustment (exp(-a)*x-1) into b. This is a
4343 // catastrophic subtraction but it really is the difference
4344 // from 1 that is of interest.
4345 // Use the internal entry point to Exp as it allows the double
4346 // range for calculating exp(-a) when a is the tiniest subnormal.
4347 a->bits^=DECNEG; // make -a
4348 decExpOp(b, a, &bset, &ignore); // b=exp(-a)
4349 a->bits^=DECNEG; // restore sign of a
4350 // now multiply by rhs and subtract 1, at the wider precision
4351 decMultiplyOp(b, b, rhs, &bset, &ignore); // b=b*rhs
4352 decAddOp(b, b, &numone, &bset, DECNEG, &ignore); // b=b-1
4353
4354 // the iteration ends when the adjustment cannot affect the
4355 // result by >=0.5 ulp (at the requested digits), which
4356 // is when its value is smaller than the accumulator by
4357 // set->digits+1 digits (or it is zero) -- this is a looser
4358 // requirement than for Exp because all that happens to the
4359 // accumulator after this is the final rounding (but note that
4360 // there must also be full precision in a, or a=0).
4361
4362 if (decNumberIsZero(b) ||
4363 (a->digits+a->exponent)>=(b->digits+b->exponent+set->digits+1)) {
4364 if (a->digits==p) break;
4365 if (decNumberIsZero(a)) {
4366 decCompareOp(&cmp, rhs, &numone, &aset, COMPARE, &ignore); // rhs=1 ?
4367 if (cmp.lsu[0]==0) a->exponent=0; // yes, exact 0
4368 else *status|=(DEC_Inexact | DEC_Rounded); // no, inexact
4369 break;
4370 }
4371 // force padding if adjustment has gone to 0 before full length
4372 if (decNumberIsZero(b)) b->exponent=a->exponent-p;
4373 }
4374
4375 // not done yet ...
4376 decAddOp(a, a, b, &aset, 0, &ignore); // a=a+b for next estimate
4377 if (pp==p) continue; // precision is at maximum
4378 // lengthen the next calculation
4379 pp=pp*2; // double precision
4380 if (pp>p) pp=p; // clamp to maximum
4381 aset.digits=pp; // working context
4382 bset.digits=pp+rhs->digits; // wider context
4383 } // Newton's iteration
4384
4385 #if DECCHECK
4386 // just a sanity check; remove the test to show always
4387 if (iterations>24)
4388 printf("Ln iterations=%d, status=%08x, p=%d, d=%d\n",
4389 iterations, *status, p, rhs->digits);
4390 #endif
4391
4392 // Copy and round the result to res
4393 residue=1; // indicate dirt to right
4394 if (ISZERO(a)) residue=0; // .. unless underflowed to 0
4395 aset.digits=set->digits; // [use default rounding]
4396 decCopyFit(res, a, &aset, &residue, status); // copy & shorten
4397 decFinish(res, set, &residue, status); // cleanup/set flags
4398 } while(0); // end protected
4399
4400 if (allocbufa!=NULL) free(allocbufa); // drop any storage used
4401 if (allocbufb!=NULL) free(allocbufb); // ..
4402 // [status is handled by caller]
4403 return res;
4404 } // decLnOp
4405
4406 /* ------------------------------------------------------------------ */
4407 /* decQuantizeOp -- force exponent to requested value */
4408 /* */
4409 /* This computes C = op(A, B), where op adjusts the coefficient */
4410 /* of C (by rounding or shifting) such that the exponent (-scale) */
4411 /* of C has the value B or matches the exponent of B. */
4412 /* The numerical value of C will equal A, except for the effects of */
4413 /* any rounding that occurred. */
4414 /* */
4415 /* res is C, the result. C may be A or B */
4416 /* lhs is A, the number to adjust */
4417 /* rhs is B, the requested exponent */
4418 /* set is the context */
4419 /* quant is 1 for quantize or 0 for rescale */
4420 /* status is the status accumulator (this can be called without */
4421 /* risk of control loss) */
4422 /* */
4423 /* C must have space for set->digits digits. */
4424 /* */
4425 /* Unless there is an error or the result is infinite, the exponent */
4426 /* after the operation is guaranteed to be that requested. */
4427 /* ------------------------------------------------------------------ */
4428 static decNumber * decQuantizeOp(decNumber *res, const decNumber *lhs,
4429 const decNumber *rhs, decContext *set,
4430 Flag quant, uInt *status) {
4431 #if DECSUBSET
4432 decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated
4433 decNumber *allocrhs=NULL; // .., rhs
4434 #endif
4435 const decNumber *inrhs=rhs; // save original rhs
4436 Int reqdigits=set->digits; // requested DIGITS
4437 Int reqexp; // requested exponent [-scale]
4438 Int residue=0; // rounding residue
4439 Int etiny=set->emin-(reqdigits-1);
4440
4441 #if DECCHECK
4442 if (decCheckOperands(res, lhs, rhs, set)) return res;
4443 #endif
4444
4445 do { // protect allocated storage
4446 #if DECSUBSET
4447 if (!set->extended) {
4448 // reduce operands and set lostDigits status, as needed
4449 if (lhs->digits>reqdigits) {
4450 alloclhs=decRoundOperand(lhs, set, status);
4451 if (alloclhs==NULL) break;
4452 lhs=alloclhs;
4453 }
4454 if (rhs->digits>reqdigits) { // [this only checks lostDigits]
4455 allocrhs=decRoundOperand(rhs, set, status);
4456 if (allocrhs==NULL) break;
4457 rhs=allocrhs;
4458 }
4459 }
4460 #endif
4461 // [following code does not require input rounding]
4462
4463 // Handle special values
4464 if (SPECIALARGS) {
4465 // NaNs get usual processing
4466 if (SPECIALARGS & (DECSNAN | DECNAN))
4467 decNaNs(res, lhs, rhs, status);
4468 // one infinity but not both is bad
4469 else if ((lhs->bits ^ rhs->bits) & DECINF)
4470 *status|=DEC_Invalid_operation;
4471 // both infinity: return lhs
4472 else decNumberCopy(res, lhs); // [nop if in place]
4473 break;
4474 }
4475
4476 // set requested exponent
4477 if (quant) reqexp=inrhs->exponent; // quantize -- match exponents
4478 else { // rescale -- use value of rhs
4479 // Original rhs must be an integer that fits and is in range,
4480 // which could be from -1999999997 to +999999999, thanks to
4481 // subnormals
4482 reqexp=decGetInt(inrhs); // [cannot fail]
4483 }
4484
4485 #if DECSUBSET
4486 if (!set->extended) etiny=set->emin; // no subnormals
4487 #endif
4488
4489 if (reqexp==BADINT // bad (rescale only) or ..
4490 || reqexp==BIGODD || reqexp==BIGEVEN // very big (ditto) or ..
4491 || (reqexp<etiny) // < lowest
4492 || (reqexp>set->emax)) { // > emax
4493 *status|=DEC_Invalid_operation;
4494 break;}
4495
4496 // the RHS has been processed, so it can be overwritten now if necessary
4497 if (ISZERO(lhs)) { // zero coefficient unchanged
4498 decNumberCopy(res, lhs); // [nop if in place]
4499 res->exponent=reqexp; // .. just set exponent
4500 #if DECSUBSET
4501 if (!set->extended) res->bits=0; // subset specification; no -0
4502 #endif
4503 }
4504 else { // non-zero lhs
4505 Int adjust=reqexp-lhs->exponent; // digit adjustment needed
4506 // if adjusted coefficient will definitely not fit, give up now
4507 if ((lhs->digits-adjust)>reqdigits) {
4508 *status|=DEC_Invalid_operation;
4509 break;
4510 }
4511
4512 if (adjust>0) { // increasing exponent
4513 // this will decrease the length of the coefficient by adjust
4514 // digits, and must round as it does so
4515 decContext workset; // work
4516 workset=*set; // clone rounding, etc.
4517 workset.digits=lhs->digits-adjust; // set requested length
4518 // [note that the latter can be <1, here]
4519 decCopyFit(res, lhs, &workset, &residue, status); // fit to result
4520 decApplyRound(res, &workset, residue, status); // .. and round
4521 residue=0; // [used]
4522 // If just rounded a 999s case, exponent will be off by one;
4523 // adjust back (after checking space), if so.
4524 if (res->exponent>reqexp) {
4525 // re-check needed, e.g., for quantize(0.9999, 0.001) under
4526 // set->digits==3
4527 if (res->digits==reqdigits) { // cannot shift by 1
4528 *status&=~(DEC_Inexact | DEC_Rounded); // [clean these]
4529 *status|=DEC_Invalid_operation;
4530 break;
4531 }
4532 res->digits=decShiftToMost(res->lsu, res->digits, 1); // shift
4533 res->exponent--; // (re)adjust the exponent.
4534 }
4535 #if DECSUBSET
4536 if (ISZERO(res) && !set->extended) res->bits=0; // subset; no -0
4537 #endif
4538 } // increase
4539 else /* adjust<=0 */ { // decreasing or = exponent
4540 // this will increase the length of the coefficient by -adjust
4541 // digits, by adding zero or more trailing zeros; this is
4542 // already checked for fit, above
4543 decNumberCopy(res, lhs); // [it will fit]
4544 // if padding needed (adjust<0), add it now...
4545 if (adjust<0) {
4546 res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
4547 res->exponent+=adjust; // adjust the exponent
4548 }
4549 } // decrease
4550 } // non-zero
4551
4552 // Check for overflow [do not use Finalize in this case, as an
4553 // overflow here is a "don't fit" situation]
4554 if (res->exponent>set->emax-res->digits+1) { // too big
4555 *status|=DEC_Invalid_operation;
4556 break;
4557 }
4558 else {
4559 decFinalize(res, set, &residue, status); // set subnormal flags
4560 *status&=~DEC_Underflow; // suppress Underflow [754r]
4561 }
4562 } while(0); // end protected
4563
4564 #if DECSUBSET
4565 if (allocrhs!=NULL) free(allocrhs); // drop any storage used
4566 if (alloclhs!=NULL) free(alloclhs); // ..
4567 #endif
4568 return res;
4569 } // decQuantizeOp
4570
4571 /* ------------------------------------------------------------------ */
4572 /* decCompareOp -- compare, min, or max two Numbers */
4573 /* */
4574 /* This computes C = A ? B and carries out one of four operations: */
4575 /* COMPARE -- returns the signum (as a number) giving the */
4576 /* result of a comparison unless one or both */
4577 /* operands is a NaN (in which case a NaN results) */
4578 /* COMPMAX -- returns the larger of the operands, using the */
4579 /* 754r maxnum operation */
4580 /* COMPMIN -- the 754r minnum operation */
4581 /* COMTOTAL -- returns the signum (as a number) giving the */
4582 /* result of a comparison using 754r total ordering */
4583 /* */
4584 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
4585 /* lhs is A */
4586 /* rhs is B */
4587 /* set is the context */
4588 /* op is the operation flag */
4589 /* status is the usual accumulator */
4590 /* */
4591 /* C must have space for one digit for COMPARE or set->digits for */
4592 /* COMPMAX and COMPMIN. */
4593 /* ------------------------------------------------------------------ */
4594 /* The emphasis here is on speed for common cases, and avoiding */
4595 /* coefficient comparison if possible. */
4596 /* ------------------------------------------------------------------ */
4597 decNumber * decCompareOp(decNumber *res, const decNumber *lhs,
4598 const decNumber *rhs, decContext *set,
4599 Flag op, uInt *status) {
4600 #if DECSUBSET
4601 decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated
4602 decNumber *allocrhs=NULL; // .., rhs
4603 #endif
4604 Int result=0; // default result value
4605 uByte merged; // work
4606
4607 #if DECCHECK
4608 if (decCheckOperands(res, lhs, rhs, set)) return res;
4609 #endif
4610
4611 do { // protect allocated storage
4612 #if DECSUBSET
4613 if (!set->extended) {
4614 // reduce operands and set lostDigits status, as needed
4615 if (lhs->digits>set->digits) {
4616 alloclhs=decRoundOperand(lhs, set, status);
4617 if (alloclhs==NULL) {result=BADINT; break;}
4618 lhs=alloclhs;
4619 }
4620 if (rhs->digits>set->digits) {
4621 allocrhs=decRoundOperand(rhs, set, status);
4622 if (allocrhs==NULL) {result=BADINT; break;}
4623 rhs=allocrhs;
4624 }
4625 }
4626 #endif
4627 // [following code does not require input rounding]
4628
4629 // If total ordering then handle differing signs 'up front'
4630 if (op == COMPTOTAL) { // total ordering
4631 if (decNumberIsNegative(lhs) & !decNumberIsNegative(rhs)) {
4632 result=-1;
4633 break;
4634 }
4635 if (!decNumberIsNegative(lhs) & decNumberIsNegative(rhs)) {
4636 result=+1;
4637 break;
4638 }
4639 }
4640
4641 // handle NaNs specially; let infinities drop through
4642 // This assumes sNaN (even just one) leads to NaN.
4643 merged=(lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
4644 if (merged) { // a NaN bit set
4645 if (op == COMPARE); // result will be NaN
4646 else if (op == COMPTOTAL) { // total ordering, always finite
4647 // signs are known to be the same; compute the ordering here
4648 // as if the signs are both positive, then invert for negatives
4649 if (!decNumberIsNaN(lhs)) result=-1;
4650 else if (!decNumberIsNaN(rhs)) result=+1;
4651 // here if both NaNs
4652 else if (decNumberIsSNaN(lhs) && decNumberIsQNaN(rhs)) result=-1;
4653 else if (decNumberIsQNaN(lhs) && decNumberIsSNaN(rhs)) result=+1;
4654 else { // both NaN or both sNaN
4655 // now it just depends on the payload
4656 result=decUnitCompare(lhs->lsu, D2U(lhs->digits),
4657 rhs->lsu, D2U(rhs->digits), 0);
4658 // [Error not possible, as these are 'aligned']
4659 } // both same NaNs
4660 if (decNumberIsNegative(lhs)) result=-result;
4661 break;
4662 } // total order
4663
4664 else if (merged & DECSNAN); // sNaN -> qNaN
4665 else { // here if MIN or MAX and one or two quiet NaNs
4666 // min or max -- 754r rules ignore single NaN
4667 if (!decNumberIsNaN(lhs) || !decNumberIsNaN(rhs)) {
4668 // just one NaN; force choice to be the non-NaN operand
4669 op=COMPMAX;
4670 if (lhs->bits & DECNAN) result=-1; // pick rhs
4671 else result=+1; // pick lhs
4672 break;
4673 }
4674 } // max or min
4675 op = COMPNAN; // use special path
4676 decNaNs(res, lhs, rhs, status);
4677 break;
4678 }
4679 result=decCompare(lhs, rhs, 0); // have numbers
4680 } while(0); // end protected
4681
4682 if (result==BADINT) *status|=DEC_Insufficient_storage; // rare
4683 else {
4684 if (op == COMPARE || op == COMPTOTAL) { // returning signum
4685 if (op == COMPTOTAL && result==0) {
4686 // operands are numerically equal or same NaN (and same sign,
4687 // tested first); if identical, leave result 0
4688 if (lhs->exponent!=rhs->exponent) {
4689 if (lhs->exponent<rhs->exponent) result=-1;
4690 else result=+1;
4691 if (decNumberIsNegative(lhs)) result=-result;
4692 } // lexp!=rexp
4693 } // total-order by exponent
4694 decNumberZero(res); // [always a valid result]
4695 if (result!=0) { // must be -1 or +1
4696 *res->lsu=1;
4697 if (result<0) res->bits=DECNEG;
4698 }
4699 }
4700 else if (op == COMPNAN); // special, drop through
4701 else { // MAX or MIN, non-NaN result
4702 Int residue=0; // rounding accumulator
4703 // choose the operand for the result
4704 const decNumber *choice;
4705 if (result==0) { // operands are numerically equal
4706 // choose according to sign then exponent (see 754r)
4707 uByte slhs=(lhs->bits & DECNEG);
4708 uByte srhs=(rhs->bits & DECNEG);
4709 #if DECSUBSET
4710 if (!set->extended) { // subset: force left-hand
4711 op=COMPMAX;
4712 result=+1;
4713 }
4714 else
4715 #endif
4716 if (slhs!=srhs) { // signs differ
4717 if (slhs) result=-1; // rhs is max
4718 else result=+1; // lhs is max
4719 }
4720 else if (slhs && srhs) { // both negative
4721 if (lhs->exponent<rhs->exponent) result=+1;
4722 else result=-1;
4723 // [if equal, use lhs, technically identical]
4724 }
4725 else { // both positive
4726 if (lhs->exponent>rhs->exponent) result=+1;
4727 else result=-1;
4728 // [ditto]
4729 }
4730 } // numerically equal
4731 // here result will be non-0
4732 if (op == COMPMIN) result=-result;// reverse if looking for MIN
4733 choice=(result>0 ? lhs : rhs); // choose
4734 // copy chosen to result, rounding if need be
4735 decCopyFit(res, choice, set, &residue, status);
4736 decFinish(res, set, &residue, status);
4737 }
4738 }
4739 #if DECSUBSET
4740 if (allocrhs!=NULL) free(allocrhs); // free any storage used
4741 if (alloclhs!=NULL) free(alloclhs); // ..
4742 #endif
4743 return res;
4744 } // decCompareOp
4745
4746 /* ------------------------------------------------------------------ */
4747 /* decCompare -- compare two decNumbers by numerical value */
4748 /* */
4749 /* This routine compares A ? B without altering them. */
4750 /* */
4751 /* Arg1 is A, a decNumber which is not a NaN */
4752 /* Arg2 is B, a decNumber which is not a NaN */
4753 /* Arg3 is 1 for a sign-independent compare, 0 otherwise */
4754 /* */
4755 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
4756 /* (the only possible failure is an allocation error) */
4757 /* ------------------------------------------------------------------ */
4758 static Int decCompare(const decNumber *lhs, const decNumber *rhs,
4759 Flag abs) {
4760 Int result; // result value
4761 Int sigr; // rhs signum
4762 Int compare; // work
4763
4764 result=1; // assume signum(lhs)
4765 if (ISZERO(lhs)) result=0;
4766 if (abs) {
4767 if (!ISZERO(rhs)) result=1; // not both 0
4768 }
4769 else { // signs matter
4770 if (result && decNumberIsNegative(lhs)) result=-1;
4771 sigr=1; // compute signum(rhs)
4772 if (ISZERO(rhs)) sigr=0;
4773 else if (decNumberIsNegative(rhs)) sigr=-1;
4774 if (result > sigr) return +1; // L > R, return 1
4775 if (result < sigr) return -1; // R < L, return -1
4776 }
4777
4778 // signums are the same
4779 if (result==0) return 0; // both 0
4780 // Both non-zero
4781 if ((lhs->bits | rhs->bits) & DECINF) { // one or more infinities
4782 if (decNumberIsInfinite(rhs)) {
4783 if (decNumberIsInfinite(lhs)) result=0;// both infinite
4784 else result=-result; // only rhs infinite
4785 }
4786 return result;
4787 }
4788
4789 // must compare the coefficients, allowing for exponents
4790 if (lhs->exponent>rhs->exponent) { // LHS exponent larger
4791 // swap sides, and sign
4792 const decNumber *temp=lhs;
4793 lhs=rhs;
4794 rhs=temp;
4795 result=-result;
4796 }
4797
4798 compare=decUnitCompare(lhs->lsu, D2U(lhs->digits),
4799 rhs->lsu, D2U(rhs->digits),
4800 rhs->exponent-lhs->exponent);
4801
4802 if (compare!=BADINT) compare*=result; // comparison succeeded
4803 return compare;
4804 } // decCompare
4805
4806 /* ------------------------------------------------------------------ */
4807 /* decUnitCompare -- compare two >=0 integers in Unit arrays */
4808 /* */
4809 /* This routine compares A ? B*10**E where A and B are unit arrays */
4810 /* A is a plain integer */
4811 /* B has an exponent of E (which must be non-negative) */
4812 /* */
4813 /* Arg1 is A first Unit (lsu) */
4814 /* Arg2 is A length in Units */
4815 /* Arg3 is B first Unit (lsu) */
4816 /* Arg4 is B length in Units */
4817 /* Arg5 is E (0 if the units are aligned) */
4818 /* */
4819 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
4820 /* (the only possible failure is an allocation error, which can */
4821 /* only occur if E!=0) */
4822 /* ------------------------------------------------------------------ */
4823 static Int decUnitCompare(const Unit *a, Int alength,
4824 const Unit *b, Int blength, Int exp) {
4825 Unit *acc; // accumulator for result
4826 Unit accbuff[SD2U(DECBUFFER+1)];// local buffer
4827 Unit *allocacc=NULL; // -> allocated acc buffer, iff allocated
4828 Int accunits, need; // units in use or needed for acc
4829 const Unit *l, *r, *u; // work
4830 Int expunits, exprem, result; // ..
4831
4832 if (exp==0) { // aligned; fastpath
4833 if (alength>blength) return 1;
4834 if (alength<blength) return -1;
4835 // same number of units in both -- need unit-by-unit compare
4836 l=a+alength-1;
4837 r=b+alength-1;
4838 for (;l>=a; l--, r--) {
4839 if (*l>*r) return 1;
4840 if (*l<*r) return -1;
4841 }
4842 return 0; // all units match
4843 } // aligned
4844
4845 // Unaligned. If one is >1 unit longer than the other, padded
4846 // approximately, then can return easily
4847 if (alength>blength+(Int)D2U(exp)) return 1;
4848 if (alength+1<blength+(Int)D2U(exp)) return -1;
4849
4850 // Need to do a real subtract. For this, a result buffer is needed
4851 // even though only the sign is of interest. Its length needs
4852 // to be the larger of alength and padded blength, +2
4853 need=blength+D2U(exp); // maximum real length of B
4854 if (need<alength) need=alength;
4855 need+=2;
4856 acc=accbuff; // assume use local buffer
4857 if (need*sizeof(Unit)>sizeof(accbuff)) {
4858 allocacc=(Unit *)malloc(need*sizeof(Unit));
4859 if (allocacc==NULL) return BADINT; // hopeless -- abandon
4860 acc=allocacc;
4861 }
4862 // Calculate units and remainder from exponent.
4863 expunits=exp/DECDPUN;
4864 exprem=exp%DECDPUN;
4865 // subtract [A+B*(-m)]
4866 accunits=decUnitAddSub(a, alength, b, blength, expunits, acc,
4867 -(Int)powers[exprem]);
4868 // [UnitAddSub result may have leading zeros, even on zero]
4869 if (accunits<0) result=-1; // negative result
4870 else { // non-negative result
4871 // check units of the result before freeing any storage
4872 for (u=acc; u<acc+accunits-1 && *u==0;) u++;
4873 result=(*u==0 ? 0 : +1);
4874 }
4875 // clean up and return the result
4876 if (allocacc!=NULL) free(allocacc); // drop any storage used
4877 return result;
4878 } // decUnitCompare
4879
4880 /* ------------------------------------------------------------------ */
4881 /* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays */
4882 /* */
4883 /* This routine performs the calculation: */
4884 /* */
4885 /* C=A+(B*M) */
4886 /* */
4887 /* Where M is in the range -DECDPUNMAX through +DECDPUNMAX. */
4888 /* */
4889 /* A may be shorter or longer than B. */
4890 /* */
4891 /* Leading zeros are not removed after a calculation. The result is */
4892 /* either the same length as the longer of A and B (adding any */
4893 /* shift), or one Unit longer than that (if a Unit carry occurred). */
4894 /* */
4895 /* A and B content are not altered unless C is also A or B. */
4896 /* C may be the same array as A or B, but only if no zero padding is */
4897 /* requested (that is, C may be B only if bshift==0). */
4898 /* C is filled from the lsu; only those units necessary to complete */
4899 /* the calculation are referenced. */
4900 /* */
4901 /* Arg1 is A first Unit (lsu) */
4902 /* Arg2 is A length in Units */
4903 /* Arg3 is B first Unit (lsu) */
4904 /* Arg4 is B length in Units */
4905 /* Arg5 is B shift in Units (>=0; pads with 0 units if positive) */
4906 /* Arg6 is C first Unit (lsu) */
4907 /* Arg7 is M, the multiplier */
4908 /* */
4909 /* returns the count of Units written to C, which will be non-zero */
4910 /* and negated if the result is negative. That is, the sign of the */
4911 /* returned Int is the sign of the result (positive for zero) and */
4912 /* the absolute value of the Int is the count of Units. */
4913 /* */
4914 /* It is the caller's responsibility to make sure that C size is */
4915 /* safe, allowing space if necessary for a one-Unit carry. */
4916 /* */
4917 /* This routine is severely performance-critical; *any* change here */
4918 /* must be measured (timed) to assure no performance degradation. */
4919 /* In particular, trickery here tends to be counter-productive, as */
4920 /* increased complexity of code hurts register optimizations on */
4921 /* register-poor architectures. Avoiding divisions is nearly */
4922 /* always a Good Idea, however. */
4923 /* */
4924 /* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark */
4925 /* (IBM Warwick, UK) for some of the ideas used in this routine. */
4926 /* ------------------------------------------------------------------ */
4927 static Int decUnitAddSub(const Unit *a, Int alength,
4928 const Unit *b, Int blength, Int bshift,
4929 Unit *c, Int m) {
4930 const Unit *alsu=a; // A lsu [need to remember it]
4931 Unit *clsu=c; // C ditto
4932 Unit *minC; // low water mark for C
4933 Unit *maxC; // high water mark for C
4934 eInt carry=0; // carry integer (could be Long)
4935 Int add; // work
4936 #if DECDPUN<=4 // myriadal, millenary, etc.
4937 Int est; // estimated quotient
4938 #endif
4939
4940 #if DECTRACE
4941 if (alength<1 || blength<1)
4942 printf("decUnitAddSub: alen blen m %d %d [%d]\n", alength, blength, m);
4943 #endif
4944
4945 maxC=c+alength; // A is usually the longer
4946 minC=c+blength; // .. and B the shorter
4947 if (bshift!=0) { // B is shifted; low As copy across
4948 minC+=bshift;
4949 // if in place [common], skip copy unless there's a gap [rare]
4950 if (a==c && bshift<=alength) {
4951 c+=bshift;
4952 a+=bshift;
4953 }
4954 else for (; c<clsu+bshift; a++, c++) { // copy needed
4955 if (a<alsu+alength) *c=*a;
4956 else *c=0;
4957 }
4958 }
4959 if (minC>maxC) { // swap
4960 Unit *hold=minC;
4961 minC=maxC;
4962 maxC=hold;
4963 }
4964
4965 // For speed, do the addition as two loops; the first where both A
4966 // and B contribute, and the second (if necessary) where only one or
4967 // other of the numbers contribute.
4968 // Carry handling is the same (i.e., duplicated) in each case.
4969 for (; c<minC; c++) {
4970 carry+=*a;
4971 a++;
4972 carry+=((eInt)*b)*m; // [special-casing m=1/-1
4973 b++; // here is not a win]
4974 // here carry is new Unit of digits; it could be +ve or -ve
4975 if ((ueInt)carry<=DECDPUNMAX) { // fastpath 0-DECDPUNMAX
4976 *c=(Unit)carry;
4977 carry=0;
4978 continue;
4979 }
4980 #if DECDPUN==4 // use divide-by-multiply
4981 if (carry>=0) {
4982 est=(((ueInt)carry>>11)*53687)>>18;
4983 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
4984 carry=est; // likely quotient [89%]
4985 if (*c<DECDPUNMAX+1) continue; // estimate was correct
4986 carry++;
4987 *c-=DECDPUNMAX+1;
4988 continue;
4989 }
4990 // negative case
4991 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
4992 est=(((ueInt)carry>>11)*53687)>>18;
4993 *c=(Unit)(carry-est*(DECDPUNMAX+1));
4994 carry=est-(DECDPUNMAX+1); // correctly negative
4995 if (*c<DECDPUNMAX+1) continue; // was OK
4996 carry++;
4997 *c-=DECDPUNMAX+1;
4998 #elif DECDPUN==3
4999 if (carry>=0) {
5000 est=(((ueInt)carry>>3)*16777)>>21;
5001 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
5002 carry=est; // likely quotient [99%]
5003 if (*c<DECDPUNMAX+1) continue; // estimate was correct
5004 carry++;
5005 *c-=DECDPUNMAX+1;
5006 continue;
5007 }
5008 // negative case
5009 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5010 est=(((ueInt)carry>>3)*16777)>>21;
5011 *c=(Unit)(carry-est*(DECDPUNMAX+1));
5012 carry=est-(DECDPUNMAX+1); // correctly negative
5013 if (*c<DECDPUNMAX+1) continue; // was OK
5014 carry++;
5015 *c-=DECDPUNMAX+1;
5016 #elif DECDPUN<=2
5017 // Can use QUOT10 as carry <= 4 digits
5018 if (carry>=0) {
5019 est=QUOT10(carry, DECDPUN);
5020 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
5021 carry=est; // quotient
5022 continue;
5023 }
5024 // negative case
5025 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5026 est=QUOT10(carry, DECDPUN);
5027 *c=(Unit)(carry-est*(DECDPUNMAX+1));
5028 carry=est-(DECDPUNMAX+1); // correctly negative
5029 #else
5030 // remainder operator is undefined if negative, so must test
5031 if ((ueInt)carry<(DECDPUNMAX+1)*2) { // fastpath carry +1
5032 *c=(Unit)(carry-(DECDPUNMAX+1)); // [helps additions]
5033 carry=1;
5034 continue;
5035 }
5036 if (carry>=0) {
5037 *c=(Unit)(carry%(DECDPUNMAX+1));
5038 carry=carry/(DECDPUNMAX+1);
5039 continue;
5040 }
5041 // negative case
5042 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5043 *c=(Unit)(carry%(DECDPUNMAX+1));
5044 carry=carry/(DECDPUNMAX+1)-(DECDPUNMAX+1);
5045 #endif
5046 } // c
5047
5048 // now may have one or other to complete
5049 // [pretest to avoid loop setup/shutdown]
5050 if (c<maxC) for (; c<maxC; c++) {
5051 if (a<alsu+alength) { // still in A
5052 carry+=*a;
5053 a++;
5054 }
5055 else { // inside B
5056 carry+=((eInt)*b)*m;
5057 b++;
5058 }
5059 // here carry is new Unit of digits; it could be +ve or -ve and
5060 // magnitude up to DECDPUNMAX squared
5061 if ((ueInt)carry<=DECDPUNMAX) { // fastpath 0-DECDPUNMAX
5062 *c=(Unit)carry;
5063 carry=0;
5064 continue;
5065 }
5066 // result for this unit is negative or >DECDPUNMAX
5067 #if DECDPUN==4 // use divide-by-multiply
5068 if (carry>=0) {
5069 est=(((ueInt)carry>>11)*53687)>>18;
5070 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
5071 carry=est; // likely quotient [79.7%]
5072 if (*c<DECDPUNMAX+1) continue; // estimate was correct
5073 carry++;
5074 *c-=DECDPUNMAX+1;
5075 continue;
5076 }
5077 // negative case
5078 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5079 est=(((ueInt)carry>>11)*53687)>>18;
5080 *c=(Unit)(carry-est*(DECDPUNMAX+1));
5081 carry=est-(DECDPUNMAX+1); // correctly negative
5082 if (*c<DECDPUNMAX+1) continue; // was OK
5083 carry++;
5084 *c-=DECDPUNMAX+1;
5085 #elif DECDPUN==3
5086 if (carry>=0) {
5087 est=(((ueInt)carry>>3)*16777)>>21;
5088 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
5089 carry=est; // likely quotient [99%]
5090 if (*c<DECDPUNMAX+1) continue; // estimate was correct
5091 carry++;
5092 *c-=DECDPUNMAX+1;
5093 continue;
5094 }
5095 // negative case
5096 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5097 est=(((ueInt)carry>>3)*16777)>>21;
5098 *c=(Unit)(carry-est*(DECDPUNMAX+1));
5099 carry=est-(DECDPUNMAX+1); // correctly negative
5100 if (*c<DECDPUNMAX+1) continue; // was OK
5101 carry++;
5102 *c-=DECDPUNMAX+1;
5103 #elif DECDPUN<=2
5104 if (carry>=0) {
5105 est=QUOT10(carry, DECDPUN);
5106 *c=(Unit)(carry-est*(DECDPUNMAX+1)); // remainder
5107 carry=est; // quotient
5108 continue;
5109 }
5110 // negative case
5111 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5112 est=QUOT10(carry, DECDPUN);
5113 *c=(Unit)(carry-est*(DECDPUNMAX+1));
5114 carry=est-(DECDPUNMAX+1); // correctly negative
5115 #else
5116 if ((ueInt)carry<(DECDPUNMAX+1)*2){ // fastpath carry 1
5117 *c=(Unit)(carry-(DECDPUNMAX+1));
5118 carry=1;
5119 continue;
5120 }
5121 // remainder operator is undefined if negative, so must test
5122 if (carry>=0) {
5123 *c=(Unit)(carry%(DECDPUNMAX+1));
5124 carry=carry/(DECDPUNMAX+1);
5125 continue;
5126 }
5127 // negative case
5128 carry=carry+(eInt)(DECDPUNMAX+1)*(DECDPUNMAX+1); // make positive
5129 *c=(Unit)(carry%(DECDPUNMAX+1));
5130 carry=carry/(DECDPUNMAX+1)-(DECDPUNMAX+1);
5131 #endif
5132 } // c
5133
5134 // OK, all A and B processed; might still have carry or borrow
5135 // return number of Units in the result, negated if a borrow
5136 if (carry==0) return c-clsu; // no carry, so no more to do
5137 if (carry>0) { // positive carry
5138 *c=(Unit)carry; // place as new unit
5139 c++; // ..
5140 return c-clsu;
5141 }
5142 // -ve carry: it's a borrow; complement needed
5143 add=1; // temporary carry...
5144 for (c=clsu; c<maxC; c++) {
5145 add=DECDPUNMAX+add-*c;
5146 if (add<=DECDPUNMAX) {
5147 *c=(Unit)add;
5148 add=0;
5149 }
5150 else {
5151 *c=0;
5152 add=1;
5153 }
5154 }
5155 // add an extra unit iff it would be non-zero
5156 #if DECTRACE
5157 printf("UAS borrow: add %d, carry %d\n", add, carry);
5158 #endif
5159 if ((add-carry-1)!=0) {
5160 *c=(Unit)(add-carry-1);
5161 c++; // interesting, include it
5162 }
5163 return clsu-c; // -ve result indicates borrowed
5164 } // decUnitAddSub
5165
5166 /* ------------------------------------------------------------------ */
5167 /* decTrim -- trim trailing zeros or normalize */
5168 /* */
5169 /* dn is the number to trim or normalize */
5170 /* all is 1 to remove all trailing zeros, 0 for just fraction ones */
5171 /* dropped returns the number of discarded trailing zeros */
5172 /* returns dn */
5173 /* */
5174 /* All fields are updated as required. This is a utility operation, */
5175 /* so special values are unchanged and no error is possible. */
5176 /* ------------------------------------------------------------------ */
5177 static decNumber * decTrim(decNumber *dn, Flag all, Int *dropped) {
5178 Int d, exp; // work
5179 uInt cut; // ..
5180 Unit *up; // -> current Unit
5181
5182 #if DECCHECK
5183 if (decCheckOperands(dn, DECUNUSED, DECUNUSED, DECUNUSED)) return dn;
5184 #endif
5185
5186 *dropped=0; // assume no zeros dropped
5187 if ((dn->bits & DECSPECIAL) // fast exit if special ..
5188 || (*dn->lsu & 0x01)) return dn; // .. or odd
5189 if (ISZERO(dn)) { // .. or 0
5190 dn->exponent=0; // (sign is preserved)
5191 return dn;
5192 }
5193
5194 // have a finite number which is even
5195 exp=dn->exponent;
5196 cut=1; // digit (1-DECDPUN) in Unit
5197 up=dn->lsu; // -> current Unit
5198 for (d=0; d<dn->digits-1; d++) { // [don't strip the final digit]
5199 // slice by powers
5200 #if DECDPUN<=4
5201 uInt quot=QUOT10(*up, cut);
5202 if ((*up-quot*powers[cut])!=0) break; // found non-0 digit
5203 #else
5204 if (*up%powers[cut]!=0) break; // found non-0 digit
5205 #endif
5206 // have a trailing 0
5207 if (!all) { // trimming
5208 // [if exp>0 then all trailing 0s are significant for trim]
5209 if (exp<=0) { // if digit might be significant
5210 if (exp==0) break; // then quit
5211 exp++; // next digit might be significant
5212 }
5213 }
5214 cut++; // next power
5215 if (cut>DECDPUN) { // need new Unit
5216 up++;
5217 cut=1;
5218 }
5219 } // d
5220 if (d==0) return dn; // none dropped
5221
5222 // effect the drop
5223 decShiftToLeast(dn->lsu, D2U(dn->digits), d);
5224 dn->exponent+=d; // maintain numerical value
5225 dn->digits-=d; // new length
5226 *dropped=d; // report the count
5227 return dn;
5228 } // decTrim
5229
5230 /* ------------------------------------------------------------------ */
5231 /* decShiftToMost -- shift digits in array towards most significant */
5232 /* */
5233 /* uar is the array */
5234 /* digits is the count of digits in use in the array */
5235 /* shift is the number of zeros to pad with (least significant); */
5236 /* it must be zero or positive */
5237 /* */
5238 /* returns the new length of the integer in the array, in digits */
5239 /* */
5240 /* No overflow is permitted (that is, the uar array must be known to */
5241 /* be large enough to hold the result, after shifting). */
5242 /* ------------------------------------------------------------------ */
5243 static Int decShiftToMost(Unit *uar, Int digits, Int shift) {
5244 Unit *target, *source, *first; // work
5245 Int cut; // odd 0's to add
5246 uInt next; // work
5247
5248 if (shift==0) return digits; // [fastpath] nothing to do
5249 if ((digits+shift)<=DECDPUN) { // [fastpath] single-unit case
5250 *uar=(Unit)(*uar*powers[shift]);
5251 return digits+shift;
5252 }
5253
5254 next=0; // all paths
5255 source=uar+D2U(digits)-1; // where msu comes from
5256 target=source+D2U(shift); // where upper part of first cut goes
5257 cut=DECDPUN-MSUDIGITS(shift); // where to slice
5258 if (cut==0) { // unit-boundary case
5259 for (; source>=uar; source--, target--) *target=*source;
5260 }
5261 else {
5262 first=uar+D2U(digits+shift)-1; // where msu of source will end up
5263 for (; source>=uar; source--, target--) {
5264 // split the source Unit and accumulate remainder for next
5265 #if DECDPUN<=4
5266 uInt quot=QUOT10(*source, cut);
5267 uInt rem=*source-quot*powers[cut];
5268 next+=quot;
5269 #else
5270 uInt rem=*source%powers[cut];
5271 next+=*source/powers[cut];
5272 #endif
5273 if (target<=first) *target=(Unit)next; // write to target iff valid
5274 next=rem*powers[DECDPUN-cut]; // save remainder for next Unit
5275 }
5276 } // shift-move
5277
5278 // propagate any partial unit to one below and clear the rest
5279 for (; target>=uar; target--) {
5280 *target=(Unit)next;
5281 next=0;
5282 }
5283 return digits+shift;
5284 } // decShiftToMost
5285
5286 /* ------------------------------------------------------------------ */
5287 /* decShiftToLeast -- shift digits in array towards least significant */
5288 /* */
5289 /* uar is the array */
5290 /* units is length of the array, in units */
5291 /* shift is the number of digits to remove from the lsu end; it */
5292 /* must be zero or positive and less than units*DECDPUN. */
5293 /* */
5294 /* returns the new length of the integer in the array, in units */
5295 /* */
5296 /* Removed digits are discarded (lost). Units not required to hold */
5297 /* the final result are unchanged. */
5298 /* ------------------------------------------------------------------ */
5299 static Int decShiftToLeast(Unit *uar, Int units, Int shift) {
5300 Unit *target, *up; // work
5301 Int cut, count; // work
5302 Int quot, rem; // for division
5303
5304 if (shift==0) return units; // [fastpath] nothing to do
5305
5306 target=uar; // both paths
5307 cut=MSUDIGITS(shift);
5308 if (cut==DECDPUN) { // unit-boundary case; easy
5309 up=uar+D2U(shift);
5310 for (; up<uar+units; target++, up++) *target=*up;
5311 return target-uar;
5312 }
5313
5314 // messier
5315 up=uar+D2U(shift-cut); // source; correct to whole Units
5316 count=units*DECDPUN-shift; // the maximum new length
5317 #if DECDPUN<=4
5318 quot=QUOT10(*up, cut);
5319 #else
5320 quot=*up/powers[cut];
5321 #endif
5322 for (; ; target++) {
5323 *target=(Unit)quot;
5324 count-=(DECDPUN-cut);
5325 if (count<=0) break;
5326 up++;
5327 quot=*up;
5328 #if DECDPUN<=4
5329 quot=QUOT10(quot, cut);
5330 rem=*up-quot*powers[cut];
5331 #else
5332 rem=quot%powers[cut];
5333 quot=quot/powers[cut];
5334 #endif
5335 *target=(Unit)(*target+rem*powers[DECDPUN-cut]);
5336 count-=cut;
5337 if (count<=0) break;
5338 }
5339 return target-uar+1;
5340 } // decShiftToLeast
5341
5342 #if DECSUBSET
5343 /* ------------------------------------------------------------------ */
5344 /* decRoundOperand -- round an operand [used for subset only] */
5345 /* */
5346 /* dn is the number to round (dn->digits is > set->digits) */
5347 /* set is the relevant context */
5348 /* status is the status accumulator */
5349 /* */
5350 /* returns an allocated decNumber with the rounded result. */
5351 /* */
5352 /* lostDigits and other status may be set by this. */
5353 /* */
5354 /* Since the input is an operand, it must not be modified. */
5355 /* Instead, return an allocated decNumber, rounded as required. */
5356 /* It is the caller's responsibility to free the allocated storage. */
5357 /* */
5358 /* If no storage is available then the result cannot be used, so NULL */
5359 /* is returned. */
5360 /* ------------------------------------------------------------------ */
5361 static decNumber *decRoundOperand(const decNumber *dn, decContext *set,
5362 uInt *status) {
5363 decNumber *res; // result structure
5364 uInt newstatus=0; // status from round
5365 Int residue=0; // rounding accumulator
5366
5367 // Allocate storage for the returned decNumber, big enough for the
5368 // length specified by the context
5369 res=(decNumber *)malloc(sizeof(decNumber)
5370 +(D2U(set->digits)-1)*sizeof(Unit));
5371 if (res==NULL) {
5372 *status|=DEC_Insufficient_storage;
5373 return NULL;
5374 }
5375 decCopyFit(res, dn, set, &residue, &newstatus);
5376 decApplyRound(res, set, residue, &newstatus);
5377
5378 // If that set Inexact then "lost digits" is raised...
5379 if (newstatus & DEC_Inexact) newstatus|=DEC_Lost_digits;
5380 *status|=newstatus;
5381 return res;
5382 } // decRoundOperand
5383 #endif
5384
5385 /* ------------------------------------------------------------------ */
5386 /* decCopyFit -- copy a number, truncating the coefficient if needed */
5387 /* */
5388 /* dest is the target decNumber */
5389 /* src is the source decNumber */
5390 /* set is the context [used for length (digits) and rounding mode] */
5391 /* residue is the residue accumulator */
5392 /* status contains the current status to be updated */
5393 /* */
5394 /* (dest==src is allowed and will be a no-op if fits) */
5395 /* All fields are updated as required. */
5396 /* ------------------------------------------------------------------ */
5397 static void decCopyFit(decNumber *dest, const decNumber *src,
5398 decContext *set, Int *residue, uInt *status) {
5399 dest->bits=src->bits;
5400 dest->exponent=src->exponent;
5401 decSetCoeff(dest, set, src->lsu, src->digits, residue, status);
5402 } // decCopyFit
5403
5404 /* ------------------------------------------------------------------ */
5405 /* decSetCoeff -- set the coefficient of a number */
5406 /* */
5407 /* dn is the number whose coefficient array is to be set. */
5408 /* It must have space for set->digits digits */
5409 /* set is the context [for size] */
5410 /* lsu -> lsu of the source coefficient [may be dn->lsu] */
5411 /* len is digits in the source coefficient [may be dn->digits] */
5412 /* residue is the residue accumulator. This has values as in */
5413 /* decApplyRound, and will be unchanged unless the */
5414 /* target size is less than len. In this case, the */
5415 /* coefficient is truncated and the residue is updated to */
5416 /* reflect the previous residue and the dropped digits. */
5417 /* status is the status accumulator, as usual */
5418 /* */
5419 /* The coefficient may already be in the number, or it can be an */
5420 /* external intermediate array. If it is in the number, lsu must == */
5421 /* dn->lsu and len must == dn->digits. */
5422 /* */
5423 /* Note that the coefficient length (len) may be < set->digits, and */
5424 /* in this case this merely copies the coefficient (or is a no-op */
5425 /* if dn->lsu==lsu). */
5426 /* */
5427 /* Note also that (only internally, from decQuantizeOp and */
5428 /* decSetSubnormal) the value of set->digits may be less than one, */
5429 /* indicating a round to left. This routine handles that case */
5430 /* correctly; caller ensures space. */
5431 /* */
5432 /* dn->digits, dn->lsu (and as required), and dn->exponent are */
5433 /* updated as necessary. dn->bits (sign) is unchanged. */
5434 /* */
5435 /* DEC_Rounded status is set if any digits are discarded. */
5436 /* DEC_Inexact status is set if any non-zero digits are discarded, or */
5437 /* incoming residue was non-0 (implies rounded) */
5438 /* ------------------------------------------------------------------ */
5439 // mapping array: maps 0-9 to canonical residues, so that a residue
5440 // can be adjusted in the range [-1, +1] and achieve correct rounding
5441 // 0 1 2 3 4 5 6 7 8 9
5442 static const uByte resmap[10]={0, 3, 3, 3, 3, 5, 7, 7, 7, 7};
5443 static void decSetCoeff(decNumber *dn, decContext *set, const Unit *lsu,
5444 Int len, Int *residue, uInt *status) {
5445 Int discard; // number of digits to discard
5446 uInt cut; // cut point in Unit
5447 const Unit *up; // work
5448 Unit *target; // ..
5449 Int count; // ..
5450 #if DECDPUN<=4
5451 uInt temp; // ..
5452 #endif
5453
5454 discard=len-set->digits; // digits to discard
5455 if (discard<=0) { // no digits are being discarded
5456 if (dn->lsu!=lsu) { // copy needed
5457 // copy the coefficient array to the result number; no shift needed
5458 count=len; // avoids D2U
5459 up=lsu;
5460 for (target=dn->lsu; count>0; target++, up++, count-=DECDPUN)
5461 *target=*up;
5462 dn->digits=len; // set the new length
5463 }
5464 // dn->exponent and residue are unchanged, record any inexactitude
5465 if (*residue!=0) *status|=(DEC_Inexact | DEC_Rounded);
5466 return;
5467 }
5468
5469 // some digits must be discarded ...
5470 dn->exponent+=discard; // maintain numerical value
5471 *status|=DEC_Rounded; // accumulate Rounded status
5472 if (*residue>1) *residue=1; // previous residue now to right, so reduce
5473
5474 if (discard>len) { // everything, +1, is being discarded
5475 // guard digit is 0
5476 // residue is all the number [NB could be all 0s]
5477 if (*residue<=0) { // not already positive
5478 count=len; // avoids D2U
5479 for (up=lsu; count>0; up++, count-=DECDPUN) if (*up!=0) { // found non-0
5480 *residue=1;
5481 break; // no need to check any others
5482 }
5483 }
5484 if (*residue!=0) *status|=DEC_Inexact; // record inexactitude
5485 *dn->lsu=0; // coefficient will now be 0
5486 dn->digits=1; // ..
5487 return;
5488 } // total discard
5489
5490 // partial discard [most common case]
5491 // here, at least the first (most significant) discarded digit exists
5492
5493 // spin up the number, noting residue during the spin, until get to
5494 // the Unit with the first discarded digit. When reach it, extract
5495 // it and remember its position
5496 count=0;
5497 for (up=lsu;; up++) {
5498 count+=DECDPUN;
5499 if (count>=discard) break; // full ones all checked
5500 if (*up!=0) *residue=1;
5501 } // up
5502
5503 // here up -> Unit with first discarded digit
5504 cut=discard-(count-DECDPUN)-1;
5505 if (cut==DECDPUN-1) { // unit-boundary case (fast)
5506 Unit half=(Unit)powers[DECDPUN]>>1;
5507 // set residue directly
5508 if (*up>=half) {
5509 if (*up>half) *residue=7;
5510 else *residue+=5; // add sticky bit
5511 }
5512 else { // <half
5513 if (*up!=0) *residue=3; // [else is 0, leave as sticky bit]
5514 }
5515 if (set->digits<=0) { // special for Quantize/Subnormal :-(
5516 *dn->lsu=0; // .. result is 0
5517 dn->digits=1; // ..
5518 }
5519 else { // shift to least
5520 count=set->digits; // now digits to end up with
5521 dn->digits=count; // set the new length
5522 up++; // move to next
5523 // on unit boundary, so shift-down copy loop is simple
5524 for (target=dn->lsu; count>0; target++, up++, count-=DECDPUN)
5525 *target=*up;
5526 }
5527 } // unit-boundary case
5528
5529 else { // discard digit is in low digit(s), and not top digit
5530 uInt discard1; // first discarded digit
5531 uInt quot, rem; // for divisions
5532 if (cut==0) quot=*up; // is at bottom of unit
5533 else /* cut>0 */ { // it's not at bottom of unit
5534 #if DECDPUN<=4
5535 quot=QUOT10(*up, cut);
5536 rem=*up-quot*powers[cut];
5537 #else
5538 rem=*up%powers[cut];
5539 quot=*up/powers[cut];
5540 #endif
5541 if (rem!=0) *residue=1;
5542 }
5543 // discard digit is now at bottom of quot
5544 #if DECDPUN<=4
5545 temp=(quot*6554)>>16; // fast /10
5546 // Vowels algorithm here not a win (9 instructions)
5547 discard1=quot-X10(temp);
5548 quot=temp;
5549 #else
5550 discard1=quot%10;
5551 quot=quot/10;
5552 #endif
5553 // here, discard1 is the guard digit, and residue is everything
5554 // else [use mapping array to accumulate residue safely]
5555 *residue+=resmap[discard1];
5556 cut++; // update cut
5557 // here: up -> Unit of the array with bottom digit
5558 // cut is the division point for each Unit
5559 // quot holds the uncut high-order digits for the current unit
5560 if (set->digits<=0) { // special for Quantize/Subnormal :-(
5561 *dn->lsu=0; // .. result is 0
5562 dn->digits=1; // ..
5563 }
5564 else { // shift to least needed
5565 count=set->digits; // now digits to end up with
5566 dn->digits=count; // set the new length
5567 // shift-copy the coefficient array to the result number
5568 for (target=dn->lsu; ; target++) {
5569 *target=(Unit)quot;
5570 count-=(DECDPUN-cut);
5571 if (count<=0) break;
5572 up++;
5573 quot=*up;
5574 #if DECDPUN<=4
5575 quot=QUOT10(quot, cut);
5576 rem=*up-quot*powers[cut];
5577 #else
5578 rem=quot%powers[cut];
5579 quot=quot/powers[cut];
5580 #endif
5581 *target=(Unit)(*target+rem*powers[DECDPUN-cut]);
5582 count-=cut;
5583 if (count<=0) break;
5584 } // shift-copy loop
5585 } // shift to least
5586 } // not unit boundary
5587
5588 if (*residue!=0) *status|=DEC_Inexact; // record inexactitude
5589 return;
5590 } // decSetCoeff
5591
5592 /* ------------------------------------------------------------------ */
5593 /* decApplyRound -- apply pending rounding to a number */
5594 /* */
5595 /* dn is the number, with space for set->digits digits */
5596 /* set is the context [for size and rounding mode] */
5597 /* residue indicates pending rounding, being any accumulated */
5598 /* guard and sticky information. It may be: */
5599 /* 6-9: rounding digit is >5 */
5600 /* 5: rounding digit is exactly half-way */
5601 /* 1-4: rounding digit is <5 and >0 */
5602 /* 0: the coefficient is exact */
5603 /* -1: as 1, but the hidden digits are subtractive, that */
5604 /* is, of the opposite sign to dn. In this case the */
5605 /* coefficient must be non-0. */
5606 /* status is the status accumulator, as usual */
5607 /* */
5608 /* This routine applies rounding while keeping the length of the */
5609 /* coefficient constant. The exponent and status are unchanged */
5610 /* except if: */
5611 /* */
5612 /* -- the coefficient was increased and is all nines (in which */
5613 /* case Overflow could occur, and is handled directly here so */
5614 /* the caller does not need to re-test for overflow) */
5615 /* */
5616 /* -- the coefficient was decreased and becomes all nines (in which */
5617 /* case Underflow could occur, and is also handled directly). */
5618 /* */
5619 /* All fields in dn are updated as required. */
5620 /* */
5621 /* ------------------------------------------------------------------ */
5622 static void decApplyRound(decNumber *dn, decContext *set, Int residue,
5623 uInt *status) {
5624 Int bump; // 1 if coefficient needs to be incremented
5625 // -1 if coefficient needs to be decremented
5626
5627 if (residue==0) return; // nothing to apply
5628
5629 bump=0; // assume a smooth ride
5630
5631 // now decide whether, and how, to round, depending on mode
5632 switch (set->round) {
5633 case DEC_ROUND_DOWN: {
5634 // no change, except if negative residue
5635 if (residue<0) bump=-1;
5636 break;} // r-d
5637
5638 case DEC_ROUND_HALF_DOWN: {
5639 if (residue>5) bump=1;
5640 break;} // r-h-d
5641
5642 case DEC_ROUND_HALF_EVEN: {
5643 if (residue>5) bump=1; // >0.5 goes up
5644 else if (residue==5) { // exactly 0.5000...
5645 // 0.5 goes up iff [new] lsd is odd
5646 if (*dn->lsu & 0x01) bump=1;
5647 }
5648 break;} // r-h-e
5649
5650 case DEC_ROUND_HALF_UP: {
5651 if (residue>=5) bump=1;
5652 break;} // r-h-u
5653
5654 case DEC_ROUND_UP: {
5655 if (residue>0) bump=1;
5656 break;} // r-u
5657
5658 case DEC_ROUND_CEILING: {
5659 // same as _UP for positive numbers, and as _DOWN for negatives
5660 // [negative residue cannot occur on 0]
5661 if (decNumberIsNegative(dn)) {
5662 if (residue<0) bump=-1;
5663 }
5664 else {
5665 if (residue>0) bump=1;
5666 }
5667 break;} // r-c
5668
5669 case DEC_ROUND_FLOOR: {
5670 // same as _UP for negative numbers, and as _DOWN for positive
5671 // [negative residue cannot occur on 0]
5672 if (!decNumberIsNegative(dn)) {
5673 if (residue<0) bump=-1;
5674 }
5675 else {
5676 if (residue>0) bump=1;
5677 }
5678 break;} // r-f
5679
5680 default: { // e.g., DEC_ROUND_MAX
5681 *status|=DEC_Invalid_context;
5682 #if DECTRACE
5683 printf("Unknown rounding mode: %d\n", set->round);
5684 #endif
5685 break;}
5686 } // switch
5687
5688 // now bump the number, up or down, if need be
5689 if (bump==0) return; // no action required
5690
5691 // Simply use decUnitAddSub unless bumping up and the number is
5692 // all nines. In this special case set to 100... explicitly
5693 // and adjust the exponent by one (as otherwise could overflow
5694 // the array)
5695 // Similarly handle all-nines result if bumping down.
5696 if (bump>0) {
5697 Unit *up; // work
5698 uInt count=dn->digits; // digits to be checked
5699 for (up=dn->lsu; ; up++) {
5700 if (count<=DECDPUN) {
5701 // this is the last Unit (the msu)
5702 if (*up!=powers[count]-1) break; // not still 9s
5703 // here if it, too, is all nines
5704 *up=(Unit)powers[count-1]; // here 999 -> 100 etc.
5705 for (up=up-1; up>=dn->lsu; up--) *up=0; // others all to 0
5706 dn->exponent++; // and bump exponent
5707 // [which, very rarely, could cause Overflow...]
5708 if ((dn->exponent+dn->digits)>set->emax+1) {
5709 decSetOverflow(dn, set, status);
5710 }
5711 return; // done
5712 }
5713 // a full unit to check, with more to come
5714 if (*up!=DECDPUNMAX) break; // not still 9s
5715 count-=DECDPUN;
5716 } // up
5717 } // bump>0
5718 else { // -1
5719 // here checking for a pre-bump of 1000... (leading 1, all
5720 // other digits zero)
5721 Unit *up, *sup; // work
5722 uInt count=dn->digits; // digits to be checked
5723 for (up=dn->lsu; ; up++) {
5724 if (count<=DECDPUN) {
5725 // this is the last Unit (the msu)
5726 if (*up!=powers[count-1]) break; // not 100..
5727 // here if have the 1000... case
5728 sup=up; // save msu pointer
5729 *up=(Unit)powers[count]-1; // here 100 in msu -> 999
5730 // others all to all-nines, too
5731 for (up=up-1; up>=dn->lsu; up--) *up=(Unit)powers[DECDPUN]-1;
5732 dn->exponent--; // and bump exponent
5733
5734 // iff the number was at the subnormal boundary (exponent=etiny)
5735 // then the exponent is now out of range, so it will in fact get
5736 // clamped to etiny and the final 9 dropped.
5737 // printf(">> emin=%d exp=%d sdig=%d\n", set->emin,
5738 // dn->exponent, set->digits);
5739 if (dn->exponent+1==set->emin-set->digits+1) {
5740 if (count==1 && dn->digits==1) *sup=0; // here 9 -> 0[.9]
5741 else {
5742 *sup=(Unit)powers[count-1]-1; // here 999.. in msu -> 99..
5743 dn->digits--;
5744 }
5745 dn->exponent++;
5746 *status|=DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5747 }
5748 return; // done
5749 }
5750
5751 // a full unit to check, with more to come
5752 if (*up!=0) break; // not still 0s
5753 count-=DECDPUN;
5754 } // up
5755
5756 } // bump<0
5757
5758 // Actual bump needed. Do it.
5759 decUnitAddSub(dn->lsu, D2U(dn->digits), uarrone, 1, 0, dn->lsu, bump);
5760 } // decApplyRound
5761
5762 #if DECSUBSET
5763 /* ------------------------------------------------------------------ */
5764 /* decFinish -- finish processing a number */
5765 /* */
5766 /* dn is the number */
5767 /* set is the context */
5768 /* residue is the rounding accumulator (as in decApplyRound) */
5769 /* status is the accumulator */
5770 /* */
5771 /* This finishes off the current number by: */
5772 /* 1. If not extended: */
5773 /* a. Converting a zero result to clean '0' */
5774 /* b. Reducing positive exponents to 0, if would fit in digits */
5775 /* 2. Checking for overflow and subnormals (always) */
5776 /* Note this is just Finalize when no subset arithmetic. */
5777 /* All fields are updated as required. */
5778 /* ------------------------------------------------------------------ */
5779 static void decFinish(decNumber *dn, decContext *set, Int *residue,
5780 uInt *status) {
5781 if (!set->extended) {
5782 if ISZERO(dn) { // value is zero
5783 dn->exponent=0; // clean exponent ..
5784 dn->bits=0; // .. and sign
5785 return; // no error possible
5786 }
5787 if (dn->exponent>=0) { // non-negative exponent
5788 // >0; reduce to integer if possible
5789 if (set->digits >= (dn->exponent+dn->digits)) {
5790 dn->digits=decShiftToMost(dn->lsu, dn->digits, dn->exponent);
5791 dn->exponent=0;
5792 }
5793 }
5794 } // !extended
5795
5796 decFinalize(dn, set, residue, status);
5797 } // decFinish
5798 #endif
5799
5800 /* ------------------------------------------------------------------ */
5801 /* decFinalize -- final check, clamp, and round of a number */
5802 /* */
5803 /* dn is the number */
5804 /* set is the context */
5805 /* residue is the rounding accumulator (as in decApplyRound) */
5806 /* status is the status accumulator */
5807 /* */
5808 /* This finishes off the current number by checking for subnormal */
5809 /* results, applying any pending rounding, checking for overflow, */
5810 /* and applying any clamping. */
5811 /* Underflow and overflow conditions are raised as appropriate. */
5812 /* All fields are updated as required. */
5813 /* ------------------------------------------------------------------ */
5814 static void decFinalize(decNumber *dn, decContext *set, Int *residue,
5815 uInt *status) {
5816 Int shift; // shift needed if clamping
5817 Int tinyexp=set->emin-dn->digits+1; // precalculate subnormal boundary
5818
5819 // Must be careful, here, when checking the exponent as the
5820 // adjusted exponent could overflow 31 bits [because it may already
5821 // be up to twice the expected].
5822
5823 // First test for subnormal. This must be done before any final
5824 // round as the result could be rounded to Nmin or 0.
5825 if (dn->exponent<=tinyexp) { // prefilter
5826 decNumber nmin;
5827 // A very nasty case here is dn == Nmin and residue<0
5828 if (dn->exponent<tinyexp) {
5829 // Go handle subnormals; this will apply round if needed.
5830 decSetSubnormal(dn, set, residue, status);
5831 return;
5832 }
5833 // Equals case: only subnormal if dn=Nmin and negative residue
5834 decNumberZero(&nmin);
5835 nmin.lsu[0]=1;
5836 nmin.exponent=set->emin;
5837 if (*residue<0 && decCompare(dn, &nmin, 1)==0) { // (signless compare)
5838 decApplyRound(dn, set, *residue, status); // might force down
5839 decSetSubnormal(dn, set, residue, status);
5840 return;
5841 }
5842 }
5843
5844 // now apply any pending round (this could raise overflow).
5845 if (*residue!=0) decApplyRound(dn, set, *residue, status);
5846
5847 // Check for overflow [redundant in the 'rare' case] or clamp
5848 if (dn->exponent<=set->emax-set->digits+1) return; // neither needed
5849
5850 // here when might have an overflow or clamp to do
5851 if (dn->exponent>set->emax-dn->digits+1) { // too big
5852 decSetOverflow(dn, set, status);
5853 return;
5854 }
5855 // here when the result is normal but in clamp range
5856 if (!set->clamp) return;
5857
5858 // here when need to apply the IEEE exponent clamp (fold-down)
5859 shift=dn->exponent-(set->emax-set->digits+1);
5860
5861 // shift coefficient (if non-zero)
5862 if (!ISZERO(dn)) {
5863 dn->digits=decShiftToMost(dn->lsu, dn->digits, shift);
5864 }
5865 dn->exponent-=shift; // adjust the exponent to match
5866 *status|=DEC_Clamped; // and record the dirty deed
5867 return;
5868 } // decFinalize
5869
5870 /* ------------------------------------------------------------------ */
5871 /* decSetOverflow -- set number to proper overflow value */
5872 /* */
5873 /* dn is the number (used for sign [only] and result) */
5874 /* set is the context [used for the rounding mode] */
5875 /* status contains the current status to be updated */
5876 /* */
5877 /* This sets the sign of a number and sets its value to either */
5878 /* Infinity or the maximum finite value, depending on the sign of */
5879 /* dn and therounding mode, following IEEE 854 rules. */
5880 /* ------------------------------------------------------------------ */
5881 static void decSetOverflow(decNumber *dn, decContext *set, uInt *status) {
5882 Flag needmax=0; // result is maximum finite value
5883 uByte sign=dn->bits&DECNEG; // clean and save sign bit
5884
5885 if (ISZERO(dn)) { // zero does not overflow magnitude
5886 Int emax=set->emax; // limit value
5887 if (set->clamp) emax-=set->digits-1; // lower if clamping
5888 if (dn->exponent>emax) { // clamp required
5889 dn->exponent=emax;
5890 *status|=DEC_Clamped;
5891 }
5892 return;
5893 }
5894
5895 decNumberZero(dn);
5896 switch (set->round) {
5897 case DEC_ROUND_DOWN: {
5898 needmax=1; // never Infinity
5899 break;} // r-d
5900 case DEC_ROUND_CEILING: {
5901 if (sign) needmax=1; // Infinity if non-negative
5902 break;} // r-c
5903 case DEC_ROUND_FLOOR: {
5904 if (!sign) needmax=1; // Infinity if negative
5905 break;} // r-f
5906 default: break; // Infinity in all other cases
5907 }
5908 if (needmax) {
5909 Unit *up; // work
5910 Int count=set->digits; // nines to add
5911 dn->digits=count;
5912 // fill in all nines to set maximum value
5913 for (up=dn->lsu; ; up++) {
5914 if (count>DECDPUN) *up=DECDPUNMAX; // unit full o'nines
5915 else { // this is the msu
5916 *up=(Unit)(powers[count]-1);
5917 break;
5918 }
5919 count-=DECDPUN; // filled those digits
5920 } // up
5921 dn->bits=sign; // sign
5922 dn->exponent=set->emax-set->digits+1;
5923 }
5924 else dn->bits=sign|DECINF; // Value is +/-Infinity
5925 *status|=DEC_Overflow | DEC_Inexact | DEC_Rounded;
5926 } // decSetOverflow
5927
5928 /* ------------------------------------------------------------------ */
5929 /* decSetSubnormal -- process value whose exponent is <Emin */
5930 /* */
5931 /* dn is the number (used as input as well as output; it may have */
5932 /* an allowed subnormal value, which may need to be rounded) */
5933 /* set is the context [used for the rounding mode] */
5934 /* residue is any pending residue */
5935 /* status contains the current status to be updated */
5936 /* */
5937 /* If subset mode, set result to zero and set Underflow flags. */
5938 /* */
5939 /* Value may be zero with a low exponent; this does not set Subnormal */
5940 /* but the exponent will be clamped to Etiny. */
5941 /* */
5942 /* Otherwise ensure exponent is not out of range, and round as */
5943 /* necessary. Underflow is set if the result is Inexact. */
5944 /* ------------------------------------------------------------------ */
5945 static void decSetSubnormal(decNumber *dn, decContext *set, Int *residue,
5946 uInt *status) {
5947 Int dnexp; // saves original exponent
5948 decContext workset; // work
5949 Int etiny, adjust; // ..
5950
5951 #if DECSUBSET
5952 // simple set to zero and 'hard underflow' for subset
5953 if (!set->extended) {
5954 decNumberZero(dn);
5955 // always full overflow
5956 *status|=DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5957 return;
5958 }
5959 #endif
5960
5961 // Full arithmetic -- allow subnormals, rounded to minimum exponent
5962 // (Etiny) if needed
5963 etiny=set->emin-(set->digits-1); // smallest allowed exponent
5964
5965 if ISZERO(dn) { // value is zero
5966 // residue can never be non-zero here
5967 #if DECCHECK
5968 if (*residue!=0) {
5969 printf("++ Subnormal 0 residue %d\n", *residue);
5970 *status|=DEC_Invalid_operation;
5971 }
5972 #endif
5973 if (dn->exponent<etiny) { // clamp required
5974 dn->exponent=etiny;
5975 *status|=DEC_Clamped;
5976 }
5977 return;
5978 }
5979
5980 *status|=DEC_Subnormal; // have a non-zero subnormal
5981 adjust=etiny-dn->exponent; // calculate digits to remove
5982 if (adjust<=0) { // not out of range; unrounded
5983 // residue can never be non-zero here, except in the Nmin-residue
5984 // case (which is a subnormal result), so can take fast-path here
5985 // it may already be inexact (from setting the coefficient)
5986 if (*status&DEC_Inexact) *status|=DEC_Underflow;
5987 return;
5988 }
5989
5990 // adjust>0, so need to rescale the result so exponent becomes Etiny
5991 // [this code is similar to that in rescale]
5992 dnexp=dn->exponent; // save exponent
5993 workset=*set; // clone rounding, etc.
5994 workset.digits=dn->digits-adjust; // set requested length
5995 workset.emin-=adjust; // and adjust emin to match
5996 // [note that the latter can be <1, here, similar to Rescale case]
5997 decSetCoeff(dn, &workset, dn->lsu, dn->digits, residue, status);
5998 decApplyRound(dn, &workset, *residue, status);
5999
6000 // Use 754R/854 default rule: Underflow is set iff Inexact
6001 // [independent of whether trapped]
6002 if (*status&DEC_Inexact) *status|=DEC_Underflow;
6003
6004 // if rounded up a 999s case, exponent will be off by one; adjust
6005 // back if so [it will fit, because it was shortened earlier]
6006 if (dn->exponent>etiny) {
6007 dn->digits=decShiftToMost(dn->lsu, dn->digits, 1);
6008 dn->exponent--; // (re)adjust the exponent.
6009 }
6010
6011 // if rounded to zero, it is by definition clamped...
6012 if (ISZERO(dn)) *status|=DEC_Clamped;
6013 } // decSetSubnormal
6014
6015 /* ------------------------------------------------------------------ */
6016 /* decCheckMath - check entry conditions for a math function */
6017 /* */
6018 /* This checks the context and the operand */
6019 /* */
6020 /* rhs is the operand to check */
6021 /* set is the context to check */
6022 /* status is unchanged if both are good */
6023 /* */
6024 /* returns non-zero if status is changed, 0 otherwise */
6025 /* */
6026 /* Restrictions enforced: */
6027 /* */
6028 /* digits, emax, and -emin in the context must be less than */
6029 /* DEC_MAX_MATH (999999), and A must be within these bounds if */
6030 /* non-zero. Invalid_operation is set in the status if a */
6031 /* restriction is violated. */
6032 /* ------------------------------------------------------------------ */
6033 static uInt decCheckMath(const decNumber *rhs, decContext *set,
6034 uInt *status) {
6035 uInt save=*status; // record
6036 if (set->digits>DEC_MAX_MATH
6037 || set->emax>DEC_MAX_MATH
6038 || -set->emin>DEC_MAX_MATH) *status|=DEC_Invalid_context;
6039 else if ((rhs->digits>DEC_MAX_MATH
6040 || rhs->exponent+rhs->digits>DEC_MAX_MATH+1
6041 || rhs->exponent+rhs->digits<2*(1-DEC_MAX_MATH))
6042 && !ISZERO(rhs)) *status|=DEC_Invalid_operation;
6043 return (*status!=save);
6044 } // decCheckMath
6045
6046 /* ------------------------------------------------------------------ */
6047 /* decGetInt -- get integer from a number */
6048 /* */
6049 /* dn is the number [which will not be altered] */
6050 /* */
6051 /* returns one of: */
6052 /* BADINT if there is a non-zero fraction */
6053 /* the converted integer */
6054 /* BIGEVEN if the integer is even and > 2*10**9 */
6055 /* BIGODD if the integer is odd and > 2*10**9 */
6056 /* */
6057 /* This checks and gets a whole number from the input decNumber. */
6058 /* The sign can be determined from dn by the caller when BIGEVEN or */
6059 /* BIGODD is returned. */
6060 /* ------------------------------------------------------------------ */
6061 static Int decGetInt(const decNumber *dn) {
6062 Int theInt; // result accumulator
6063 const Unit *up; // work
6064 Int got; // digits (real or not) processed
6065 Int ilength=dn->digits+dn->exponent; // integral length
6066 Flag neg=decNumberIsNegative(dn); // 1 if -ve
6067
6068 // The number must be an integer that fits in 10 digits
6069 // Assert, here, that 10 is enough for any rescale Etiny
6070 #if DEC_MAX_EMAX > 999999999
6071 #error GetInt may need updating [for Emax]
6072 #endif
6073 #if DEC_MIN_EMIN < -999999999
6074 #error GetInt may need updating [for Emin]
6075 #endif
6076 if (ISZERO(dn)) return 0; // zeros are OK, with any exponent
6077
6078 up=dn->lsu; // ready for lsu
6079 theInt=0; // ready to accumulate
6080 if (dn->exponent>=0) { // relatively easy
6081 // no fractional part [usual]; allow for positive exponent
6082 got=dn->exponent;
6083 }
6084 else { // -ve exponent; some fractional part to check and discard
6085 Int count=-dn->exponent; // digits to discard
6086 // spin up whole units until reach the Unit with the unit digit
6087 for (; count>=DECDPUN; up++) {
6088 if (*up!=0) return BADINT; // non-zero Unit to discard
6089 count-=DECDPUN;
6090 }
6091 if (count==0) got=0; // [a multiple of DECDPUN]
6092 else { // [not multiple of DECDPUN]
6093 Int rem; // work
6094 // slice off fraction digits and check for non-zero
6095 #if DECDPUN<=4
6096 theInt=QUOT10(*up, count);
6097 rem=*up-theInt*powers[count];
6098 #else
6099 rem=*up%powers[count]; // slice off discards
6100 theInt=*up/powers[count];
6101 #endif
6102 if (rem!=0) return BADINT; // non-zero fraction
6103 // it looks good
6104 got=DECDPUN-count; // number of digits so far
6105 up++; // ready for next
6106 }
6107 }
6108 // now it's known there's no fractional part
6109
6110 // tricky code now, to accumulate up to 9.3 digits
6111 if (got==0) {theInt=*up; got+=DECDPUN; up++;} // ensure lsu is there
6112
6113 if (ilength<11) {
6114 Int save=theInt;
6115 // collect any remaining unit(s)
6116 for (; got<ilength; up++) {
6117 theInt+=*up*powers[got];
6118 got+=DECDPUN;
6119 }
6120 if (ilength==10) { // need to check for wrap
6121 if (theInt/(Int)powers[got-DECDPUN]!=(Int)*(up-1)) ilength=11;
6122 // [that test also disallows the BADINT result case]
6123 else if (neg && theInt>1999999997) ilength=11;
6124 else if (!neg && theInt>999999999) ilength=11;
6125 if (ilength==11) theInt=save; // restore correct low bit
6126 }
6127 }
6128
6129 if (ilength>10) { // too big
6130 if (theInt&1) return BIGODD; // bottom bit 1
6131 return BIGEVEN; // bottom bit 0
6132 }
6133
6134 if (neg) theInt=-theInt; // apply sign
6135 return theInt;
6136 } // decGetInt
6137
6138 /* ------------------------------------------------------------------ */
6139 /* decPutInt -- put integer into a number */
6140 /* */
6141 /* res is the target number, with enough space for the biggest */
6142 /* integer that the second argument will be */
6143 /* in is the input integer */
6144 /* */
6145 /* Returns res, an integral value; no error is possible. */
6146 /* ------------------------------------------------------------------ */
6147 static decNumber *decPutInt(decNumber *res, Int in) {
6148 Unit *up; // work pointer
6149
6150 decNumberZero(res); // clean
6151 if (in<=0) { // handle sign, 0, and BADINT
6152 if (in==0) return res; // easy
6153 else if (in!=BADINT) {
6154 res->bits=DECNEG; // set sign
6155 in=-in; // invert and drop through
6156 }
6157 else { // BADINT: invert would fail
6158 decContext set;
6159 decContextDefault(&set, DEC_INIT_DECIMAL64); // 16 digits
6160 decNumberFromString(res, "-2147483648", &set); // ugh
6161 return res;
6162 }
6163 }
6164 // in is now positive
6165 for (up=res->lsu; in>0; up++) {
6166 *up=(Unit)(in%(DECDPUNMAX+1));
6167 in=in/(DECDPUNMAX+1);
6168 }
6169 res->digits=decGetDigits(res->lsu, up-res->lsu);
6170 return res;
6171 } // decPutInt
6172
6173 /* ------------------------------------------------------------------ */
6174 /* decBiStr -- compare string with pairwise options */
6175 /* */
6176 /* targ is the string to compare */
6177 /* str1 is one of the strings to compare against (length may be 0) */
6178 /* str2 is the other; it must be the same length as str1 */
6179 /* */
6180 /* returns 1 if strings compare equal, (that is, it is the same */
6181 /* length as str1 and str2, and each character of targ is in either */
6182 /* str1 or str2 in the corresponding position), or 0 otherwise */
6183 /* */
6184 /* This is used for generic caseless compare, including the awkward */
6185 /* case of the Turkish dotted and dotless Is. Use as (for example): */
6186 /* if (decBiStr(test, "mike", "MIKE")) ... */
6187 /* ------------------------------------------------------------------ */
6188 static Flag decBiStr(const char *targ, const char *str1, const char *str2) {
6189 for (;;targ++, str1++, str2++) {
6190 if (*targ!=*str1 && *targ!=*str2) return 0;
6191 // *targ has a match in one (or both, if terminator)
6192 if (*targ=='\0') break;
6193 } // forever
6194 return 1;
6195 } // decBiStr
6196
6197 /* ------------------------------------------------------------------ */
6198 /* decNaNs -- handle NaN operand or operands */
6199 /* */
6200 /* res is the result number */
6201 /* lhs is the first operand */
6202 /* rhs is the second operand, or NULL if none */
6203 /* status contains the current status */
6204 /* returns res in case convenient */
6205 /* */
6206 /* Called when one or both operands is a NaN, and propagates the */
6207 /* appropriate result to res. When an sNaN is found, it is changed */
6208 /* to a qNaN and Invalid operation is set. */
6209 /* ------------------------------------------------------------------ */
6210 static decNumber * decNaNs(decNumber *res, const decNumber *lhs,
6211 const decNumber *rhs, uInt *status) {
6212 // This decision tree ends up with LHS being the source pointer,
6213 // and status updated if need be
6214 if (lhs->bits & DECSNAN)
6215 *status|=DEC_Invalid_operation | DEC_sNaN;
6216 else if (rhs==NULL);
6217 else if (rhs->bits & DECSNAN) {
6218 lhs=rhs;
6219 *status|=DEC_Invalid_operation | DEC_sNaN;
6220 }
6221 else if (lhs->bits & DECNAN);
6222 else lhs=rhs;
6223
6224 decNumberCopy(res, lhs);
6225 res->bits&=~DECSNAN; // convert any sNaN to NaN, while
6226 res->bits|=DECNAN; // .. preserving sign
6227 res->exponent=0; // clean exponent
6228 // [coefficient was copied]
6229 return res;
6230 } // decNaNs
6231
6232 /* ------------------------------------------------------------------ */
6233 /* decStatus -- apply non-zero status */
6234 /* */
6235 /* dn is the number to set if error */
6236 /* status contains the current status (not yet in context) */
6237 /* set is the context */
6238 /* */
6239 /* If the status is an error status, the number is set to a NaN, */
6240 /* unless the error was an overflow, divide-by-zero, or underflow, */
6241 /* in which case the number will have already been set. */
6242 /* */
6243 /* The context status is then updated with the new status. Note that */
6244 /* this may raise a signal, so control may never return from this */
6245 /* routine (hence resources must be recovered before it is called). */
6246 /* ------------------------------------------------------------------ */
6247 static void decStatus(decNumber *dn, uInt status, decContext *set) {
6248 if (status & DEC_NaNs) { // error status -> NaN
6249 // if cause was an sNaN, clear and propagate [NaN is already set up]
6250 if (status & DEC_sNaN) status&=~DEC_sNaN;
6251 else {
6252 decNumberZero(dn); // other error: clean throughout
6253 dn->bits=DECNAN; // and make a quiet NaN
6254 }
6255 }
6256 decContextSetStatus(set, status); // [may not return]
6257 return;
6258 } // decStatus
6259
6260 /* ------------------------------------------------------------------ */
6261 /* decGetDigits -- count digits in a Units array */
6262 /* */
6263 /* uar is the Unit array holding the number (this is often an */
6264 /* accumulator of some sort) */
6265 /* len is the length of the array in units */
6266 /* */
6267 /* returns the number of (significant) digits in the array */
6268 /* */
6269 /* All leading zeros are excluded, except the last if the array has */
6270 /* only zero Units. */
6271 /* ------------------------------------------------------------------ */
6272 // This may be called twice during some operations.
6273 static Int decGetDigits(Unit *uar, Int len) {
6274 Unit *up=uar+(len-1); // -> msu
6275 Int digits=(len-1)*DECDPUN+1; // possible digits excluding msu
6276 // (at least 1 in final msu)
6277 #if DECDPUN>4
6278 uInt const *pow; // work
6279 #endif
6280
6281 for (; up>=uar; up--) {
6282 if (*up==0) { // unit is all 0s
6283 if (digits==1) break; // a zero has one digit
6284 digits-=DECDPUN; // adjust for 0 unit
6285 continue;}
6286 // found the first (most significant) non-zero Unit
6287 #if DECDPUN>1 // not done yet
6288 if (*up<10) break; // is 1-9
6289 digits++;
6290 #if DECDPUN>2 // not done yet
6291 if (*up<100) break; // is 10-99
6292 digits++;
6293 #if DECDPUN>3 // not done yet
6294 if (*up<1000) break; // is 100-999
6295 digits++;
6296 #if DECDPUN>4 // count the rest ...
6297 for (pow=&powers[4]; *up>=*pow; pow++) digits++;
6298 #endif
6299 #endif
6300 #endif
6301 #endif
6302 break;
6303 } // up
6304 return digits;
6305 } // decGetDigits
6306
6307
6308 #if DECTRACE | DECCHECK
6309 /* ------------------------------------------------------------------ */
6310 /* decNumberShow -- display a number [debug aid] */
6311 /* dn is the number to show */
6312 /* */
6313 /* Shows: sign, exponent, coefficient (msu first), digits */
6314 /* or: sign, special-value */
6315 /* ------------------------------------------------------------------ */
6316 // this is public so other modules can use it
6317 void decNumberShow(const decNumber *dn) {
6318 const Unit *up; // work
6319 uInt u, d; // ..
6320 Int cut; // ..
6321 char isign='+'; // main sign
6322 if (dn==NULL) {
6323 printf("NULL\n");
6324 return;}
6325 if (decNumberIsNegative(dn)) isign='-';
6326 printf(" >> %c ", isign);
6327 if (dn->bits&DECSPECIAL) { // Is a special value
6328 if (decNumberIsInfinite(dn)) printf("Infinity");
6329 else { // a NaN
6330 if (dn->bits&DECSNAN) printf("sNaN"); // signalling NaN
6331 else printf("NaN");
6332 }
6333 // if coefficient and exponent are 0, no more to do
6334 if (dn->exponent==0 && dn->digits==1 && *dn->lsu==0) {
6335 printf("\n");
6336 return;}
6337 // drop through to report other information
6338 printf(" ");
6339 }
6340
6341 // now carefully display the coefficient
6342 up=dn->lsu+D2U(dn->digits)-1; // msu
6343 printf("%d", *up);
6344 for (up=up-1; up>=dn->lsu; up--) {
6345 u=*up;
6346 printf(":");
6347 for (cut=DECDPUN-1; cut>=0; cut--) {
6348 d=u/powers[cut];
6349 u-=d*powers[cut];
6350 printf("%d", d);
6351 } // cut
6352 } // up
6353 if (dn->exponent!=0) {
6354 char esign='+';
6355 if (dn->exponent<0) esign='-';
6356 printf(" E%c%d", esign, abs(dn->exponent));
6357 }
6358 printf(" [%d]\n", dn->digits);
6359 } // decNumberShow
6360 #endif
6361
6362 #if DECTRACE || DECCHECK
6363 /* ------------------------------------------------------------------ */
6364 /* decDumpAr -- display a unit array [debug aid] */
6365 /* name is a single-character tag name */
6366 /* ar is the array to display */
6367 /* len is the length of the array in Units */
6368 /* ------------------------------------------------------------------ */
6369 static void decDumpAr(char name, const Unit *ar, Int len) {
6370 Int i;
6371 #if DECDPUN==9
6372 char *spec="%09d ";
6373 #elif DECDPUN==8
6374 char *spec="%08d ";
6375 #elif DECDPUN==7
6376 char *spec="%07d ";
6377 #elif DECDPUN==6
6378 char *spec="%06d ";
6379 #elif DECDPUN==5
6380 char *spec="%05d ";
6381 #elif DECDPUN==4
6382 char *spec="%04d ";
6383 #elif DECDPUN==3
6384 char *spec="%03d ";
6385 #elif DECDPUN==2
6386 char *spec="%02d ";
6387 #else
6388 char *spec="%d ";
6389 #endif
6390 printf(" :%c: ", name);
6391 for (i=len-1; i>=0; i--) {
6392 if (i==len-1) printf("%d ", ar[i]);
6393 else printf(spec, ar[i]);
6394 }
6395 printf("\n");
6396 return;}
6397 #endif
6398
6399 #if DECCHECK
6400 /* ------------------------------------------------------------------ */
6401 /* decCheckOperands -- check operand(s) to a routine */
6402 /* res is the result structure (not checked; it will be set to */
6403 /* quiet NaN if error found (and it is not NULL)) */
6404 /* lhs is the first operand (may be DECUNUSED) */
6405 /* rhs is the second (may be DECUNUSED) */
6406 /* set is the context (may be DECUNUSED) */
6407 /* returns 0 if both operands, and the context are clean, or 1 */
6408 /* otherwise (in which case the context will show an error, */
6409 /* unless NULL). Note that res is not cleaned; caller should */
6410 /* handle this so res=NULL case is safe. */
6411 /* The caller is expected to abandon immediately if 1 is returned. */
6412 /* ------------------------------------------------------------------ */
6413 static Flag decCheckOperands(decNumber *res, const decNumber *lhs,
6414 const decNumber *rhs, decContext *set) {
6415 Flag bad=0;
6416 if (set==NULL) { // oops; hopeless
6417 #if DECTRACE
6418 printf("Context is NULL.\n");
6419 #endif
6420 bad=1;
6421 return 1;}
6422 else if (set!=DECUNUSED
6423 && (set->digits<1 || set->round<0 || set->round>=DEC_ROUND_MAX)) {
6424 bad=1;
6425 #if DECTRACE
6426 printf("Bad context [digits=%d round=%d].\n", set->digits, set->round);
6427 #endif
6428 }
6429 else {
6430 if (res==NULL) {
6431 bad=1;
6432 #if DECTRACE
6433 printf("Bad result [is NULL].\n");
6434 #endif
6435 }
6436 if (!bad && lhs!=DECUNUSED) bad=(decCheckNumber(lhs, set));
6437 if (!bad && rhs!=DECUNUSED) bad=(decCheckNumber(rhs, set));
6438 }
6439 if (bad) {
6440 if (set!=DECUNUSED) decContextSetStatus(set, DEC_Invalid_operation);
6441 if (res!=DECUNUSED && res!=NULL) {
6442 decNumberZero(res);
6443 res->bits=DECNAN; // qNaN
6444 }
6445 }
6446 return bad;
6447 } // decCheckOperands
6448
6449 /* ------------------------------------------------------------------ */
6450 /* decCheckNumber -- check a number */
6451 /* dn is the number to check */
6452 /* set is the context (may be DECUNUSED) */
6453 /* returns 0 if the number is clean, or 1 otherwise */
6454 /* */
6455 /* The number is considered valid if it could be a result from some */
6456 /* operation in some valid context (not necessarily the current one). */
6457 /* ------------------------------------------------------------------ */
6458 static Flag decCheckNumber(const decNumber *dn, decContext *set) {
6459 const Unit *up; // work
6460 uInt maxuint; // ..
6461 Int ae, d, digits; // ..
6462 Int emin, emax; // ..
6463
6464 if (dn==NULL) { // hopeless
6465 #if DECTRACE
6466 printf("Reference to decNumber is NULL.\n");
6467 #endif
6468 return 1;}
6469
6470 // check special values
6471 if (dn->bits & DECSPECIAL) {
6472 if (dn->exponent!=0) {
6473 #if DECTRACE
6474 printf("Exponent %d (not 0) for a special value.\n", dn->exponent);
6475 #endif
6476 return 1;}
6477
6478 // 2003.09.08: NaNs may now have coefficients, so next tests Inf only
6479 if (decNumberIsInfinite(dn)) {
6480 if (dn->digits!=1) {
6481 #if DECTRACE
6482 printf("Digits %d (not 1) for an infinity.\n", dn->digits);
6483 #endif
6484 return 1;}
6485 if (*dn->lsu!=0) {
6486 #if DECTRACE
6487 printf("LSU %d (not 0) for an infinity.\n", *dn->lsu);
6488 #endif
6489 return 1;}
6490 } // Inf
6491 // 2002.12.26: negative NaNs can now appear through proposed IEEE
6492 // concrete formats (decimal64, etc.), though they are
6493 // never visible in strings.
6494 return 0;
6495
6496 // if ((dn->bits & DECINF) || (dn->bits & DECNEG)==0) return 0;
6497 // #if DECTRACE
6498 // printf("Negative NaN in number.\n");
6499 // #endif
6500 // return 1;
6501 }
6502
6503 // check the coefficient
6504 if (dn->digits<1 || dn->digits>DECNUMMAXP) {
6505 #if DECTRACE
6506 printf("Digits %d in number.\n", dn->digits);
6507 #endif
6508 return 1;}
6509
6510 d=dn->digits;
6511
6512 for (up=dn->lsu; d>0; up++) {
6513 if (d>DECDPUN) maxuint=DECDPUNMAX;
6514 else { // reached the msu
6515 maxuint=powers[d]-1;
6516 if (dn->digits>1 && *up<powers[d-1]) {
6517 #if DECTRACE
6518 printf("Leading 0 in number.\n");
6519 decNumberShow(dn);
6520 #endif
6521 return 1;}
6522 }
6523 if (*up>maxuint) {
6524 #if DECTRACE
6525 printf("Bad Unit [%08x] in %d-digit number at offset %d [maxuint %d].\n",
6526 *up, dn->digits, up-dn->lsu, maxuint);
6527 #endif
6528 return 1;}
6529 d-=DECDPUN;
6530 }
6531
6532 // check the exponent. Note that input operands can have exponents
6533 // which are out of the set->emin/set->emax and set->digits range
6534 // (just as they can have more digits than set->digits).
6535 ae=dn->exponent+dn->digits-1; // adjusted exponent
6536 emax=DECNUMMAXE;
6537 emin=DECNUMMINE;
6538 digits=DECNUMMAXP;
6539 if (ae<emin-(digits-1)) {
6540 #if DECTRACE
6541 printf("Adjusted exponent underflow [%d].\n", ae);
6542 decNumberShow(dn);
6543 #endif
6544 return 1;}
6545 if (ae>+emax) {
6546 #if DECTRACE
6547 printf("Adjusted exponent overflow [%d].\n", ae);
6548 decNumberShow(dn);
6549 #endif
6550 return 1;}
6551
6552 return 0; // it's OK
6553 } // decCheckNumber
6554
6555
6556 /* ------------------------------------------------------------------ */
6557 /* decCheckInexact -- check a normal finite inexact result has digits */
6558 /* dn is the number to check */
6559 /* set is the context (for status and precision) */
6560 /* sets Invalid operation, etc., if some digits are missing */
6561 /* [this check is not made for DECSUBSET compilation] */
6562 /* ------------------------------------------------------------------ */
6563 static void decCheckInexact(const decNumber *dn, decContext *set) {
6564 #if !DECSUBSET
6565 if ((set->status & (DEC_Inexact|DEC_Subnormal))==DEC_Inexact
6566 && (set->digits!=dn->digits) && !(dn->bits & DECSPECIAL)) {
6567 decContextSetStatus(set, DEC_Invalid_operation);
6568 #if DECTRACE
6569 printf("Insufficient digits [%d] on normal Inexact result.\n", dn->digits);
6570 decNumberShow(dn);
6571 #endif
6572 }
6573 #endif
6574 return;
6575 } // decCheckInexact
6576 #endif
6577
6578 #if DECALLOC
6579 #undef malloc
6580 #undef free
6581 /* ------------------------------------------------------------------ */
6582 /* decMalloc -- accountable allocation routine */
6583 /* n is the number of bytes to allocate */
6584 /* */
6585 /* Semantics is the same as the stdlib malloc routine, but bytes */
6586 /* allocated are accounted for globally, and corruption fences are */
6587 /* added before and after the 'actual' storage. */
6588 /* ------------------------------------------------------------------ */
6589 /* This routine allocates storage with an extra twelve bytes; 8 are */
6590 /* at the start and hold: */
6591 /* 0-3 the original length requested */
6592 /* 4-7 buffer corruption detection fence (DECFENCE, x4) */
6593 /* The 4 bytes at the end also hold a corruption fence (DECFENCE, x4) */
6594 /* ------------------------------------------------------------------ */
6595 static void *decMalloc(size_t n) {
6596 uInt size=n+12; // true size
6597 void *alloc; // -> allocated storage
6598 uInt *j; // work
6599 uByte *b, *b0; // ..
6600
6601 alloc=malloc(size); // -> allocated storage
6602 if (alloc==NULL) return NULL; // out of strorage
6603 b0=(uByte *)alloc; // as bytes
6604 decAllocBytes+=n; // account for storage
6605 j=(uInt *)alloc; // -> first four bytes
6606 *j=n; // save n
6607 // printf(" allo ++ dAB: %d (%d)\n", decAllocBytes, n);
6608 for (b=b0+4; b<b0+8; b++) *b=DECFENCE;
6609 for (b=b0+n+8; b<b0+n+12; b++) *b=DECFENCE;
6610 return b0+8; // -> play area
6611 } // decMalloc
6612
6613 /* ------------------------------------------------------------------ */
6614 /* decFree -- accountable free routine */
6615 /* alloc is the storage to free */
6616 /* */
6617 /* Semantics is the same as the stdlib malloc routine, except that */
6618 /* the global storage accounting is updated and the fences are */
6619 /* checked to ensure that no routine has written 'out of bounds'. */
6620 /* ------------------------------------------------------------------ */
6621 /* This routine first checks that the fences have not been corrupted. */
6622 /* It then frees the storage using the 'truw' storage address (that */
6623 /* is, offset by 8). */
6624 /* ------------------------------------------------------------------ */
6625 static void decFree(void *alloc) {
6626 uInt *j, n; // pointer, original length
6627 uByte *b, *b0; // work
6628
6629 if (alloc==NULL) return; // allowed; it's a nop
6630 b0=(uByte *)alloc; // as bytes
6631 b0-=8; // -> true start of storage
6632 j=(uInt *)b0; // -> first four bytes
6633 n=*j; // lift
6634 for (b=b0+4; b<b0+8; b++) if (*b!=DECFENCE)
6635 printf("=== Corrupt byte [%02x] at offset %d from %d ===\n", *b,
6636 b-b0-8, (Int)b0);
6637 for (b=b0+n+8; b<b0+n+12; b++) if (*b!=DECFENCE)
6638 printf("=== Corrupt byte [%02x] at offset +%d from %d, n=%d ===\n", *b,
6639 b-b0-8, (Int)b0, n);
6640 free(b0); // drop the storage
6641 decAllocBytes-=n; // account for storage
6642 // printf(" free -- dAB: %d (%d)\n", decAllocBytes, -n);
6643 } // decFree
6644 #define malloc(a) decMalloc(a)
6645 #define free(a) decFree(a)
6646 #endif
6647