/* Common base code for the decNumber C Library.
- Copyright (C) 2007 Free Software Foundation, Inc.
+ Copyright (C) 2007, 2009 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
+ Software Foundation; either version 3, 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. */
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
+<http://www.gnu.org/licenses/>. */
/* ------------------------------------------------------------------ */
/* 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). */
+/* 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 */
#define DIVIDE 0x80000000 /* Divide operations [as flags] */
#define REMAINDER 0x40000000 /* .. */
#define DIVIDEINT 0x20000000 /* .. */
-#define REMNEAR 0x10000000 /* .. */
+#define REMNEAR 0x10000000 /* .. */
/* Private functions (local, used only by routines in this module) */
static decFloat *decDivide(decFloat *, const decFloat *,
/* decCanonical -- copy a decFloat, making canonical */
/* */
/* result gets the canonicalized df */
-/* df is the decFloat to copy and make canonical */
+/* df is the decFloat to copy and make canonical */
/* returns result */
/* */
/* This is exposed via decFloatCanonical for Double and Quad only. */
uoff-=32;
dpd|=encode<<(10-uoff); /* get pending bits */
}
- dpd&=0x3ff; /* clear uninteresting 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 */
+ encode|=canon<<(uoff-10); /* insert the canonical form */
DFWORD(result, inword)=encode; /* .. and save */
continue;
}
/* decDivide -- divide operations */
/* */
/* result gets the result of dividing dfl by dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
-/* op is the operation selector */
+/* 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 DIVBASE ((uInt)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,
decFloat quotient; /* for remainders */
bcdnum num; /* for final conversion */
uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
- uInt div[DIVOPLEN]; /* divisor 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 */
+ 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 *ua, *ud, *uq; /* .. */
+ uByte *ub; /* .. */
+ uInt uiwork; /* for macros */
uInt divtop; /* top unit of div adjusted for estimating */
#if DIVCOUNT
static uInt maxcount=0; /* worst-seen subtractions count */
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 */
+ return decInfinity(result, result); /* x/0 -> signed Infinity */
}
num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
if (DFISZERO(dfl)) { /* 0/x (x!=0) */
DFWORD(result, 0)|=num.sign; /* add sign */
return result;
}
- if (!(op&DIVIDE)) { /* a remainder */
+ 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 */
#endif
/* set msu and lsu pointers */
- msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
+ 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--;
/* 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 */
+ 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, */
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 */
+ if (*ud==*ua) { /* all compared equal */
*lsuq+=1; /* increment result */
msua=lsua; /* collapse acc units */
*msua=0; /* .. to a zero */
}
else if (divunits==1) {
mul=(uLong)*msua * DIVBASE + *(msua-1);
- mul/=*msud; /* no more to the right */
+ mul/=*msud; /* no more to the right */
}
else {
- mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) + (*(msua-1)<<2);
+ mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
+ + (*(msua-1)<<2);
mul/=divtop; /* [divtop already allows for sticky bits] */
}
multiplier=(Int)mul;
/* 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 */
+ uInt top, mid, rem; /* work */
if (*uq==0) { /* no split needed */
- UINTAT(ub)=0; /* clear 9 BCD8s */
- UINTAT(ub+4)=0; /* .. */
+ UBFROMUI(ub, 0); /* clear 9 BCD8s */
+ UBFROMUI(ub+4, 0); /* .. */
*(ub+8)=0; /* .. */
continue;
}
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]);
+ UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
+ UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
+ UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
} /* BCD conversion loop */
- ub--; /* -> lsu */
+ ub--; /* -> lsu */
/* complete the bcdnum; quodigits is correct, so the position of */
/* the first non-zero is known */
num.msd--; /* use the 0 .. */
num.lsd=num.msd; /* .. at the new MSD place */
}
- if (reround!=0) { /* discarding non-zero */
+ 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 */
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 (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
for (; *ub==9; ub--) *ub=0; /* at most 3 more */
*ub+=1;
if (ub<num.msd) num.msd--; /* carried */
/* */
/* num gets the result of multiplying dfl and dfr */
/* bcdacc .. with the coefficient in this array */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* */
/* This effects the multiplication of two decFloats, both known to be */
/* 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 */
+/* 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. */
/* */
/* during lazy carry splitting because the initial quotient estimate */
/* (est) can exceed 32 bits. */
-#define MULTBASE BILLION /* the base used for multiply */
+#define MULTBASE ((uInt)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 */
uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
uInt *ui, *uj; /* work */
- uByte *ub; /* .. */
+ uByte *ub; /* .. */
+ uInt uiwork; /* for macros */
#if DECUSE64
- uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
- uLong *pl; /* work -> lazy accumulator */
+ 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 .. */
/* zero the accumulator */
#if MULACCLEN==4
accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
- #else /* use a loop */
+ #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] */
/* */
/* 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 */
+ /* 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 */
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) {
+ if (*pl>=MULTBASE) {
/* *pl holds a binary number which needs to be split */
hop=(uInt)(*pl>>MULSHIFTA);
est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
/* 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 */
+ /* extremely slow on 32-bit machines. This algorithm splits X */
/* using: */
/* */
/* magic=2**(A+B)/1E+9; // 'magic number' */
/* */
/* 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 */
+ /* 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 */
printf("\n");
#endif
- for (pa=acc;; pa++) { /* each low uInt */
+ for (pa=acc;; pa++) { /* each low uInt */
uInt hi, lo; /* words of exact multiply result */
uInt hop, estlo; /* work */
#if QUAD
- uInt esthi; /* .. */
+ uInt esthi; /* .. */
#endif
lo=*pa;
- hi=*(pa+MULACCLEN); /* top 32 bits */
+ hi=*(pa+MULACCLEN); /* top 32 bits */
/* hi and lo now hold a binary number which needs to be split */
#if DOUBLE
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 mulsplit9 1000000 /* divisor */
#define mulsplit6 1000 /* divisor */
/* The splitting is done by simple divides and remainders, */
/* assuming the compiler will optimize these where useful */
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]);
+ UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
+ UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
+ UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
}
else { /* *pa==0 */
- UINTAT(ub)=0; /* clear 9 BCD8s */
- UINTAT(ub+4)=0; /* .. */
+ UBFROMUI(ub, 0); /* clear 9 BCD8s */
+ UBFROMUI(ub+4, 0); /* .. */
*(ub+8)=0; /* .. */
}
if (pa==acc) break;
/* decFloatAbs -- absolute value, heeding NaNs, etc. */
/* */
/* result gets the canonicalized df with sign 0 */
-/* df is the decFloat to abs */
+/* df is the decFloat to abs */
/* set is the context */
/* returns result */
/* */
/* decFloatAdd -- add two decFloats */
/* */
/* result gets the result of adding dfl and dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* */
/* ------------------------------------------------------------------ */
+#if QUAD
+/* Table for testing MSDs for fastpath elimination; returns the MSD of */
+/* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
+/* Infinities return -32 and NaNs return -128 so that summing the two */
+/* MSDs also allows rapid tests for the Specials (see code below). */
+const Int DECTESTMSD[64]={
+ 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
+ 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
+#else
+/* The table for testing MSDs is shared between the modules */
+extern const Int DECTESTMSD[64];
+#endif
+
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; /* .. */
+ Int bexpl, bexpr; /* left and right biased exponents */
+ uByte *ub, *us, *ut; /* work */
+ uInt uiwork; /* for macros */
+ #if QUAD
+ uShort uswork; /* .. */
+ #endif
uInt sourhil, sourhir; /* top words from source decFloats */
- /* [valid only until specials */
- /* handled or exponents decoded] */
+ /* [valid only through end of */
+ /* fastpath code -- before swap] */
uInt diffsign; /* non-zero if signs differ */
uInt carry; /* carry: 0 or 1 before add loop */
- Int overlap; /* coefficient overlap (if full) */
+ Int overlap; /* coefficient overlap (if full) */
+ Int summ; /* sum of the MSDs */
/* the following buffers hold coefficients with various alignments */
/* (see commentary and diagrams below) */
uByte acc[4+2+DECPMAX*3+8];
uByte *umsd, *ulsd; /* local MSD and LSD pointers */
#if DECLITEND
- #define CARRYPAT 0x01000000 /* carry=1 pattern */
+ #define CARRYPAT 0x01000000 /* carry=1 pattern */
#else
- #define CARRYPAT 0x00000001 /* carry=1 pattern */
+ #define CARRYPAT 0x00000001 /* carry=1 pattern */
#endif
/* Start decoding the arguments */
- /* the initial exponents are placed into the opposite Ints to */
+ /* 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) */
+ /* whichever path is taken (thereby minimising worst-case path). */
+ /* The calculated exponents will be nonsense when the arguments are */
+ /* Special, but are not used in that path */
sourhil=DFWORD(dfl, 0); /* LHS top word */
- expr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
+ summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */
+ bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
+ bexpr+=GETECON(dfl); /* .. + continuation */
+
sourhir=DFWORD(dfr, 0); /* RHS top word */
- expl=DECCOMBEXP[sourhir>>26];
+ summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */
+ bexpl=DECCOMBEXP[sourhir>>26];
+ bexpl+=GETECON(dfr);
+
+ /* here bexpr has biased exponent from lhs, and vice versa */
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 */
- }
+ /* now determine whether to take a fast path or the full-function */
+ /* slow path. The slow path must be taken when: */
+ /* -- both numbers are finite, and: */
+ /* the exponents are different, or */
+ /* the signs are different, or */
+ /* the sum of the MSDs is >8 (hence might overflow) */
+ /* specialness and the sum of the MSDs can be tested at once using */
+ /* the summ value just calculated, so the test for specials is no */
+ /* longer on the worst-case path (as of 3.60) */
+
+ if (summ<=8) { /* MSD+MSD is good, or there is a special */
+ if (summ<0) { /* there is a special */
+ /* Inf+Inf would give -64; Inf+finite is -32 or higher */
+ if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */
+ /* two infinities with different signs is invalid */
+ if (summ==-64 && diffsign) return decInvalid(result, set);
+ if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
+ return decInfinity(result, dfr); /* RHS must be Inf */
+ }
+ /* Here when both arguments are finite; fast path is possible */
+ /* (currently only for aligned and same-sign) */
+ if (bexpr==bexpl && !diffsign) {
+ uInt tac[DECLETS+1]; /* base-1000 coefficient */
+ uInt encode; /* work */
+
+ /* Get one coefficient as base-1000 and add the other */
+ GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */
+ ADDCOEFFTHOU(dfr, tac);
+ /* here the sum of the MSDs (plus any carry) will be <10 due to */
+ /* the fastpath test earlier */
+
+ /* construct the result; low word is the same for both formats */
+ encode =BIN2DPD[tac[0]];
+ encode|=BIN2DPD[tac[1]]<<10;
+ encode|=BIN2DPD[tac[2]]<<20;
+ encode|=BIN2DPD[tac[3]]<<30;
+ DFWORD(result, (DECBYTES/4)-1)=encode;
+
+ /* collect next two declets (all that remains, for Double) */
+ encode =BIN2DPD[tac[3]]>>2;
+ encode|=BIN2DPD[tac[4]]<<8;
- /* Here when both arguments are finite */
+ #if QUAD
+ /* complete and lay out middling words */
+ encode|=BIN2DPD[tac[5]]<<18;
+ encode|=BIN2DPD[tac[6]]<<28;
+ DFWORD(result, 2)=encode;
+
+ encode =BIN2DPD[tac[6]]>>4;
+ encode|=BIN2DPD[tac[7]]<<6;
+ encode|=BIN2DPD[tac[8]]<<16;
+ encode|=BIN2DPD[tac[9]]<<26;
+ DFWORD(result, 1)=encode;
+
+ /* and final two declets */
+ encode =BIN2DPD[tac[9]]>>6;
+ encode|=BIN2DPD[tac[10]]<<4;
+ #endif
+
+ /* add exponent continuation and sign (from either argument) */
+ encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
+
+ /* create lookup index = MSD + top two bits of biased exponent <<4 */
+ tac[DECLETS]|=(bexpl>>DECECONL)<<4;
+ encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
+ DFWORD(result, 0)=encode; /* complete */
+
+ /* decFloatShow(result, ">"); */
+ return result;
+ } /* fast path OK */
+ /* drop through to slow path */
+ } /* low sum or Special(s) */
- /* 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 */
+ /* Slow path required -- arguments are finite and might overflow, */
+ /* or require alignment, or might have different signs */
/* now swap either exponents or argument pointers */
- if (expl<=expr) {
+ if (bexpl<=bexpr) {
/* original left is bigger */
- Int expswap=expl;
- expl=expr;
- expr=expswap;
+ Int bexpswap=bexpl;
+ bexpl=bexpr;
+ bexpr=bexpswap;
/* printf("left bigger\n"); */
}
else {
dfr=dfswap;
/* printf("right bigger\n"); */
}
- /* [here dfl and expl refer to the datum with the larger exponent, */
+ /* [here dfl and bexpl 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 */
#if DOUBLE
#define COFF 4 /* offset into acc */
#elif QUAD
- USHORTAT(acc+4)=0; /* prefix 00 */
+ UBFROMUS(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.exponent=bexpl-DECBIAS;
tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
decShowNum(&tum, "dflx");}
#endif
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);
+ UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
+ UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
+ UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
+ UBFROMUI(acc+16, 0x09090909-UBTOUI(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);
+ UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
+ UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
+ UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
+ UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
+ UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
#endif
} /* diffsign */
/* 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);
+ overlap=DECPMAX-(bexpl-bexpr);
#if DECTRACE
- printf("exps: %ld %ld\n", (LI)expl, (LI)expr);
+ printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
#endif
/* safe because the lhs is non-zero]. */
gap=-overlap;
if (gap>DECPMAX) {
- expr+=gap-1;
+ bexpr+=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 */
+ ut=acc+COFF+DECPMAX; /* start of gap */
+ for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */
if (overlap<-DECPMAX) { /* gap was > DECPMAX */
*ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
}
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;
+ if (overlap==DECPMAX) { /* aligned */
+ ub=buf+COFF; /* where msd will go */
+ #if QUAD
+ UBFROMUS(buf+4, 0); /* clear quad's 00 */
+ #endif
+ GETCOEFF(dfr, ub); /* decode from decFloat */
}
+ else { /* unaligned */
+ 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) */
+ UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */
+ UBFROMUI(buf+8, 0);
+ if (ub>buf+12) {
+ ut=buf+12; /* start any remaining */
+ for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
+ }
+ GETCOEFF(dfr, ub); /* decode from decFloat */
+
+ /* now move tail of rhs across to main acc; again use direct */
+ /* copies for 8 digits-worth */
+ UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX));
+ UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
+ if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
+ us=buf+COFF+DECPMAX+8; /* source */
+ ut=acc+COFF+DECPMAX+8; /* target */
+ for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
+ }
+ } /* unaligned */
ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
- /* now do the add of the non-tail; this is all nicely aligned, */
+ /* 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 */
+ /* zero 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] */
+ /* [byte-by-byte add, here, is about 15% slower total effect 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) */
+ ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */
+ us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */
#if !DECLITEND
- for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+ for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */
/* bcd8 add */
- carry+=*uj; /* rhs + carry */
+ carry+=UBTOUI(us); /* rhs + carry */
if (carry==0) continue; /* no-op */
- carry+=*ui; /* lhs */
+ carry+=UBTOUI(ut); /* lhs */
/* Big-endian BCD adjust (uses internal carry) */
carry+=0x76f6f6f6; /* note top nibble not all bits */
- *ui=(carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4); /* BCD adjust */
+ /* apply BCD adjust and save */
+ UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
carry>>=31; /* true carry was at far left */
} /* add loop */
#else
- for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+ for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */
/* bcd8 add */
- carry+=*uj; /* rhs + carry */
+ carry+=UBTOUI(us); /* rhs + carry */
if (carry==0) continue; /* no-op [common if unaligned] */
- carry+=*ui; /* lhs */
+ carry+=UBTOUI(ut); /* 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+=(carry & 0x00800000)>>15;
carry+=(carry & 0x00008000)>>15;
carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
- *ui=carry & 0x0f0f0f0f; /* clear debris and save */
+ UBFROMUI(ut, 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);
*(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);
+ UBFROMUI(umsd, 0x09090909-UBTOUI(umsd));
+ UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4));
+ UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8));
+ UBFROMUI(umsd+12, 0x09090909-UBTOUI(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);
+ UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
+ UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
+ UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
+ UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
+ UBFROMUI(umsd+32, 0x09090909-UBTOUI(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);
+ UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT));
+ UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(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);
+ UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8));
+ UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(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 */
+ ut=umsd+BNEXTY; /* -> continue */
+ for (;;ut+=4) {
+ UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
+ if (ut>=ulsd-3) break; /* all done */
}
}
}
umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
if (ulsd>umsd) { /* more to check */
umsd++; /* to align after checked area */
- for (; U
INTAT(umsd)==0 && umsd+3<ulsd;) umsd+=4;
+ for (; U
BTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
for (; *umsd==0 && umsd<ulsd;) umsd++;
}
- if (*umsd==0) { /* must be true zero (and diffsign) */
+ if (*umsd==0) { /* must be true zero (and diffsign) */
num.sign=0; /* assume + */
if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
}
#endif
} /* same sign */
- num.msd=umsd; /* set MSD .. */
- num.lsd=ulsd; /* .. and LSD */
- num.exponent=expr; /* set exponent to smaller */
+ num.msd=umsd; /* set MSD .. */
+ num.lsd=ulsd; /* .. and LSD */
+ num.exponent=bexpr-DECBIAS; /* set exponent to smaller, unbiassed */
#if DECTRACE
decFloatShow(dfl, "dfl");
/* decFloatAnd -- logical digitwise AND of two decFloats */
/* */
/* result gets the result of ANDing dfl and dfr */
-/* dfl is the first decFloat (lhs) */
+/* 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 */
+/* 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,
/* decFloatCanonical -- copy a decFloat, making canonical */
/* */
/* result gets the canonicalized df */
-/* df is the decFloat to copy and make canonical */
+/* df is the decFloat to copy and make canonical */
/* returns result */
/* */
/* This works on specials, too; no error or exception is possible. */
/* ------------------------------------------------------------------ */
/* decFloatClass -- return the class of a decFloat */
/* */
-/* df is the decFloat to test */
+/* df is the decFloat to test */
/* returns the decClass that df falls into */
/* ------------------------------------------------------------------ */
enum decClass decFloatClass(const decFloat *df) {
/* ------------------------------------------------------------------ */
/* decFloatClassString -- return the class of a decFloat as a string */
/* */
-/* df is the decFloat to test */
+/* df is the decFloat to test */
/* returns a constant string describing the class df falls into */
/* ------------------------------------------------------------------ */
const char *decFloatClassString(const decFloat *df) {
} /* decFloatClassString */
/* ------------------------------------------------------------------ */
-/* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
+/* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
/* */
/* result gets the result of comparing dfl and dfr */
-/* dfl is the first decFloat (lhs) */
+/* 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) */
/* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
/* */
/* result gets the result of comparing dfl and dfr */
-/* dfl is the first decFloat (lhs) */
+/* 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) */
/* decFloatCompareTotal -- compare two decFloats with total ordering */
/* */
/* result gets the result of comparing dfl and dfr */
-/* dfl is the first decFloat (lhs) */
+/* 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 */
+ Int comp; /* work */
+ uInt uiwork; /* for macros */
+ #if QUAD
+ uShort uswork; /* .. */
+ #endif
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 */
+ nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
if (DFISSIGNED(dfl)) nanl=-nanl;
nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
if (DFISSIGNED(dfr)) nanr=-nanr;
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 */
+ 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
- ub = bufl; /* avoid type-pun violation */
- USHORTAT(ub)=0;
- uc = bufr; /* avoid type-pun violation */
- USHORTAT(uc)=0;
+ UBFROMUS(bufl, 0);
+ UBFROMUS(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 */
+ uInt ui=UBTOUI(ub);
+ if (ui==UBTOUI(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;
/* decFloatCompareTotalMag -- compare magnitudes with total ordering */
/* */
/* result gets the result of comparing abs(dfl) and abs(dfr) */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* returns result, which may be -1, 0, or 1 */
/* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
/* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
/* */
-/* result gets the copy of dfl with sign bit inverted */
+/* result gets the copy of dfl with sign bit inverted */
/* dfl is the decFloat to copy */
/* returns result */
/* */
} /* decFloatCopyNegate */
/* ------------------------------------------------------------------ */
-/* decFloatCopySign -- copy a decFloat with the sign of another */
+/* 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) */
+/* 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 */
/* */
/* 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; \
+ if (dpd==0) return 1; \
return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
uInt decFloatDigits(const decFloat *df) {
} /* [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 */
+ dpdlenchk(1, sourlo>>30); /* very rare */
dpdlendun(2, sourlo>>20);
} /* [cannot drop through] */
dpdlenchk(3, sourlo>>10);
} /* [cannot drop through] */
sourlo=DFWORD(df, 3);
if (sourlo&0xfff00000) { /* in one of first two */
- dpdlenchk(7, sourlo>>30); /* very rare */
+ dpdlenchk(7, sourlo>>30); /* very rare */
dpdlendun(8, sourlo>>20);
} /* [cannot drop through] */
dpdlenchk(9, sourlo>>10);
/* decFloatDivide -- divide a decFloat by another */
/* */
/* result gets the result of dividing dfl by dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* decFloatDivideInteger -- integer divide a decFloat by another */
/* */
/* result gets the result of dividing dfl by dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* 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) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
-/* dff is the final decFloat (fhs) */
+/* 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 */
+ #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +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 */
+ uByte coe[ROUNDUP4(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 */
+ 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 */
+ uInt carry; /* +1 for ten's complement and during add */
+ uByte *ub, *uh, *ul; /* work */
+ uInt uiwork; /* for macros */
/* handle all the special values [any special operand leads to a */
/* special result] */
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] */
+ /* has the higher exponent and lo points to the other [don't care, */
+ /* if the same]. One coefficient will be in acc, the other in coe. */
if (mul.exponent>=fin.exponent) {
hi=&mul;
lo=&fin;
}
/* 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;
+ /* and make testing for zero trivial (tests are safe because acc */
+ /* and coe are rounded up to uInts) */
+ for (; UBTOUI(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 (; UBTOUI(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 */
+ if (*hi->msd==0) { /* 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 */
+ if (*lo->msd==0) { /* lo is zero */
lo->sign=0;
if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
} /* diffsign && lo=0 */
return decFinalize(result, lo, set); /* may need clamping */
} /* numfl is zero */
/* [here, both are minimal length and hi is non-zero] */
+ /* (if lo is zero then padding with zeros may be needed, below) */
/* 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 */
+ /* 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;
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->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
}
/* 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 */
+ /* starting from the far right; this could be either hi or lo, and */
+ /* will be aligned */
ub=acc+FMALEN-1; /* where lsd of result will go */
ul=lo->lsd; /* lsd of rhs */
/* 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] */
+
+ Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
+ Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
+
+ if (hilen+padding-lolen > DECPMAX+2) { /* can reduce lo to single */
+ /* make sure it is virtually at least DECPMAX from hi->msd, at */
+ /* least to right of hi->lsd (in case of destructive subtract), */
+ /* and separated by at least two digits from either of those */
+ /* (the 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 to: */
+ /* hi: 1000000000000000000 E-16 */
+ /* lo: .................1000000000000000 E-16 */
+ Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
+
/* 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 */
+ lo->lsd=lo->msd; /* to single digit [maybe 0] */
+ lo->exponent=newexp; /* new lowest exponent */
+ padding=hi->exponent-lo->exponent; /* recalculate */
+ ul=lo->lsd; /* .. and repoint */
+ }
+
+ /* padding is still > 0, but will fit in acc (less leading carry slot) */
#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);
+ if (hilen+padding+1>FMALEN)
+ printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (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-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
+ UBFROMUI(ub-3, UBTOUI(ul-3)); /* [cannot overlap] */
+ }
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... */
+ /* dow do the add from hi->lsd to the left */
+ /* [bytewise, because either operand can run out at any time] */
+ /* carry was set up depending on ten's complement above */
+ /* first assume both operands have some digits */
for (;; ub--) {
- uInt hid, lod;
- if (uh<hi->msd) {
+ if (uh<hi->msd || ul<lo->msd) break;
+ *ub=(uByte)(carry+(*uh--)+(*ul--));
+ carry=0;
+ if (*ub<10) continue;
+ *ub-=10;
+ carry=1;
+ } /* both loop */
+
+ if (ul<lo->msd) { /* to left of lo */
+ for (;; ub--) {
+ if (uh<hi->msd) break;
+ *ub=(uByte)(carry+(*uh--)); /* [+0] */
+ carry=0;
+ if (*ub<10) continue;
+ *ub-=10;
+ carry=1;
+ } /* hi loop */
+ }
+ else { /* to left of hi */
+ for (;; ub--) {
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=(uByte)(carry+hipad+(*ul--));
+ carry=0;
+ if (*ub<10) continue;
*ub-=10;
carry=1;
- }
- } /* addition loop */
+ } /* lo loop */
+ }
/* addition complete -- now handle carry, borrow, etc. */
/* use lo to set up the num (its exponent is already correct, and */
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->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(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;
/* 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 (; UBTOUI(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 + */
} /* subtraction gave positive result */
} /* diffsign */
+ #if DECCHECK
+ /* assert no left underrun */
+ if (lo->msd<acc) {
+ printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
+ }
+ #endif
+
return decFinalize(result, lo, set); /* round, check, and lay out */
} /* decFloatFMA */
/* ------------------------------------------------------------------ */
-/* decFloatFromInt -- initialise a decFloat from an Int */
+/* decFloatFromInt -- initialise a decFloat from an Int */
/* */
/* result gets the converted Int */
/* n is the Int to convert */
/* decFloatInvert -- logical digitwise INVERT of a decFloat */
/* */
/* result gets the result of INVERTing df */
-/* df is the decFloat to invert */
+/* df is the decFloat to invert */
/* set is the context */
/* returns result, which will be canonical with sign=0 */
/* */
/* ------------------------------------------------------------------ */
/* decFloatIs -- decFloat tests (IsSigned, etc.) */
/* */
-/* df is the decFloat to test */
-/* returns 0 or 1 in an int32_t */
+/* df is the decFloat to test */
+/* returns 0 or 1 in a uInt */
/* */
/* 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) */
+/* is a little cleaner (and they can be referred to here by the */
+/* generic names) */
/* ------------------------------------------------------------------ */
uInt decFloatIsCanonical(const decFloat *df) {
if (DFISSPECIAL(df)) {
} /* decFloatIs... */
/* ------------------------------------------------------------------ */
-/* decFloatLogB -- return adjusted exponent, by 754r rules */
+/* decFloatLogB -- return adjusted exponent, by 754 rules */
/* */
/* result gets the adjusted exponent as an integer, or a NaN etc. */
-/* df is the decFloat to be examined */
+/* df is the decFloat to be examined */
/* set is the context */
/* returns result */
/* */
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
if (DFISINF(df)) {
DFWORD(result, 0)=0; /* need +ve */
- return decInfinity(result, result); /* canonical +Infinity */
+ return decInfinity(result, result); /* canonical +Infinity */
}
if (DFISZERO(df)) {
- set->status|=DEC_Division_by_zero; /* as per 754r */
+ set->status|=DEC_Division_by_zero; /* as per 754 */
DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
- return decInfinity(result, result); /* canonical -Infinity */
+ return decInfinity(result, result); /* canonical -Infinity */
}
ae=GETEXPUN(df) /* get unbiased exponent .. */
+decFloatDigits(df)-1; /* .. and make adjusted exponent */
} /* decFloatLogB */
/* ------------------------------------------------------------------ */
-/* decFloatMax -- return maxnum of two operands */
+/* decFloatMax -- return maxnum of two operands */
/* */
/* result gets the chosen decFloat */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* decFloatMaxMag -- return maxnummag of two operands */
/* */
/* result gets the chosen decFloat */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
} /* decFloatMaxMag */
/* ------------------------------------------------------------------ */
-/* decFloatMin -- return minnum of two operands */
+/* decFloatMin -- return minnum of two operands */
/* */
/* result gets the chosen decFloat */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* decFloatMinMag -- return minnummag of two operands */
/* */
/* result gets the chosen decFloat */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* ------------------------------------------------------------------ */
/* decFloatMinus -- negate value, heeding NaNs, etc. */
/* */
-/* result gets the canonicalized 0-df */
-/* df is the decFloat to minus */
+/* result gets the canonicalized 0-df */
+/* df is the decFloat to minus */
/* set is the context */
/* returns result */
/* */
/* ------------------------------------------------------------------ */
/* decFloatMultiply -- multiply two decFloats */
/* */
-/* result gets the result of multiplying dfl and dfr: */
-/* dfl is the first decFloat (lhs) */
+/* 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 */
const decFloat *dfl, const decFloat *dfr,
decContext *set) {
bcdnum num; /* for final conversion */
- uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
+ uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
/* NaNs are handled as usual */
/* set is the context */
/* returns result */
/* */
-/* This is 754r nextdown; Invalid is the only status possible (from */
+/* This is 754 nextdown; Invalid is the only status possible (from */
/* an sNaN). */
/* ------------------------------------------------------------------ */
decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
/* 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 */
+ 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 */
+ saveround=set->round; /* save mode */
set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
- savestat=set->status; /* save status */
+ 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 */
+ set->round=saveround; /* .. and mode */
return result;
} /* decFloatNextMinus */
/* set is the context */
/* returns result */
/* */
-/* This is 754r nextup; Invalid is the only status possible (from */
+/* This is 754 nextup; Invalid is the only status possible (from */
/* an sNaN). */
/* ------------------------------------------------------------------ */
decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
/* 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 */
+ 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 */
+ 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 */
+ set->round=saveround; /* .. and mode */
return result;
} /* decFloatNextPlus */
/* set is the context */
/* returns result */
/* */
-/* This is 754r nextafter; status may be set unless the result is a */
-/* normal number. */
+/* This is 754-1985 nextafter, as modified during revision (dropped */
+/* from 754-2008); status may be set unless the result is a normal */
+/* number. */
/* ------------------------------------------------------------------ */
decFloat * decFloatNextToward(decFloat *result,
const decFloat *dfl, const decFloat *dfr,
}
saveround=set->round; /* save mode */
set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
- deltatop=0; /* positive delta */
+ deltatop=0; /* positive delta */
}
else { /* lhs>rhs, do NextMinus, see above for commentary */
if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
return result;
}
saveround=set->round; /* save mode */
- set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
+ set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
deltatop=DECFLOAT_Sign; /* negative delta */
}
- savestat=set->status; /* save status */
+ 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 */
+ 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 */
+ 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) */
+/* 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 */
+/* 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,
/* ------------------------------------------------------------------ */
/* decFloatPlus -- add value to 0, heeding NaNs, etc. */
/* */
-/* result gets the canonicalized 0+df */
-/* df is the decFloat to plus */
+/* 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 */
+/* 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. */
/* ------------------------------------------------------------------ */
/* decFloatQuantize -- quantize a decFloat */
/* */
/* result gets the result of quantizing dfl to match dfr */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs), which sets the exponent */
/* set is the context */
/* returns result */
decContext *set) {
Int explb, exprb; /* left and right biased exponents */
uByte *ulsd; /* local LSD pointer */
- uInt *ui; /* work */
- uByte *ub; /* .. */
+ uByte *ub, *uc; /* work */
Int drop; /* .. */
uInt dpd; /* .. */
- uInt encode; /* encoding accumulator */
+ uInt encode; /* encoding accumulator */
uInt sourhil, sourhir; /* top words from source decFloats */
+ uInt uiwork; /* for macros */
+ #if QUAD
+ uShort uswork; /* .. */
+ #endif
/* the following buffer holds the coefficient for manipulation */
- uByte buf[4+DECPMAX*3]; /* + space for zeros to left or right */
+ uByte buf[4+DECPMAX*3+2*QUAD]; /* + space for zeros to left or right */
#if DECTRACE
- bcdnum num; /* for trace displays */
+ bcdnum num; /* for trace displays */
#endif
/* Start decoding the arguments */
decShowNum(&num, "dfl");
#endif
- if (drop>0) { /* [most common case] */
+ 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) */
/* 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;
+ UBFROMUI(BUFOFF-4, 0);
/* Three cases here: */
/* 1. new LSD is in coefficient (almost always) */
/* [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 */
+ if (drop<DECPMAX) { /* NB lengths not addresses */
roundat=BUFOFF+DECPMAX-drop;
reround=*roundat;
for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
/* 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;
+ for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
/* now at most 3 digits left to non-9 (usually just the one) */
for (; *ub==9; ub--) *ub=0;
*ub+=1;
/* 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;
+ UBFROMUI(BUFOFF-8, 0); /* must be at least 5 */
+ for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0);
}
} /* need round (drop>0) */
#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 */
+ /* note that here zeros to the right are added by fours, so in */
+ /* the Quad case this could write 36 zeros if the coefficient has */
+ /* fewer than three significant digits (hence the +2*QUAD for buf) */
+ for (uc=BUFOFF+DECPMAX;; uc+=4) {
+ UBFROMUI(uc, 0);
+ if (UBTOUI(uc-DECPMAX)!=0) { /* could be bad */
/* if all four digits should be zero, definitely bad */
- if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
+ if (uc<=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)
+ if ((UBTOUI(uc-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 */
+ if (uc>=BUFOFF+DECPMAX+(-drop)-4) break; /* done */
}
ulsd=BUFOFF+DECPMAX+(-drop)-1;
} /* pad and check leading zeros */
/* decFloatReduce -- reduce finite coefficient to minimum length */
/* */
/* result gets the reduced decFloat */
-/* df is the source decFloat */
+/* df is the source decFloat */
/* set is the context */
/* returns result, which will be canonical */
/* */
/* decFloatRemainder -- integer divide and return remainder */
/* */
/* result gets the remainder of dividing dfl by dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
/* decFloatRemainderNear -- integer divide to nearest and remainder */
/* */
/* result gets the remainder of dividing dfl by dfr: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
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 */
+ 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] */
num.lsd=num.msd+DECPMAX-1;
num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
num.exponent=GETEXPUN(dfl);
- savestat=set->status; /* record */
+ savestat=set->status; /* record */
decFinalize(result, &num, set);
- set->status=savestat; /* restore */
+ set->status=savestat; /* restore */
return result;
} /* decFloatRotate */
/* ------------------------------------------------------------------ */
/* decFloatSameQuantum -- test decFloats for same quantum */
/* */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* returns 1 if the operands have the same quantum, 0 otherwise */
/* */
} /* decFloatSameQuantum */
/* ------------------------------------------------------------------ */
-/* decFloatScaleB -- multiply by a power of 10, as per 754r */
+/* decFloatScaleB -- multiply by a power of 10, as per 754 */
/* */
/* result gets the result of the operation */
-/* dfl is the first decFloat (lhs) */
-/* dfr is the second decFloat (rhs), am integer (with q=0) */
+/* dfl is the first decFloat (lhs) */
+/* dfr is the second decFloat (rhs), am integer (with q=0) */
/* set is the context */
/* returns result */
/* */
digits=decFloatDigits(dfr); /* calculate digits */
#if DOUBLE
- if (digits>3) return decInvalid(result, set); /* definitely out of range */
+ 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 */
+ 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 (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
if (DFISSIGNED(dfr)) expr=-expr;
/* dfl is finite and expr is valid */
- *result=*dfl; /* copy to target */
+ *result=*dfl; /* copy to target */
return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
} /* decFloatScaleB */
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 */
+ Int shift; /* dfr as an Int */
+ uByte buf[DECPMAX*2]; /* coefficient + padding */
+ uInt digits, savestat; /* work */
bcdnum num; /* .. */
+ uInt uiwork; /* for macros */
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 (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 */
+ 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 */
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; /* .. */
+ UBFROMUI(buf+DECPMAX, 0); /* 8 will handle most cases */
+ UBFROMUI(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 */
+ savestat=set->status; /* record */
decFinalize(result, &num, set);
- set->status=savestat; /* restore */
+ set->status=savestat; /* restore */
return result;
} /* decFloatShift */
/* ------------------------------------------------------------------ */
-/* decFloatSubtract -- subtract a decFloat from another */
+/* decFloatSubtract -- subtract a decFloat from another */
/* */
/* result gets the result of subtracting dfr from dfl: */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* set is the context */
/* returns result */
} /* decFloatSubtract */
/* ------------------------------------------------------------------ */
-/* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
+/* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
/* */
-/* df is the decFloat to round */
+/* 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 */
return (Int)decToInt32(df, set, round, 1, 0);}
/* ------------------------------------------------------------------ */
-/* decFloatToIntegral -- round to integral value (two flavours) */
+/* decFloatToIntegral -- round to integral value (two flavours) */
/* */
/* result gets the result */
-/* df is the decFloat to round */
+/* df is the decFloat to round */
/* set is the context */
-/* round is the rounding mode to use */
+/* round is the rounding mode to use */
/* returns result */
/* */
/* No exceptions, even Inexact, are raised except for sNaN input, or */
/* decFloatXor -- logical digitwise XOR of two decFloats */
/* */
/* result gets the result of XORing dfl and dfr */
-/* dfl is the first decFloat (lhs) */
+/* 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 */
+/* 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,
/* decInfinity -- set canonical Infinity with sign from a decFloat */
/* */
/* result gets a canonical Infinity */
-/* df is source decFloat (only the sign is used) */
+/* df is source decFloat (only the sign is used) */
/* returns result */
/* */
-/* df may be the same as 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 */
+ decFloatZero(result); /* clear everything */
DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
return result;
} /* decInfinity */
/* */
/* result gets the result of handling dfl and dfr, one or both of */
/* which is a NaN */
-/* dfl is the first decFloat (lhs) */
+/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) -- may be NULL for a single- */
/* operand operation */
/* set is the context */
} /* decNaNs */
/* ------------------------------------------------------------------ */
-/* decNumCompare -- numeric comparison of two decFloats */
+/* 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 */
+/* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
/* */
-/* No error is possible; status and mode are unchanged. */
+/* 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 */
+ uInt uiwork; /* for macros */
/* 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 */
sigr=-sigl; /* sign to return if abs(RHS) wins */
if (DFISINF(dfl)) {
- if (DFISINF(dfr)) return 0; /* both infinite & same sign */
+ if (DFISINF(dfr)) return 0; /* both infinite & same sign */
return sigl; /* inf > n */
}
if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
/* decode the coefficients */
/* (shift both right two if Quad to make a multiple of four) */
#if QUAD
- ub=bufl; /* avoid type-pun violation */
- UINTAT(ub)=0;
- uc=bufr; /* avoid type-pun violation */
- UINTAT(uc)=0;
+ UBFROMUI(bufl, 0);
+ UBFROMUI(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 */
+ uInt ui=UBTOUI(ub);
+ if (ui==UBTOUI(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 */
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; /* .. */
+ UBFROMUI(bufl+DECPMAX+QUAD*2, 0); /* add eight zeros */
+ UBFROMUI(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 *up; /* work */
uByte *upend=bufl+DECPMAX+QUAD*2+shift;
- for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+ for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
/* [pads up to 36 in all for Quad] */
for (;; ub+=4) {
- if (UINTAT(ub)!=0) return sigl;
+ if (UBTOUI(ub)!=0) return sigl;
if (ub+4>bufl+shift-4) break;
}
}
/* 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 */
+ uInt ui=UBTOUI(ub);
+ if (ui!=UBTOUI(uc)) { /* mismatch found */
for (;; uc++, ub++) { /* check from left [little-endian?] */
if (*ub>*uc) return sigl; /* difference found */
if (*ub<*uc) return sigr; /* .. */
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; /* .. */
+ UBFROMUI(bufr+DECPMAX+QUAD*2, 0); /* add eight zeros */
+ UBFROMUI(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 *up; /* work */
uByte *upend=bufr+DECPMAX+QUAD*2-shift;
- for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+ for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
/* [pads up to 36 in all for Quad] */
for (;; uc+=4) {
- if (UINTAT(uc)!=0) return sigr;
+ if (UBTOUI(uc)!=0) return sigr;
if (uc+4>bufr-shift-4) break;
}
}
/* 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 */
+ uInt ui=UBTOUI(ub);
+ if (ui!=UBTOUI(uc)) { /* mismatch found */
for (;; ub++, uc++) { /* check from left [little-endian?] */
if (*ub>*uc) return sigl; /* difference found */
if (*ub<*uc) return sigr; /* .. */
/* ------------------------------------------------------------------ */
/* decToInt32 -- local routine to effect ToInteger conversions */
/* */
-/* df is the decFloat to convert */
+/* 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 */
+/* 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 */
/* */
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 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 */
+ 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 */
/* decToIntegral -- local routine to effect ToIntegral value */
/* */
/* result gets the result */
-/* df is the decFloat to round */
+/* 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 */
+/* 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,
decFloat zero; /* work */
/* Start decoding the argument */
- sourhi=DFWORD(df, 0); /* top word */
+ sourhi=DFWORD(df, 0); /* top word */
exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
if (EXPISSPECIAL(exp)) { /* is special? */
if (exp>=0) return decCanonical(result, df); /* already integral */
- saveround=set->round; /* save rounding mode .. */
+ 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 .. */
+ set->round=saveround; /* restore rounding mode .. */
if (!exact) set->status=savestatus; /* .. and status, unless exact */
return result;
} /* decToIntegral */
OSDN Git Service
