| /* Decimal number arithmetic module for the decNumber C Library. |
| Copyright (C) 2005-2014 Free Software Foundation, Inc. |
| Contributed by IBM Corporation. Author Mike Cowlishaw. |
| |
| This file is part of GCC. |
| |
| GCC is free software; you can redistribute it and/or modify it under |
| the terms of the GNU General Public License as published by the Free |
| Software Foundation; either version 3, or (at your option) any later |
| version. |
| |
| GCC is distributed in the hope that it will be useful, but WITHOUT ANY |
| WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| for more details. |
| |
| Under Section 7 of GPL version 3, you are granted additional |
| permissions described in the GCC Runtime Library Exception, version |
| 3.1, as published by the Free Software Foundation. |
| |
| You should have received a copy of the GNU General Public License and |
| a copy of the GCC Runtime Library Exception along with this program; |
| see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| <http://www.gnu.org/licenses/>. */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* Decimal Number arithmetic module */ |
| /* ------------------------------------------------------------------ */ |
| /* This module comprises the routines for arbitrary-precision General */ |
| /* Decimal Arithmetic as defined in the specification which may be */ |
| /* found on the General Decimal Arithmetic pages. It implements both */ |
| /* the full ('extended') arithmetic and the simpler ('subset') */ |
| /* arithmetic. */ |
| /* */ |
| /* Usage notes: */ |
| /* */ |
| /* 1. This code is ANSI C89 except: */ |
| /* */ |
| /* a) C99 line comments (double forward slash) are used. (Most C */ |
| /* compilers accept these. If yours does not, a simple script */ |
| /* can be used to convert them to ANSI C comments.) */ |
| /* */ |
| /* b) Types from C99 stdint.h are used. If you do not have this */ |
| /* header file, see the User's Guide section of the decNumber */ |
| /* documentation; this lists the necessary definitions. */ |
| /* */ |
| /* c) If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and */ |
| /* uint64_t types may be used. To avoid these, set DECUSE64=0 */ |
| /* and DECDPUN<=4 (see documentation). */ |
| /* */ |
| /* The code also conforms to C99 restrictions; in particular, */ |
| /* strict aliasing rules are observed. */ |
| /* */ |
| /* 2. The decNumber format which this library uses is optimized for */ |
| /* efficient processing of relatively short numbers; in particular */ |
| /* it allows the use of fixed sized structures and minimizes copy */ |
| /* and move operations. It does, however, support arbitrary */ |
| /* precision (up to 999,999,999 digits) and arbitrary exponent */ |
| /* range (Emax in the range 0 through 999,999,999 and Emin in the */ |
| /* range -999,999,999 through 0). Mathematical functions (for */ |
| /* example decNumberExp) as identified below are restricted more */ |
| /* tightly: digits, emax, and -emin in the context must be <= */ |
| /* DEC_MAX_MATH (999999), and their operand(s) must be within */ |
| /* these bounds. */ |
| /* */ |
| /* 3. Logical functions are further restricted; their operands must */ |
| /* be finite, positive, have an exponent of zero, and all digits */ |
| /* must be either 0 or 1. The result will only contain digits */ |
| /* which are 0 or 1 (and will have exponent=0 and a sign of 0). */ |
| /* */ |
| /* 4. Operands to operator functions are never modified unless they */ |
| /* are also specified to be the result number (which is always */ |
| /* permitted). Other than that case, operands must not overlap. */ |
| /* */ |
| /* 5. Error handling: the type of the error is ORed into the status */ |
| /* flags in the current context (decContext structure). The */ |
| /* SIGFPE signal is then raised if the corresponding trap-enabler */ |
| /* flag in the decContext is set (is 1). */ |
| /* */ |
| /* It is the responsibility of the caller to clear the status */ |
| /* flags as required. */ |
| /* */ |
| /* The result of any routine which returns a number will always */ |
| /* be a valid number (which may be a special value, such as an */ |
| /* Infinity or NaN). */ |
| /* */ |
| /* 6. The decNumber format is not an exchangeable concrete */ |
| /* representation as it comprises fields which may be machine- */ |
| /* dependent (packed or unpacked, or special length, for example). */ |
| /* Canonical conversions to and from strings are provided; other */ |
| /* conversions are available in separate modules. */ |
| /* */ |
| /* 7. Normally, input operands are assumed to be valid. Set DECCHECK */ |
| /* to 1 for extended operand checking (including NULL operands). */ |
| /* Results are undefined if a badly-formed structure (or a NULL */ |
| /* pointer to a structure) is provided, though with DECCHECK */ |
| /* enabled the operator routines are protected against exceptions. */ |
| /* (Except if the result pointer is NULL, which is unrecoverable.) */ |
| /* */ |
| /* However, the routines will never cause exceptions if they are */ |
| /* given well-formed operands, even if the value of the operands */ |
| /* is inappropriate for the operation and DECCHECK is not set. */ |
| /* (Except for SIGFPE, as and where documented.) */ |
| /* */ |
| /* 8. Subset arithmetic is available only if DECSUBSET is set to 1. */ |
| /* ------------------------------------------------------------------ */ |
| /* Implementation notes for maintenance of this module: */ |
| /* */ |
| /* 1. Storage leak protection: Routines which use malloc are not */ |
| /* permitted to use return for fastpath or error exits (i.e., */ |
| /* they follow strict structured programming conventions). */ |
| /* Instead they have a do{}while(0); construct surrounding the */ |
| /* code which is protected -- break may be used to exit this. */ |
| /* Other routines can safely use the return statement inline. */ |
| /* */ |
| /* Storage leak accounting can be enabled using DECALLOC. */ |
| /* */ |
| /* 2. All loops use the for(;;) construct. Any do construct does */ |
| /* not loop; it is for allocation protection as just described. */ |
| /* */ |
| /* 3. Setting status in the context must always be the very last */ |
| /* action in a routine, as non-0 status may raise a trap and hence */ |
| /* the call to set status may not return (if the handler uses long */ |
| /* jump). Therefore all cleanup must be done first. In general, */ |
| /* to achieve this status is accumulated and is only applied just */ |
| /* before return by calling decContextSetStatus (via decStatus). */ |
| /* */ |
| /* Routines which allocate storage cannot, in general, use the */ |
| /* 'top level' routines which could cause a non-returning */ |
| /* transfer of control. The decXxxxOp routines are safe (do not */ |
| /* call decStatus even if traps are set in the context) and should */ |
| /* be used instead (they are also a little faster). */ |
| /* */ |
| /* 4. Exponent checking is minimized by allowing the exponent to */ |
| /* grow outside its limits during calculations, provided that */ |
| /* the decFinalize function is called later. Multiplication and */ |
| /* division, and intermediate calculations in exponentiation, */ |
| /* require more careful checks because of the risk of 31-bit */ |
| /* overflow (the most negative valid exponent is -1999999997, for */ |
| /* a 999999999-digit number with adjusted exponent of -999999999). */ |
| /* */ |
| /* 5. Rounding is deferred until finalization of results, with any */ |
| /* 'off to the right' data being represented as a single digit */ |
| /* residue (in the range -1 through 9). This avoids any double- */ |
| /* rounding when more than one shortening takes place (for */ |
| /* example, when a result is subnormal). */ |
| /* */ |
| /* 6. The digits count is allowed to rise to a multiple of DECDPUN */ |
| /* during many operations, so whole Units are handled and exact */ |
| /* accounting of digits is not needed. The correct digits value */ |
| /* is found by decGetDigits, which accounts for leading zeros. */ |
| /* This must be called before any rounding if the number of digits */ |
| /* is not known exactly. */ |
| /* */ |
| /* 7. The multiply-by-reciprocal 'trick' is used for partitioning */ |
| /* numbers up to four digits, using appropriate constants. This */ |
| /* is not useful for longer numbers because overflow of 32 bits */ |
| /* would lead to 4 multiplies, which is almost as expensive as */ |
| /* a divide (unless a floating-point or 64-bit multiply is */ |
| /* assumed to be available). */ |
| /* */ |
| /* 8. Unusual abbreviations that may be used in the commentary: */ |
| /* lhs -- left hand side (operand, of an operation) */ |
| /* lsd -- least significant digit (of coefficient) */ |
| /* lsu -- least significant Unit (of coefficient) */ |
| /* msd -- most significant digit (of coefficient) */ |
| /* msi -- most significant item (in an array) */ |
| /* msu -- most significant Unit (of coefficient) */ |
| /* rhs -- right hand side (operand, of an operation) */ |
| /* +ve -- positive */ |
| /* -ve -- negative */ |
| /* ** -- raise to the power */ |
| /* ------------------------------------------------------------------ */ |
| |
| #include <stdlib.h> /* for malloc, free, etc. */ |
| #include <stdio.h> /* for printf [if needed] */ |
| #include <string.h> /* for strcpy */ |
| #include <ctype.h> /* for lower */ |
| #include "dconfig.h" /* for GCC definitions */ |
| #include "decNumber.h" /* base number library */ |
| #include "decNumberLocal.h" /* decNumber local types, etc. */ |
| |
| /* Constants */ |
| /* Public lookup table used by the D2U macro */ |
| const uByte d2utable[DECMAXD2U+1]=D2UTABLE; |
| |
| #define DECVERB 1 /* set to 1 for verbose DECCHECK */ |
| #define powers DECPOWERS /* old internal name */ |
| |
| /* Local constants */ |
| #define DIVIDE 0x80 /* Divide operators */ |
| #define REMAINDER 0x40 /* .. */ |
| #define DIVIDEINT 0x20 /* .. */ |
| #define REMNEAR 0x10 /* .. */ |
| #define COMPARE 0x01 /* Compare operators */ |
| #define COMPMAX 0x02 /* .. */ |
| #define COMPMIN 0x03 /* .. */ |
| #define COMPTOTAL 0x04 /* .. */ |
| #define COMPNAN 0x05 /* .. [NaN processing] */ |
| #define COMPSIG 0x06 /* .. [signaling COMPARE] */ |
| #define COMPMAXMAG 0x07 /* .. */ |
| #define COMPMINMAG 0x08 /* .. */ |
| |
| #define DEC_sNaN 0x40000000 /* local status: sNaN signal */ |
| #define BADINT (Int)0x80000000 /* most-negative Int; error indicator */ |
| /* Next two indicate an integer >= 10**6, and its parity (bottom bit) */ |
| #define BIGEVEN (Int)0x80000002 |
| #define BIGODD (Int)0x80000003 |
| |
| static Unit uarrone[1]={1}; /* Unit array of 1, used for incrementing */ |
| |
| /* Granularity-dependent code */ |
| #if DECDPUN<=4 |
| #define eInt Int /* extended integer */ |
| #define ueInt uInt /* unsigned extended integer */ |
| /* Constant multipliers for divide-by-power-of five using reciprocal */ |
| /* multiply, after removing powers of 2 by shifting, and final shift */ |
| /* of 17 [we only need up to **4] */ |
| static const uInt multies[]={131073, 26215, 5243, 1049, 210}; |
| /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */ |
| #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17) |
| #else |
| /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */ |
| #if !DECUSE64 |
| #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4 |
| #endif |
| #define eInt Long /* extended integer */ |
| #define ueInt uLong /* unsigned extended integer */ |
| #endif |
| |
| /* Local routines */ |
| static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *, |
| decContext *, uByte, uInt *); |
| static Flag decBiStr(const char *, const char *, const char *); |
| static uInt decCheckMath(const decNumber *, decContext *, uInt *); |
| static void decApplyRound(decNumber *, decContext *, Int, uInt *); |
| static Int decCompare(const decNumber *lhs, const decNumber *rhs, Flag); |
| static decNumber * decCompareOp(decNumber *, const decNumber *, |
| const decNumber *, decContext *, |
| Flag, uInt *); |
| static void decCopyFit(decNumber *, const decNumber *, decContext *, |
| Int *, uInt *); |
| static decNumber * decDecap(decNumber *, Int); |
| static decNumber * decDivideOp(decNumber *, const decNumber *, |
| const decNumber *, decContext *, Flag, uInt *); |
| static decNumber * decExpOp(decNumber *, const decNumber *, |
| decContext *, uInt *); |
| static void decFinalize(decNumber *, decContext *, Int *, uInt *); |
| static Int decGetDigits(Unit *, Int); |
| static Int decGetInt(const decNumber *); |
| static decNumber * decLnOp(decNumber *, const decNumber *, |
| decContext *, uInt *); |
| static decNumber * decMultiplyOp(decNumber *, const decNumber *, |
| const decNumber *, decContext *, |
| uInt *); |
| static decNumber * decNaNs(decNumber *, const decNumber *, |
| const decNumber *, decContext *, uInt *); |
| static decNumber * decQuantizeOp(decNumber *, const decNumber *, |
| const decNumber *, decContext *, Flag, |
| uInt *); |
| static void decReverse(Unit *, Unit *); |
| static void decSetCoeff(decNumber *, decContext *, const Unit *, |
| Int, Int *, uInt *); |
| static void decSetMaxValue(decNumber *, decContext *); |
| static void decSetOverflow(decNumber *, decContext *, uInt *); |
| static void decSetSubnormal(decNumber *, decContext *, Int *, uInt *); |
| static Int decShiftToLeast(Unit *, Int, Int); |
| static Int decShiftToMost(Unit *, Int, Int); |
| static void decStatus(decNumber *, uInt, decContext *); |
| static void decToString(const decNumber *, char[], Flag); |
| static decNumber * decTrim(decNumber *, decContext *, Flag, Flag, Int *); |
| static Int decUnitAddSub(const Unit *, Int, const Unit *, Int, Int, |
| Unit *, Int); |
| static Int decUnitCompare(const Unit *, Int, const Unit *, Int, Int); |
| |
| #if !DECSUBSET |
| /* decFinish == decFinalize when no subset arithmetic needed */ |
| #define decFinish(a,b,c,d) decFinalize(a,b,c,d) |
| #else |
| static void decFinish(decNumber *, decContext *, Int *, uInt *); |
| static decNumber * decRoundOperand(const decNumber *, decContext *, uInt *); |
| #endif |
| |
| /* Local macros */ |
| /* masked special-values bits */ |
| #define SPECIALARG (rhs->bits & DECSPECIAL) |
| #define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL) |
| |
| /* Diagnostic macros, etc. */ |
| #if DECALLOC |
| /* Handle malloc/free accounting. If enabled, our accountable routines */ |
| /* are used; otherwise the code just goes straight to the system malloc */ |
| /* and free routines. */ |
| #define malloc(a) decMalloc(a) |
| #define free(a) decFree(a) |
| #define DECFENCE 0x5a /* corruption detector */ |
| /* 'Our' malloc and free: */ |
| static void *decMalloc(size_t); |
| static void decFree(void *); |
| uInt decAllocBytes=0; /* count of bytes allocated */ |
| /* Note that DECALLOC code only checks for storage buffer overflow. */ |
| /* To check for memory leaks, the decAllocBytes variable must be */ |
| /* checked to be 0 at appropriate times (e.g., after the test */ |
| /* harness completes a set of tests). This checking may be unreliable */ |
| /* if the testing is done in a multi-thread environment. */ |
| #endif |
| |
| #if DECCHECK |
| /* Optional checking routines. Enabling these means that decNumber */ |
| /* and decContext operands to operator routines are checked for */ |
| /* correctness. This roughly doubles the execution time of the */ |
| /* fastest routines (and adds 600+ bytes), so should not normally be */ |
| /* used in 'production'. */ |
| /* decCheckInexact is used to check that inexact results have a full */ |
| /* complement of digits (where appropriate -- this is not the case */ |
| /* for Quantize, for example) */ |
| #define DECUNRESU ((decNumber *)(void *)0xffffffff) |
| #define DECUNUSED ((const decNumber *)(void *)0xffffffff) |
| #define DECUNCONT ((decContext *)(void *)(0xffffffff)) |
| static Flag decCheckOperands(decNumber *, const decNumber *, |
| const decNumber *, decContext *); |
| static Flag decCheckNumber(const decNumber *); |
| static void decCheckInexact(const decNumber *, decContext *); |
| #endif |
| |
| #if DECTRACE || DECCHECK |
| /* Optional trace/debugging routines (may or may not be used) */ |
| void decNumberShow(const decNumber *); /* displays the components of a number */ |
| static void decDumpAr(char, const Unit *, Int); |
| #endif |
| |
| /* ================================================================== */ |
| /* Conversions */ |
| /* ================================================================== */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* from-int32 -- conversion from Int or uInt */ |
| /* */ |
| /* dn is the decNumber to receive the integer */ |
| /* in or uin is the integer to be converted */ |
| /* returns dn */ |
| /* */ |
| /* No error is possible. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberFromInt32(decNumber *dn, Int in) { |
| uInt unsig; |
| if (in>=0) unsig=in; |
| else { /* negative (possibly BADINT) */ |
| if (in==BADINT) unsig=(uInt)1073741824*2; /* special case */ |
| else unsig=-in; /* invert */ |
| } |
| /* in is now positive */ |
| decNumberFromUInt32(dn, unsig); |
| if (in<0) dn->bits=DECNEG; /* sign needed */ |
| return dn; |
| } /* decNumberFromInt32 */ |
| |
| decNumber * decNumberFromUInt32(decNumber *dn, uInt uin) { |
| Unit *up; /* work pointer */ |
| decNumberZero(dn); /* clean */ |
| if (uin==0) return dn; /* [or decGetDigits bad call] */ |
| for (up=dn->lsu; uin>0; up++) { |
| *up=(Unit)(uin%(DECDPUNMAX+1)); |
| uin=uin/(DECDPUNMAX+1); |
| } |
| dn->digits=decGetDigits(dn->lsu, up-dn->lsu); |
| return dn; |
| } /* decNumberFromUInt32 */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* to-int32 -- conversion to Int or uInt */ |
| /* */ |
| /* dn is the decNumber to convert */ |
| /* set is the context for reporting errors */ |
| /* returns the converted decNumber, or 0 if Invalid is set */ |
| /* */ |
| /* Invalid is set if the decNumber does not have exponent==0 or if */ |
| /* it is a NaN, Infinite, or out-of-range. */ |
| /* ------------------------------------------------------------------ */ |
| Int decNumberToInt32(const decNumber *dn, decContext *set) { |
| #if DECCHECK |
| if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0; |
| #endif |
| |
| /* special or too many digits, or bad exponent */ |
| if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0) ; /* bad */ |
| else { /* is a finite integer with 10 or fewer digits */ |
| Int d; /* work */ |
| const Unit *up; /* .. */ |
| uInt hi=0, lo; /* .. */ |
| up=dn->lsu; /* -> lsu */ |
| lo=*up; /* get 1 to 9 digits */ |
| #if DECDPUN>1 /* split to higher */ |
| hi=lo/10; |
| lo=lo%10; |
| #endif |
| up++; |
| /* collect remaining Units, if any, into hi */ |
| for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1]; |
| /* now low has the lsd, hi the remainder */ |
| if (hi>214748364 || (hi==214748364 && lo>7)) { /* out of range? */ |
| /* most-negative is a reprieve */ |
| if (dn->bits&DECNEG && hi==214748364 && lo==8) return 0x80000000; |
| /* bad -- drop through */ |
| } |
| else { /* in-range always */ |
| Int i=X10(hi)+lo; |
| if (dn->bits&DECNEG) return -i; |
| return i; |
| } |
| } /* integer */ |
| decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */ |
| return 0; |
| } /* decNumberToInt32 */ |
| |
| uInt decNumberToUInt32(const decNumber *dn, decContext *set) { |
| #if DECCHECK |
| if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0; |
| #endif |
| /* special or too many digits, or bad exponent, or negative (<0) */ |
| if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0 |
| || (dn->bits&DECNEG && !ISZERO(dn))); /* bad */ |
| else { /* is a finite integer with 10 or fewer digits */ |
| Int d; /* work */ |
| const Unit *up; /* .. */ |
| uInt hi=0, lo; /* .. */ |
| up=dn->lsu; /* -> lsu */ |
| lo=*up; /* get 1 to 9 digits */ |
| #if DECDPUN>1 /* split to higher */ |
| hi=lo/10; |
| lo=lo%10; |
| #endif |
| up++; |
| /* collect remaining Units, if any, into hi */ |
| for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1]; |
| |
| /* now low has the lsd, hi the remainder */ |
| if (hi>429496729 || (hi==429496729 && lo>5)) ; /* no reprieve possible */ |
| else return X10(hi)+lo; |
| } /* integer */ |
| decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */ |
| return 0; |
| } /* decNumberToUInt32 */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* to-scientific-string -- conversion to numeric string */ |
| /* to-engineering-string -- conversion to numeric string */ |
| /* */ |
| /* decNumberToString(dn, string); */ |
| /* decNumberToEngString(dn, string); */ |
| /* */ |
| /* dn is the decNumber to convert */ |
| /* string is the string where the result will be laid out */ |
| /* */ |
| /* string must be at least dn->digits+14 characters long */ |
| /* */ |
| /* No error is possible, and no status can be set. */ |
| /* ------------------------------------------------------------------ */ |
| char * decNumberToString(const decNumber *dn, char *string){ |
| decToString(dn, string, 0); |
| return string; |
| } /* DecNumberToString */ |
| |
| char * decNumberToEngString(const decNumber *dn, char *string){ |
| decToString(dn, string, 1); |
| return string; |
| } /* DecNumberToEngString */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* to-number -- conversion from numeric string */ |
| /* */ |
| /* decNumberFromString -- convert string to decNumber */ |
| /* dn -- the number structure to fill */ |
| /* chars[] -- the string to convert ('\0' terminated) */ |
| /* set -- the context used for processing any error, */ |
| /* determining the maximum precision available */ |
| /* (set.digits), determining the maximum and minimum */ |
| /* exponent (set.emax and set.emin), determining if */ |
| /* extended values are allowed, and checking the */ |
| /* rounding mode if overflow occurs or rounding is */ |
| /* needed. */ |
| /* */ |
| /* The length of the coefficient and the size of the exponent are */ |
| /* checked by this routine, so the correct error (Underflow or */ |
| /* Overflow) can be reported or rounding applied, as necessary. */ |
| /* */ |
| /* If bad syntax is detected, the result will be a quiet NaN. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberFromString(decNumber *dn, const char chars[], |
| decContext *set) { |
| Int exponent=0; /* working exponent [assume 0] */ |
| uByte bits=0; /* working flags [assume +ve] */ |
| Unit *res; /* where result will be built */ |
| Unit resbuff[SD2U(DECBUFFER+9)];/* local buffer in case need temporary */ |
| /* [+9 allows for ln() constants] */ |
| Unit *allocres=NULL; /* -> allocated result, iff allocated */ |
| Int d=0; /* count of digits found in decimal part */ |
| const char *dotchar=NULL; /* where dot was found */ |
| const char *cfirst=chars; /* -> first character of decimal part */ |
| const char *last=NULL; /* -> last digit of decimal part */ |
| const char *c; /* work */ |
| Unit *up; /* .. */ |
| #if DECDPUN>1 |
| Int cut, out; /* .. */ |
| #endif |
| Int residue; /* rounding residue */ |
| uInt status=0; /* error code */ |
| |
| #if DECCHECK |
| if (decCheckOperands(DECUNRESU, DECUNUSED, DECUNUSED, set)) |
| return decNumberZero(dn); |
| #endif |
| |
| do { /* status & malloc protection */ |
| for (c=chars;; c++) { /* -> input character */ |
| if (*c>='0' && *c<='9') { /* test for Arabic digit */ |
| last=c; |
| d++; /* count of real digits */ |
| continue; /* still in decimal part */ |
| } |
| if (*c=='.' && dotchar==NULL) { /* first '.' */ |
| dotchar=c; /* record offset into decimal part */ |
| if (c==cfirst) cfirst++; /* first digit must follow */ |
| continue;} |
| if (c==chars) { /* first in string... */ |
| if (*c=='-') { /* valid - sign */ |
| cfirst++; |
| bits=DECNEG; |
| continue;} |
| if (*c=='+') { /* valid + sign */ |
| cfirst++; |
| continue;} |
| } |
| /* *c is not a digit, or a valid +, -, or '.' */ |
| break; |
| } /* c */ |
| |
| if (last==NULL) { /* no digits yet */ |
| status=DEC_Conversion_syntax;/* assume the worst */ |
| if (*c=='\0') break; /* and no more to come... */ |
| #if DECSUBSET |
| /* if subset then infinities and NaNs are not allowed */ |
| if (!set->extended) break; /* hopeless */ |
| #endif |
| /* Infinities and NaNs are possible, here */ |
| if (dotchar!=NULL) break; /* .. unless had a dot */ |
| decNumberZero(dn); /* be optimistic */ |
| if (decBiStr(c, "infinity", "INFINITY") |
| || decBiStr(c, "inf", "INF")) { |
| dn->bits=bits | DECINF; |
| status=0; /* is OK */ |
| break; /* all done */ |
| } |
| /* a NaN expected */ |
| /* 2003.09.10 NaNs are now permitted to have a sign */ |
| dn->bits=bits | DECNAN; /* assume simple NaN */ |
| if (*c=='s' || *c=='S') { /* looks like an sNaN */ |
| c++; |
| dn->bits=bits | DECSNAN; |
| } |
| if (*c!='n' && *c!='N') break; /* check caseless "NaN" */ |
| c++; |
| if (*c!='a' && *c!='A') break; /* .. */ |
| c++; |
| if (*c!='n' && *c!='N') break; /* .. */ |
| c++; |
| /* now either nothing, or nnnn payload, expected */ |
| /* -> start of integer and skip leading 0s [including plain 0] */ |
| for (cfirst=c; *cfirst=='0';) cfirst++; |
| if (*cfirst=='\0') { /* "NaN" or "sNaN", maybe with all 0s */ |
| status=0; /* it's good */ |
| break; /* .. */ |
| } |
| /* something other than 0s; setup last and d as usual [no dots] */ |
| for (c=cfirst;; c++, d++) { |
| if (*c<'0' || *c>'9') break; /* test for Arabic digit */ |
| last=c; |
| } |
| if (*c!='\0') break; /* not all digits */ |
| if (d>set->digits-1) { |
| /* [NB: payload in a decNumber can be full length unless */ |
| /* clamped, in which case can only be digits-1] */ |
| if (set->clamp) break; |
| if (d>set->digits) break; |
| } /* too many digits? */ |
| /* good; drop through to convert the integer to coefficient */ |
| status=0; /* syntax is OK */ |
| bits=dn->bits; /* for copy-back */ |
| } /* last==NULL */ |
| |
| else if (*c!='\0') { /* more to process... */ |
| /* had some digits; exponent is only valid sequence now */ |
| Flag nege; /* 1=negative exponent */ |
| const char *firstexp; /* -> first significant exponent digit */ |
| status=DEC_Conversion_syntax;/* assume the worst */ |
| if (*c!='e' && *c!='E') break; |
| /* Found 'e' or 'E' -- now process explicit exponent */ |
| /* 1998.07.11: sign no longer required */ |
| nege=0; |
| c++; /* to (possible) sign */ |
| if (*c=='-') {nege=1; c++;} |
| else if (*c=='+') c++; |
| if (*c=='\0') break; |
| |
| for (; *c=='0' && *(c+1)!='\0';) c++; /* strip insignificant zeros */ |
| firstexp=c; /* save exponent digit place */ |
| for (; ;c++) { |
| if (*c<'0' || *c>'9') break; /* not a digit */ |
| exponent=X10(exponent)+(Int)*c-(Int)'0'; |
| } /* c */ |
| /* if not now on a '\0', *c must not be a digit */ |
| if (*c!='\0') break; |
| |
| /* (this next test must be after the syntax checks) */ |
| /* if it was too long the exponent may have wrapped, so check */ |
| /* carefully and set it to a certain overflow if wrap possible */ |
| if (c>=firstexp+9+1) { |
| if (c>firstexp+9+1 || *firstexp>'1') exponent=DECNUMMAXE*2; |
| /* [up to 1999999999 is OK, for example 1E-1000000998] */ |
| } |
| if (nege) exponent=-exponent; /* was negative */ |
| status=0; /* is OK */ |
| } /* stuff after digits */ |
| |
| /* Here when whole string has been inspected; syntax is good */ |
| /* cfirst->first digit (never dot), last->last digit (ditto) */ |
| |
| /* strip leading zeros/dot [leave final 0 if all 0's] */ |
| if (*cfirst=='0') { /* [cfirst has stepped over .] */ |
| for (c=cfirst; c<last; c++, cfirst++) { |
| if (*c=='.') continue; /* ignore dots */ |
| if (*c!='0') break; /* non-zero found */ |
| d--; /* 0 stripped */ |
| } /* c */ |
| #if DECSUBSET |
| /* make a rapid exit for easy zeros if !extended */ |
| if (*cfirst=='0' && !set->extended) { |
| decNumberZero(dn); /* clean result */ |
| break; /* [could be return] */ |
| } |
| #endif |
| } /* at least one leading 0 */ |
| |
| /* Handle decimal point... */ |
| if (dotchar!=NULL && dotchar<last) /* non-trailing '.' found? */ |
| exponent-=(last-dotchar); /* adjust exponent */ |
| /* [we can now ignore the .] */ |
| |
| /* OK, the digits string is good. Assemble in the decNumber, or in */ |
| /* a temporary units array if rounding is needed */ |
| if (d<=set->digits) res=dn->lsu; /* fits into supplied decNumber */ |
| else { /* rounding needed */ |
| Int needbytes=D2U(d)*sizeof(Unit);/* bytes needed */ |
| res=resbuff; /* assume use local buffer */ |
| if (needbytes>(Int)sizeof(resbuff)) { /* too big for local */ |
| allocres=(Unit *)malloc(needbytes); |
| if (allocres==NULL) {status|=DEC_Insufficient_storage; break;} |
| res=allocres; |
| } |
| } |
| /* res now -> number lsu, buffer, or allocated storage for Unit array */ |
| |
| /* Place the coefficient into the selected Unit array */ |
| /* [this is often 70% of the cost of this function when DECDPUN>1] */ |
| #if DECDPUN>1 |
| out=0; /* accumulator */ |
| up=res+D2U(d)-1; /* -> msu */ |
| cut=d-(up-res)*DECDPUN; /* digits in top unit */ |
| for (c=cfirst;; c++) { /* along the digits */ |
| if (*c=='.') continue; /* ignore '.' [don't decrement cut] */ |
| out=X10(out)+(Int)*c-(Int)'0'; |
| if (c==last) break; /* done [never get to trailing '.'] */ |
| cut--; |
| if (cut>0) continue; /* more for this unit */ |
| *up=(Unit)out; /* write unit */ |
| up--; /* prepare for unit below.. */ |
| cut=DECDPUN; /* .. */ |
| out=0; /* .. */ |
| } /* c */ |
| *up=(Unit)out; /* write lsu */ |
| |
| #else |
| /* DECDPUN==1 */ |
| up=res; /* -> lsu */ |
| for (c=last; c>=cfirst; c--) { /* over each character, from least */ |
| if (*c=='.') continue; /* ignore . [don't step up] */ |
| *up=(Unit)((Int)*c-(Int)'0'); |
| up++; |
| } /* c */ |
| #endif |
| |
| dn->bits=bits; |
| dn->exponent=exponent; |
| dn->digits=d; |
| |
| /* if not in number (too long) shorten into the number */ |
| if (d>set->digits) { |
| residue=0; |
| decSetCoeff(dn, set, res, d, &residue, &status); |
| /* always check for overflow or subnormal and round as needed */ |
| decFinalize(dn, set, &residue, &status); |
| } |
| else { /* no rounding, but may still have overflow or subnormal */ |
| /* [these tests are just for performance; finalize repeats them] */ |
| if ((dn->exponent-1<set->emin-dn->digits) |
| || (dn->exponent-1>set->emax-set->digits)) { |
| residue=0; |
| decFinalize(dn, set, &residue, &status); |
| } |
| } |
| /* decNumberShow(dn); */ |
| } while(0); /* [for break] */ |
| |
| free(allocres); /* drop any storage used */ |
| if (status!=0) decStatus(dn, status, set); |
| return dn; |
| } /* decNumberFromString */ |
| |
| /* ================================================================== */ |
| /* Operators */ |
| /* ================================================================== */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberAbs -- absolute value operator */ |
| /* */ |
| /* This computes C = abs(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* See also decNumberCopyAbs for a quiet bitwise version of this. */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* This has the same effect as decNumberPlus unless A is negative, */ |
| /* in which case it has the same effect as decNumberMinus. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberAbs(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decNumber dzero; /* for 0 */ |
| uInt status=0; /* accumulator */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| decNumberZero(&dzero); /* set 0 */ |
| dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ |
| decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberAbs */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberAdd -- add two Numbers */ |
| /* */ |
| /* This computes C = A + B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X+X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* This just calls the routine shared with Subtract */ |
| decNumber * decNumberAdd(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decAddOp(res, lhs, rhs, set, 0, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberAdd */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberAnd -- AND two Numbers, digitwise */ |
| /* */ |
| /* This computes C = A & B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X&X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context (used for result length and error report) */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Logical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberAnd(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| const Unit *ua, *ub; /* -> operands */ |
| const Unit *msua, *msub; /* -> operand msus */ |
| Unit *uc, *msuc; /* -> result and its msu */ |
| Int msudigs; /* digits in res msu */ |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs) |
| || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| |
| /* operands are valid */ |
| ua=lhs->lsu; /* bottom-up */ |
| ub=rhs->lsu; /* .. */ |
| uc=res->lsu; /* .. */ |
| msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */ |
| msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */ |
| msuc=uc+D2U(set->digits)-1; /* -> msu of result */ |
| msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ |
| for (; uc<=msuc; ua++, ub++, uc++) { /* Unit loop */ |
| Unit a, b; /* extract units */ |
| if (ua>msua) a=0; |
| else a=*ua; |
| if (ub>msub) b=0; |
| else b=*ub; |
| *uc=0; /* can now write back */ |
| if (a|b) { /* maybe 1 bits to examine */ |
| Int i, j; |
| *uc=0; /* can now write back */ |
| /* This loop could be unrolled and/or use BIN2BCD tables */ |
| for (i=0; i<DECDPUN; i++) { |
| if (a&b&1) *uc=*uc+(Unit)powers[i]; /* effect AND */ |
| j=a%10; |
| a=a/10; |
| j|=b%10; |
| b=b/10; |
| if (j>1) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| if (uc==msuc && i==msudigs-1) break; /* just did final digit */ |
| } /* each digit */ |
| } /* both OK */ |
| } /* each unit */ |
| /* [here uc-1 is the msu of the result] */ |
| res->digits=decGetDigits(res->lsu, uc-res->lsu); |
| res->exponent=0; /* integer */ |
| res->bits=0; /* sign=0 */ |
| return res; /* [no status to set] */ |
| } /* decNumberAnd */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberCompare -- compare two Numbers */ |
| /* */ |
| /* This computes C = A ? B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for one digit (or NaN). */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberCompare(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPARE, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberCompare */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberCompareSignal -- compare, signalling on all NaNs */ |
| /* */ |
| /* This computes C = A ? B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for one digit (or NaN). */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberCompareSignal(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPSIG, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberCompareSignal */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberCompareTotal -- compare two Numbers, using total ordering */ |
| /* */ |
| /* This computes C = A ? B, under total ordering */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for one digit; the result will always be one of */ |
| /* -1, 0, or 1. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberCompareTotal */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberCompareTotalMag -- compare, total ordering of magnitudes */ |
| /* */ |
| /* This computes C = |A| ? |B|, under total ordering */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for one digit; the result will always be one of */ |
| /* -1, 0, or 1. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberCompareTotalMag(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| uInt needbytes; /* for space calculations */ |
| decNumber bufa[D2N(DECBUFFER+1)];/* +1 in case DECBUFFER=0 */ |
| decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ |
| decNumber bufb[D2N(DECBUFFER+1)]; |
| decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ |
| decNumber *a, *b; /* temporary pointers */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| do { /* protect allocated storage */ |
| /* if either is negative, take a copy and absolute */ |
| if (decNumberIsNegative(lhs)) { /* lhs<0 */ |
| a=bufa; |
| needbytes=sizeof(decNumber)+(D2U(lhs->digits)-1)*sizeof(Unit); |
| if (needbytes>sizeof(bufa)) { /* need malloc space */ |
| allocbufa=(decNumber *)malloc(needbytes); |
| if (allocbufa==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| a=allocbufa; /* use the allocated space */ |
| } |
| decNumberCopy(a, lhs); /* copy content */ |
| a->bits&=~DECNEG; /* .. and clear the sign */ |
| lhs=a; /* use copy from here on */ |
| } |
| if (decNumberIsNegative(rhs)) { /* rhs<0 */ |
| b=bufb; |
| needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit); |
| if (needbytes>sizeof(bufb)) { /* need malloc space */ |
| allocbufb=(decNumber *)malloc(needbytes); |
| if (allocbufb==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| b=allocbufb; /* use the allocated space */ |
| } |
| decNumberCopy(b, rhs); /* copy content */ |
| b->bits&=~DECNEG; /* .. and clear the sign */ |
| rhs=b; /* use copy from here on */ |
| } |
| decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status); |
| } while(0); /* end protected */ |
| |
| free(allocbufa); /* drop any storage used */ |
| free(allocbufb); /* .. */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberCompareTotalMag */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberDivide -- divide one number by another */ |
| /* */ |
| /* This computes C = A / B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X/X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberDivide(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decDivideOp(res, lhs, rhs, set, DIVIDE, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberDivide */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberDivideInteger -- divide and return integer quotient */ |
| /* */ |
| /* This computes C = A # B, where # is the integer divide operator */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X#X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decDivideOp(res, lhs, rhs, set, DIVIDEINT, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberDivideInteger */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberExp -- exponentiation */ |
| /* */ |
| /* This computes C = exp(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context; note that rounding mode has no effect */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Mathematical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* */ |
| /* Finite results will always be full precision and Inexact, except */ |
| /* when A is a zero or -Infinity (giving 1 or 0 respectively). */ |
| /* */ |
| /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */ |
| /* almost always be correctly rounded, but may be up to 1 ulp in */ |
| /* error in rare cases. */ |
| /* ------------------------------------------------------------------ */ |
| /* This is a wrapper for decExpOp which can handle the slightly wider */ |
| /* (double) range needed by Ln (which has to be able to calculate */ |
| /* exp(-a) where a can be the tiniest number (Ntiny). */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberExp(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| uInt status=0; /* accumulator */ |
| #if DECSUBSET |
| decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ |
| #endif |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* Check restrictions; these restrictions ensure that if h=8 (see */ |
| /* decExpOp) then the result will either overflow or underflow to 0. */ |
| /* Other math functions restrict the input range, too, for inverses. */ |
| /* If not violated then carry out the operation. */ |
| if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */ |
| #if DECSUBSET |
| if (!set->extended) { |
| /* reduce operand and set lostDigits status, as needed */ |
| if (rhs->digits>set->digits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| rhs=allocrhs; |
| } |
| } |
| #endif |
| decExpOp(res, rhs, set, &status); |
| } while(0); /* end protected */ |
| |
| #if DECSUBSET |
| free(allocrhs); /* drop any storage used */ |
| #endif |
| /* apply significant status */ |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberExp */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberFMA -- fused multiply add */ |
| /* */ |
| /* This computes D = (A * B) + C with only one rounding */ |
| /* */ |
| /* res is D, the result. D may be A or B or C (e.g., X=FMA(X,X,X)) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* fhs is C [far hand side] */ |
| /* set is the context */ |
| /* */ |
| /* Mathematical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberFMA(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, const decNumber *fhs, |
| decContext *set) { |
| uInt status=0; /* accumulator */ |
| decContext dcmul; /* context for the multiplication */ |
| uInt needbytes; /* for space calculations */ |
| decNumber bufa[D2N(DECBUFFER*2+1)]; |
| decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ |
| decNumber *acc; /* accumulator pointer */ |
| decNumber dzero; /* work */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| if (decCheckOperands(res, fhs, DECUNUSED, set)) return res; |
| #endif |
| |
| do { /* protect allocated storage */ |
| #if DECSUBSET |
| if (!set->extended) { /* [undefined if subset] */ |
| status|=DEC_Invalid_operation; |
| break;} |
| #endif |
| /* Check math restrictions [these ensure no overflow or underflow] */ |
| if ((!decNumberIsSpecial(lhs) && decCheckMath(lhs, set, &status)) |
| || (!decNumberIsSpecial(rhs) && decCheckMath(rhs, set, &status)) |
| || (!decNumberIsSpecial(fhs) && decCheckMath(fhs, set, &status))) break; |
| /* set up context for multiply */ |
| dcmul=*set; |
| dcmul.digits=lhs->digits+rhs->digits; /* just enough */ |
| /* [The above may be an over-estimate for subset arithmetic, but that's OK] */ |
| dcmul.emax=DEC_MAX_EMAX; /* effectively unbounded .. */ |
| dcmul.emin=DEC_MIN_EMIN; /* [thanks to Math restrictions] */ |
| /* set up decNumber space to receive the result of the multiply */ |
| acc=bufa; /* may fit */ |
| needbytes=sizeof(decNumber)+(D2U(dcmul.digits)-1)*sizeof(Unit); |
| if (needbytes>sizeof(bufa)) { /* need malloc space */ |
| allocbufa=(decNumber *)malloc(needbytes); |
| if (allocbufa==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| acc=allocbufa; /* use the allocated space */ |
| } |
| /* multiply with extended range and necessary precision */ |
| /*printf("emin=%ld\n", dcmul.emin); */ |
| decMultiplyOp(acc, lhs, rhs, &dcmul, &status); |
| /* Only Invalid operation (from sNaN or Inf * 0) is possible in */ |
| /* status; if either is seen than ignore fhs (in case it is */ |
| /* another sNaN) and set acc to NaN unless we had an sNaN */ |
| /* [decMultiplyOp leaves that to caller] */ |
| /* Note sNaN has to go through addOp to shorten payload if */ |
| /* necessary */ |
| if ((status&DEC_Invalid_operation)!=0) { |
| if (!(status&DEC_sNaN)) { /* but be true invalid */ |
| decNumberZero(res); /* acc not yet set */ |
| res->bits=DECNAN; |
| break; |
| } |
| decNumberZero(&dzero); /* make 0 (any non-NaN would do) */ |
| fhs=&dzero; /* use that */ |
| } |
| #if DECCHECK |
| else { /* multiply was OK */ |
| if (status!=0) printf("Status=%08lx after FMA multiply\n", (LI)status); |
| } |
| #endif |
| /* add the third operand and result -> res, and all is done */ |
| decAddOp(res, acc, fhs, set, 0, &status); |
| } while(0); /* end protected */ |
| |
| free(allocbufa); /* drop any storage used */ |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberFMA */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberInvert -- invert a Number, digitwise */ |
| /* */ |
| /* This computes C = ~A */ |
| /* */ |
| /* res is C, the result. C may be A (e.g., X=~X) */ |
| /* rhs is A */ |
| /* set is the context (used for result length and error report) */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Logical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberInvert(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| const Unit *ua, *msua; /* -> operand and its msu */ |
| Unit *uc, *msuc; /* -> result and its msu */ |
| Int msudigs; /* digits in res msu */ |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| if (rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| /* operand is valid */ |
| ua=rhs->lsu; /* bottom-up */ |
| uc=res->lsu; /* .. */ |
| msua=ua+D2U(rhs->digits)-1; /* -> msu of rhs */ |
| msuc=uc+D2U(set->digits)-1; /* -> msu of result */ |
| msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ |
| for (; uc<=msuc; ua++, uc++) { /* Unit loop */ |
| Unit a; /* extract unit */ |
| Int i, j; /* work */ |
| if (ua>msua) a=0; |
| else a=*ua; |
| *uc=0; /* can now write back */ |
| /* always need to examine all bits in rhs */ |
| /* This loop could be unrolled and/or use BIN2BCD tables */ |
| for (i=0; i<DECDPUN; i++) { |
| if ((~a)&1) *uc=*uc+(Unit)powers[i]; /* effect INVERT */ |
| j=a%10; |
| a=a/10; |
| if (j>1) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| if (uc==msuc && i==msudigs-1) break; /* just did final digit */ |
| } /* each digit */ |
| } /* each unit */ |
| /* [here uc-1 is the msu of the result] */ |
| res->digits=decGetDigits(res->lsu, uc-res->lsu); |
| res->exponent=0; /* integer */ |
| res->bits=0; /* sign=0 */ |
| return res; /* [no status to set] */ |
| } /* decNumberInvert */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberLn -- natural logarithm */ |
| /* */ |
| /* This computes C = ln(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context; note that rounding mode has no effect */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Notable cases: */ |
| /* A<0 -> Invalid */ |
| /* A=0 -> -Infinity (Exact) */ |
| /* A=+Infinity -> +Infinity (Exact) */ |
| /* A=1 exactly -> 0 (Exact) */ |
| /* */ |
| /* Mathematical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* */ |
| /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */ |
| /* almost always be correctly rounded, but may be up to 1 ulp in */ |
| /* error in rare cases. */ |
| /* ------------------------------------------------------------------ */ |
| /* This is a wrapper for decLnOp which can handle the slightly wider */ |
| /* (+11) range needed by Ln, Log10, etc. (which may have to be able */ |
| /* to calculate at p+e+2). */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberLn(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| uInt status=0; /* accumulator */ |
| #if DECSUBSET |
| decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ |
| #endif |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* Check restrictions; this is a math function; if not violated */ |
| /* then carry out the operation. */ |
| if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */ |
| #if DECSUBSET |
| if (!set->extended) { |
| /* reduce operand and set lostDigits status, as needed */ |
| if (rhs->digits>set->digits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| rhs=allocrhs; |
| } |
| /* special check in subset for rhs=0 */ |
| if (ISZERO(rhs)) { /* +/- zeros -> error */ |
| status|=DEC_Invalid_operation; |
| break;} |
| } /* extended=0 */ |
| #endif |
| decLnOp(res, rhs, set, &status); |
| } while(0); /* end protected */ |
| |
| #if DECSUBSET |
| free(allocrhs); /* drop any storage used */ |
| #endif |
| /* apply significant status */ |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberLn */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberLogB - get adjusted exponent, by 754 rules */ |
| /* */ |
| /* This computes C = adjustedexponent(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context, used only for digits and status */ |
| /* */ |
| /* C must have space for 10 digits (A might have 10**9 digits and */ |
| /* an exponent of +999999999, or one digit and an exponent of */ |
| /* -1999999999). */ |
| /* */ |
| /* This returns the adjusted exponent of A after (in theory) padding */ |
| /* with zeros on the right to set->digits digits while keeping the */ |
| /* same value. The exponent is not limited by emin/emax. */ |
| /* */ |
| /* Notable cases: */ |
| /* A<0 -> Use |A| */ |
| /* A=0 -> -Infinity (Division by zero) */ |
| /* A=Infinite -> +Infinity (Exact) */ |
| /* A=1 exactly -> 0 (Exact) */ |
| /* NaNs are propagated as usual */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberLogB(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| uInt status=0; /* accumulator */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* NaNs as usual; Infinities return +Infinity; 0->oops */ |
| if (decNumberIsNaN(rhs)) decNaNs(res, rhs, NULL, set, &status); |
| else if (decNumberIsInfinite(rhs)) decNumberCopyAbs(res, rhs); |
| else if (decNumberIsZero(rhs)) { |
| decNumberZero(res); /* prepare for Infinity */ |
| res->bits=DECNEG|DECINF; /* -Infinity */ |
| status|=DEC_Division_by_zero; /* as per 754 */ |
| } |
| else { /* finite non-zero */ |
| Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */ |
| decNumberFromInt32(res, ae); /* lay it out */ |
| } |
| |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberLogB */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberLog10 -- logarithm in base 10 */ |
| /* */ |
| /* This computes C = log10(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context; note that rounding mode has no effect */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Notable cases: */ |
| /* A<0 -> Invalid */ |
| /* A=0 -> -Infinity (Exact) */ |
| /* A=+Infinity -> +Infinity (Exact) */ |
| /* A=10**n (if n is an integer) -> n (Exact) */ |
| /* */ |
| /* Mathematical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* */ |
| /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */ |
| /* almost always be correctly rounded, but may be up to 1 ulp in */ |
| /* error in rare cases. */ |
| /* ------------------------------------------------------------------ */ |
| /* This calculates ln(A)/ln(10) using appropriate precision. For */ |
| /* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the */ |
| /* requested digits and t is the number of digits in the exponent */ |
| /* (maximum 6). For ln(10) it is p + 3; this is often handled by the */ |
| /* fastpath in decLnOp. The final division is done to the requested */ |
| /* precision. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberLog10(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| uInt status=0, ignore=0; /* status accumulators */ |
| uInt needbytes; /* for space calculations */ |
| Int p; /* working precision */ |
| Int t; /* digits in exponent of A */ |
| |
| /* buffers for a and b working decimals */ |
| /* (adjustment calculator, same size) */ |
| decNumber bufa[D2N(DECBUFFER+2)]; |
| decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ |
| decNumber *a=bufa; /* temporary a */ |
| decNumber bufb[D2N(DECBUFFER+2)]; |
| decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ |
| decNumber *b=bufb; /* temporary b */ |
| decNumber bufw[D2N(10)]; /* working 2-10 digit number */ |
| decNumber *w=bufw; /* .. */ |
| #if DECSUBSET |
| decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ |
| #endif |
| |
| decContext aset; /* working context */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* Check restrictions; this is a math function; if not violated */ |
| /* then carry out the operation. */ |
| if (!decCheckMath(rhs, set, &status)) do { /* protect malloc */ |
| #if DECSUBSET |
| if (!set->extended) { |
| /* reduce operand and set lostDigits status, as needed */ |
| if (rhs->digits>set->digits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| rhs=allocrhs; |
| } |
| /* special check in subset for rhs=0 */ |
| if (ISZERO(rhs)) { /* +/- zeros -> error */ |
| status|=DEC_Invalid_operation; |
| break;} |
| } /* extended=0 */ |
| #endif |
| |
| decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */ |
| |
| /* handle exact powers of 10; only check if +ve finite */ |
| if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) { |
| Int residue=0; /* (no residue) */ |
| uInt copystat=0; /* clean status */ |
| |
| /* round to a single digit... */ |
| aset.digits=1; |
| decCopyFit(w, rhs, &aset, &residue, ©stat); /* copy & shorten */ |
| /* if exact and the digit is 1, rhs is a power of 10 */ |
| if (!(copystat&DEC_Inexact) && w->lsu[0]==1) { |
| /* the exponent, conveniently, is the power of 10; making */ |
| /* this the result needs a little care as it might not fit, */ |
| /* so first convert it into the working number, and then move */ |
| /* to res */ |
| decNumberFromInt32(w, w->exponent); |
| residue=0; |
| decCopyFit(res, w, set, &residue, &status); /* copy & round */ |
| decFinish(res, set, &residue, &status); /* cleanup/set flags */ |
| break; |
| } /* not a power of 10 */ |
| } /* not a candidate for exact */ |
| |
| /* simplify the information-content calculation to use 'total */ |
| /* number of digits in a, including exponent' as compared to the */ |
| /* requested digits, as increasing this will only rarely cost an */ |
| /* iteration in ln(a) anyway */ |
| t=6; /* it can never be >6 */ |
| |
| /* allocate space when needed... */ |
| p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3; |
| needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit); |
| if (needbytes>sizeof(bufa)) { /* need malloc space */ |
| allocbufa=(decNumber *)malloc(needbytes); |
| if (allocbufa==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| a=allocbufa; /* use the allocated space */ |
| } |
| aset.digits=p; /* as calculated */ |
| aset.emax=DEC_MAX_MATH; /* usual bounds */ |
| aset.emin=-DEC_MAX_MATH; /* .. */ |
| aset.clamp=0; /* and no concrete format */ |
| decLnOp(a, rhs, &aset, &status); /* a=ln(rhs) */ |
| |
| /* skip the division if the result so far is infinite, NaN, or */ |
| /* zero, or there was an error; note NaN from sNaN needs copy */ |
| if (status&DEC_NaNs && !(status&DEC_sNaN)) break; |
| if (a->bits&DECSPECIAL || ISZERO(a)) { |
| decNumberCopy(res, a); /* [will fit] */ |
| break;} |
| |
| /* for ln(10) an extra 3 digits of precision are needed */ |
| p=set->digits+3; |
| needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit); |
| if (needbytes>sizeof(bufb)) { /* need malloc space */ |
| allocbufb=(decNumber *)malloc(needbytes); |
| if (allocbufb==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| b=allocbufb; /* use the allocated space */ |
| } |
| decNumberZero(w); /* set up 10... */ |
| #if DECDPUN==1 |
| w->lsu[1]=1; w->lsu[0]=0; /* .. */ |
| #else |
| w->lsu[0]=10; /* .. */ |
| #endif |
| w->digits=2; /* .. */ |
| |
| aset.digits=p; |
| decLnOp(b, w, &aset, &ignore); /* b=ln(10) */ |
| |
| aset.digits=set->digits; /* for final divide */ |
| decDivideOp(res, a, b, &aset, DIVIDE, &status); /* into result */ |
| } while(0); /* [for break] */ |
| |
| free(allocbufa); /* drop any storage used */ |
| free(allocbufb); /* .. */ |
| #if DECSUBSET |
| free(allocrhs); /* .. */ |
| #endif |
| /* apply significant status */ |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberLog10 */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMax -- compare two Numbers and return the maximum */ |
| /* */ |
| /* This computes C = A ? B, returning the maximum by 754 rules */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMax(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPMAX, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMax */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMaxMag -- compare and return the maximum by magnitude */ |
| /* */ |
| /* This computes C = A ? B, returning the maximum by 754 rules */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMaxMag(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPMAXMAG, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMaxMag */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMin -- compare two Numbers and return the minimum */ |
| /* */ |
| /* This computes C = A ? B, returning the minimum by 754 rules */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMin(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPMIN, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMin */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMinMag -- compare and return the minimum by magnitude */ |
| /* */ |
| /* This computes C = A ? B, returning the minimum by 754 rules */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X?X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMinMag(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decCompareOp(res, lhs, rhs, set, COMPMINMAG, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMinMag */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMinus -- prefix minus operator */ |
| /* */ |
| /* This computes C = 0 - A */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* See also decNumberCopyNegate for a quiet bitwise version of this. */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* Simply use AddOp for the subtract, which will do the necessary. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMinus(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decNumber dzero; |
| uInt status=0; /* accumulator */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| decNumberZero(&dzero); /* make 0 */ |
| dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ |
| decAddOp(res, &dzero, rhs, set, DECNEG, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMinus */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberNextMinus -- next towards -Infinity */ |
| /* */ |
| /* This computes C = A - infinitesimal, rounded towards -Infinity */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* This is a generalization of 754 NextDown. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberNextMinus(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decNumber dtiny; /* constant */ |
| decContext workset=*set; /* work */ |
| uInt status=0; /* accumulator */ |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* +Infinity is the special case */ |
| if ((rhs->bits&(DECINF|DECNEG))==DECINF) { |
| decSetMaxValue(res, set); /* is +ve */ |
| /* there is no status to set */ |
| return res; |
| } |
| decNumberZero(&dtiny); /* start with 0 */ |
| dtiny.lsu[0]=1; /* make number that is .. */ |
| dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ |
| workset.round=DEC_ROUND_FLOOR; |
| decAddOp(res, rhs, &dtiny, &workset, DECNEG, &status); |
| status&=DEC_Invalid_operation|DEC_sNaN; /* only sNaN Invalid please */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberNextMinus */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberNextPlus -- next towards +Infinity */ |
| /* */ |
| /* This computes C = A + infinitesimal, rounded towards +Infinity */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* This is a generalization of 754 NextUp. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberNextPlus(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decNumber dtiny; /* constant */ |
| decContext workset=*set; /* work */ |
| uInt status=0; /* accumulator */ |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| /* -Infinity is the special case */ |
| if ((rhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) { |
| decSetMaxValue(res, set); |
| res->bits=DECNEG; /* negative */ |
| /* there is no status to set */ |
| return res; |
| } |
| decNumberZero(&dtiny); /* start with 0 */ |
| dtiny.lsu[0]=1; /* make number that is .. */ |
| dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ |
| workset.round=DEC_ROUND_CEILING; |
| decAddOp(res, rhs, &dtiny, &workset, 0, &status); |
| status&=DEC_Invalid_operation|DEC_sNaN; /* only sNaN Invalid please */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberNextPlus */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberNextToward -- next towards rhs */ |
| /* */ |
| /* This computes C = A +/- infinitesimal, rounded towards */ |
| /* +/-Infinity in the direction of B, as per 754-1985 nextafter */ |
| /* modified during revision but dropped from 754-2008. */ |
| /* */ |
| /* res is C, the result. C may be A or B. */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* This is a generalization of 754-1985 NextAfter. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberNextToward(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| decNumber dtiny; /* constant */ |
| decContext workset=*set; /* work */ |
| Int result; /* .. */ |
| uInt status=0; /* accumulator */ |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { |
| decNaNs(res, lhs, rhs, set, &status); |
| } |
| else { /* Is numeric, so no chance of sNaN Invalid, etc. */ |
| result=decCompare(lhs, rhs, 0); /* sign matters */ |
| if (result==BADINT) status|=DEC_Insufficient_storage; /* rare */ |
| else { /* valid compare */ |
| if (result==0) decNumberCopySign(res, lhs, rhs); /* easy */ |
| else { /* differ: need NextPlus or NextMinus */ |
| uByte sub; /* add or subtract */ |
| if (result<0) { /* lhs<rhs, do nextplus */ |
| /* -Infinity is the special case */ |
| if ((lhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) { |
| decSetMaxValue(res, set); |
| res->bits=DECNEG; /* negative */ |
| return res; /* there is no status to set */ |
| } |
| workset.round=DEC_ROUND_CEILING; |
| sub=0; /* add, please */ |
| } /* plus */ |
| else { /* lhs>rhs, do nextminus */ |
| /* +Infinity is the special case */ |
| if ((lhs->bits&(DECINF|DECNEG))==DECINF) { |
| decSetMaxValue(res, set); |
| return res; /* there is no status to set */ |
| } |
| workset.round=DEC_ROUND_FLOOR; |
| sub=DECNEG; /* subtract, please */ |
| } /* minus */ |
| decNumberZero(&dtiny); /* start with 0 */ |
| dtiny.lsu[0]=1; /* make number that is .. */ |
| dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ |
| decAddOp(res, lhs, &dtiny, &workset, sub, &status); /* + or - */ |
| /* turn off exceptions if the result is a normal number */ |
| /* (including Nmin), otherwise let all status through */ |
| if (decNumberIsNormal(res, set)) status=0; |
| } /* unequal */ |
| } /* compare OK */ |
| } /* numeric */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberNextToward */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberOr -- OR two Numbers, digitwise */ |
| /* */ |
| /* This computes C = A | B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X|X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context (used for result length and error report) */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Logical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberOr(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| const Unit *ua, *ub; /* -> operands */ |
| const Unit *msua, *msub; /* -> operand msus */ |
| Unit *uc, *msuc; /* -> result and its msu */ |
| Int msudigs; /* digits in res msu */ |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs) |
| || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| /* operands are valid */ |
| ua=lhs->lsu; /* bottom-up */ |
| ub=rhs->lsu; /* .. */ |
| uc=res->lsu; /* .. */ |
| msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */ |
| msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */ |
| msuc=uc+D2U(set->digits)-1; /* -> msu of result */ |
| msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ |
| for (; uc<=msuc; ua++, ub++, uc++) { /* Unit loop */ |
| Unit a, b; /* extract units */ |
| if (ua>msua) a=0; |
| else a=*ua; |
| if (ub>msub) b=0; |
| else b=*ub; |
| *uc=0; /* can now write back */ |
| if (a|b) { /* maybe 1 bits to examine */ |
| Int i, j; |
| /* This loop could be unrolled and/or use BIN2BCD tables */ |
| for (i=0; i<DECDPUN; i++) { |
| if ((a|b)&1) *uc=*uc+(Unit)powers[i]; /* effect OR */ |
| j=a%10; |
| a=a/10; |
| j|=b%10; |
| b=b/10; |
| if (j>1) { |
| decStatus(res, DEC_Invalid_operation, set); |
| return res; |
| } |
| if (uc==msuc && i==msudigs-1) break; /* just did final digit */ |
| } /* each digit */ |
| } /* non-zero */ |
| } /* each unit */ |
| /* [here uc-1 is the msu of the result] */ |
| res->digits=decGetDigits(res->lsu, uc-res->lsu); |
| res->exponent=0; /* integer */ |
| res->bits=0; /* sign=0 */ |
| return res; /* [no status to set] */ |
| } /* decNumberOr */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberPlus -- prefix plus operator */ |
| /* */ |
| /* This computes C = 0 + A */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* See also decNumberCopy for a quiet bitwise version of this. */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* This simply uses AddOp; Add will take fast path after preparing A. */ |
| /* Performance is a concern here, as this routine is often used to */ |
| /* check operands and apply rounding and overflow/underflow testing. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberPlus(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decNumber dzero; |
| uInt status=0; /* accumulator */ |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| decNumberZero(&dzero); /* make 0 */ |
| dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ |
| decAddOp(res, &dzero, rhs, set, 0, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberPlus */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberMultiply -- multiply two Numbers */ |
| /* */ |
| /* This computes C = A x B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X+X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decMultiplyOp(res, lhs, rhs, set, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberMultiply */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberPower -- raise a number to a power */ |
| /* */ |
| /* This computes C = A ** B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X**X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Mathematical function restrictions apply (see above); a NaN is */ |
| /* returned with Invalid_operation if a restriction is violated. */ |
| /* */ |
| /* However, if 1999999997<=B<=999999999 and B is an integer then the */ |
| /* restrictions on A and the context are relaxed to the usual bounds, */ |
| /* for compatibility with the earlier (integer power only) version */ |
| /* of this function. */ |
| /* */ |
| /* When B is an integer, the result may be exact, even if rounded. */ |
| /* */ |
| /* The final result is rounded according to the context; it will */ |
| /* almost always be correctly rounded, but may be up to 1 ulp in */ |
| /* error in rare cases. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberPower(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| #if DECSUBSET |
| decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */ |
| decNumber *allocrhs=NULL; /* .., rhs */ |
| #endif |
| decNumber *allocdac=NULL; /* -> allocated acc buffer, iff used */ |
| decNumber *allocinv=NULL; /* -> allocated 1/x buffer, iff used */ |
| Int reqdigits=set->digits; /* requested DIGITS */ |
| Int n; /* rhs in binary */ |
| Flag rhsint=0; /* 1 if rhs is an integer */ |
| Flag useint=0; /* 1 if can use integer calculation */ |
| Flag isoddint=0; /* 1 if rhs is an integer and odd */ |
| Int i; /* work */ |
| #if DECSUBSET |
| Int dropped; /* .. */ |
| #endif |
| uInt needbytes; /* buffer size needed */ |
| Flag seenbit; /* seen a bit while powering */ |
| Int residue=0; /* rounding residue */ |
| uInt status=0; /* accumulators */ |
| uByte bits=0; /* result sign if errors */ |
| decContext aset; /* working context */ |
| decNumber dnOne; /* work value 1... */ |
| /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */ |
| decNumber dacbuff[D2N(DECBUFFER+9)]; |
| decNumber *dac=dacbuff; /* -> result accumulator */ |
| /* same again for possible 1/lhs calculation */ |
| decNumber invbuff[D2N(DECBUFFER+9)]; |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| do { /* protect allocated storage */ |
| #if DECSUBSET |
| if (!set->extended) { /* reduce operands and set status, as needed */ |
| if (lhs->digits>reqdigits) { |
| alloclhs=decRoundOperand(lhs, set, &status); |
| if (alloclhs==NULL) break; |
| lhs=alloclhs; |
| } |
| if (rhs->digits>reqdigits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| rhs=allocrhs; |
| } |
| } |
| #endif |
| /* [following code does not require input rounding] */ |
| |
| /* handle NaNs and rhs Infinity (lhs infinity is harder) */ |
| if (SPECIALARGS) { |
| if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { /* NaNs */ |
| decNaNs(res, lhs, rhs, set, &status); |
| break;} |
| if (decNumberIsInfinite(rhs)) { /* rhs Infinity */ |
| Flag rhsneg=rhs->bits&DECNEG; /* save rhs sign */ |
| if (decNumberIsNegative(lhs) /* lhs<0 */ |
| && !decNumberIsZero(lhs)) /* .. */ |
| status|=DEC_Invalid_operation; |
| else { /* lhs >=0 */ |
| decNumberZero(&dnOne); /* set up 1 */ |
| dnOne.lsu[0]=1; |
| decNumberCompare(dac, lhs, &dnOne, set); /* lhs ? 1 */ |
| decNumberZero(res); /* prepare for 0/1/Infinity */ |
| if (decNumberIsNegative(dac)) { /* lhs<1 */ |
| if (rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */ |
| } |
| else if (dac->lsu[0]==0) { /* lhs=1 */ |
| /* 1**Infinity is inexact, so return fully-padded 1.0000 */ |
| Int shift=set->digits-1; |
| *res->lsu=1; /* was 0, make int 1 */ |
| res->digits=decShiftToMost(res->lsu, 1, shift); |
| res->exponent=-shift; /* make 1.0000... */ |
| status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */ |
| } |
| else { /* lhs>1 */ |
| if (!rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */ |
| } |
| } /* lhs>=0 */ |
| break;} |
| /* [lhs infinity drops through] */ |
| } /* specials */ |
| |
| /* Original rhs may be an integer that fits and is in range */ |
| n=decGetInt(rhs); |
| if (n!=BADINT) { /* it is an integer */ |
| rhsint=1; /* record the fact for 1**n */ |
| isoddint=(Flag)n&1; /* [works even if big] */ |
| if (n!=BIGEVEN && n!=BIGODD) /* can use integer path? */ |
| useint=1; /* looks good */ |
| } |
| |
| if (decNumberIsNegative(lhs) /* -x .. */ |
| && isoddint) bits=DECNEG; /* .. to an odd power */ |
| |
| /* handle LHS infinity */ |
| if (decNumberIsInfinite(lhs)) { /* [NaNs already handled] */ |
| uByte rbits=rhs->bits; /* save */ |
| decNumberZero(res); /* prepare */ |
| if (n==0) *res->lsu=1; /* [-]Inf**0 => 1 */ |
| else { |
| /* -Inf**nonint -> error */ |
| if (!rhsint && decNumberIsNegative(lhs)) { |
| status|=DEC_Invalid_operation; /* -Inf**nonint is error */ |
| break;} |
| if (!(rbits & DECNEG)) bits|=DECINF; /* was not a **-n */ |
| /* [otherwise will be 0 or -0] */ |
| res->bits=bits; |
| } |
| break;} |
| |
| /* similarly handle LHS zero */ |
| if (decNumberIsZero(lhs)) { |
| if (n==0) { /* 0**0 => Error */ |
| #if DECSUBSET |
| if (!set->extended) { /* [unless subset] */ |
| decNumberZero(res); |
| *res->lsu=1; /* return 1 */ |
| break;} |
| #endif |
| status|=DEC_Invalid_operation; |
| } |
| else { /* 0**x */ |
| uByte rbits=rhs->bits; /* save */ |
| if (rbits & DECNEG) { /* was a 0**(-n) */ |
| #if DECSUBSET |
| if (!set->extended) { /* [bad if subset] */ |
| status|=DEC_Invalid_operation; |
| break;} |
| #endif |
| bits|=DECINF; |
| } |
| decNumberZero(res); /* prepare */ |
| /* [otherwise will be 0 or -0] */ |
| res->bits=bits; |
| } |
| break;} |
| |
| /* here both lhs and rhs are finite; rhs==0 is handled in the */ |
| /* integer path. Next handle the non-integer cases */ |
| if (!useint) { /* non-integral rhs */ |
| /* any -ve lhs is bad, as is either operand or context out of */ |
| /* bounds */ |
| if (decNumberIsNegative(lhs)) { |
| status|=DEC_Invalid_operation; |
| break;} |
| if (decCheckMath(lhs, set, &status) |
| || decCheckMath(rhs, set, &status)) break; /* variable status */ |
| |
| decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */ |
| aset.emax=DEC_MAX_MATH; /* usual bounds */ |
| aset.emin=-DEC_MAX_MATH; /* .. */ |
| aset.clamp=0; /* and no concrete format */ |
| |
| /* calculate the result using exp(ln(lhs)*rhs), which can */ |
| /* all be done into the accumulator, dac. The precision needed */ |
| /* is enough to contain the full information in the lhs (which */ |
| /* is the total digits, including exponent), or the requested */ |
| /* precision, if larger, + 4; 6 is used for the exponent */ |
| /* maximum length, and this is also used when it is shorter */ |
| /* than the requested digits as it greatly reduces the >0.5 ulp */ |
| /* cases at little cost (because Ln doubles digits each */ |
| /* iteration so a few extra digits rarely causes an extra */ |
| /* iteration) */ |
| aset.digits=MAXI(lhs->digits, set->digits)+6+4; |
| } /* non-integer rhs */ |
| |
| else { /* rhs is in-range integer */ |
| if (n==0) { /* x**0 = 1 */ |
| /* (0**0 was handled above) */ |
| decNumberZero(res); /* result=1 */ |
| *res->lsu=1; /* .. */ |
| break;} |
| /* rhs is a non-zero integer */ |
| if (n<0) n=-n; /* use abs(n) */ |
| |
| aset=*set; /* clone the context */ |
| aset.round=DEC_ROUND_HALF_EVEN; /* internally use balanced */ |
| /* calculate the working DIGITS */ |
| aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2; |
| #if DECSUBSET |
| if (!set->extended) aset.digits--; /* use classic precision */ |
| #endif |
| /* it's an error if this is more than can be handled */ |
| if (aset.digits>DECNUMMAXP) {status|=DEC_Invalid_operation; break;} |
| } /* integer path */ |
| |
| /* aset.digits is the count of digits for the accumulator needed */ |
| /* if accumulator is too long for local storage, then allocate */ |
| needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit); |
| /* [needbytes also used below if 1/lhs needed] */ |
| if (needbytes>sizeof(dacbuff)) { |
| allocdac=(decNumber *)malloc(needbytes); |
| if (allocdac==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| dac=allocdac; /* use the allocated space */ |
| } |
| /* here, aset is set up and accumulator is ready for use */ |
| |
| if (!useint) { /* non-integral rhs */ |
| /* x ** y; special-case x=1 here as it will otherwise always */ |
| /* reduce to integer 1; decLnOp has a fastpath which detects */ |
| /* the case of x=1 */ |
| decLnOp(dac, lhs, &aset, &status); /* dac=ln(lhs) */ |
| /* [no error possible, as lhs 0 already handled] */ |
| if (ISZERO(dac)) { /* x==1, 1.0, etc. */ |
| /* need to return fully-padded 1.0000 etc., but rhsint->1 */ |
| *dac->lsu=1; /* was 0, make int 1 */ |
| if (!rhsint) { /* add padding */ |
| Int shift=set->digits-1; |
| dac->digits=decShiftToMost(dac->lsu, 1, shift); |
| dac->exponent=-shift; /* make 1.0000... */ |
| status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */ |
| } |
| } |
| else { |
| decMultiplyOp(dac, dac, rhs, &aset, &status); /* dac=dac*rhs */ |
| decExpOp(dac, dac, &aset, &status); /* dac=exp(dac) */ |
| } |
| /* and drop through for final rounding */ |
| } /* non-integer rhs */ |
| |
| else { /* carry on with integer */ |
| decNumberZero(dac); /* acc=1 */ |
| *dac->lsu=1; /* .. */ |
| |
| /* if a negative power the constant 1 is needed, and if not subset */ |
| /* invert the lhs now rather than inverting the result later */ |
| if (decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */ |
| decNumber *inv=invbuff; /* asssume use fixed buffer */ |
| decNumberCopy(&dnOne, dac); /* dnOne=1; [needed now or later] */ |
| #if DECSUBSET |
| if (set->extended) { /* need to calculate 1/lhs */ |
| #endif |
| /* divide lhs into 1, putting result in dac [dac=1/dac] */ |
| decDivideOp(dac, &dnOne, lhs, &aset, DIVIDE, &status); |
| /* now locate or allocate space for the inverted lhs */ |
| if (needbytes>sizeof(invbuff)) { |
| allocinv=(decNumber *)malloc(needbytes); |
| if (allocinv==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| inv=allocinv; /* use the allocated space */ |
| } |
| /* [inv now points to big-enough buffer or allocated storage] */ |
| decNumberCopy(inv, dac); /* copy the 1/lhs */ |
| decNumberCopy(dac, &dnOne); /* restore acc=1 */ |
| lhs=inv; /* .. and go forward with new lhs */ |
| #if DECSUBSET |
| } |
| #endif |
| } |
| |
| /* Raise-to-the-power loop... */ |
| seenbit=0; /* set once a 1-bit is encountered */ |
| for (i=1;;i++){ /* for each bit [top bit ignored] */ |
| /* abandon if had overflow or terminal underflow */ |
| if (status & (DEC_Overflow|DEC_Underflow)) { /* interesting? */ |
| if (status&DEC_Overflow || ISZERO(dac)) break; |
| } |
| /* [the following two lines revealed an optimizer bug in a C++ */ |
| /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */ |
| n=n<<1; /* move next bit to testable position */ |
| if (n<0) { /* top bit is set */ |
| seenbit=1; /* OK, significant bit seen */ |
| decMultiplyOp(dac, dac, lhs, &aset, &status); /* dac=dac*x */ |
| } |
| if (i==31) break; /* that was the last bit */ |
| if (!seenbit) continue; /* no need to square 1 */ |
| decMultiplyOp(dac, dac, dac, &aset, &status); /* dac=dac*dac [square] */ |
| } /*i*/ /* 32 bits */ |
| |
| /* complete internal overflow or underflow processing */ |
| if (status & (DEC_Overflow|DEC_Underflow)) { |
| #if DECSUBSET |
| /* If subset, and power was negative, reverse the kind of -erflow */ |
| /* [1/x not yet done] */ |
| if (!set->extended && decNumberIsNegative(rhs)) { |
| if (status & DEC_Overflow) |
| status^=DEC_Overflow | DEC_Underflow | DEC_Subnormal; |
| else { /* trickier -- Underflow may or may not be set */ |
| status&=~(DEC_Underflow | DEC_Subnormal); /* [one or both] */ |
| status|=DEC_Overflow; |
| } |
| } |
| #endif |
| dac->bits=(dac->bits & ~DECNEG) | bits; /* force correct sign */ |
| /* round subnormals [to set.digits rather than aset.digits] */ |
| /* or set overflow result similarly as required */ |
| decFinalize(dac, set, &residue, &status); |
| decNumberCopy(res, dac); /* copy to result (is now OK length) */ |
| break; |
| } |
| |
| #if DECSUBSET |
| if (!set->extended && /* subset math */ |
| decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */ |
| /* so divide result into 1 [dac=1/dac] */ |
| decDivideOp(dac, &dnOne, dac, &aset, DIVIDE, &status); |
| } |
| #endif |
| } /* rhs integer path */ |
| |
| /* reduce result to the requested length and copy to result */ |
| decCopyFit(res, dac, set, &residue, &status); |
| decFinish(res, set, &residue, &status); /* final cleanup */ |
| #if DECSUBSET |
| if (!set->extended) decTrim(res, set, 0, 1, &dropped); /* trailing zeros */ |
| #endif |
| } while(0); /* end protected */ |
| |
| free(allocdac); /* drop any storage used */ |
| free(allocinv); /* .. */ |
| #if DECSUBSET |
| free(alloclhs); /* .. */ |
| free(allocrhs); /* .. */ |
| #endif |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberPower */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberQuantize -- force exponent to requested value */ |
| /* */ |
| /* This computes C = op(A, B), where op adjusts the coefficient */ |
| /* of C (by rounding or shifting) such that the exponent (-scale) */ |
| /* of C has exponent of B. The numerical value of C will equal A, */ |
| /* except for the effects of any rounding that occurred. */ |
| /* */ |
| /* res is C, the result. C may be A or B */ |
| /* lhs is A, the number to adjust */ |
| /* rhs is B, the number with exponent to match */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Unless there is an error or the result is infinite, the exponent */ |
| /* after the operation is guaranteed to be equal to that of B. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberQuantize(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decQuantizeOp(res, lhs, rhs, set, 1, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberQuantize */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberReduce -- remove trailing zeros */ |
| /* */ |
| /* This computes C = 0 + A, and normalizes the result */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* Previously known as Normalize */ |
| decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| return decNumberReduce(res, rhs, set); |
| } /* decNumberNormalize */ |
| |
| decNumber * decNumberReduce(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| #if DECSUBSET |
| decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ |
| #endif |
| uInt status=0; /* as usual */ |
| Int residue=0; /* as usual */ |
| Int dropped; /* work */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| do { /* protect allocated storage */ |
| #if DECSUBSET |
| if (!set->extended) { |
| /* reduce operand and set lostDigits status, as needed */ |
| if (rhs->digits>set->digits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| rhs=allocrhs; |
| } |
| } |
| #endif |
| /* [following code does not require input rounding] */ |
| |
| /* Infinities copy through; NaNs need usual treatment */ |
| if (decNumberIsNaN(rhs)) { |
| decNaNs(res, rhs, NULL, set, &status); |
| break; |
| } |
| |
| /* reduce result to the requested length and copy to result */ |
| decCopyFit(res, rhs, set, &residue, &status); /* copy & round */ |
| decFinish(res, set, &residue, &status); /* cleanup/set flags */ |
| decTrim(res, set, 1, 0, &dropped); /* normalize in place */ |
| /* [may clamp] */ |
| } while(0); /* end protected */ |
| |
| #if DECSUBSET |
| free(allocrhs); /* .. */ |
| #endif |
| if (status!=0) decStatus(res, status, set);/* then report status */ |
| return res; |
| } /* decNumberReduce */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberRescale -- force exponent to requested value */ |
| /* */ |
| /* This computes C = op(A, B), where op adjusts the coefficient */ |
| /* of C (by rounding or shifting) such that the exponent (-scale) */ |
| /* of C has the value B. The numerical value of C will equal A, */ |
| /* except for the effects of any rounding that occurred. */ |
| /* */ |
| /* res is C, the result. C may be A or B */ |
| /* lhs is A, the number to adjust */ |
| /* rhs is B, the requested exponent */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* Unless there is an error or the result is infinite, the exponent */ |
| /* after the operation is guaranteed to be equal to B. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberRescale(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decQuantizeOp(res, lhs, rhs, set, 0, &status); |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberRescale */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberRemainder -- divide and return remainder */ |
| /* */ |
| /* This computes C = A % B */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X%X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decDivideOp(res, lhs, rhs, set, REMAINDER, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberRemainder */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberRemainderNear -- divide and return remainder from nearest */ |
| /* */ |
| /* This computes C = A % B, where % is the IEEE remainder operator */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X%X) */ |
| /* lhs is A */ |
| /* rhs is B */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| decDivideOp(res, lhs, rhs, set, REMNEAR, &status); |
| if (status!=0) decStatus(res, status, set); |
| #if DECCHECK |
| decCheckInexact(res, set); |
| #endif |
| return res; |
| } /* decNumberRemainderNear */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberRotate -- rotate the coefficient of a Number left/right */ |
| /* */ |
| /* This computes C = A rot B (in base ten and rotating set->digits */ |
| /* digits). */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=XrotX) */ |
| /* lhs is A */ |
| /* rhs is B, the number of digits to rotate (-ve to right) */ |
| /* set is the context */ |
| /* */ |
| /* The digits of the coefficient of A are rotated to the left (if B */ |
| /* is positive) or to the right (if B is negative) without adjusting */ |
| /* the exponent or the sign of A. If lhs->digits is less than */ |
| /* set->digits the coefficient is padded with zeros on the left */ |
| /* before the rotate. Any leading zeros in the result are removed */ |
| /* as usual. */ |
| /* */ |
| /* B must be an integer (q=0) and in the range -set->digits through */ |
| /* +set->digits. */ |
| /* C must have space for set->digits digits. */ |
| /* NaNs are propagated as usual. Infinities are unaffected (but */ |
| /* B must be valid). No status is set unless B is invalid or an */ |
| /* operand is an sNaN. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberRotate(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| Int rotate; /* rhs as an Int */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| /* NaNs propagate as normal */ |
| if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) |
| decNaNs(res, lhs, rhs, set, &status); |
| /* rhs must be an integer */ |
| else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) |
| status=DEC_Invalid_operation; |
| else { /* both numeric, rhs is an integer */ |
| rotate=decGetInt(rhs); /* [cannot fail] */ |
| if (rotate==BADINT /* something bad .. */ |
| || rotate==BIGODD || rotate==BIGEVEN /* .. very big .. */ |
| || abs(rotate)>set->digits) /* .. or out of range */ |
| status=DEC_Invalid_operation; |
| else { /* rhs is OK */ |
| decNumberCopy(res, lhs); |
| /* convert -ve rotate to equivalent positive rotation */ |
| if (rotate<0) rotate=set->digits+rotate; |
| if (rotate!=0 && rotate!=set->digits /* zero or full rotation */ |
| && !decNumberIsInfinite(res)) { /* lhs was infinite */ |
| /* left-rotate to do; 0 < rotate < set->digits */ |
| uInt units, shift; /* work */ |
| uInt msudigits; /* digits in result msu */ |
| Unit *msu=res->lsu+D2U(res->digits)-1; /* current msu */ |
| Unit *msumax=res->lsu+D2U(set->digits)-1; /* rotation msu */ |
| for (msu++; msu<=msumax; msu++) *msu=0; /* ensure high units=0 */ |
| res->digits=set->digits; /* now full-length */ |
| msudigits=MSUDIGITS(res->digits); /* actual digits in msu */ |
| |
| /* rotation here is done in-place, in three steps */ |
| /* 1. shift all to least up to one unit to unit-align final */ |
| /* lsd [any digits shifted out are rotated to the left, */ |
| /* abutted to the original msd (which may require split)] */ |
| /* */ |
| /* [if there are no whole units left to rotate, the */ |
| /* rotation is now complete] */ |
| /* */ |
| /* 2. shift to least, from below the split point only, so that */ |
| /* the final msd is in the right place in its Unit [any */ |
| /* digits shifted out will fit exactly in the current msu, */ |
| /* left aligned, no split required] */ |
| /* */ |
| /* 3. rotate all the units by reversing left part, right */ |
| /* part, and then whole */ |
| /* */ |
| /* example: rotate right 8 digits (2 units + 2), DECDPUN=3. */ |
| /* */ |
| /* start: 00a bcd efg hij klm npq */ |
| /* */ |
| /* 1a 000 0ab cde fgh|ijk lmn [pq saved] */ |
| /* 1b 00p qab cde fgh|ijk lmn */ |
| /* */ |
| /* 2a 00p qab cde fgh|00i jkl [mn saved] */ |
| /* 2b mnp qab cde fgh|00i jkl */ |
| /* */ |
| /* 3a fgh cde qab mnp|00i jkl */ |
| /* 3b fgh cde qab mnp|jkl 00i */ |
| /* 3c 00i jkl mnp qab cde fgh */ |
| |
| /* Step 1: amount to shift is the partial right-rotate count */ |
| rotate=set->digits-rotate; /* make it right-rotate */ |
| units=rotate/DECDPUN; /* whole units to rotate */ |
| shift=rotate%DECDPUN; /* left-over digits count */ |
| if (shift>0) { /* not an exact number of units */ |
| uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */ |
| decShiftToLeast(res->lsu, D2U(res->digits), shift); |
| if (shift>msudigits) { /* msumax-1 needs >0 digits */ |
| uInt rem=save%powers[shift-msudigits];/* split save */ |
| *msumax=(Unit)(save/powers[shift-msudigits]); /* and insert */ |
| *(msumax-1)=*(msumax-1) |
| +(Unit)(rem*powers[DECDPUN-(shift-msudigits)]); /* .. */ |
| } |
| else { /* all fits in msumax */ |
| *msumax=*msumax+(Unit)(save*powers[msudigits-shift]); /* [maybe *1] */ |
| } |
| } /* digits shift needed */ |
| |
| /* If whole units to rotate... */ |
| if (units>0) { /* some to do */ |
| /* Step 2: the units to touch are the whole ones in rotate, */ |
| /* if any, and the shift is DECDPUN-msudigits (which may be */ |
| /* 0, again) */ |
| shift=DECDPUN-msudigits; |
| if (shift>0) { /* not an exact number of units */ |
| uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */ |
| decShiftToLeast(res->lsu, units, shift); |
| *msumax=*msumax+(Unit)(save*powers[msudigits]); |
| } /* partial shift needed */ |
| |
| /* Step 3: rotate the units array using triple reverse */ |
| /* (reversing is easy and fast) */ |
| decReverse(res->lsu+units, msumax); /* left part */ |
| decReverse(res->lsu, res->lsu+units-1); /* right part */ |
| decReverse(res->lsu, msumax); /* whole */ |
| } /* whole units to rotate */ |
| /* the rotation may have left an undetermined number of zeros */ |
| /* on the left, so true length needs to be calculated */ |
| res->digits=decGetDigits(res->lsu, msumax-res->lsu+1); |
| } /* rotate needed */ |
| } /* rhs OK */ |
| } /* numerics */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberRotate */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberSameQuantum -- test for equal exponents */ |
| /* */ |
| /* res is the result number, which will contain either 0 or 1 */ |
| /* lhs is a number to test */ |
| /* rhs is the second (usually a pattern) */ |
| /* */ |
| /* No errors are possible and no context is needed. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs) { |
| Unit ret=0; /* return value */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, DECUNCONT)) return res; |
| #endif |
| |
| if (SPECIALARGS) { |
| if (decNumberIsNaN(lhs) && decNumberIsNaN(rhs)) ret=1; |
| else if (decNumberIsInfinite(lhs) && decNumberIsInfinite(rhs)) ret=1; |
| /* [anything else with a special gives 0] */ |
| } |
| else if (lhs->exponent==rhs->exponent) ret=1; |
| |
| decNumberZero(res); /* OK to overwrite an operand now */ |
| *res->lsu=ret; |
| return res; |
| } /* decNumberSameQuantum */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberScaleB -- multiply by a power of 10 */ |
| /* */ |
| /* This computes C = A x 10**B where B is an integer (q=0) with */ |
| /* maximum magnitude 2*(emax+digits) */ |
| /* */ |
| /* res is C, the result. C may be A or B */ |
| /* lhs is A, the number to adjust */ |
| /* rhs is B, the requested power of ten to use */ |
| /* set is the context */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* */ |
| /* The result may underflow or overflow. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberScaleB(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| Int reqexp; /* requested exponent change [B] */ |
| uInt status=0; /* accumulator */ |
| Int residue; /* work */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| /* Handle special values except lhs infinite */ |
| if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) |
| decNaNs(res, lhs, rhs, set, &status); |
| /* rhs must be an integer */ |
| else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) |
| status=DEC_Invalid_operation; |
| else { |
| /* lhs is a number; rhs is a finite with q==0 */ |
| reqexp=decGetInt(rhs); /* [cannot fail] */ |
| if (reqexp==BADINT /* something bad .. */ |
| || reqexp==BIGODD || reqexp==BIGEVEN /* .. very big .. */ |
| || abs(reqexp)>(2*(set->digits+set->emax))) /* .. or out of range */ |
| status=DEC_Invalid_operation; |
| else { /* rhs is OK */ |
| decNumberCopy(res, lhs); /* all done if infinite lhs */ |
| if (!decNumberIsInfinite(res)) { /* prepare to scale */ |
| res->exponent+=reqexp; /* adjust the exponent */ |
| residue=0; |
| decFinalize(res, set, &residue, &status); /* .. and check */ |
| } /* finite LHS */ |
| } /* rhs OK */ |
| } /* rhs finite */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberScaleB */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberShift -- shift the coefficient of a Number left or right */ |
| /* */ |
| /* This computes C = A << B or C = A >> -B (in base ten). */ |
| /* */ |
| /* res is C, the result. C may be A and/or B (e.g., X=X<<X) */ |
| /* lhs is A */ |
| /* rhs is B, the number of digits to shift (-ve to right) */ |
| /* set is the context */ |
| /* */ |
| /* The digits of the coefficient of A are shifted to the left (if B */ |
| /* is positive) or to the right (if B is negative) without adjusting */ |
| /* the exponent or the sign of A. */ |
| /* */ |
| /* B must be an integer (q=0) and in the range -set->digits through */ |
| /* +set->digits. */ |
| /* C must have space for set->digits digits. */ |
| /* NaNs are propagated as usual. Infinities are unaffected (but */ |
| /* B must be valid). No status is set unless B is invalid or an */ |
| /* operand is an sNaN. */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberShift(decNumber *res, const decNumber *lhs, |
| const decNumber *rhs, decContext *set) { |
| uInt status=0; /* accumulator */ |
| Int shift; /* rhs as an Int */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, lhs, rhs, set)) return res; |
| #endif |
| |
| /* NaNs propagate as normal */ |
| if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) |
| decNaNs(res, lhs, rhs, set, &status); |
| /* rhs must be an integer */ |
| else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) |
| status=DEC_Invalid_operation; |
| else { /* both numeric, rhs is an integer */ |
| shift=decGetInt(rhs); /* [cannot fail] */ |
| if (shift==BADINT /* something bad .. */ |
| || shift==BIGODD || shift==BIGEVEN /* .. very big .. */ |
| || abs(shift)>set->digits) /* .. or out of range */ |
| status=DEC_Invalid_operation; |
| else { /* rhs is OK */ |
| decNumberCopy(res, lhs); |
| if (shift!=0 && !decNumberIsInfinite(res)) { /* something to do */ |
| if (shift>0) { /* to left */ |
| if (shift==set->digits) { /* removing all */ |
| *res->lsu=0; /* so place 0 */ |
| res->digits=1; /* .. */ |
| } |
| else { /* */ |
| /* first remove leading digits if necessary */ |
| if (res->digits+shift>set->digits) { |
| decDecap(res, res->digits+shift-set->digits); |
| /* that updated res->digits; may have gone to 1 (for a */ |
| /* single digit or for zero */ |
| } |
| if (res->digits>1 || *res->lsu) /* if non-zero.. */ |
| res->digits=decShiftToMost(res->lsu, res->digits, shift); |
| } /* partial left */ |
| } /* left */ |
| else { /* to right */ |
| if (-shift>=res->digits) { /* discarding all */ |
| *res->lsu=0; /* so place 0 */ |
| res->digits=1; /* .. */ |
| } |
| else { |
| decShiftToLeast(res->lsu, D2U(res->digits), -shift); |
| res->digits-=(-shift); |
| } |
| } /* to right */ |
| } /* non-0 non-Inf shift */ |
| } /* rhs OK */ |
| } /* numerics */ |
| if (status!=0) decStatus(res, status, set); |
| return res; |
| } /* decNumberShift */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decNumberSquareRoot -- square root operator */ |
| /* */ |
| /* This computes C = squareroot(A) */ |
| /* */ |
| /* res is C, the result. C may be A */ |
| /* rhs is A */ |
| /* set is the context; note that rounding mode has no effect */ |
| /* */ |
| /* C must have space for set->digits digits. */ |
| /* ------------------------------------------------------------------ */ |
| /* This uses the following varying-precision algorithm in: */ |
| /* */ |
| /* Properly Rounded Variable Precision Square Root, T. E. Hull and */ |
| /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */ |
| /* pp229-237, ACM, September 1985. */ |
| /* */ |
| /* The square-root is calculated using Newton's method, after which */ |
| /* a check is made to ensure the result is correctly rounded. */ |
| /* */ |
| /* % [Reformatted original Numerical Turing source code follows.] */ |
| /* function sqrt(x : real) : real */ |
| /* % sqrt(x) returns the properly rounded approximation to the square */ |
| /* % root of x, in the precision of the calling environment, or it */ |
| /* % fails if x < 0. */ |
| /* % t e hull and a abrham, august, 1984 */ |
| /* if x <= 0 then */ |
| /* if x < 0 then */ |
| /* assert false */ |
| /* else */ |
| /* result 0 */ |
| /* end if */ |
| /* end if */ |
| /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */ |
| /* var e := getexp(x) % exponent part of x */ |
| /* var approx : real */ |
| /* if e mod 2 = 0 then */ |
| /* approx := .259 + .819 * f % approx to root of f */ |
| /* else */ |
| /* f := f/l0 % adjustments */ |
| /* e := e + 1 % for odd */ |
| /* approx := .0819 + 2.59 * f % exponent */ |
| /* end if */ |
| /* */ |
| /* var p:= 3 */ |
| /* const maxp := currentprecision + 2 */ |
| /* loop */ |
| /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */ |
| /* precision p */ |
| /* approx := .5 * (approx + f/approx) */ |
| /* exit when p = maxp */ |
| /* end loop */ |
| /* */ |
| /* % approx is now within 1 ulp of the properly rounded square root */ |
| /* % of f; to ensure proper rounding, compare squares of (approx - */ |
| /* % l/2 ulp) and (approx + l/2 ulp) with f. */ |
| /* p := currentprecision */ |
| /* begin */ |
| /* precision p + 2 */ |
| /* const approxsubhalf := approx - setexp(.5, -p) */ |
| /* if mulru(approxsubhalf, approxsubhalf) > f then */ |
| /* approx := approx - setexp(.l, -p + 1) */ |
| /* else */ |
| /* const approxaddhalf := approx + setexp(.5, -p) */ |
| /* if mulrd(approxaddhalf, approxaddhalf) < f then */ |
| /* approx := approx + setexp(.l, -p + 1) */ |
| /* end if */ |
| /* end if */ |
| /* end */ |
| /* result setexp(approx, e div 2) % fix exponent */ |
| /* end sqrt */ |
| /* ------------------------------------------------------------------ */ |
| decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs, |
| decContext *set) { |
| decContext workset, approxset; /* work contexts */ |
| decNumber dzero; /* used for constant zero */ |
| Int maxp; /* largest working precision */ |
| Int workp; /* working precision */ |
| Int residue=0; /* rounding residue */ |
| uInt status=0, ignore=0; /* status accumulators */ |
| uInt rstatus; /* .. */ |
| Int exp; /* working exponent */ |
| Int ideal; /* ideal (preferred) exponent */ |
| Int needbytes; /* work */ |
| Int dropped; /* .. */ |
| |
| #if DECSUBSET |
| decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ |
| #endif |
| /* buffer for f [needs +1 in case DECBUFFER 0] */ |
| decNumber buff[D2N(DECBUFFER+1)]; |
| /* buffer for a [needs +2 to match likely maxp] */ |
| decNumber bufa[D2N(DECBUFFER+2)]; |
| /* buffer for temporary, b [must be same size as a] */ |
| decNumber bufb[D2N(DECBUFFER+2)]; |
| decNumber *allocbuff=NULL; /* -> allocated buff, iff allocated */ |
| decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ |
| decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ |
| decNumber *f=buff; /* reduced fraction */ |
| decNumber *a=bufa; /* approximation to result */ |
| decNumber *b=bufb; /* intermediate result */ |
| /* buffer for temporary variable, up to 3 digits */ |
| decNumber buft[D2N(3)]; |
| decNumber *t=buft; /* up-to-3-digit constant or work */ |
| |
| #if DECCHECK |
| if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; |
| #endif |
| |
| do { /* protect allocated storage */ |
| #if DECSUBSET |
| if (!set->extended) { |
| /* reduce operand and set lostDigits status, as needed */ |
| if (rhs->digits>set->digits) { |
| allocrhs=decRoundOperand(rhs, set, &status); |
| if (allocrhs==NULL) break; |
| /* [Note: 'f' allocation below could reuse this buffer if */ |
| /* used, but as this is rare they are kept separate for clarity.] */ |
| rhs=allocrhs; |
| } |
| } |
| #endif |
| /* [following code does not require input rounding] */ |
| |
| /* handle infinities and NaNs */ |
| if (SPECIALARG) { |
| if (decNumberIsInfinite(rhs)) { /* an infinity */ |
| if (decNumberIsNegative(rhs)) status|=DEC_Invalid_operation; |
| else decNumberCopy(res, rhs); /* +Infinity */ |
| } |
| else decNaNs(res, rhs, NULL, set, &status); /* a NaN */ |
| break; |
| } |
| |
| /* calculate the ideal (preferred) exponent [floor(exp/2)] */ |
| /* [It would be nicer to write: ideal=rhs->exponent>>1, but this */ |
| /* generates a compiler warning. Generated code is the same.] */ |
| ideal=(rhs->exponent&~1)/2; /* target */ |
| |
| /* handle zeros */ |
| if (ISZERO(rhs)) { |
| decNumberCopy(res, rhs); /* could be 0 or -0 */ |
| res->exponent=ideal; /* use the ideal [safe] */ |
| /* use decFinish to clamp any out-of-range exponent, etc. */ |
| decFinish(res, set, &residue, &status); |
| break; |
| } |
| |
| /* any other -x is an oops */ |
| if (decNumberIsNegative(rhs)) { |
| status|=DEC_Invalid_operation; |
| break; |
| } |
| |
| /* space is needed for three working variables */ |
| /* f -- the same precision as the RHS, reduced to 0.01->0.99... */ |
| /* a -- Hull's approximation -- precision, when assigned, is */ |
| /* currentprecision+1 or the input argument precision, */ |
| /* whichever is larger (+2 for use as temporary) */ |
| /* b -- intermediate temporary result (same size as a) */ |
| /* if any is too long for local storage, then allocate */ |
| workp=MAXI(set->digits+1, rhs->digits); /* actual rounding precision */ |
| workp=MAXI(workp, 7); /* at least 7 for low cases */ |
| maxp=workp+2; /* largest working precision */ |
| |
| needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit); |
| if (needbytes>(Int)sizeof(buff)) { |
| allocbuff=(decNumber *)malloc(needbytes); |
| if (allocbuff==NULL) { /* hopeless -- abandon */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| f=allocbuff; /* use the allocated space */ |
| } |
| /* a and b both need to be able to hold a maxp-length number */ |
| needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit); |
| if (needbytes>(Int)sizeof(bufa)) { /* [same applies to b] */ |
| allocbufa=(decNumber *)malloc(needbytes); |
| allocbufb=(decNumber *)malloc(needbytes); |
| if (allocbufa==NULL || allocbufb==NULL) { /* hopeless */ |
| status|=DEC_Insufficient_storage; |
| break;} |
| a=allocbufa; /* use the allocated spaces */ |
| b=allocbufb; /* .. */ |
| } |
| |
| /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */ |
| decNumberCopy(f, rhs); |
| exp=f->exponent+f->digits; /* adjusted to Hull rules */ |
| f->exponent=-(f->digits); /* to range */ |
| |
| /* set up working context */ |
| decContextDefault(&workset, DEC_INIT_DECIMAL64); |
| workset.emax=DEC_MAX_EMAX; |
| workset.emin=DEC_MIN_EMIN; |
| |
| /* [Until further notice, no error is possible and status bits */ |
| /* (Rounded, etc.) should be ignored, not accumulated.] */ |
| |
| /* Calculate initial approximation, and allow for odd exponent */ |
| workset.digits=workp; /* p for initial calculation */ |
| t->bits=0; t->digits=3; |
| a->bits=0; a->digits=3; |
| if ((exp & 1)==0) { /* even exponent */
|