Google Git
Sign in
gnu / gcc / a8404c07e7fca388c02c39077865f7d5fa928430 / . / libdecnumber / decNumber.c
blob: 22ac359a6b085bfd4a5c74a06195653e77b645d2 [file] [log] [blame]
/* Decimal number arithmetic module for the decNumber C Library.
Copyright (C) 2005-2022 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, &copystat); /* 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; /* assume 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 */