w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

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