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