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