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