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