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