aboutsummaryrefslogtreecommitdiff
path: root/libdecnumber/decBasic.c
diff options
context:
space:
mode:
authorDaniel Jacobowitz <drow@false.org>2007-10-11 18:15:02 +0000
committerDaniel Jacobowitz <drow@false.org>2007-10-11 18:15:02 +0000
commitf5bc1778c3150f6e9364a1af81f3ba2658f9e68c (patch)
tree70a9ff23817fe9ed2dadc01f6fd3a179bb7d18b9 /libdecnumber/decBasic.c
parent607ae575a784d7d02956031883ae407faf06fd89 (diff)
downloadgdb-f5bc1778c3150f6e9364a1af81f3ba2658f9e68c.zip
gdb-f5bc1778c3150f6e9364a1af81f3ba2658f9e68c.tar.gz
gdb-f5bc1778c3150f6e9364a1af81f3ba2658f9e68c.tar.bz2
* src-release (GDB_SUPPORT_DIRS): Add libdecnumber.
* libdecnumber: New directory, imported from GCC.
Diffstat (limited to 'libdecnumber/decBasic.c')
-rw-r--r--libdecnumber/decBasic.c3769
1 files changed, 3769 insertions, 0 deletions
diff --git a/libdecnumber/decBasic.c b/libdecnumber/decBasic.c
new file mode 100644
index 0000000..9ce277d
--- /dev/null
+++ b/libdecnumber/decBasic.c
@@ -0,0 +1,3769 @@
+/* Common base code for the decNumber C Library.
+ Copyright (C) 2007 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 2, or (at your option) any later
+ version.
+
+ In addition to the permissions in the GNU General Public License,
+ the Free Software Foundation gives you unlimited permission to link
+ the compiled version of this file into combinations with other
+ programs, and to distribute those combinations without any
+ restriction coming from the use of this file. (The General Public
+ License restrictions do apply in other respects; for example, they
+ cover modification of the file, and distribution when not linked
+ into a combine executable.)
+
+ 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.
+
+ You should have received a copy of the GNU General Public License
+ along with GCC; see the file COPYING. If not, write to the Free
+ Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301, USA. */
+
+/* ------------------------------------------------------------------ */
+/* decBasic.c -- common base code for Basic decimal types */
+/* ------------------------------------------------------------------ */
+/* This module comprises code that is shared between decDouble and */
+/* decQuad (but not decSingle). The main arithmetic operations are */
+/* here (Add, Subtract, Multiply, FMA, and Division operators). */
+/* */
+/* Unlike decNumber, parameterization takes place at compile time */
+/* rather than at runtime. The parameters are set in the decDouble.c */
+/* (etc.) files, which then include this one to produce the compiled */
+/* code. The functions here, therefore, are code shared between */
+/* multiple formats. */
+/* */
+/* This must be included after decCommon.c. */
+/* ------------------------------------------------------------------ */
+/* Names here refer to decFloat rather than to decDouble, etc., and */
+/* the functions are in strict alphabetical order. */
+
+/* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
+/* decCommon.c */
+#if !defined(QUAD)
+ #error decBasic.c must be included after decCommon.c
+#endif
+#if SINGLE
+ #error Routines in decBasic.c are for decDouble and decQuad only
+#endif
+
+/* Private constants */
+#define DIVIDE 0x80000000 /* Divide operations [as flags] */
+#define REMAINDER 0x40000000 /* .. */
+#define DIVIDEINT 0x20000000 /* .. */
+#define REMNEAR 0x10000000 /* .. */
+
+/* Private functions (local, used only by routines in this module) */
+static decFloat *decDivide(decFloat *, const decFloat *,
+ const decFloat *, decContext *, uInt);
+static decFloat *decCanonical(decFloat *, const decFloat *);
+static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
+ const decFloat *);
+static decFloat *decInfinity(decFloat *, const decFloat *);
+static decFloat *decInvalid(decFloat *, decContext *);
+static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
+ decContext *);
+static Int decNumCompare(const decFloat *, const decFloat *, Flag);
+static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
+ enum rounding, Flag);
+static uInt decToInt32(const decFloat *, decContext *, enum rounding,
+ Flag, Flag);
+
+/* ------------------------------------------------------------------ */
+/* decCanonical -- copy a decFloat, making canonical */
+/* */
+/* result gets the canonicalized df */
+/* df is the decFloat to copy and make canonical */
+/* returns result */
+/* */
+/* This is exposed via decFloatCanonical for Double and Quad only. */
+/* This works on specials, too; no error or exception is possible. */
+/* ------------------------------------------------------------------ */
+static decFloat * decCanonical(decFloat *result, const decFloat *df) {
+ uInt encode, precode, dpd; /* work */
+ uInt inword, uoff, canon; /* .. */
+ Int n; /* counter (down) */
+ if (df!=result) *result=*df; /* effect copy if needed */
+ if (DFISSPECIAL(result)) {
+ if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
+ /* is a NaN */
+ DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */
+ if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */
+ /* drop through to check payload */
+ }
+ /* return quickly if the coefficient continuation is canonical */
+ { /* declare block */
+ #if DOUBLE
+ uInt sourhi=DFWORD(df, 0);
+ uInt sourlo=DFWORD(df, 1);
+ if (CANONDPDOFF(sourhi, 8)
+ && CANONDPDTWO(sourhi, sourlo, 30)
+ && CANONDPDOFF(sourlo, 20)
+ && CANONDPDOFF(sourlo, 10)
+ && CANONDPDOFF(sourlo, 0)) return result;
+ #elif QUAD
+ uInt sourhi=DFWORD(df, 0);
+ uInt sourmh=DFWORD(df, 1);
+ uInt sourml=DFWORD(df, 2);
+ uInt sourlo=DFWORD(df, 3);
+ if (CANONDPDOFF(sourhi, 4)
+ && CANONDPDTWO(sourhi, sourmh, 26)
+ && CANONDPDOFF(sourmh, 16)
+ && CANONDPDOFF(sourmh, 6)
+ && CANONDPDTWO(sourmh, sourml, 28)
+ && CANONDPDOFF(sourml, 18)
+ && CANONDPDOFF(sourml, 8)
+ && CANONDPDTWO(sourml, sourlo, 30)
+ && CANONDPDOFF(sourlo, 20)
+ && CANONDPDOFF(sourlo, 10)
+ && CANONDPDOFF(sourlo, 0)) return result;
+ #endif
+ } /* block */
+
+ /* Loop to repair a non-canonical coefficent, as needed */
+ inword=DECWORDS-1; /* current input word */
+ uoff=0; /* bit offset of declet */
+ encode=DFWORD(result, inword);
+ for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */
+ dpd=encode>>uoff;
+ uoff+=10;
+ if (uoff>32) { /* crossed uInt boundary */
+ inword--;
+ encode=DFWORD(result, inword);
+ uoff-=32;
+ dpd|=encode<<(10-uoff); /* get pending bits */
+ }
+ dpd&=0x3ff; /* clear uninteresting bits */
+ if (dpd<0x16e) continue; /* must be canonical */
+ canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */
+ if (canon==dpd) continue; /* have canonical declet */
+ /* need to replace declet */
+ if (uoff>=10) { /* all within current word */
+ encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
+ encode|=canon<<(uoff-10); /* insert the canonical form */
+ DFWORD(result, inword)=encode; /* .. and save */
+ continue;
+ }
+ /* straddled words */
+ precode=DFWORD(result, inword+1); /* get previous */
+ precode&=0xffffffff>>(10-uoff); /* clear top bits */
+ DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
+ encode&=0xffffffff<<uoff; /* clear bottom bits */
+ encode|=canon>>(10-uoff); /* insert canonical */
+ DFWORD(result, inword)=encode; /* .. and save */
+ } /* n */
+ return result;
+ } /* decCanonical */
+
+/* ------------------------------------------------------------------ */
+/* decDivide -- divide operations */
+/* */
+/* result gets the result of dividing dfl by dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* op is the operation selector */
+/* returns result */
+/* */
+/* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
+/* ------------------------------------------------------------------ */
+#define DIVCOUNT 0 /* 1 to instrument subtractions counter */
+#define DIVBASE BILLION /* the base used for divide */
+#define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
+#define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */
+static decFloat * decDivide(decFloat *result, const decFloat *dfl,
+ const decFloat *dfr, decContext *set, uInt op) {
+ decFloat quotient; /* for remainders */
+ bcdnum num; /* for final conversion */
+ uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
+ uInt div[DIVOPLEN]; /* divisor in base-billion .. */
+ uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */
+ uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
+ uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */
+ Int divunits, accunits; /* lengths */
+ Int quodigits; /* digits in quotient */
+ uInt *lsua, *lsuq; /* -> current acc and quo lsus */
+ Int length, multiplier; /* work */
+ uInt carry, sign; /* .. */
+ uInt *ua, *ud, *uq; /* .. */
+ uByte *ub; /* .. */
+ uInt divtop; /* top unit of div adjusted for estimating */
+ #if DIVCOUNT
+ static uInt maxcount=0; /* worst-seen subtractions count */
+ uInt divcount=0; /* subtractions count [this divide] */
+ #endif
+
+ /* calculate sign */
+ num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
+
+ if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
+ /* NaNs are handled as usual */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* one or two infinities */
+ if (DFISINF(dfl)) {
+ if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
+ if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
+ /* Infinity/x is infinite and quiet, even if x=0 */
+ DFWORD(result, 0)=num.sign;
+ return decInfinity(result, result);
+ }
+ /* must be x/Infinity -- remainders are lhs */
+ if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
+ /* divides: return zero with correct sign and exponent depending */
+ /* on op (Etiny for divide, 0 for divideInt) */
+ decFloatZero(result);
+ if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
+ else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */
+ return result;
+ }
+ /* next, handle zero operands (x/0 and 0/x) */
+ if (DFISZERO(dfr)) { /* x/0 */
+ if (DFISZERO(dfl)) { /* 0/0 is undefined */
+ decFloatZero(result);
+ DFWORD(result, 0)=DECFLOAT_qNaN;
+ set->status|=DEC_Division_undefined;
+ return result;
+ }
+ if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
+ set->status|=DEC_Division_by_zero;
+ DFWORD(result, 0)=num.sign;
+ return decInfinity(result, result); /* x/0 -> signed Infinity */
+ }
+ num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
+ if (DFISZERO(dfl)) { /* 0/x (x!=0) */
+ /* if divide, result is 0 with ideal exponent; divideInt has */
+ /* exponent=0, remainders give zero with lower exponent */
+ if (op&DIVIDEINT) {
+ decFloatZero(result);
+ DFWORD(result, 0)|=num.sign; /* add sign */
+ return result;
+ }
+ if (!(op&DIVIDE)) { /* a remainder */
+ /* exponent is the minimum of the operands */
+ num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
+ /* if the result is zero the sign shall be sign of dfl */
+ num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
+ }
+ bcdacc[0]=0;
+ num.msd=bcdacc; /* -> 0 */
+ num.lsd=bcdacc; /* .. */
+ return decFinalize(result, &num, set); /* [divide may clamp exponent] */
+ } /* 0/x */
+ /* [here, both operands are known to be finite and non-zero] */
+
+ /* extract the operand coefficents into 'units' which are */
+ /* base-billion; the lhs is high-aligned in acc and the msu of both */
+ /* acc and div is at the right-hand end of array (offset length-1); */
+ /* the quotient can need one more unit than the operands as digits */
+ /* in it are not necessarily aligned neatly; further, the quotient */
+ /* may not start accumulating until after the end of the initial */
+ /* operand in acc if that is small (e.g., 1) so the accumulator */
+ /* must have at least that number of units extra (at the ls end) */
+ GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
+ GETCOEFFBILL(dfr, div);
+ /* zero the low uInts of acc */
+ acc[0]=0;
+ acc[1]=0;
+ acc[2]=0;
+ acc[3]=0;
+ #if DOUBLE
+ #if DIVOPLEN!=2
+ #error Unexpected Double DIVOPLEN
+ #endif
+ #elif QUAD
+ acc[4]=0;
+ acc[5]=0;
+ acc[6]=0;
+ acc[7]=0;
+ #if DIVOPLEN!=4
+ #error Unexpected Quad DIVOPLEN
+ #endif
+ #endif
+
+ /* set msu and lsu pointers */
+ msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
+ msuq=quo+DIVOPLEN;
+ /*[loop for div will terminate because operands are non-zero] */
+ for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
+ /* the initial least-significant unit of acc is set so acc appears */
+ /* to have the same length as div. */
+ /* This moves one position towards the least possible for each */
+ /* iteration */
+ divunits=(Int)(msud-div+1); /* precalculate */
+ lsua=msua-divunits+1; /* initial working lsu of acc */
+ lsuq=msuq; /* and of quo */
+
+ /* set up the estimator for the multiplier; this is the msu of div, */
+ /* plus two bits from the unit below (if any) rounded up by one if */
+ /* there are any non-zero bits or units below that [the extra two */
+ /* bits makes for a much better estimate when the top unit is small] */
+ divtop=*msud<<2;
+ if (divunits>1) {
+ uInt *um=msud-1;
+ uInt d=*um;
+ if (d>=750000000) {divtop+=3; d-=750000000;}
+ else if (d>=500000000) {divtop+=2; d-=500000000;}
+ else if (d>=250000000) {divtop++; d-=250000000;}
+ if (d) divtop++;
+ else for (um--; um>=div; um--) if (*um) {
+ divtop++;
+ break;
+ }
+ } /* >1 unit */
+
+ #if DECTRACE
+ {Int i;
+ printf("----- div=");
+ for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
+ printf("\n");}
+ #endif
+
+ /* now collect up to DECPMAX+1 digits in the quotient (this may */
+ /* need OPLEN+1 uInts if unaligned) */
+ quodigits=0; /* no digits yet */
+ for (;; lsua--) { /* outer loop -- each input position */
+ #if DECCHECK
+ if (lsua<acc) {
+ printf("Acc underrun...\n");
+ break;
+ }
+ #endif
+ #if DECTRACE
+ printf("Outer: quodigits=%ld acc=", (LI)quodigits);
+ for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
+ printf("\n");
+ #endif
+ *lsuq=0; /* default unit result is 0 */
+ for (;;) { /* inner loop -- calculate quotient unit */
+ /* strip leading zero units from acc (either there initially or */
+ /* from subtraction below); this may strip all if exactly 0 */
+ for (; *msua==0 && msua>=lsua;) msua--;
+ accunits=(Int)(msua-lsua+1); /* [maybe 0] */
+ /* subtraction is only necessary and possible if there are as */
+ /* least as many units remaining in acc for this iteration as */
+ /* there are in div */
+ if (accunits<divunits) {
+ if (accunits==0) msua++; /* restore */
+ break;
+ }
+
+ /* If acc is longer than div then subtraction is definitely */
+ /* possible (as msu of both is non-zero), but if they are the */
+ /* same length a comparison is needed. */
+ /* If a subtraction is needed then a good estimate of the */
+ /* multiplier for the subtraction is also needed in order to */
+ /* minimise the iterations of this inner loop because the */
+ /* subtractions needed dominate division performance. */
+ if (accunits==divunits) {
+ /* compare the high divunits of acc and div: */
+ /* acc<div: this quotient unit is unchanged; subtraction */
+ /* will be possible on the next iteration */
+ /* acc==div: quotient gains 1, set acc=0 */
+ /* acc>div: subtraction necessary at this position */
+ for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
+ /* [now at first mismatch or lsu] */
+ if (*ud>*ua) break; /* next time... */
+ if (*ud==*ua) { /* all compared equal */
+ *lsuq+=1; /* increment result */
+ msua=lsua; /* collapse acc units */
+ *msua=0; /* .. to a zero */
+ break;
+ }
+
+ /* subtraction necessary; estimate multiplier [see above] */
+ /* if both *msud and *msua are small it is cost-effective to */
+ /* bring in part of the following units (if any) to get a */
+ /* better estimate (assume some other non-zero in div) */
+ #define DIVLO 1000000U
+ #define DIVHI (DIVBASE/DIVLO)
+ #if DECUSE64
+ if (divunits>1) {
+ /* there cannot be a *(msud-2) for DECDOUBLE so next is */
+ /* an exact calculation unless DECQUAD (which needs to */
+ /* assume bits out there if divunits>2) */
+ uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
+ uLong div=(uLong)*msud * DIVBASE + *(msud-1);
+ #if QUAD
+ if (divunits>2) div++;
+ #endif
+ mul/=div;
+ multiplier=(Int)mul;
+ }
+ else multiplier=*msua/(*msud);
+ #else
+ if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
+ multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
+ /(*msud*DIVHI + *(msud-1)/DIVLO +1);
+ }
+ else multiplier=(*msua<<2)/divtop;
+ #endif
+ }
+ else { /* accunits>divunits */
+ /* msud is one unit 'lower' than msua, so estimate differently */
+ #if DECUSE64
+ uLong mul;
+ /* as before, bring in extra digits if possible */
+ if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
+ mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
+ + *(msua-2)/DIVLO;
+ mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
+ }
+ else if (divunits==1) {
+ mul=(uLong)*msua * DIVBASE + *(msua-1);
+ mul/=*msud; /* no more to the right */
+ }
+ else {
+ mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) + (*(msua-1)<<2);
+ mul/=divtop; /* [divtop already allows for sticky bits] */
+ }
+ multiplier=(Int)mul;
+ #else
+ multiplier=*msua * ((DIVBASE<<2)/divtop);
+ #endif
+ }
+ if (multiplier==0) multiplier=1; /* marginal case */
+ *lsuq+=multiplier;
+
+ #if DIVCOUNT
+ /* printf("Multiplier: %ld\n", (LI)multiplier); */
+ divcount++;
+ #endif
+
+ /* Carry out the subtraction acc-(div*multiplier); for each */
+ /* unit in div, do the multiply, split to units (see */
+ /* decFloatMultiply for the algorithm), and subtract from acc */
+ #define DIVMAGIC 2305843009U /* 2**61/10**9 */
+ #define DIVSHIFTA 29
+ #define DIVSHIFTB 32
+ carry=0;
+ for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
+ uInt lo, hop;
+ #if DECUSE64
+ uLong sub=(uLong)multiplier*(*ud)+carry;
+ if (sub<DIVBASE) {
+ carry=0;
+ lo=(uInt)sub;
+ }
+ else {
+ hop=(uInt)(sub>>DIVSHIFTA);
+ carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
+ /* the estimate is now in hi; now calculate sub-hi*10**9 */
+ /* to get the remainder (which will be <DIVBASE)) */
+ lo=(uInt)sub;
+ lo-=carry*DIVBASE; /* low word of result */
+ if (lo>=DIVBASE) {
+ lo-=DIVBASE; /* correct by +1 */
+ carry++;
+ }
+ }
+ #else /* 32-bit */
+ uInt hi;
+ /* calculate multiplier*(*ud) into hi and lo */
+ LONGMUL32HI(hi, *ud, multiplier); /* get the high word */
+ lo=multiplier*(*ud); /* .. and the low */
+ lo+=carry; /* add the old hi */
+ carry=hi+(lo<carry); /* .. with any carry */
+ if (carry || lo>=DIVBASE) { /* split is needed */
+ hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */
+ LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */
+ /* [DIVSHIFTB is 32, so carry can be used directly] */
+ /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
+ /* happily the top word of the result is irrelevant because it */
+ /* will always be zero so this needs only one multiplication */
+ lo-=(carry*DIVBASE);
+ /* the correction here will be at most +1; do it */
+ if (lo>=DIVBASE) {
+ lo-=DIVBASE;
+ carry++;
+ }
+ }
+ #endif
+ if (lo>*ua) { /* borrow needed */
+ *ua+=DIVBASE;
+ carry++;
+ }
+ *ua-=lo;
+ } /* ud loop */
+ if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */
+ } /* inner loop */
+
+ /* the outer loop terminates when there is either an exact result */
+ /* or enough digits; first update the quotient digit count and */
+ /* pointer (if any significant digits) */
+ #if DECTRACE
+ if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
+ #endif
+ if (quodigits) {
+ quodigits+=9; /* had leading unit earlier */
+ lsuq--;
+ if (quodigits>DECPMAX+1) break; /* have enough */
+ }
+ else if (*lsuq) { /* first quotient digits */
+ const uInt *pow;
+ for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
+ lsuq--;
+ /* [cannot have >DECPMAX+1 on first unit] */
+ }
+
+ if (*msua!=0) continue; /* not an exact result */
+ /* acc is zero iff used all of original units and zero down to lsua */
+ /* (must also continue to original lsu for correct quotient length) */
+ if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
+ for (; msua>lsua && *msua==0;) msua--;
+ if (*msua==0 && msua==lsua) break;
+ } /* outer loop */
+
+ /* all of the original operand in acc has been covered at this point */
+ /* quotient now has at least DECPMAX+2 digits */
+ /* *msua is now non-0 if inexact and sticky bits */
+ /* lsuq is one below the last uint of the quotient */
+ lsuq++; /* set -> true lsu of quo */
+ if (*msua) *lsuq|=1; /* apply sticky bit */
+
+ /* quo now holds the (unrounded) quotient in base-billion; one */
+ /* base-billion 'digit' per uInt. */
+ #if DECTRACE
+ printf("DivQuo:");
+ for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
+ printf("\n");
+ #endif
+
+ /* Now convert to BCD for rounding and cleanup, starting from the */
+ /* most significant end [offset by one into bcdacc to leave room */
+ /* for a possible carry digit if rounding for REMNEAR is needed] */
+ for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
+ uInt top, mid, rem; /* work */
+ if (*uq==0) { /* no split needed */
+ UINTAT(ub)=0; /* clear 9 BCD8s */
+ UINTAT(ub+4)=0; /* .. */
+ *(ub+8)=0; /* .. */
+ continue;
+ }
+ /* *uq is non-zero -- split the base-billion digit into */
+ /* hi, mid, and low three-digits */
+ #define divsplit9 1000000 /* divisor */
+ #define divsplit6 1000 /* divisor */
+ /* The splitting is done by simple divides and remainders, */
+ /* assuming the compiler will optimize these [GCC does] */
+ top=*uq/divsplit9;
+ rem=*uq%divsplit9;
+ mid=rem/divsplit6;
+ rem=rem%divsplit6;
+ /* lay out the nine BCD digits (plus one unwanted byte) */
+ UINTAT(ub) =UINTAT(&BIN2BCD8[top*4]);
+ UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
+ UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
+ } /* BCD conversion loop */
+ ub--; /* -> lsu */
+
+ /* complete the bcdnum; quodigits is correct, so the position of */
+ /* the first non-zero is known */
+ num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
+ num.lsd=ub;
+
+ /* make exponent adjustments, etc */
+ if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */
+ num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
+ /* if the result was exact then there may be up to 8 extra */
+ /* trailing zeros in the overflowed quotient final unit */
+ if (*msua==0) {
+ for (; *ub==0;) ub--; /* drop zeros */
+ num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */
+ num.lsd=ub;
+ }
+ } /* adjustment needed */
+
+ #if DIVCOUNT
+ if (divcount>maxcount) { /* new high-water nark */
+ maxcount=divcount;
+ printf("DivNewMaxCount: %ld\n", (LI)maxcount);
+ }
+ #endif
+
+ if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
+
+ /* Is DIVIDEINT or a remainder; there is more to do -- first form */
+ /* the integer (this is done 'after the fact', unlike as in */
+ /* decNumber, so as not to tax DIVIDE) */
+
+ /* The first non-zero digit will be in the first 9 digits, known */
+ /* from quodigits and num.msd, so there is always space for DECPMAX */
+ /* digits */
+
+ length=(Int)(num.lsd-num.msd+1);
+ /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
+
+ if (length+num.exponent>DECPMAX) { /* cannot fit */
+ decFloatZero(result);
+ DFWORD(result, 0)=DECFLOAT_qNaN;
+ set->status|=DEC_Division_impossible;
+ return result;
+ }
+
+ if (num.exponent>=0) { /* already an int, or need pad zeros */
+ for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
+ num.lsd+=num.exponent;
+ }
+ else { /* too long: round or truncate needed */
+ Int drop=-num.exponent;
+ if (!(op&REMNEAR)) { /* simple truncate */
+ num.lsd-=drop;
+ if (num.lsd<num.msd) { /* truncated all */
+ num.lsd=num.msd; /* make 0 */
+ *num.lsd=0; /* .. [sign still relevant] */
+ }
+ }
+ else { /* round to nearest even [sigh] */
+ /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
+ /* (this is a special case of Quantize -- q.v. for commentary) */
+ uByte *roundat; /* -> re-round digit */
+ uByte reround; /* reround value */
+ *(num.msd-1)=0; /* in case of left carry, or make 0 */
+ if (drop<length) roundat=num.lsd-drop+1;
+ else if (drop==length) roundat=num.msd;
+ else roundat=num.msd-1; /* [-> 0] */
+ reround=*roundat;
+ for (ub=roundat+1; ub<=num.lsd; ub++) {
+ if (*ub!=0) {
+ reround=DECSTICKYTAB[reround];
+ break;
+ }
+ } /* check stickies */
+ if (roundat>num.msd) num.lsd=roundat-1;
+ else {
+ num.msd--; /* use the 0 .. */
+ num.lsd=num.msd; /* .. at the new MSD place */
+ }
+ if (reround!=0) { /* discarding non-zero */
+ uInt bump=0;
+ /* rounding is DEC_ROUND_HALF_EVEN always */
+ if (reround>5) bump=1; /* >0.5 goes up */
+ else if (reround==5) /* exactly 0.5000 .. */
+ bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */
+ if (bump!=0) { /* need increment */
+ /* increment the coefficient; this might end up with 1000... */
+ ub=num.lsd;
+ for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
+ for (; *ub==9; ub--) *ub=0; /* at most 3 more */
+ *ub+=1;
+ if (ub<num.msd) num.msd--; /* carried */
+ } /* bump needed */
+ } /* reround!=0 */
+ } /* remnear */
+ } /* round or truncate needed */
+ num.exponent=0; /* all paths */
+ /*decShowNum(&num, "int"); */
+
+ if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
+
+ /* Have a remainder to calculate */
+ decFinalize(&quotient, &num, set); /* lay out the integer so far */
+ DFWORD(&quotient, 0)^=DECFLOAT_Sign; /* negate it */
+ sign=DFWORD(dfl, 0); /* save sign of dfl */
+ decFloatFMA(result, &quotient, dfr, dfl, set);
+ if (!DFISZERO(result)) return result;
+ /* if the result is zero the sign shall be sign of dfl */
+ DFWORD(&quotient, 0)=sign; /* construct decFloat of sign */
+ return decFloatCopySign(result, result, &quotient);
+ } /* decDivide */
+
+/* ------------------------------------------------------------------ */
+/* decFiniteMultiply -- multiply two finite decFloats */
+/* */
+/* num gets the result of multiplying dfl and dfr */
+/* bcdacc .. with the coefficient in this array */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* */
+/* This effects the multiplication of two decFloats, both known to be */
+/* finite, leaving the result in a bcdnum ready for decFinalize (for */
+/* use in Multiply) or in a following addition (FMA). */
+/* */
+/* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
+/* No error is possible and no status is set. */
+/* ------------------------------------------------------------------ */
+/* This routine has two separate implementations of the core */
+/* multiplication; both using base-billion. One uses only 32-bit */
+/* variables (Ints and uInts) or smaller; the other uses uLongs (for */
+/* multiplication and addition only). Both implementations cover */
+/* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
+/* comparisons. In any one compilation only one implementation for */
+/* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
+/* version is forced. */
+/* */
+/* Historical note: an earlier version of this code also supported the */
+/* 256-bit format and has been preserved. That is somewhat trickier */
+/* during lazy carry splitting because the initial quotient estimate */
+/* (est) can exceed 32 bits. */
+
+#define MULTBASE BILLION /* the base used for multiply */
+#define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
+#define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */
+#define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
+
+/* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
+#if DECEMAXD>9
+ #error Exponent may overflow when doubled for Multiply
+#endif
+#if MULACCLEN!=(MULACCLEN/4)*4
+ /* This assumption is used below only for initialization */
+ #error MULACCLEN is not a multiple of 4
+#endif
+
+static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
+ const decFloat *dfl, const decFloat *dfr) {
+ uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
+ uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
+ uInt *ui, *uj; /* work */
+ uByte *ub; /* .. */
+
+ #if DECUSE64
+ uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
+ uLong *pl; /* work -> lazy accumulator */
+ uInt acc[MULACCLEN]; /* coefficent in base-billion .. */
+ #else
+ uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */
+ #endif
+ uInt *pa; /* work -> accumulator */
+ /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
+
+ /* Calculate sign and exponent */
+ num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
+ num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
+
+ /* Extract the coefficients and prepare the accumulator */
+ /* the coefficients of the operands are decoded into base-billion */
+ /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
+ /* appropriate size. */
+ GETCOEFFBILL(dfl, bufl);
+ GETCOEFFBILL(dfr, bufr);
+ #if DECTRACE && 0
+ printf("CoeffbL:");
+ for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
+ printf("\n");
+ printf("CoeffbR:");
+ for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
+ printf("\n");
+ #endif
+
+ /* start the 64-bit/32-bit differing paths... */
+#if DECUSE64
+
+ /* zero the accumulator */
+ #if MULACCLEN==4
+ accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
+ #else /* use a loop */
+ /* MULACCLEN is a multiple of four, asserted above */
+ for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
+ *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
+ } /* pl */
+ #endif
+
+ /* Effect the multiplication */
+ /* The multiplcation proceeds using MFC's lazy-carry resolution */
+ /* algorithm from decNumber. First, the multiplication is */
+ /* effected, allowing accumulation of the partial products (which */
+ /* are in base-billion at each column position) into 64 bits */
+ /* without resolving back to base=billion after each addition. */
+ /* These 64-bit numbers (which may contain up to 19 decimal digits) */
+ /* are then split using the Clark & Cowlishaw algorithm (see below). */
+ /* [Testing for 0 in the inner loop is not really a 'win'] */
+ for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
+ if (*ui==0) continue; /* product cannot affect result */
+ pl=accl+(ui-bufr); /* where to add the lhs */
+ for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
+ /* if (*uj==0) continue; // product cannot affect result */
+ *pl+=((uLong)*ui)*(*uj);
+ } /* uj */
+ } /* ui */
+
+ /* The 64-bit carries must now be resolved; this means that a */
+ /* quotient/remainder has to be calculated for base-billion (1E+9). */
+ /* For this, Clark & Cowlishaw's quotient estimation approach (also */
+ /* used in decNumber) is needed, because 64-bit divide is generally */
+ /* extremely slow on 32-bit machines, and may be slower than this */
+ /* approach even on 64-bit machines. This algorithm splits X */
+ /* using: */
+ /* */
+ /* magic=2**(A+B)/1E+9; // 'magic number' */
+ /* hop=X/2**A; // high order part of X (by shift) */
+ /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
+ /* */
+ /* A and B are quite constrained; hop and magic must fit in 32 bits, */
+ /* and 2**(A+B) must be as large as possible (which is 2**61 if */
+ /* magic is to fit). Further, maxX increases with the length of */
+ /* the operands (and hence the number of partial products */
+ /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
+ /* */
+ /* It can be shown that when OPLEN is 2 then the maximum error in */
+ /* the estimated quotient is <1, but for larger maximum x the */
+ /* maximum error is above 1 so a correction that is >1 may be */
+ /* needed. Values of A and B are chosen to satisfy the constraints */
+ /* just mentioned while minimizing the maximum error (and hence the */
+ /* maximum correction), as shown in the following table: */
+ /* */
+ /* Type OPLEN A B maxX maxError maxCorrection */
+ /* --------------------------------------------------------- */
+ /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
+ /* QUAD 4 30 31 <4*10**18 1.17 2 */
+ /* */
+ /* In the OPLEN==2 case there is most choice, but the value for B */
+ /* of 32 has a big advantage as then the calculation of the */
+ /* estimate requires no shifting; the compiler can extract the high */
+ /* word directly after multiplying magic*hop. */
+ #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
+ #if DOUBLE
+ #define MULSHIFTA 29
+ #define MULSHIFTB 32
+ #elif QUAD
+ #define MULSHIFTA 30
+ #define MULSHIFTB 31
+ #else
+ #error Unexpected type
+ #endif
+
+ #if DECTRACE
+ printf("MulAccl:");
+ for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
+ printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
+ printf("\n");
+ #endif
+
+ for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
+ uInt lo, hop; /* work */
+ uInt est; /* cannot exceed 4E+9 */
+ if (*pl>MULTBASE) {
+ /* *pl holds a binary number which needs to be split */
+ hop=(uInt)(*pl>>MULSHIFTA);
+ est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
+ /* the estimate is now in est; now calculate hi:lo-est*10**9; */
+ /* happily the top word of the result is irrelevant because it */
+ /* will always be zero so this needs only one multiplication */
+ lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */
+ /* If QUAD, the correction here could be +2 */
+ if (lo>=MULTBASE) {
+ lo-=MULTBASE; /* correct by +1 */
+ est++;
+ #if QUAD
+ /* may need to correct by +2 */
+ if (lo>=MULTBASE) {
+ lo-=MULTBASE;
+ est++;
+ }
+ #endif
+ }
+ /* finally place lo as the new coefficient 'digit' and add est to */
+ /* the next place up [this is safe because this path is never */
+ /* taken on the final iteration as *pl will fit] */
+ *pa=lo;
+ *(pl+1)+=est;
+ } /* *pl needed split */
+ else { /* *pl<MULTBASE */
+ *pa=(uInt)*pl; /* just copy across */
+ }
+ } /* pl loop */
+
+#else /* 32-bit */
+ for (pa=acc;; pa+=4) { /* zero the accumulator */
+ *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */
+ if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */
+ } /* pa */
+
+ /* Effect the multiplication */
+ /* uLongs are not available (and in particular, there is no uLong */
+ /* divide) but it is still possible to use MFC's lazy-carry */
+ /* resolution algorithm from decNumber. First, the multiplication */
+ /* is effected, allowing accumulation of the partial products */
+ /* (which are in base-billion at each column position) into 64 bits */
+ /* [with the high-order 32 bits in each position being held at */
+ /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
+ /* These 64-bit numbers (which may contain up to 19 decimal digits) */
+ /* are then split using the Clark & Cowlishaw algorithm (see */
+ /* below). */
+ for (ui=bufr;; ui++) { /* over each item in rhs */
+ uInt hi, lo; /* words of exact multiply result */
+ pa=acc+(ui-bufr); /* where to add the lhs */
+ for (uj=bufl;; uj++, pa++) { /* over each item in lhs */
+ LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */
+ lo=(*ui)*(*uj); /* .. */
+ *pa+=lo; /* accumulate low bits and .. */
+ *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */
+ if (uj==bufl+MULOPLEN-1) break;
+ }
+ if (ui==bufr+MULOPLEN-1) break;
+ }
+
+ /* The 64-bit carries must now be resolved; this means that a */
+ /* quotient/remainder has to be calculated for base-billion (1E+9). */
+ /* For this, Clark & Cowlishaw's quotient estimation approach (also */
+ /* used in decNumber) is needed, because 64-bit divide is generally */
+ /* extremely slow on 32-bit machines. This algorithm splits X */
+ /* using: */
+ /* */
+ /* magic=2**(A+B)/1E+9; // 'magic number' */
+ /* hop=X/2**A; // high order part of X (by shift) */
+ /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
+ /* */
+ /* A and B are quite constrained; hop and magic must fit in 32 bits, */
+ /* and 2**(A+B) must be as large as possible (which is 2**61 if */
+ /* magic is to fit). Further, maxX increases with the length of */
+ /* the operands (and hence the number of partial products */
+ /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
+ /* */
+ /* It can be shown that when OPLEN is 2 then the maximum error in */
+ /* the estimated quotient is <1, but for larger maximum x the */
+ /* maximum error is above 1 so a correction that is >1 may be */
+ /* needed. Values of A and B are chosen to satisfy the constraints */
+ /* just mentioned while minimizing the maximum error (and hence the */
+ /* maximum correction), as shown in the following table: */
+ /* */
+ /* Type OPLEN A B maxX maxError maxCorrection */
+ /* --------------------------------------------------------- */
+ /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
+ /* QUAD 4 30 31 <4*10**18 1.17 2 */
+ /* */
+ /* In the OPLEN==2 case there is most choice, but the value for B */
+ /* of 32 has a big advantage as then the calculation of the */
+ /* estimate requires no shifting; the high word is simply */
+ /* calculated from multiplying magic*hop. */
+ #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
+ #if DOUBLE
+ #define MULSHIFTA 29
+ #define MULSHIFTB 32
+ #elif QUAD
+ #define MULSHIFTA 30
+ #define MULSHIFTB 31
+ #else
+ #error Unexpected type
+ #endif
+
+ #if DECTRACE
+ printf("MulHiLo:");
+ for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
+ printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
+ printf("\n");
+ #endif
+
+ for (pa=acc;; pa++) { /* each low uInt */
+ uInt hi, lo; /* words of exact multiply result */
+ uInt hop, estlo; /* work */
+ #if QUAD
+ uInt esthi; /* .. */
+ #endif
+
+ lo=*pa;
+ hi=*(pa+MULACCLEN); /* top 32 bits */
+ /* hi and lo now hold a binary number which needs to be split */
+
+ #if DOUBLE
+ hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */
+ LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
+ /* [MULSHIFTB is 32, so estlo can be used directly] */
+ /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
+ /* happily the top word of the result is irrelevant because it */
+ /* will always be zero so this needs only one multiplication */
+ lo-=(estlo*MULTBASE);
+ /* esthi=0; // high word is ignored below */
+ /* the correction here will be at most +1; do it */
+ if (lo>=MULTBASE) {
+ lo-=MULTBASE;
+ estlo++;
+ }
+ #elif QUAD
+ hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */
+ LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
+ estlo=hop*MULMAGIC; /* .. so low word needed */
+ estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
+ /* esthi=0; // high word is ignored below */
+ lo-=(estlo*MULTBASE); /* as above */
+ /* the correction here could be +1 or +2 */
+ if (lo>=MULTBASE) {
+ lo-=MULTBASE;
+ estlo++;
+ }
+ if (lo>=MULTBASE) {
+ lo-=MULTBASE;
+ estlo++;
+ }
+ #else
+ #error Unexpected type
+ #endif
+
+ /* finally place lo as the new accumulator digit and add est to */
+ /* the next place up; this latter add could cause a carry of 1 */
+ /* to the high word of the next place */
+ *pa=lo;
+ *(pa+1)+=estlo;
+ /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
+ /* *(pa+1+MULACCLEN)+=esthi; */
+ if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
+ if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */
+ } /* pa loop */
+#endif
+
+ /* At this point, whether using the 64-bit or the 32-bit paths, the */
+ /* accumulator now holds the (unrounded) result in base-billion; */
+ /* one base-billion 'digit' per uInt. */
+ #if DECTRACE
+ printf("MultAcc:");
+ for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
+ printf("\n");
+ #endif
+
+ /* Now convert to BCD for rounding and cleanup, starting from the */
+ /* most significant end */
+ pa=acc+MULACCLEN-1;
+ if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
+ else { /* >=1 word of leading zeros */
+ num->msd=bcdacc; /* known leading zeros are gone */
+ pa--; /* skip first word .. */
+ for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
+ }
+ for (ub=bcdacc;; pa--, ub+=9) {
+ if (*pa!=0) { /* split(s) needed */
+ uInt top, mid, rem; /* work */
+ /* *pa is non-zero -- split the base-billion acc digit into */
+ /* hi, mid, and low three-digits */
+ #define mulsplit9 1000000 /* divisor */
+ #define mulsplit6 1000 /* divisor */
+ /* The splitting is done by simple divides and remainders, */
+ /* assuming the compiler will optimize these where useful */
+ /* [GCC does] */
+ top=*pa/mulsplit9;
+ rem=*pa%mulsplit9;
+ mid=rem/mulsplit6;
+ rem=rem%mulsplit6;
+ /* lay out the nine BCD digits (plus one unwanted byte) */
+ UINTAT(ub) =UINTAT(&BIN2BCD8[top*4]);
+ UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
+ UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
+ }
+ else { /* *pa==0 */
+ UINTAT(ub)=0; /* clear 9 BCD8s */
+ UINTAT(ub+4)=0; /* .. */
+ *(ub+8)=0; /* .. */
+ }
+ if (pa==acc) break;
+ } /* BCD conversion loop */
+
+ num->lsd=ub+8; /* complete the bcdnum .. */
+
+ #if DECTRACE
+ decShowNum(num, "postmult");
+ decFloatShow(dfl, "dfl");
+ decFloatShow(dfr, "dfr");
+ #endif
+ return;
+ } /* decFiniteMultiply */
+
+/* ------------------------------------------------------------------ */
+/* decFloatAbs -- absolute value, heeding NaNs, etc. */
+/* */
+/* result gets the canonicalized df with sign 0 */
+/* df is the decFloat to abs */
+/* set is the context */
+/* returns result */
+/* */
+/* This has the same effect as decFloatPlus unless df is negative, */
+/* in which case it has the same effect as decFloatMinus. The */
+/* effect is also the same as decFloatCopyAbs except that NaNs are */
+/* handled normally (the sign of a NaN is not affected, and an sNaN */
+/* will signal) and the result will be canonical. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatAbs(decFloat *result, const decFloat *df,
+ decContext *set) {
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
+ decCanonical(result, df); /* copy and check */
+ DFBYTE(result, 0)&=~0x80; /* zero sign bit */
+ return result;
+ } /* decFloatAbs */
+
+/* ------------------------------------------------------------------ */
+/* decFloatAdd -- add two decFloats */
+/* */
+/* result gets the result of adding dfl and dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatAdd(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ bcdnum num; /* for final conversion */
+ Int expl, expr; /* left and right exponents */
+ uInt *ui, *uj; /* work */
+ uByte *ub; /* .. */
+
+ uInt sourhil, sourhir; /* top words from source decFloats */
+ /* [valid only until specials */
+ /* handled or exponents decoded] */
+ uInt diffsign; /* non-zero if signs differ */
+ uInt carry; /* carry: 0 or 1 before add loop */
+ Int overlap; /* coefficient overlap (if full) */
+ /* the following buffers hold coefficients with various alignments */
+ /* (see commentary and diagrams below) */
+ uByte acc[4+2+DECPMAX*3+8];
+ uByte buf[4+2+DECPMAX*2];
+ uByte *umsd, *ulsd; /* local MSD and LSD pointers */
+
+ #if DECLITEND
+ #define CARRYPAT 0x01000000 /* carry=1 pattern */
+ #else
+ #define CARRYPAT 0x00000001 /* carry=1 pattern */
+ #endif
+
+ /* Start decoding the arguments */
+ /* the initial exponents are placed into the opposite Ints to */
+ /* that which might be expected; there are two sets of data to */
+ /* keep track of (each decFloat and the corresponding exponent), */
+ /* and this scheme means that at the swap point (after comparing */
+ /* exponents) only one pair of words needs to be swapped */
+ /* whichever path is taken (thereby minimising worst-case path) */
+ sourhil=DFWORD(dfl, 0); /* LHS top word */
+ expr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
+ sourhir=DFWORD(dfr, 0); /* RHS top word */
+ expl=DECCOMBEXP[sourhir>>26];
+
+ diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
+
+ if (EXPISSPECIAL(expl | expr)) { /* either is special? */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* one or two infinities */
+ /* two infinities with different signs is invalid */
+ if (diffsign && DFISINF(dfl) && DFISINF(dfr))
+ return decInvalid(result, set);
+ if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
+ return decInfinity(result, dfr); /* RHS must be Infinite */
+ }
+
+ /* Here when both arguments are finite */
+
+ /* complete exponent gathering (keeping swapped) */
+ expr+=GETECON(dfl)-DECBIAS; /* .. + continuation and unbias */
+ expl+=GETECON(dfr)-DECBIAS;
+ /* here expr has exponent from lhs, and vice versa */
+
+ /* now swap either exponents or argument pointers */
+ if (expl<=expr) {
+ /* original left is bigger */
+ Int expswap=expl;
+ expl=expr;
+ expr=expswap;
+ /* printf("left bigger\n"); */
+ }
+ else {
+ const decFloat *dfswap=dfl;
+ dfl=dfr;
+ dfr=dfswap;
+ /* printf("right bigger\n"); */
+ }
+ /* [here dfl and expl refer to the datum with the larger exponent, */
+ /* of if the exponents are equal then the original LHS argument] */
+
+ /* if lhs is zero then result will be the rhs (now known to have */
+ /* the smaller exponent), which also may need to be tested for zero */
+ /* for the weird IEEE 754 sign rules */
+ if (DFISZERO(dfl)) {
+ decCanonical(result, dfr); /* clean copy */
+ /* "When the sum of two operands with opposite signs is */
+ /* exactly zero, the sign of that sum shall be '+' in all */
+ /* rounding modes except round toward -Infinity, in which */
+ /* mode that sign shall be '-'." */
+ if (diffsign && DFISZERO(result)) {
+ DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */
+ if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
+ }
+ return result;
+ } /* numfl is zero */
+ /* [here, LHS is non-zero; code below assumes that] */
+
+ /* Coefficients layout during the calculations to follow: */
+ /* */
+ /* Overlap case: */
+ /* +------------------------------------------------+ */
+ /* acc: |0000| coeffa | tail B | | */
+ /* +------------------------------------------------+ */
+ /* buf: |0000| pad0s | coeffb | | */
+ /* +------------------------------------------------+ */
+ /* */
+ /* Touching coefficients or gap: */
+ /* +------------------------------------------------+ */
+ /* acc: |0000| coeffa | gap | coeffb | */
+ /* +------------------------------------------------+ */
+ /* [buf not used or needed; gap clamped to Pmax] */
+
+ /* lay out lhs coefficient into accumulator; this starts at acc+4 */
+ /* for decDouble or acc+6 for decQuad so the LSD is word- */
+ /* aligned; the top word gap is there only in case a carry digit */
+ /* is prefixed after the add -- it does not need to be zeroed */
+ #if DOUBLE
+ #define COFF 4 /* offset into acc */
+ #elif QUAD
+ USHORTAT(acc+4)=0; /* prefix 00 */
+ #define COFF 6 /* offset into acc */
+ #endif
+
+ GETCOEFF(dfl, acc+COFF); /* decode from decFloat */
+ ulsd=acc+COFF+DECPMAX-1;
+ umsd=acc+4; /* [having this here avoids */
+ /* weird GCC optimizer failure] */
+ #if DECTRACE
+ {bcdnum tum;
+ tum.msd=umsd;
+ tum.lsd=ulsd;
+ tum.exponent=expl;
+ tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
+ decShowNum(&tum, "dflx");}
+ #endif
+
+ /* if signs differ, take ten's complement of lhs (here the */
+ /* coefficient is subtracted from all-nines; the 1 is added during */
+ /* the later add cycle -- zeros to the right do not matter because */
+ /* the complement of zero is zero); these are fixed-length inverts */
+ /* where the lsd is known to be at a 4-byte boundary (so no borrow */
+ /* possible) */
+ carry=0; /* assume no carry */
+ if (diffsign) {
+ carry=CARRYPAT; /* for +1 during add */
+ UINTAT(acc+ 4)=0x09090909-UINTAT(acc+ 4);
+ UINTAT(acc+ 8)=0x09090909-UINTAT(acc+ 8);
+ UINTAT(acc+12)=0x09090909-UINTAT(acc+12);
+ UINTAT(acc+16)=0x09090909-UINTAT(acc+16);
+ #if QUAD
+ UINTAT(acc+20)=0x09090909-UINTAT(acc+20);
+ UINTAT(acc+24)=0x09090909-UINTAT(acc+24);
+ UINTAT(acc+28)=0x09090909-UINTAT(acc+28);
+ UINTAT(acc+32)=0x09090909-UINTAT(acc+32);
+ UINTAT(acc+36)=0x09090909-UINTAT(acc+36);
+ #endif
+ } /* diffsign */
+
+ /* now process the rhs coefficient; if it cannot overlap lhs then */
+ /* it can be put straight into acc (with an appropriate gap, if */
+ /* needed) because no actual addition will be needed (except */
+ /* possibly to complete ten's complement) */
+ overlap=DECPMAX-(expl-expr);
+ #if DECTRACE
+ printf("exps: %ld %ld\n", (LI)expl, (LI)expr);
+ printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
+ #endif
+
+ if (overlap<=0) { /* no overlap possible */
+ uInt gap; /* local work */
+ /* since a full addition is not needed, a ten's complement */
+ /* calculation started above may need to be completed */
+ if (carry) {
+ for (ub=ulsd; *ub==9; ub--) *ub=0;
+ *ub+=1;
+ carry=0; /* taken care of */
+ }
+ /* up to DECPMAX-1 digits of the final result can extend down */
+ /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
+ /* rhs will be simply sticky bits. In this case the gap is */
+ /* clamped to DECPMAX and the exponent adjusted to suit [this is */
+ /* safe because the lhs is non-zero]. */
+ gap=-overlap;
+ if (gap>DECPMAX) {
+ expr+=gap-1;
+ gap=DECPMAX;
+ }
+ ub=ulsd+gap+1; /* where MSD will go */
+ /* Fill the gap with 0s; note that there is no addition to do */
+ ui=&UINTAT(acc+COFF+DECPMAX); /* start of gap */
+ for (; ui<&UINTAT(ub); ui++) *ui=0; /* mind the gap */
+ if (overlap<-DECPMAX) { /* gap was > DECPMAX */
+ *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
+ }
+ else { /* need full coefficient */
+ GETCOEFF(dfr, ub); /* decode from decFloat */
+ ub+=DECPMAX-1; /* new LSD... */
+ }
+ ulsd=ub; /* save new LSD */
+ } /* no overlap possible */
+
+ else { /* overlap>0 */
+ /* coefficients overlap (perhaps completely, although also */
+ /* perhaps only where zeros) */
+ ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */
+ /* Fill the prefix gap with 0s; 8 will cover most common */
+ /* unalignments, so start with direct assignments (a loop is */
+ /* then used for any remaining -- the loop (and the one in a */
+ /* moment) is not then on the critical path because the number */
+ /* of additions is reduced by (at least) two in this case) */
+ UINTAT(buf+4)=0; /* [clears decQuad 00 too] */
+ UINTAT(buf+8)=0;
+ if (ub>buf+12) {
+ ui=&UINTAT(buf+12); /* start of any remaining */
+ for (; ui<&UINTAT(ub); ui++) *ui=0; /* fill them */
+ }
+ GETCOEFF(dfr, ub); /* decode from decFloat */
+
+ /* now move tail of rhs across to main acc; again use direct */
+ /* assignment for 8 digits-worth */
+ UINTAT(acc+COFF+DECPMAX)=UINTAT(buf+COFF+DECPMAX);
+ UINTAT(acc+COFF+DECPMAX+4)=UINTAT(buf+COFF+DECPMAX+4);
+ if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
+ uj=&UINTAT(buf+COFF+DECPMAX+8); /* source */
+ ui=&UINTAT(acc+COFF+DECPMAX+8); /* target */
+ for (; uj<&UINTAT(ub+DECPMAX); ui++, uj++) *ui=*uj;
+ }
+
+ ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
+
+ /* now do the add of the non-tail; this is all nicely aligned, */
+ /* and is over a multiple of four digits (because for Quad two */
+ /* two 0 digits were added on the left); words in both acc and */
+ /* buf (buf especially) will often be zero */
+ /* [byte-by-byte add, here, is about 15% slower than the by-fours] */
+
+ /* Now effect the add; this is harder on a little-endian */
+ /* machine as the inter-digit carry cannot use the usual BCD */
+ /* addition trick because the bytes are loaded in the wrong order */
+ /* [this loop could be unrolled, but probably scarcely worth it] */
+
+ ui=&UINTAT(acc+COFF+DECPMAX-4); /* target LSW (acc) */
+ uj=&UINTAT(buf+COFF+DECPMAX-4); /* source LSW (buf, to add to acc) */
+
+ #if !DECLITEND
+ for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+ /* bcd8 add */
+ carry+=*uj; /* rhs + carry */
+ if (carry==0) continue; /* no-op */
+ carry+=*ui; /* lhs */
+ /* Big-endian BCD adjust (uses internal carry) */
+ carry+=0x76f6f6f6; /* note top nibble not all bits */
+ *ui=(carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4); /* BCD adjust */
+ carry>>=31; /* true carry was at far left */
+ } /* add loop */
+ #else
+ for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+ /* bcd8 add */
+ carry+=*uj; /* rhs + carry */
+ if (carry==0) continue; /* no-op [common if unaligned] */
+ carry+=*ui; /* lhs */
+ /* Little-endian BCD adjust; inter-digit carry must be manual */
+ /* because the lsb from the array will be in the most-significant */
+ /* byte of carry */
+ carry+=0x76767676; /* note no inter-byte carries */
+ carry+=(carry & 0x80000000)>>15;
+ carry+=(carry & 0x00800000)>>15;
+ carry+=(carry & 0x00008000)>>15;
+ carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
+ *ui=carry & 0x0f0f0f0f; /* clear debris and save */
+ /* here, final carry-out bit is at 0x00000080; move it ready */
+ /* for next word-add (i.e., to 0x01000000) */
+ carry=(carry & 0x00000080)<<17;
+ } /* add loop */
+ #endif
+ #if DECTRACE
+ {bcdnum tum;
+ printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
+ tum.msd=umsd; /* acc+4; */
+ tum.lsd=ulsd;
+ tum.exponent=0;
+ tum.sign=0;
+ decShowNum(&tum, "dfadd");}
+ #endif
+ } /* overlap possible */
+
+ /* ordering here is a little strange in order to have slowest path */
+ /* first in GCC asm listing */
+ if (diffsign) { /* subtraction */
+ if (!carry) { /* no carry out means RHS<LHS */
+ /* borrowed -- take ten's complement */
+ /* sign is lhs sign */
+ num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
+
+ /* invert the coefficient first by fours, then add one; space */
+ /* at the end of the buffer ensures the by-fours is always */
+ /* safe, but lsd+1 must be cleared to prevent a borrow */
+ /* if big-endian */
+ #if !DECLITEND
+ *(ulsd+1)=0;
+ #endif
+ /* there are always at least four coefficient words */
+ UINTAT(umsd) =0x09090909-UINTAT(umsd);
+ UINTAT(umsd+4) =0x09090909-UINTAT(umsd+4);
+ UINTAT(umsd+8) =0x09090909-UINTAT(umsd+8);
+ UINTAT(umsd+12)=0x09090909-UINTAT(umsd+12);
+ #if DOUBLE
+ #define BNEXT 16
+ #elif QUAD
+ UINTAT(umsd+16)=0x09090909-UINTAT(umsd+16);
+ UINTAT(umsd+20)=0x09090909-UINTAT(umsd+20);
+ UINTAT(umsd+24)=0x09090909-UINTAT(umsd+24);
+ UINTAT(umsd+28)=0x09090909-UINTAT(umsd+28);
+ UINTAT(umsd+32)=0x09090909-UINTAT(umsd+32);
+ #define BNEXT 36
+ #endif
+ if (ulsd>=umsd+BNEXT) { /* unaligned */
+ /* eight will handle most unaligments for Double; 16 for Quad */
+ UINTAT(umsd+BNEXT)=0x09090909-UINTAT(umsd+BNEXT);
+ UINTAT(umsd+BNEXT+4)=0x09090909-UINTAT(umsd+BNEXT+4);
+ #if DOUBLE
+ #define BNEXTY (BNEXT+8)
+ #elif QUAD
+ UINTAT(umsd+BNEXT+8)=0x09090909-UINTAT(umsd+BNEXT+8);
+ UINTAT(umsd+BNEXT+12)=0x09090909-UINTAT(umsd+BNEXT+12);
+ #define BNEXTY (BNEXT+16)
+ #endif
+ if (ulsd>=umsd+BNEXTY) { /* very unaligned */
+ ui=&UINTAT(umsd+BNEXTY); /* -> continue */
+ for (;;ui++) {
+ *ui=0x09090909-*ui; /* invert four digits */
+ if (ui>=&UINTAT(ulsd-3)) break; /* all done */
+ }
+ }
+ }
+ /* complete the ten's complement by adding 1 */
+ for (ub=ulsd; *ub==9; ub--) *ub=0;
+ *ub+=1;
+ } /* borrowed */
+
+ else { /* carry out means RHS>=LHS */
+ num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
+ /* all done except for the special IEEE 754 exact-zero-result */
+ /* rule (see above); while testing for zero, strip leading */
+ /* zeros (which will save decFinalize doing it) (this is in */
+ /* diffsign path, so carry impossible and true umsd is */
+ /* acc+COFF) */
+
+ /* Check the initial coefficient area using the fast macro; */
+ /* this will often be all that needs to be done (as on the */
+ /* worst-case path when the subtraction was aligned and */
+ /* full-length) */
+ if (ISCOEFFZERO(acc+COFF)) {
+ umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
+ if (ulsd>umsd) { /* more to check */
+ umsd++; /* to align after checked area */
+ for (; UINTAT(umsd)==0 && umsd+3<ulsd;) umsd+=4;
+ for (; *umsd==0 && umsd<ulsd;) umsd++;
+ }
+ if (*umsd==0) { /* must be true zero (and diffsign) */
+ num.sign=0; /* assume + */
+ if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
+ }
+ }
+ /* [else was not zero, might still have leading zeros] */
+ } /* subtraction gave positive result */
+ } /* diffsign */
+
+ else { /* same-sign addition */
+ num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
+ #if DOUBLE
+ if (carry) { /* only possible with decDouble */
+ *(acc+3)=1; /* [Quad has leading 00] */
+ umsd=acc+3;
+ }
+ #endif
+ } /* same sign */
+
+ num.msd=umsd; /* set MSD .. */
+ num.lsd=ulsd; /* .. and LSD */
+ num.exponent=expr; /* set exponent to smaller */
+
+ #if DECTRACE
+ decFloatShow(dfl, "dfl");
+ decFloatShow(dfr, "dfr");
+ decShowNum(&num, "postadd");
+ #endif
+ return decFinalize(result, &num, set); /* round, check, and lay out */
+ } /* decFloatAdd */
+
+/* ------------------------------------------------------------------ */
+/* decFloatAnd -- logical digitwise AND of two decFloats */
+/* */
+/* result gets the result of ANDing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result, which will be canonical with sign=0 */
+/* */
+/* The operands must be positive, finite with exponent q=0, and */
+/* comprise just zeros and ones; if not, Invalid operation results. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatAnd(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
+ || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
+ /* the operands are positive finite integers (q=0) with just 0s and 1s */
+ #if DOUBLE
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
+ #elif QUAD
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
+ DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
+ DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
+ #endif
+ return result;
+ } /* decFloatAnd */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCanonical -- copy a decFloat, making canonical */
+/* */
+/* result gets the canonicalized df */
+/* df is the decFloat to copy and make canonical */
+/* returns result */
+/* */
+/* This works on specials, too; no error or exception is possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
+ return decCanonical(result, df);
+ } /* decFloatCanonical */
+
+/* ------------------------------------------------------------------ */
+/* decFloatClass -- return the class of a decFloat */
+/* */
+/* df is the decFloat to test */
+/* returns the decClass that df falls into */
+/* ------------------------------------------------------------------ */
+enum decClass decFloatClass(const decFloat *df) {
+ Int exp; /* exponent */
+ if (DFISSPECIAL(df)) {
+ if (DFISQNAN(df)) return DEC_CLASS_QNAN;
+ if (DFISSNAN(df)) return DEC_CLASS_SNAN;
+ /* must be an infinity */
+ if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
+ return DEC_CLASS_POS_INF;
+ }
+ if (DFISZERO(df)) { /* quite common */
+ if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
+ return DEC_CLASS_POS_ZERO;
+ }
+ /* is finite and non-zero; similar code to decFloatIsNormal, here */
+ /* [this could be speeded up slightly by in-lining decFloatDigits] */
+ exp=GETEXPUN(df) /* get unbiased exponent .. */
+ +decFloatDigits(df)-1; /* .. and make adjusted exponent */
+ if (exp>=DECEMIN) { /* is normal */
+ if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
+ return DEC_CLASS_POS_NORMAL;
+ }
+ /* is subnormal */
+ if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
+ return DEC_CLASS_POS_SUBNORMAL;
+ } /* decFloatClass */
+
+/* ------------------------------------------------------------------ */
+/* decFloatClassString -- return the class of a decFloat as a string */
+/* */
+/* df is the decFloat to test */
+/* returns a constant string describing the class df falls into */
+/* ------------------------------------------------------------------ */
+const char *decFloatClassString(const decFloat *df) {
+ enum decClass eclass=decFloatClass(df);
+ if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
+ if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
+ if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
+ if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
+ if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
+ if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
+ if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
+ if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
+ if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
+ if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
+ return DEC_ClassString_UN; /* Unknown */
+ } /* decFloatClassString */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
+/* */
+/* result gets the result of comparing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result, which may be -1, 0, 1, or NaN (Unordered) */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCompare(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp; /* work */
+ /* NaNs are handled as usual */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* numeric comparison needed */
+ comp=decNumCompare(dfl, dfr, 0);
+ decFloatZero(result);
+ if (comp==0) return result;
+ DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
+ if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
+ return result;
+ } /* decFloatCompare */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
+/* */
+/* result gets the result of comparing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result, which may be -1, 0, 1, or NaN (Unordered) */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCompareSignal(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp; /* work */
+ /* NaNs are handled as usual, except that all NaNs signal */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) {
+ set->status|=DEC_Invalid_operation;
+ return decNaNs(result, dfl, dfr, set);
+ }
+ /* numeric comparison needed */
+ comp=decNumCompare(dfl, dfr, 0);
+ decFloatZero(result);
+ if (comp==0) return result;
+ DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
+ if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
+ return result;
+ } /* decFloatCompareSignal */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCompareTotal -- compare two decFloats with total ordering */
+/* */
+/* result gets the result of comparing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* returns result, which may be -1, 0, or 1 */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCompareTotal(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr) {
+ Int comp; /* work */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) {
+ Int nanl, nanr; /* work */
+ /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
+ nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
+ if (DFISSIGNED(dfl)) nanl=-nanl;
+ nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
+ if (DFISSIGNED(dfr)) nanr=-nanr;
+ if (nanl>nanr) comp=+1;
+ else if (nanl<nanr) comp=-1;
+ else { /* NaNs are the same type and sign .. must compare payload */
+ /* buffers need +2 for QUAD */
+ uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */
+ uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */
+ uByte *ub, *uc; /* work */
+ Int sigl; /* signum of LHS */
+ sigl=(DFISSIGNED(dfl) ? -1 : +1);
+
+ /* decode the coefficients */
+ /* (shift both right two if Quad to make a multiple of four) */
+ #if QUAD
+ USHORTAT(bufl)=0;
+ USHORTAT(bufr)=0;
+ #endif
+ GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
+ GETCOEFF(dfr, bufr+QUAD*2); /* .. */
+ /* all multiples of four, here */
+ comp=0; /* assume equal */
+ for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
+ if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
+ /* about to find a winner; go by bytes in case little-endian */
+ for (;; ub++, uc++) {
+ if (*ub==*uc) continue;
+ if (*ub>*uc) comp=sigl; /* difference found */
+ else comp=-sigl; /* .. */
+ break;
+ }
+ }
+ } /* same NaN type and sign */
+ }
+ else {
+ /* numeric comparison needed */
+ comp=decNumCompare(dfl, dfr, 1); /* total ordering */
+ }
+ decFloatZero(result);
+ if (comp==0) return result;
+ DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
+ if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
+ return result;
+ } /* decFloatCompareTotal */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCompareTotalMag -- compare magnitudes with total ordering */
+/* */
+/* result gets the result of comparing abs(dfl) and abs(dfr) */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* returns result, which may be -1, 0, or 1 */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCompareTotalMag(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr) {
+ decFloat a, b; /* for copy if needed */
+ /* copy and redirect signed operand(s) */
+ if (DFISSIGNED(dfl)) {
+ decFloatCopyAbs(&a, dfl);
+ dfl=&a;
+ }
+ if (DFISSIGNED(dfr)) {
+ decFloatCopyAbs(&b, dfr);
+ dfr=&b;
+ }
+ return decFloatCompareTotal(result, dfl, dfr);
+ } /* decFloatCompareTotalMag */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCopy -- copy a decFloat as-is */
+/* */
+/* result gets the copy of dfl */
+/* dfl is the decFloat to copy */
+/* returns result */
+/* */
+/* This is a bitwise operation; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
+ if (dfl!=result) *result=*dfl; /* copy needed */
+ return result;
+ } /* decFloatCopy */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
+/* */
+/* result gets the copy of dfl with sign bit 0 */
+/* dfl is the decFloat to copy */
+/* returns result */
+/* */
+/* This is a bitwise operation; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
+ if (dfl!=result) *result=*dfl; /* copy needed */
+ DFBYTE(result, 0)&=~0x80; /* zero sign bit */
+ return result;
+ } /* decFloatCopyAbs */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
+/* */
+/* result gets the copy of dfl with sign bit inverted */
+/* dfl is the decFloat to copy */
+/* returns result */
+/* */
+/* This is a bitwise operation; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
+ if (dfl!=result) *result=*dfl; /* copy needed */
+ DFBYTE(result, 0)^=0x80; /* invert sign bit */
+ return result;
+ } /* decFloatCopyNegate */
+
+/* ------------------------------------------------------------------ */
+/* decFloatCopySign -- copy a decFloat with the sign of another */
+/* */
+/* result gets the result of copying dfl with the sign of dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* returns result */
+/* */
+/* This is a bitwise operation; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatCopySign(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr) {
+ uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
+ if (dfl!=result) *result=*dfl; /* copy needed */
+ DFBYTE(result, 0)&=~0x80; /* clear sign .. */
+ DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
+ return result;
+ } /* decFloatCopySign */
+
+/* ------------------------------------------------------------------ */
+/* decFloatDigits -- return the number of digits in a decFloat */
+/* */
+/* df is the decFloat to investigate */
+/* returns the number of significant digits in the decFloat; a */
+/* zero coefficient returns 1 as does an infinity (a NaN returns */
+/* the number of digits in the payload) */
+/* ------------------------------------------------------------------ */
+/* private macro to extract a declet according to provided formula */
+/* (form), and if it is non-zero then return the calculated digits */
+/* depending on the declet number (n), where n=0 for the most */
+/* significant declet; uses uInt dpd for work */
+#define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
+ if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
+/* next one is used when it is known that the declet must be */
+/* non-zero, or is the final zero declet */
+#define dpdlendun(n, form) {dpd=(form)&0x3ff; \
+ if (dpd==0) return 1; \
+ return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
+
+uInt decFloatDigits(const decFloat *df) {
+ uInt dpd; /* work */
+ uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
+ #if QUAD
+ uInt sourmh, sourml;
+ #endif
+ uInt sourlo;
+
+ if (DFISINF(df)) return 1;
+ /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
+ /* then the coefficient is full-length */
+ if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
+
+ #if DOUBLE
+ if (sourhi&0x0003ffff) { /* ends in first */
+ dpdlenchk(0, sourhi>>8);
+ sourlo=DFWORD(df, 1);
+ dpdlendun(1, (sourhi<<2) | (sourlo>>30));
+ } /* [cannot drop through] */
+ sourlo=DFWORD(df, 1); /* sourhi not involved now */
+ if (sourlo&0xfff00000) { /* in one of first two */
+ dpdlenchk(1, sourlo>>30); /* very rare */
+ dpdlendun(2, sourlo>>20);
+ } /* [cannot drop through] */
+ dpdlenchk(3, sourlo>>10);
+ dpdlendun(4, sourlo);
+ /* [cannot drop through] */
+
+ #elif QUAD
+ if (sourhi&0x00003fff) { /* ends in first */
+ dpdlenchk(0, sourhi>>4);
+ sourmh=DFWORD(df, 1);
+ dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
+ } /* [cannot drop through] */
+ sourmh=DFWORD(df, 1);
+ if (sourmh) {
+ dpdlenchk(1, sourmh>>26);
+ dpdlenchk(2, sourmh>>16);
+ dpdlenchk(3, sourmh>>6);
+ sourml=DFWORD(df, 2);
+ dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
+ } /* [cannot drop through] */
+ sourml=DFWORD(df, 2);
+ if (sourml) {
+ dpdlenchk(4, sourml>>28);
+ dpdlenchk(5, sourml>>18);
+ dpdlenchk(6, sourml>>8);
+ sourlo=DFWORD(df, 3);
+ dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
+ } /* [cannot drop through] */
+ sourlo=DFWORD(df, 3);
+ if (sourlo&0xfff00000) { /* in one of first two */
+ dpdlenchk(7, sourlo>>30); /* very rare */
+ dpdlendun(8, sourlo>>20);
+ } /* [cannot drop through] */
+ dpdlenchk(9, sourlo>>10);
+ dpdlendun(10, sourlo);
+ /* [cannot drop through] */
+ #endif
+ } /* decFloatDigits */
+
+/* ------------------------------------------------------------------ */
+/* decFloatDivide -- divide a decFloat by another */
+/* */
+/* result gets the result of dividing dfl by dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+/* This is just a wrapper. */
+decFloat * decFloatDivide(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ return decDivide(result, dfl, dfr, set, DIVIDE);
+ } /* decFloatDivide */
+
+/* ------------------------------------------------------------------ */
+/* decFloatDivideInteger -- integer divide a decFloat by another */
+/* */
+/* result gets the result of dividing dfl by dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatDivideInteger(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ return decDivide(result, dfl, dfr, set, DIVIDEINT);
+ } /* decFloatDivideInteger */
+
+/* ------------------------------------------------------------------ */
+/* decFloatFMA -- multiply and add three decFloats, fused */
+/* */
+/* result gets the result of (dfl*dfr)+dff with a single rounding */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* dff is the final decFloat (fhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
+ const decFloat *dfr, const decFloat *dff,
+ decContext *set) {
+ /* The accumulator has the bytes needed for FiniteMultiply, plus */
+ /* one byte to the left in case of carry, plus DECPMAX+2 to the */
+ /* right for the final addition (up to full fhs + round & sticky) */
+ #define FMALEN (1+ (DECPMAX9*18) +DECPMAX+2)
+ uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
+ /* .. and for final result */
+ bcdnum mul; /* for multiplication result */
+ bcdnum fin; /* for final operand, expanded */
+ uByte coe[DECPMAX]; /* dff coefficient in BCD */
+ bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
+ uInt diffsign; /* non-zero if signs differ */
+ uInt hipad; /* pad digit for hi if needed */
+ Int padding; /* excess exponent */
+ uInt carry; /* +1 for ten's complement and during add */
+ uByte *ub, *uh, *ul; /* work */
+
+ /* handle all the special values [any special operand leads to a */
+ /* special result] */
+ if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
+ decFloat proxy; /* multiplication result proxy */
+ /* NaNs are handled as usual, giving priority to sNaNs */
+ if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
+ /* One or more of the three is infinite */
+ /* infinity times zero is bad */
+ decFloatZero(&proxy);
+ if (DFISINF(dfl)) {
+ if (DFISZERO(dfr)) return decInvalid(result, set);
+ decInfinity(&proxy, &proxy);
+ }
+ else if (DFISINF(dfr)) {
+ if (DFISZERO(dfl)) return decInvalid(result, set);
+ decInfinity(&proxy, &proxy);
+ }
+ /* compute sign of multiplication and place in proxy */
+ DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
+ if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
+ /* dff is Infinite */
+ if (!DFISINF(&proxy)) return decInfinity(result, dff);
+ /* both sides of addition are infinite; different sign is bad */
+ if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
+ return decInvalid(result, set);
+ return decFloatCopy(result, &proxy);
+ }
+
+ /* Here when all operands are finite */
+
+ /* First multiply dfl*dfr */
+ decFiniteMultiply(&mul, acc+1, dfl, dfr);
+ /* The multiply is complete, exact and unbounded, and described in */
+ /* mul with the coefficient held in acc[1...] */
+
+ /* now add in dff; the algorithm is essentially the same as */
+ /* decFloatAdd, but the code is different because the code there */
+ /* is highly optimized for adding two numbers of the same size */
+ fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
+ fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
+ diffsign=mul.sign^fin.sign; /* note if signs differ */
+ fin.msd=coe;
+ fin.lsd=coe+DECPMAX-1;
+ GETCOEFF(dff, coe); /* extract the coefficient */
+
+ /* now set hi and lo so that hi points to whichever of mul and fin */
+ /* has the higher exponent and lo point to the other [don't care if */
+ /* the same] */
+ if (mul.exponent>=fin.exponent) {
+ hi=&mul;
+ lo=&fin;
+ }
+ else {
+ hi=&fin;
+ lo=&mul;
+ }
+
+ /* remove leading zeros on both operands; this will save time later */
+ /* and make testing for zero trivial */
+ for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
+ for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
+ for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
+ for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
+
+ /* if hi is zero then result will be lo (which has the smaller */
+ /* exponent), which also may need to be tested for zero for the */
+ /* weird IEEE 754 sign rules */
+ if (*hi->msd==0 && hi->msd==hi->lsd) { /* hi is zero */
+ /* "When the sum of two operands with opposite signs is */
+ /* exactly zero, the sign of that sum shall be '+' in all */
+ /* rounding modes except round toward -Infinity, in which */
+ /* mode that sign shall be '-'." */
+ if (diffsign) {
+ if (*lo->msd==0 && lo->msd==lo->lsd) { /* lo is zero */
+ lo->sign=0;
+ if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
+ } /* diffsign && lo=0 */
+ } /* diffsign */
+ return decFinalize(result, lo, set); /* may need clamping */
+ } /* numfl is zero */
+ /* [here, both are minimal length and hi is non-zero] */
+
+ /* if signs differ, take the ten's complement of hi (zeros to the */
+ /* right do not matter because the complement of zero is zero); */
+ /* the +1 is done later, as part of the addition, inserted at the */
+ /* correct digit */
+ hipad=0;
+ carry=0;
+ if (diffsign) {
+ hipad=9;
+ carry=1;
+ /* exactly the correct number of digits must be inverted */
+ for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UINTAT(uh)=0x09090909-UINTAT(uh);
+ for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
+ }
+
+ /* ready to add; note that hi has no leading zeros so gap */
+ /* calculation does not have to be as pessimistic as in decFloatAdd */
+ /* (this is much more like the arbitrary-precision algorithm in */
+ /* Rexx and decNumber) */
+
+ /* padding is the number of zeros that would need to be added to hi */
+ /* for its lsd to be aligned with the lsd of lo */
+ padding=hi->exponent-lo->exponent;
+ /* printf("FMA pad %ld\n", (LI)padding); */
+
+ /* the result of the addition will be built into the accumulator, */
+ /* starting from the far right; this could be either hi or lo */
+ ub=acc+FMALEN-1; /* where lsd of result will go */
+ ul=lo->lsd; /* lsd of rhs */
+
+ if (padding!=0) { /* unaligned */
+ /* if the msd of lo is more than DECPMAX+2 digits to the right of */
+ /* the original msd of hi then it can be reduced to a single */
+ /* digit at the right place, as it stays clear of hi digits */
+ /* [it must be DECPMAX+2 because during a subtraction the msd */
+ /* could become 0 after a borrow from 1.000 to 0.9999...] */
+ Int hilen=(Int)(hi->lsd-hi->msd+1); /* lengths */
+ Int lolen=(Int)(lo->lsd-lo->msd+1); /* .. */
+ Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
+ Int reduce=newexp-lo->exponent;
+ if (reduce>0) { /* [= case gives reduce=0 nop] */
+ /* printf("FMA reduce: %ld\n", (LI)reduce); */
+ if (reduce>=lolen) { /* eating all */
+ lo->lsd=lo->msd; /* reduce to single digit */
+ lo->exponent=newexp; /* [known to be non-zero] */
+ }
+ else { /* < */
+ uByte *up=lo->lsd;
+ lo->lsd=lo->lsd-reduce;
+ if (*lo->lsd==0) /* could need sticky bit */
+ for (; up>lo->lsd; up--) { /* search discarded digits */
+ if (*up!=0) { /* found one... */
+ *lo->lsd=1; /* set sticky bit */
+ break;
+ }
+ }
+ lo->exponent+=reduce;
+ }
+ padding=hi->exponent-lo->exponent; /* recalculate */
+ ul=lo->lsd; /* .. */
+ } /* maybe reduce */
+ /* padding is now <= DECPMAX+2 but still > 0; tricky DOUBLE case */
+ /* is when hi is a 1 that will become a 0.9999... by subtraction: */
+ /* hi: 1 E+16 */
+ /* lo: .................1000000000000000 E-16 */
+ /* which for the addition pads and reduces to: */
+ /* hi: 1000000000000000000 E-2 */
+ /* lo: .................1 E-2 */
+ #if DECCHECK
+ if (padding>DECPMAX+2) printf("FMA excess padding: %ld\n", (LI)padding);
+ if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
+ /* printf("FMA padding: %ld\n", (LI)padding); */
+ #endif
+ /* padding digits can now be set in the result; one or more of */
+ /* these will come from lo; others will be zeros in the gap */
+ for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
+ for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
+ }
+
+ /* addition now complete to the right of the rightmost digit of hi */
+ uh=hi->lsd;
+
+ /* carry was set up depending on ten's complement above; do the add... */
+ for (;; ub--) {
+ uInt hid, lod;
+ if (uh<hi->msd) {
+ if (ul<lo->msd) break;
+ hid=hipad;
+ }
+ else hid=*uh--;
+ if (ul<lo->msd) lod=0;
+ else lod=*ul--;
+ *ub=(uByte)(carry+hid+lod);
+ if (*ub<10) carry=0;
+ else {
+ *ub-=10;
+ carry=1;
+ }
+ } /* addition loop */
+
+ /* addition complete -- now handle carry, borrow, etc. */
+ /* use lo to set up the num (its exponent is already correct, and */
+ /* sign usually is) */
+ lo->msd=ub+1;
+ lo->lsd=acc+FMALEN-1;
+ /* decShowNum(lo, "lo"); */
+ if (!diffsign) { /* same-sign addition */
+ if (carry) { /* carry out */
+ *ub=1; /* place the 1 .. */
+ lo->msd--; /* .. and update */
+ }
+ } /* same sign */
+ else { /* signs differed (subtraction) */
+ if (!carry) { /* no carry out means hi<lo */
+ /* borrowed -- take ten's complement of the right digits */
+ lo->sign=hi->sign; /* sign is lhs sign */
+ for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UINTAT(ul)=0x09090909-UINTAT(ul);
+ for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
+ /* complete the ten's complement by adding 1 [cannot overrun] */
+ for (ul--; *ul==9; ul--) *ul=0;
+ *ul+=1;
+ } /* borrowed */
+ else { /* carry out means hi>=lo */
+ /* sign to use is lo->sign */
+ /* all done except for the special IEEE 754 exact-zero-result */
+ /* rule (see above); while testing for zero, strip leading */
+ /* zeros (which will save decFinalize doing it) */
+ for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
+ for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
+ if (*lo->msd==0) { /* must be true zero (and diffsign) */
+ lo->sign=0; /* assume + */
+ if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
+ }
+ /* [else was not zero, might still have leading zeros] */
+ } /* subtraction gave positive result */
+ } /* diffsign */
+
+ return decFinalize(result, lo, set); /* round, check, and lay out */
+ } /* decFloatFMA */
+
+/* ------------------------------------------------------------------ */
+/* decFloatFromInt -- initialise a decFloat from an Int */
+/* */
+/* result gets the converted Int */
+/* n is the Int to convert */
+/* returns result */
+/* */
+/* The result is Exact; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatFromInt32(decFloat *result, Int n) {
+ uInt u=(uInt)n; /* copy as bits */
+ uInt encode; /* work */
+ DFWORD(result, 0)=ZEROWORD; /* always */
+ #if QUAD
+ DFWORD(result, 1)=0;
+ DFWORD(result, 2)=0;
+ #endif
+ if (n<0) { /* handle -n with care */
+ /* [This can be done without the test, but is then slightly slower] */
+ u=(~u)+1;
+ DFWORD(result, 0)|=DECFLOAT_Sign;
+ }
+ /* Since the maximum value of u now is 2**31, only the low word of */
+ /* result is affected */
+ encode=BIN2DPD[u%1000];
+ u/=1000;
+ encode|=BIN2DPD[u%1000]<<10;
+ u/=1000;
+ encode|=BIN2DPD[u%1000]<<20;
+ u/=1000; /* now 0, 1, or 2 */
+ encode|=u<<30;
+ DFWORD(result, DECWORDS-1)=encode;
+ return result;
+ } /* decFloatFromInt32 */
+
+/* ------------------------------------------------------------------ */
+/* decFloatFromUInt -- initialise a decFloat from a uInt */
+/* */
+/* result gets the converted uInt */
+/* n is the uInt to convert */
+/* returns result */
+/* */
+/* The result is Exact; no errors or exceptions are possible. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
+ uInt encode; /* work */
+ DFWORD(result, 0)=ZEROWORD; /* always */
+ #if QUAD
+ DFWORD(result, 1)=0;
+ DFWORD(result, 2)=0;
+ #endif
+ encode=BIN2DPD[u%1000];
+ u/=1000;
+ encode|=BIN2DPD[u%1000]<<10;
+ u/=1000;
+ encode|=BIN2DPD[u%1000]<<20;
+ u/=1000; /* now 0 -> 4 */
+ encode|=u<<30;
+ DFWORD(result, DECWORDS-1)=encode;
+ DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
+ return result;
+ } /* decFloatFromUInt32 */
+
+/* ------------------------------------------------------------------ */
+/* decFloatInvert -- logical digitwise INVERT of a decFloat */
+/* */
+/* result gets the result of INVERTing df */
+/* df is the decFloat to invert */
+/* set is the context */
+/* returns result, which will be canonical with sign=0 */
+/* */
+/* The operand must be positive, finite with exponent q=0, and */
+/* comprise just zeros and ones; if not, Invalid operation results. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatInvert(decFloat *result, const decFloat *df,
+ decContext *set) {
+ uInt sourhi=DFWORD(df, 0); /* top word of dfs */
+
+ if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
+ /* the operand is a finite integer (q=0) */
+ #if DOUBLE
+ DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
+ DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
+ #elif QUAD
+ DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
+ DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
+ DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
+ DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
+ #endif
+ return result;
+ } /* decFloatInvert */
+
+/* ------------------------------------------------------------------ */
+/* decFloatIs -- decFloat tests (IsSigned, etc.) */
+/* */
+/* df is the decFloat to test */
+/* returns 0 or 1 in an int32_t */
+/* */
+/* Many of these could be macros, but having them as real functions */
+/* is a bit cleaner (and they can be referred to here by the generic */
+/* names) */
+/* ------------------------------------------------------------------ */
+uInt decFloatIsCanonical(const decFloat *df) {
+ if (DFISSPECIAL(df)) {
+ if (DFISINF(df)) {
+ if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
+ if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
+ return 1;
+ }
+ /* is a NaN */
+ if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
+ if (DFISCCZERO(df)) return 1; /* coefficient continuation */
+ /* drop through to check payload */
+ }
+ { /* declare block */
+ #if DOUBLE
+ uInt sourhi=DFWORD(df, 0);
+ uInt sourlo=DFWORD(df, 1);
+ if (CANONDPDOFF(sourhi, 8)
+ && CANONDPDTWO(sourhi, sourlo, 30)
+ && CANONDPDOFF(sourlo, 20)
+ && CANONDPDOFF(sourlo, 10)
+ && CANONDPDOFF(sourlo, 0)) return 1;
+ #elif QUAD
+ uInt sourhi=DFWORD(df, 0);
+ uInt sourmh=DFWORD(df, 1);
+ uInt sourml=DFWORD(df, 2);
+ uInt sourlo=DFWORD(df, 3);
+ if (CANONDPDOFF(sourhi, 4)
+ && CANONDPDTWO(sourhi, sourmh, 26)
+ && CANONDPDOFF(sourmh, 16)
+ && CANONDPDOFF(sourmh, 6)
+ && CANONDPDTWO(sourmh, sourml, 28)
+ && CANONDPDOFF(sourml, 18)
+ && CANONDPDOFF(sourml, 8)
+ && CANONDPDTWO(sourml, sourlo, 30)
+ && CANONDPDOFF(sourlo, 20)
+ && CANONDPDOFF(sourlo, 10)
+ && CANONDPDOFF(sourlo, 0)) return 1;
+ #endif
+ } /* block */
+ return 0; /* a declet is non-canonical */
+ }
+
+uInt decFloatIsFinite(const decFloat *df) {
+ return !DFISSPECIAL(df);
+ }
+uInt decFloatIsInfinite(const decFloat *df) {
+ return DFISINF(df);
+ }
+uInt decFloatIsInteger(const decFloat *df) {
+ return DFISINT(df);
+ }
+uInt decFloatIsNaN(const decFloat *df) {
+ return DFISNAN(df);
+ }
+uInt decFloatIsNormal(const decFloat *df) {
+ Int exp; /* exponent */
+ if (DFISSPECIAL(df)) return 0;
+ if (DFISZERO(df)) return 0;
+ /* is finite and non-zero */
+ exp=GETEXPUN(df) /* get unbiased exponent .. */
+ +decFloatDigits(df)-1; /* .. and make adjusted exponent */
+ return (exp>=DECEMIN); /* < DECEMIN is subnormal */
+ }
+uInt decFloatIsSignaling(const decFloat *df) {
+ return DFISSNAN(df);
+ }
+uInt decFloatIsSignalling(const decFloat *df) {
+ return DFISSNAN(df);
+ }
+uInt decFloatIsSigned(const decFloat *df) {
+ return DFISSIGNED(df);
+ }
+uInt decFloatIsSubnormal(const decFloat *df) {
+ if (DFISSPECIAL(df)) return 0;
+ /* is finite */
+ if (decFloatIsNormal(df)) return 0;
+ /* it is <Nmin, but could be zero */
+ if (DFISZERO(df)) return 0;
+ return 1; /* is subnormal */
+ }
+uInt decFloatIsZero(const decFloat *df) {
+ return DFISZERO(df);
+ } /* decFloatIs... */
+
+/* ------------------------------------------------------------------ */
+/* decFloatLogB -- return adjusted exponent, by 754r rules */
+/* */
+/* result gets the adjusted exponent as an integer, or a NaN etc. */
+/* df is the decFloat to be examined */
+/* set is the context */
+/* returns result */
+/* */
+/* 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 */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatLogB(decFloat *result, const decFloat *df,
+ decContext *set) {
+ Int ae; /* adjusted exponent */
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
+ if (DFISINF(df)) {
+ DFWORD(result, 0)=0; /* need +ve */
+ return decInfinity(result, result); /* canonical +Infinity */
+ }
+ if (DFISZERO(df)) {
+ set->status|=DEC_Division_by_zero; /* as per 754r */
+ DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
+ return decInfinity(result, result); /* canonical -Infinity */
+ }
+ ae=GETEXPUN(df) /* get unbiased exponent .. */
+ +decFloatDigits(df)-1; /* .. and make adjusted exponent */
+ /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
+ /* it is worth using a special case of decFloatFromInt32 */
+ DFWORD(result, 0)=ZEROWORD; /* always */
+ if (ae<0) {
+ DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
+ ae=-ae;
+ }
+ #if DOUBLE
+ DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
+ #elif QUAD
+ DFWORD(result, 1)=0;
+ DFWORD(result, 2)=0;
+ DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
+ DFWORD(result, 3)|=BIN2DPD[ae%1000];
+ #endif
+ return result;
+ } /* decFloatLogB */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMax -- return maxnum of two operands */
+/* */
+/* result gets the chosen decFloat */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* If just one operand is a quiet NaN it is ignored. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMax(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp;
+ if (DFISNAN(dfl)) {
+ /* sNaN or both NaNs leads to normal NaN processing */
+ if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
+ return decCanonical(result, dfr); /* RHS is numeric */
+ }
+ if (DFISNAN(dfr)) {
+ /* sNaN leads to normal NaN processing (both NaN handled above) */
+ if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ return decCanonical(result, dfl); /* LHS is numeric */
+ }
+ /* Both operands are numeric; numeric comparison needed -- use */
+ /* total order for a well-defined choice (and +0 > -0) */
+ comp=decNumCompare(dfl, dfr, 1);
+ if (comp>=0) return decCanonical(result, dfl);
+ return decCanonical(result, dfr);
+ } /* decFloatMax */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMaxMag -- return maxnummag of two operands */
+/* */
+/* result gets the chosen decFloat */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* Returns according to the magnitude comparisons if both numeric and */
+/* unequal, otherwise returns maxnum */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMaxMag(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp;
+ decFloat absl, absr;
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
+
+ decFloatCopyAbs(&absl, dfl);
+ decFloatCopyAbs(&absr, dfr);
+ comp=decNumCompare(&absl, &absr, 0);
+ if (comp>0) return decCanonical(result, dfl);
+ if (comp<0) return decCanonical(result, dfr);
+ return decFloatMax(result, dfl, dfr, set);
+ } /* decFloatMaxMag */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMin -- return minnum of two operands */
+/* */
+/* result gets the chosen decFloat */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* If just one operand is a quiet NaN it is ignored. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMin(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp;
+ if (DFISNAN(dfl)) {
+ /* sNaN or both NaNs leads to normal NaN processing */
+ if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
+ return decCanonical(result, dfr); /* RHS is numeric */
+ }
+ if (DFISNAN(dfr)) {
+ /* sNaN leads to normal NaN processing (both NaN handled above) */
+ if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ return decCanonical(result, dfl); /* LHS is numeric */
+ }
+ /* Both operands are numeric; numeric comparison needed -- use */
+ /* total order for a well-defined choice (and +0 > -0) */
+ comp=decNumCompare(dfl, dfr, 1);
+ if (comp<=0) return decCanonical(result, dfl);
+ return decCanonical(result, dfr);
+ } /* decFloatMin */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMinMag -- return minnummag of two operands */
+/* */
+/* result gets the chosen decFloat */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* Returns according to the magnitude comparisons if both numeric and */
+/* unequal, otherwise returns minnum */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMinMag(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int comp;
+ decFloat absl, absr;
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
+
+ decFloatCopyAbs(&absl, dfl);
+ decFloatCopyAbs(&absr, dfr);
+ comp=decNumCompare(&absl, &absr, 0);
+ if (comp<0) return decCanonical(result, dfl);
+ if (comp>0) return decCanonical(result, dfr);
+ return decFloatMin(result, dfl, dfr, set);
+ } /* decFloatMinMag */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMinus -- negate value, heeding NaNs, etc. */
+/* */
+/* result gets the canonicalized 0-df */
+/* df is the decFloat to minus */
+/* set is the context */
+/* returns result */
+/* */
+/* This has the same effect as 0-df where the exponent of the zero is */
+/* the same as that of df (if df is finite). */
+/* The effect is also the same as decFloatCopyNegate except that NaNs */
+/* are handled normally (the sign of a NaN is not affected, and an */
+/* sNaN will signal), the result is canonical, and zero gets sign 0. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMinus(decFloat *result, const decFloat *df,
+ decContext *set) {
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
+ decCanonical(result, df); /* copy and check */
+ if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
+ else DFBYTE(result, 0)^=0x80; /* flip sign bit */
+ return result;
+ } /* decFloatMinus */
+
+/* ------------------------------------------------------------------ */
+/* decFloatMultiply -- multiply two decFloats */
+/* */
+/* result gets the result of multiplying dfl and dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatMultiply(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ bcdnum num; /* for final conversion */
+ uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
+
+ if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
+ /* NaNs are handled as usual */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* infinity times zero is bad */
+ if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
+ if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
+ /* both infinite; return canonical infinity with computed sign */
+ DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
+ return decInfinity(result, result);
+ }
+
+ /* Here when both operands are finite */
+ decFiniteMultiply(&num, bcdacc, dfl, dfr);
+ return decFinalize(result, &num, set); /* round, check, and lay out */
+ } /* decFloatMultiply */
+
+/* ------------------------------------------------------------------ */
+/* decFloatNextMinus -- next towards -Infinity */
+/* */
+/* result gets the next lesser decFloat */
+/* dfl is the decFloat to start with */
+/* set is the context */
+/* returns result */
+/* */
+/* This is 754r nextdown; Invalid is the only status possible (from */
+/* an sNaN). */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
+ decContext *set) {
+ decFloat delta; /* tiny increment */
+ uInt savestat; /* saves status */
+ enum rounding saveround; /* .. and mode */
+
+ /* +Infinity is the special case */
+ if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
+ DFSETNMAX(result);
+ return result; /* [no status to set] */
+ }
+ /* other cases are effected by sutracting a tiny delta -- this */
+ /* should be done in a wider format as the delta is unrepresentable */
+ /* here (but can be done with normal add if the sign of zero is */
+ /* treated carefully, because no Inexactitude is interesting); */
+ /* rounding to -Infinity then pushes the result to next below */
+ decFloatZero(&delta); /* set up tiny delta */
+ DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
+ DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
+ /* set up for the directional round */
+ saveround=set->round; /* save mode */
+ set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
+ savestat=set->status; /* save status */
+ decFloatAdd(result, dfl, &delta, set);
+ /* Add rules mess up the sign when going from +Ntiny to 0 */
+ if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
+ set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
+ set->status|=savestat; /* restore pending flags */
+ set->round=saveround; /* .. and mode */
+ return result;
+ } /* decFloatNextMinus */
+
+/* ------------------------------------------------------------------ */
+/* decFloatNextPlus -- next towards +Infinity */
+/* */
+/* result gets the next larger decFloat */
+/* dfl is the decFloat to start with */
+/* set is the context */
+/* returns result */
+/* */
+/* This is 754r nextup; Invalid is the only status possible (from */
+/* an sNaN). */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
+ decContext *set) {
+ uInt savestat; /* saves status */
+ enum rounding saveround; /* .. and mode */
+ decFloat delta; /* tiny increment */
+
+ /* -Infinity is the special case */
+ if (DFISINF(dfl) && DFISSIGNED(dfl)) {
+ DFSETNMAX(result);
+ DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
+ return result; /* [no status to set] */
+ }
+ /* other cases are effected by sutracting a tiny delta -- this */
+ /* should be done in a wider format as the delta is unrepresentable */
+ /* here (but can be done with normal add if the sign of zero is */
+ /* treated carefully, because no Inexactitude is interesting); */
+ /* rounding to +Infinity then pushes the result to next above */
+ decFloatZero(&delta); /* set up tiny delta */
+ DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
+ DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
+ /* set up for the directional round */
+ saveround=set->round; /* save mode */
+ set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
+ savestat=set->status; /* save status */
+ decFloatAdd(result, dfl, &delta, set);
+ /* Add rules mess up the sign when going from -Ntiny to -0 */
+ if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
+ set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
+ set->status|=savestat; /* restore pending flags */
+ set->round=saveround; /* .. and mode */
+ return result;
+ } /* decFloatNextPlus */
+
+/* ------------------------------------------------------------------ */
+/* decFloatNextToward -- next towards a decFloat */
+/* */
+/* result gets the next decFloat */
+/* dfl is the decFloat to start with */
+/* dfr is the decFloat to move toward */
+/* set is the context */
+/* returns result */
+/* */
+/* This is 754r nextafter; status may be set unless the result is a */
+/* normal number. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatNextToward(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ decFloat delta; /* tiny increment or decrement */
+ decFloat pointone; /* 1e-1 */
+ uInt savestat; /* saves status */
+ enum rounding saveround; /* .. and mode */
+ uInt deltatop; /* top word for delta */
+ Int comp; /* work */
+
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* Both are numeric, so Invalid no longer a possibility */
+ comp=decNumCompare(dfl, dfr, 0);
+ if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
+ /* unequal; do NextPlus or NextMinus but with different status rules */
+
+ if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
+ if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
+ DFSETNMAX(result);
+ DFWORD(result, 0)|=DECFLOAT_Sign;
+ return result;
+ }
+ saveround=set->round; /* save mode */
+ set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
+ deltatop=0; /* positive delta */
+ }
+ else { /* lhs>rhs, do NextMinus, see above for commentary */
+ if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
+ DFSETNMAX(result);
+ return result;
+ }
+ saveround=set->round; /* save mode */
+ set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
+ deltatop=DECFLOAT_Sign; /* negative delta */
+ }
+ savestat=set->status; /* save status */
+ /* Here, Inexact is needed where appropriate (and hence Underflow, */
+ /* etc.). Therefore the tiny delta which is otherwise */
+ /* unrepresentable (see NextPlus and NextMinus) is constructed */
+ /* using the multiplication of FMA. */
+ decFloatZero(&delta); /* set up tiny delta */
+ DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
+ DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
+ decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
+ decFloatFMA(result, &delta, &pointone, dfl, set);
+ /* [Delta is truly tiny, so no need to correct sign of zero] */
+ /* use new status unless the result is normal */
+ if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
+ set->round=saveround; /* restore mode */
+ return result;
+ } /* decFloatNextToward */
+
+/* ------------------------------------------------------------------ */
+/* decFloatOr -- logical digitwise OR of two decFloats */
+/* */
+/* result gets the result of ORing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result, which will be canonical with sign=0 */
+/* */
+/* The operands must be positive, finite with exponent q=0, and */
+/* comprise just zeros and ones; if not, Invalid operation results. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatOr(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
+ || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
+ /* the operands are positive finite integers (q=0) with just 0s and 1s */
+ #if DOUBLE
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
+ #elif QUAD
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
+ DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
+ DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
+ #endif
+ return result;
+ } /* decFloatOr */
+
+/* ------------------------------------------------------------------ */
+/* decFloatPlus -- add value to 0, heeding NaNs, etc. */
+/* */
+/* result gets the canonicalized 0+df */
+/* df is the decFloat to plus */
+/* set is the context */
+/* returns result */
+/* */
+/* This has the same effect as 0+df where the exponent of the zero is */
+/* the same as that of df (if df is finite). */
+/* The effect is also the same as decFloatCopy except that NaNs */
+/* are handled normally (the sign of a NaN is not affected, and an */
+/* sNaN will signal), the result is canonical, and zero gets sign 0. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatPlus(decFloat *result, const decFloat *df,
+ decContext *set) {
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
+ decCanonical(result, df); /* copy and check */
+ if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
+ return result;
+ } /* decFloatPlus */
+
+/* ------------------------------------------------------------------ */
+/* decFloatQuantize -- quantize a decFloat */
+/* */
+/* result gets the result of quantizing dfl to match dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs), which sets the exponent */
+/* set is the context */
+/* returns result */
+/* */
+/* Unless there is an error or the result is infinite, the exponent */
+/* of result is guaranteed to be the same as that of dfr. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatQuantize(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int explb, exprb; /* left and right biased exponents */
+ uByte *ulsd; /* local LSD pointer */
+ uInt *ui; /* work */
+ uByte *ub; /* .. */
+ Int drop; /* .. */
+ uInt dpd; /* .. */
+ uInt encode; /* encoding accumulator */
+ uInt sourhil, sourhir; /* top words from source decFloats */
+ /* the following buffer holds the coefficient for manipulation */
+ uByte buf[4+DECPMAX*3]; /* + space for zeros to left or right */
+ #if DECTRACE
+ bcdnum num; /* for trace displays */
+ #endif
+
+ /* Start decoding the arguments */
+ sourhil=DFWORD(dfl, 0); /* LHS top word */
+ explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
+ sourhir=DFWORD(dfr, 0); /* RHS top word */
+ exprb=DECCOMBEXP[sourhir>>26];
+
+ if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
+ /* NaNs are handled as usual */
+ if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ /* one infinity but not both is bad */
+ if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
+ /* both infinite; return canonical infinity with sign of LHS */
+ return decInfinity(result, dfl);
+ }
+
+ /* Here when both arguments are finite */
+ /* complete extraction of the exponents [no need to unbias] */
+ explb+=GETECON(dfl); /* + continuation */
+ exprb+=GETECON(dfr); /* .. */
+
+ /* calculate the number of digits to drop from the coefficient */
+ drop=exprb-explb; /* 0 if nothing to do */
+ if (drop==0) return decCanonical(result, dfl); /* return canonical */
+
+ /* the coefficient is needed; lay it out into buf, offset so zeros */
+ /* can be added before or after as needed -- an extra heading is */
+ /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
+ /* fours */
+ #define BUFOFF (buf+4+DECPMAX)
+ GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
+ /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
+
+ #if DECTRACE
+ num.msd=BUFOFF;
+ num.lsd=BUFOFF+DECPMAX-1;
+ num.exponent=explb-DECBIAS;
+ num.sign=sourhil & DECFLOAT_Sign;
+ decShowNum(&num, "dfl");
+ #endif
+
+ if (drop>0) { /* [most common case] */
+ /* (this code is very similar to that in decFloatFinalize, but */
+ /* has many differences so is duplicated here -- so any changes */
+ /* may need to be made there, too) */
+ uByte *roundat; /* -> re-round digit */
+ uByte reround; /* reround value */
+ /* printf("Rounding; drop=%ld\n", (LI)drop); */
+
+ /* there is at least one zero needed to the left, in all but one */
+ /* exceptional (all-nines) case, so place four zeros now; this is */
+ /* needed almost always and makes rounding all-nines by fours safe */
+ UINTAT(BUFOFF-4)=0;
+
+ /* Three cases here: */
+ /* 1. new LSD is in coefficient (almost always) */
+ /* 2. new LSD is digit to left of coefficient (so MSD is */
+ /* round-for-reround digit) */
+ /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
+ /* Note that leading zeros can safely be treated as useful digits */
+
+ /* [duplicate check-stickies code to save a test] */
+ /* [by-digit check for stickies as runs of zeros are rare] */
+ if (drop<DECPMAX) { /* NB lengths not addresses */
+ roundat=BUFOFF+DECPMAX-drop;
+ reround=*roundat;
+ for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
+ if (*ub!=0) { /* non-zero to be discarded */
+ reround=DECSTICKYTAB[reround]; /* apply sticky bit */
+ break; /* [remainder don't-care] */
+ }
+ } /* check stickies */
+ ulsd=roundat-1; /* set LSD */
+ }
+ else { /* edge case */
+ if (drop==DECPMAX) {
+ roundat=BUFOFF;
+ reround=*roundat;
+ }
+ else {
+ roundat=BUFOFF-1;
+ reround=0;
+ }
+ for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
+ if (*ub!=0) { /* non-zero to be discarded */
+ reround=DECSTICKYTAB[reround]; /* apply sticky bit */
+ break; /* [remainder don't-care] */
+ }
+ } /* check stickies */
+ *BUFOFF=0; /* make a coefficient of 0 */
+ ulsd=BUFOFF; /* .. at the MSD place */
+ }
+
+ if (reround!=0) { /* discarding non-zero */
+ uInt bump=0;
+ set->status|=DEC_Inexact;
+
+ /* next decide whether to increment the coefficient */
+ if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
+ if (reround>5) bump=1; /* >0.5 goes up */
+ else if (reround==5) /* exactly 0.5000 .. */
+ bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */
+ } /* r-h-e */
+ else switch (set->round) {
+ case DEC_ROUND_DOWN: {
+ /* no change */
+ break;} /* r-d */
+ case DEC_ROUND_HALF_DOWN: {
+ if (reround>5) bump=1;
+ break;} /* r-h-d */
+ case DEC_ROUND_HALF_UP: {
+ if (reround>=5) bump=1;
+ break;} /* r-h-u */
+ case DEC_ROUND_UP: {
+ if (reround>0) bump=1;
+ break;} /* r-u */
+ case DEC_ROUND_CEILING: {
+ /* same as _UP for positive numbers, and as _DOWN for negatives */
+ if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
+ break;} /* r-c */
+ case DEC_ROUND_FLOOR: {
+ /* same as _UP for negative numbers, and as _DOWN for positive */
+ /* [negative reround cannot occur on 0] */
+ if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
+ break;} /* r-f */
+ case DEC_ROUND_05UP: {
+ if (reround>0) { /* anything out there is 'sticky' */
+ /* bump iff lsd=0 or 5; this cannot carry so it could be */
+ /* effected immediately with no bump -- but the code */
+ /* is clearer if this is done the same way as the others */
+ if (*ulsd==0 || *ulsd==5) bump=1;
+ }
+ break;} /* r-r */
+ default: { /* e.g., DEC_ROUND_MAX */
+ set->status|=DEC_Invalid_context;
+ #if DECCHECK
+ printf("Unknown rounding mode: %ld\n", (LI)set->round);
+ #endif
+ break;}
+ } /* switch (not r-h-e) */
+ /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
+
+ if (bump!=0) { /* need increment */
+ /* increment the coefficient; this could give 1000... (after */
+ /* the all nines case) */
+ ub=ulsd;
+ for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
+ /* now at most 3 digits left to non-9 (usually just the one) */
+ for (; *ub==9; ub--) *ub=0;
+ *ub+=1;
+ /* [the all-nines case will have carried one digit to the */
+ /* left of the original MSD -- just where it is needed] */
+ } /* bump needed */
+ } /* inexact rounding */
+
+ /* now clear zeros to the left so exactly DECPMAX digits will be */
+ /* available in the coefficent -- the first word to the left was */
+ /* cleared earlier for safe carry; now add any more needed */
+ if (drop>4) {
+ UINTAT(BUFOFF-8)=0; /* must be at least 5 */
+ for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
+ }
+ } /* need round (drop>0) */
+
+ else { /* drop<0; padding with -drop digits is needed */
+ /* This is the case where an error can occur if the padded */
+ /* coefficient will not fit; checking for this can be done in the */
+ /* same loop as padding for zeros if the no-hope and zero cases */
+ /* are checked first */
+ if (-drop>DECPMAX-1) { /* cannot fit unless 0 */
+ if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
+ /* a zero can have any exponent; just drop through and use it */
+ ulsd=BUFOFF+DECPMAX-1;
+ }
+ else { /* padding will fit (but may still be too long) */
+ /* final-word mask depends on endianess */
+ #if DECLITEND
+ static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
+ #else
+ static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
+ #endif
+ for (ui=&UINTAT(BUFOFF+DECPMAX);; ui++) {
+ *ui=0;
+ if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
+ /* if all four digits should be zero, definitely bad */
+ if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
+ return decInvalid(result, set);
+ /* must be a 1- to 3-digit sequence; check more carefully */
+ if ((UINTAT(&UBYTEAT(ui)-DECPMAX)&dmask[(-drop)%4])!=0)
+ return decInvalid(result, set);
+ break; /* no need for loop end test */
+ }
+ if (ui>=&UINTAT(BUFOFF+DECPMAX+(-drop)-4)) break; /* done */
+ }
+ ulsd=BUFOFF+DECPMAX+(-drop)-1;
+ } /* pad and check leading zeros */
+ } /* drop<0 */
+
+ #if DECTRACE
+ num.msd=ulsd-DECPMAX+1;
+ num.lsd=ulsd;
+ num.exponent=explb-DECBIAS;
+ num.sign=sourhil & DECFLOAT_Sign;
+ decShowNum(&num, "res");
+ #endif
+
+ /*------------------------------------------------------------------*/
+ /* At this point the result is DECPMAX digits, ending at ulsd, so */
+ /* fits the encoding exactly; there is no possibility of error */
+ /*------------------------------------------------------------------*/
+ encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
+ encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */
+ /* the exponent continuation can be extracted from the original RHS */
+ encode|=sourhir & ECONMASK;
+ encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */
+
+ /* finally encode the coefficient */
+ /* private macro to encode a declet; this version can be used */
+ /* because all coefficient digits exist */
+ #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
+ dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
+
+ #if DOUBLE
+ getDPD3q(dpd, 4); encode|=dpd<<8;
+ getDPD3q(dpd, 3); encode|=dpd>>2;
+ DFWORD(result, 0)=encode;
+ encode=dpd<<30;
+ getDPD3q(dpd, 2); encode|=dpd<<20;
+ getDPD3q(dpd, 1); encode|=dpd<<10;
+ getDPD3q(dpd, 0); encode|=dpd;
+ DFWORD(result, 1)=encode;
+
+ #elif QUAD
+ getDPD3q(dpd,10); encode|=dpd<<4;
+ getDPD3q(dpd, 9); encode|=dpd>>6;
+ DFWORD(result, 0)=encode;
+ encode=dpd<<26;
+ getDPD3q(dpd, 8); encode|=dpd<<16;
+ getDPD3q(dpd, 7); encode|=dpd<<6;
+ getDPD3q(dpd, 6); encode|=dpd>>4;
+ DFWORD(result, 1)=encode;
+ encode=dpd<<28;
+ getDPD3q(dpd, 5); encode|=dpd<<18;
+ getDPD3q(dpd, 4); encode|=dpd<<8;
+ getDPD3q(dpd, 3); encode|=dpd>>2;
+ DFWORD(result, 2)=encode;
+ encode=dpd<<30;
+ getDPD3q(dpd, 2); encode|=dpd<<20;
+ getDPD3q(dpd, 1); encode|=dpd<<10;
+ getDPD3q(dpd, 0); encode|=dpd;
+ DFWORD(result, 3)=encode;
+ #endif
+ return result;
+ } /* decFloatQuantize */
+
+/* ------------------------------------------------------------------ */
+/* decFloatReduce -- reduce finite coefficient to minimum length */
+/* */
+/* result gets the reduced decFloat */
+/* df is the source decFloat */
+/* set is the context */
+/* returns result, which will be canonical */
+/* */
+/* This removes all possible trailing zeros from the coefficient; */
+/* some may remain when the number is very close to Nmax. */
+/* Special values are unchanged and no status is set unless df=sNaN. */
+/* Reduced zero has an exponent q=0. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatReduce(decFloat *result, const decFloat *df,
+ decContext *set) {
+ bcdnum num; /* work */
+ uByte buf[DECPMAX], *ub; /* coefficient and pointer */
+ if (df!=result) *result=*df; /* copy, if needed */
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */
+ /* zeros and infinites propagate too */
+ if (DFISINF(df)) return decInfinity(result, df); /* canonical */
+ if (DFISZERO(df)) {
+ uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
+ decFloatZero(result);
+ DFWORD(result, 0)|=sign;
+ return result; /* exponent dropped, sign OK */
+ }
+ /* non-zero finite */
+ GETCOEFF(df, buf);
+ ub=buf+DECPMAX-1; /* -> lsd */
+ if (*ub) return result; /* no trailing zeros */
+ for (ub--; *ub==0;) ub--; /* terminates because non-zero */
+ /* *ub is the first non-zero from the right */
+ num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
+ num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
+ num.msd=buf;
+ num.lsd=ub;
+ return decFinalize(result, &num, set);
+ } /* decFloatReduce */
+
+/* ------------------------------------------------------------------ */
+/* decFloatRemainder -- integer divide and return remainder */
+/* */
+/* result gets the remainder of dividing dfl by dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatRemainder(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ return decDivide(result, dfl, dfr, set, REMAINDER);
+ } /* decFloatRemainder */
+
+/* ------------------------------------------------------------------ */
+/* decFloatRemainderNear -- integer divide to nearest and remainder */
+/* */
+/* result gets the remainder of dividing dfl by dfr: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* This is the IEEE remainder, where the nearest integer is used. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatRemainderNear(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ return decDivide(result, dfl, dfr, set, REMNEAR);
+ } /* decFloatRemainderNear */
+
+/* ------------------------------------------------------------------ */
+/* decFloatRotate -- rotate the coefficient of a decFloat left/right */
+/* */
+/* result gets the result of rotating dfl */
+/* dfl is the source decFloat to rotate */
+/* dfr is the count of digits to rotate, an integer (with q=0) */
+/* set is the context */
+/* returns result */
+/* */
+/* The digits of the coefficient of dfl are rotated to the left (if */
+/* dfr is positive) or to the right (if dfr is negative) without */
+/* adjusting the exponent or the sign of dfl. */
+/* */
+/* dfr must be in the range -DECPMAX through +DECPMAX. */
+/* NaNs are propagated as usual. An infinite dfl is unaffected (but */
+/* dfr must be valid). No status is set unless dfr is invalid or an */
+/* operand is an sNaN. The result is canonical. */
+/* ------------------------------------------------------------------ */
+#define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
+decFloat * decFloatRotate(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int rotate; /* dfr as an Int */
+ uByte buf[DECPMAX+PHALF]; /* coefficient + half */
+ uInt digits, savestat; /* work */
+ bcdnum num; /* .. */
+ uByte *ub; /* .. */
+
+ if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ if (!DFISINT(dfr)) return decInvalid(result, set);
+ digits=decFloatDigits(dfr); /* calculate digits */
+ if (digits>2) return decInvalid(result, set); /* definitely out of range */
+ rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
+ if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
+ /* [from here on no error or status change is possible] */
+ if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
+ /* handle no-rotate cases */
+ if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
+ /* a real rotate is needed: 0 < rotate < DECPMAX */
+ /* reduce the rotation to no more than half to reduce copying later */
+ /* (for QUAD in fact half + 2 digits) */
+ if (DFISSIGNED(dfr)) rotate=-rotate;
+ if (abs(rotate)>PHALF) {
+ if (rotate<0) rotate=DECPMAX+rotate;
+ else rotate=rotate-DECPMAX;
+ }
+ /* now lay out the coefficient, leaving room to the right or the */
+ /* left depending on the direction of rotation */
+ ub=buf;
+ if (rotate<0) ub+=PHALF; /* rotate right, so space to left */
+ GETCOEFF(dfl, ub);
+ /* copy half the digits to left or right, and set num.msd */
+ if (rotate<0) {
+ memcpy(buf, buf+DECPMAX, PHALF);
+ num.msd=buf+PHALF+rotate;
+ }
+ else {
+ memcpy(buf+DECPMAX, buf, PHALF);
+ num.msd=buf+rotate;
+ }
+ /* fill in rest of num */
+ num.lsd=num.msd+DECPMAX-1;
+ num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
+ num.exponent=GETEXPUN(dfl);
+ savestat=set->status; /* record */
+ decFinalize(result, &num, set);
+ set->status=savestat; /* restore */
+ return result;
+ } /* decFloatRotate */
+
+/* ------------------------------------------------------------------ */
+/* decFloatSameQuantum -- test decFloats for same quantum */
+/* */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* returns 1 if the operands have the same quantum, 0 otherwise */
+/* */
+/* No error is possible and no status results. */
+/* ------------------------------------------------------------------ */
+uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
+ if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
+ if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
+ if (DFISINF(dfl) && DFISINF(dfr)) return 1;
+ return 0; /* any other special mixture gives false */
+ }
+ if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
+ return 0;
+ } /* decFloatSameQuantum */
+
+/* ------------------------------------------------------------------ */
+/* decFloatScaleB -- multiply by a power of 10, as per 754r */
+/* */
+/* result gets the result of the operation */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs), am integer (with q=0) */
+/* set is the context */
+/* returns result */
+/* */
+/* This computes result=dfl x 10**dfr where dfr is an integer in the */
+/* range +/-2*(emax+pmax), typically resulting from LogB. */
+/* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
+/* as usual. */
+/* ------------------------------------------------------------------ */
+#define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
+decFloat * decFloatScaleB(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ uInt digits; /* work */
+ Int expr; /* dfr as an Int */
+
+ if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ if (!DFISINT(dfr)) return decInvalid(result, set);
+ digits=decFloatDigits(dfr); /* calculate digits */
+
+ #if DOUBLE
+ if (digits>3) return decInvalid(result, set); /* definitely out of range */
+ expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */
+ #elif QUAD
+ if (digits>5) return decInvalid(result, set); /* definitely out of range */
+ expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */
+ +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */
+ #endif
+ if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */
+ /* [from now on no error possible] */
+ if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
+ if (DFISSIGNED(dfr)) expr=-expr;
+ /* dfl is finite and expr is valid */
+ *result=*dfl; /* copy to target */
+ return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
+ } /* decFloatScaleB */
+
+/* ------------------------------------------------------------------ */
+/* decFloatShift -- shift the coefficient of a decFloat left or right */
+/* */
+/* result gets the result of shifting dfl */
+/* dfl is the source decFloat to shift */
+/* dfr is the count of digits to shift, an integer (with q=0) */
+/* set is the context */
+/* returns result */
+/* */
+/* The digits of the coefficient of dfl are shifted to the left (if */
+/* dfr is positive) or to the right (if dfr is negative) without */
+/* adjusting the exponent or the sign of dfl. */
+/* */
+/* dfr must be in the range -DECPMAX through +DECPMAX. */
+/* NaNs are propagated as usual. An infinite dfl is unaffected (but */
+/* dfr must be valid). No status is set unless dfr is invalid or an */
+/* operand is an sNaN. The result is canonical. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatShift(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ Int shift; /* dfr as an Int */
+ uByte buf[DECPMAX*2]; /* coefficient + padding */
+ uInt digits, savestat; /* work */
+ bcdnum num; /* .. */
+
+ if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+ if (!DFISINT(dfr)) return decInvalid(result, set);
+ digits=decFloatDigits(dfr); /* calculate digits */
+ if (digits>2) return decInvalid(result, set); /* definitely out of range */
+ shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
+ if (shift>DECPMAX) return decInvalid(result, set); /* too big */
+ /* [from here on no error or status change is possible] */
+
+ if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
+ /* handle no-shift and all-shift (clear to zero) cases */
+ if (shift==0) return decCanonical(result, dfl);
+ if (shift==DECPMAX) { /* zero with sign */
+ uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
+ decFloatZero(result); /* make +0 */
+ DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
+ /* [cannot safely use CopySign] */
+ return result;
+ }
+ /* a real shift is needed: 0 < shift < DECPMAX */
+ num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
+ num.exponent=GETEXPUN(dfl);
+ num.msd=buf;
+ GETCOEFF(dfl, buf);
+ if (DFISSIGNED(dfr)) { /* shift right */
+ /* edge cases are taken care of, so this is easy */
+ num.lsd=buf+DECPMAX-shift-1;
+ }
+ else { /* shift left -- zero padding needed to right */
+ UINTAT(buf+DECPMAX)=0; /* 8 will handle most cases */
+ UINTAT(buf+DECPMAX+4)=0; /* .. */
+ if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
+ num.msd+=shift;
+ num.lsd=num.msd+DECPMAX-1;
+ }
+ savestat=set->status; /* record */
+ decFinalize(result, &num, set);
+ set->status=savestat; /* restore */
+ return result;
+ } /* decFloatShift */
+
+/* ------------------------------------------------------------------ */
+/* decFloatSubtract -- subtract a decFloat from another */
+/* */
+/* result gets the result of subtracting dfr from dfl: */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result */
+/* */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatSubtract(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ decFloat temp;
+ /* NaNs must propagate without sign change */
+ if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
+ temp=*dfr; /* make a copy */
+ DFBYTE(&temp, 0)^=0x80; /* flip sign */
+ return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
+ } /* decFloatSubtract */
+
+/* ------------------------------------------------------------------ */
+/* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
+/* */
+/* df is the decFloat to round */
+/* set is the context */
+/* round is the rounding mode to use */
+/* returns a uInt or an Int, rounded according to the name */
+/* */
+/* Invalid will always be signaled if df is a NaN, is Infinite, or is */
+/* outside the range of the target; Inexact will not be signaled for */
+/* simple rounding unless 'Exact' appears in the name. */
+/* ------------------------------------------------------------------ */
+uInt decFloatToUInt32(const decFloat *df, decContext *set,
+ enum rounding round) {
+ return decToInt32(df, set, round, 0, 1);}
+
+uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
+ enum rounding round) {
+ return decToInt32(df, set, round, 1, 1);}
+
+Int decFloatToInt32(const decFloat *df, decContext *set,
+ enum rounding round) {
+ return (Int)decToInt32(df, set, round, 0, 0);}
+
+Int decFloatToInt32Exact(const decFloat *df, decContext *set,
+ enum rounding round) {
+ return (Int)decToInt32(df, set, round, 1, 0);}
+
+/* ------------------------------------------------------------------ */
+/* decFloatToIntegral -- round to integral value (two flavours) */
+/* */
+/* result gets the result */
+/* df is the decFloat to round */
+/* set is the context */
+/* round is the rounding mode to use */
+/* returns result */
+/* */
+/* No exceptions, even Inexact, are raised except for sNaN input, or */
+/* if 'Exact' appears in the name. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
+ decContext *set, enum rounding round) {
+ return decToIntegral(result, df, set, round, 0);}
+
+decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
+ decContext *set) {
+ return decToIntegral(result, df, set, set->round, 1);}
+
+/* ------------------------------------------------------------------ */
+/* decFloatXor -- logical digitwise XOR of two decFloats */
+/* */
+/* result gets the result of XORing dfl and dfr */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) */
+/* set is the context */
+/* returns result, which will be canonical with sign=0 */
+/* */
+/* The operands must be positive, finite with exponent q=0, and */
+/* comprise just zeros and ones; if not, Invalid operation results. */
+/* ------------------------------------------------------------------ */
+decFloat * decFloatXor(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
+ || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
+ /* the operands are positive finite integers (q=0) with just 0s and 1s */
+ #if DOUBLE
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
+ #elif QUAD
+ DFWORD(result, 0)=ZEROWORD
+ |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
+ DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
+ DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
+ DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
+ #endif
+ return result;
+ } /* decFloatXor */
+
+/* ------------------------------------------------------------------ */
+/* decInvalid -- set Invalid_operation result */
+/* */
+/* result gets a canonical NaN */
+/* set is the context */
+/* returns result */
+/* */
+/* status has Invalid_operation added */
+/* ------------------------------------------------------------------ */
+static decFloat *decInvalid(decFloat *result, decContext *set) {
+ decFloatZero(result);
+ DFWORD(result, 0)=DECFLOAT_qNaN;
+ set->status|=DEC_Invalid_operation;
+ return result;
+ } /* decInvalid */
+
+/* ------------------------------------------------------------------ */
+/* decInfinity -- set canonical Infinity with sign from a decFloat */
+/* */
+/* result gets a canonical Infinity */
+/* df is source decFloat (only the sign is used) */
+/* returns result */
+/* */
+/* df may be the same as result */
+/* ------------------------------------------------------------------ */
+static decFloat *decInfinity(decFloat *result, const decFloat *df) {
+ uInt sign=DFWORD(df, 0); /* save source signword */
+ decFloatZero(result); /* clear everything */
+ DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
+ return result;
+ } /* decInfinity */
+
+/* ------------------------------------------------------------------ */
+/* decNaNs -- handle NaN argument(s) */
+/* */
+/* result gets the result of handling dfl and dfr, one or both of */
+/* which is a NaN */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs) -- may be NULL for a single- */
+/* operand operation */
+/* set is the context */
+/* returns result */
+/* */
+/* Called when one or both operands is a NaN, and propagates the */
+/* appropriate result to res. When an sNaN is found, it is changed */
+/* to a qNaN and Invalid operation is set. */
+/* ------------------------------------------------------------------ */
+static decFloat *decNaNs(decFloat *result,
+ const decFloat *dfl, const decFloat *dfr,
+ decContext *set) {
+ /* handle sNaNs first */
+ if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
+ if (DFISSNAN(dfl)) {
+ decCanonical(result, dfl); /* propagate canonical sNaN */
+ DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
+ set->status|=DEC_Invalid_operation;
+ return result;
+ }
+ /* one or both is a quiet NaN */
+ if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */
+ return decCanonical(result, dfl); /* propagate canonical qNaN */
+ } /* decNaNs */
+
+/* ------------------------------------------------------------------ */
+/* decNumCompare -- numeric comparison of two decFloats */
+/* */
+/* dfl is the left-hand decFloat, which is not a NaN */
+/* dfr is the right-hand decFloat, which is not a NaN */
+/* tot is 1 for total order compare, 0 for simple numeric */
+/* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
+/* */
+/* No error is possible; status and mode are unchanged. */
+/* ------------------------------------------------------------------ */
+static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
+ Int sigl, sigr; /* LHS and RHS non-0 signums */
+ Int shift; /* shift needed to align operands */
+ uByte *ub, *uc; /* work */
+ /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
+ uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */
+ uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */
+
+ sigl=1;
+ if (DFISSIGNED(dfl)) {
+ if (!DFISSIGNED(dfr)) { /* -LHS +RHS */
+ if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
+ return -1; /* RHS wins */
+ }
+ sigl=-1;
+ }
+ if (DFISSIGNED(dfr)) {
+ if (!DFISSIGNED(dfl)) { /* +LHS -RHS */
+ if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
+ return +1; /* LHS wins */
+ }
+ }
+
+ /* signs are the same; operand(s) could be zero */
+ sigr=-sigl; /* sign to return if abs(RHS) wins */
+
+ if (DFISINF(dfl)) {
+ if (DFISINF(dfr)) return 0; /* both infinite & same sign */
+ return sigl; /* inf > n */
+ }
+ if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
+
+ /* here, both are same sign and finite; calculate their offset */
+ shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */
+ /* [bias can be ignored -- the absolute exponent is not relevant] */
+
+ if (DFISZERO(dfl)) {
+ if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */
+ /* both are zero, return 0 if both same exponent or numeric compare */
+ if (shift==0 || !tot) return 0;
+ if (shift>0) return sigl;
+ return sigr; /* [shift<0] */
+ }
+ else { /* LHS!=0 */
+ if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */
+ }
+ /* both are known to be non-zero at this point */
+
+ /* if the exponents are so different that the coefficients do not */
+ /* overlap (by even one digit) then a full comparison is not needed */
+ if (abs(shift)>=DECPMAX) { /* no overlap */
+ /* coefficients are known to be non-zero */
+ if (shift>0) return sigl;
+ return sigr; /* [shift<0] */
+ }
+
+ /* decode the coefficients */
+ /* (shift both right two if Quad to make a multiple of four) */
+ #if QUAD
+ UINTAT(bufl)=0;
+ UINTAT(bufr)=0;
+ #endif
+ GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
+ GETCOEFF(dfr, bufr+QUAD*2); /* .. */
+ if (shift==0) { /* aligned; common and easy */
+ /* all multiples of four, here */
+ for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
+ if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
+ /* about to find a winner; go by bytes in case little-endian */
+ for (;; ub++, uc++) {
+ if (*ub>*uc) return sigl; /* difference found */
+ if (*ub<*uc) return sigr; /* .. */
+ }
+ }
+ } /* aligned */
+ else if (shift>0) { /* lhs to left */
+ ub=bufl; /* RHS pointer */
+ /* pad bufl so right-aligned; most shifts will fit in 8 */
+ UINTAT(bufl+DECPMAX+QUAD*2)=0; /* add eight zeros */
+ UINTAT(bufl+DECPMAX+QUAD*2+4)=0; /* .. */
+ if (shift>8) {
+ /* more than eight; fill the rest, and also worth doing the */
+ /* lead-in by fours */
+ uByte *up; /* work */
+ uByte *upend=bufl+DECPMAX+QUAD*2+shift;
+ for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+ /* [pads up to 36 in all for Quad] */
+ for (;; ub+=4) {
+ if (UINTAT(ub)!=0) return sigl;
+ if (ub+4>bufl+shift-4) break;
+ }
+ }
+ /* check remaining leading digits */
+ for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
+ /* now start the overlapped part; bufl has been padded, so the */
+ /* comparison can go for the full length of bufr, which is a */
+ /* multiple of 4 bytes */
+ for (uc=bufr; ; uc+=4, ub+=4) {
+ if (UINTAT(uc)!=UINTAT(ub)) { /* mismatch found */
+ for (;; uc++, ub++) { /* check from left [little-endian?] */
+ if (*ub>*uc) return sigl; /* difference found */
+ if (*ub<*uc) return sigr; /* .. */
+ }
+ } /* mismatch */
+ if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
+ }
+ } /* shift>0 */
+
+ else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
+ uc=bufr; /* RHS pointer */
+ /* pad bufr so right-aligned; most shifts will fit in 8 */
+ UINTAT(bufr+DECPMAX+QUAD*2)=0; /* add eight zeros */
+ UINTAT(bufr+DECPMAX+QUAD*2+4)=0; /* .. */
+ if (shift<-8) {
+ /* more than eight; fill the rest, and also worth doing the */
+ /* lead-in by fours */
+ uByte *up; /* work */
+ uByte *upend=bufr+DECPMAX+QUAD*2-shift;
+ for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+ /* [pads up to 36 in all for Quad] */
+ for (;; uc+=4) {
+ if (UINTAT(uc)!=0) return sigr;
+ if (uc+4>bufr-shift-4) break;
+ }
+ }
+ /* check remaining leading digits */
+ for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
+ /* now start the overlapped part; bufr has been padded, so the */
+ /* comparison can go for the full length of bufl, which is a */
+ /* multiple of 4 bytes */
+ for (ub=bufl; ; ub+=4, uc+=4) {
+ if (UINTAT(ub)!=UINTAT(uc)) { /* mismatch found */
+ for (;; ub++, uc++) { /* check from left [little-endian?] */
+ if (*ub>*uc) return sigl; /* difference found */
+ if (*ub<*uc) return sigr; /* .. */
+ }
+ } /* mismatch */
+ if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
+ }
+ } /* shift<0 */
+
+ /* Here when compare equal */
+ if (!tot) return 0; /* numerically equal */
+ /* total ordering .. exponent matters */
+ if (shift>0) return sigl; /* total order by exponent */
+ if (shift<0) return sigr; /* .. */
+ return 0;
+ } /* decNumCompare */
+
+/* ------------------------------------------------------------------ */
+/* decToInt32 -- local routine to effect ToInteger conversions */
+/* */
+/* df is the decFloat to convert */
+/* set is the context */
+/* rmode is the rounding mode to use */
+/* exact is 1 if Inexact should be signalled */
+/* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
+/* returns 32-bit result as a uInt */
+/* */
+/* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
+/* these cases 0 is returned. */
+/* ------------------------------------------------------------------ */
+static uInt decToInt32(const decFloat *df, decContext *set,
+ enum rounding rmode, Flag exact, Flag unsign) {
+ Int exp; /* exponent */
+ uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */
+ uInt hi, lo; /* .. penultimate, least, etc. */
+ decFloat zero, result; /* work */
+ Int i; /* .. */
+
+ /* Start decoding the argument */
+ sourhi=DFWORD(df, 0); /* top word */
+ exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
+ if (EXPISSPECIAL(exp)) { /* is special? */
+ set->status|=DEC_Invalid_operation; /* signal */
+ return 0;
+ }
+
+ /* Here when the argument is finite */
+ if (GETEXPUN(df)==0) result=*df; /* already a true integer */
+ else { /* need to round to integer */
+ enum rounding saveround; /* saver */
+ uInt savestatus; /* .. */
+ saveround=set->round; /* save rounding mode .. */
+ savestatus=set->status; /* .. and status */
+ set->round=rmode; /* set mode */
+ decFloatZero(&zero); /* make 0E+0 */
+ set->status=0; /* clear */
+ decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
+ set->round=saveround; /* restore rounding mode .. */
+ if (exact) set->status|=savestatus; /* include Inexact */
+ else set->status=savestatus; /* .. or just original status */
+ }
+
+ /* only the last four declets of the coefficient can contain */
+ /* non-zero; check for others (and also NaN or Infinity from the */
+ /* Quantize) first (see DFISZERO for explanation): */
+ /* decFloatShow(&result, "sofar"); */
+ #if DOUBLE
+ if ((DFWORD(&result, 0)&0x1c03ff00)!=0
+ || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
+ #elif QUAD
+ if ((DFWORD(&result, 2)&0xffffff00)!=0
+ || DFWORD(&result, 1)!=0
+ || (DFWORD(&result, 0)&0x1c003fff)!=0
+ || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
+ #endif
+ set->status|=DEC_Invalid_operation; /* Invalid or out of range */
+ return 0;
+ }
+ /* get last twelve digits of the coefficent into hi & ho, base */
+ /* 10**9 (see GETCOEFFBILL): */
+ sourlo=DFWORD(&result, DECWORDS-1);
+ lo=DPD2BIN0[sourlo&0x3ff]
+ +DPD2BINK[(sourlo>>10)&0x3ff]
+ +DPD2BINM[(sourlo>>20)&0x3ff];
+ sourpen=DFWORD(&result, DECWORDS-2);
+ hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
+
+ /* according to request, check range carefully */
+ if (unsign) {
+ if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
+ set->status|=DEC_Invalid_operation; /* out of range */
+ return 0;
+ }
+ return hi*BILLION+lo;
+ }
+ /* signed */
+ if (hi>2 || (hi==2 && lo>147483647)) {
+ /* handle the usual edge case */
+ if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
+ set->status|=DEC_Invalid_operation; /* truly out of range */
+ return 0;
+ }
+ i=hi*BILLION+lo;
+ if (DFISSIGNED(&result)) i=-i;
+ return (uInt)i;
+ } /* decToInt32 */
+
+/* ------------------------------------------------------------------ */
+/* decToIntegral -- local routine to effect ToIntegral value */
+/* */
+/* result gets the result */
+/* df is the decFloat to round */
+/* set is the context */
+/* rmode is the rounding mode to use */
+/* exact is 1 if Inexact should be signalled */
+/* returns result */
+/* ------------------------------------------------------------------ */
+static decFloat * decToIntegral(decFloat *result, const decFloat *df,
+ decContext *set, enum rounding rmode,
+ Flag exact) {
+ Int exp; /* exponent */
+ uInt sourhi; /* top word from source decFloat */
+ enum rounding saveround; /* saver */
+ uInt savestatus; /* .. */
+ decFloat zero; /* work */
+
+ /* Start decoding the argument */
+ sourhi=DFWORD(df, 0); /* top word */
+ exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
+
+ if (EXPISSPECIAL(exp)) { /* is special? */
+ /* NaNs are handled as usual */
+ if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
+ /* must be infinite; return canonical infinity with sign of df */
+ return decInfinity(result, df);
+ }
+
+ /* Here when the argument is finite */
+ /* complete extraction of the exponent */
+ exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */
+
+ if (exp>=0) return decCanonical(result, df); /* already integral */
+
+ saveround=set->round; /* save rounding mode .. */
+ savestatus=set->status; /* .. and status */
+ set->round=rmode; /* set mode */
+ decFloatZero(&zero); /* make 0E+0 */
+ decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
+ set->round=saveround; /* restore rounding mode .. */
+ if (!exact) set->status=savestatus; /* .. and status, unless exact */
+ return result;
+ } /* decToIntegral */
generated by cgit v1.1 at 2024-11-23 05:42:06 +0000