diff options
Diffstat (limited to 'gcc-4.4.3/libdecnumber/decBasic.c')
| -rw-r--r-- | gcc-4.4.3/libdecnumber/decBasic.c | 3768 |
1 files changed, 0 insertions, 3768 deletions
diff --git a/gcc-4.4.3/libdecnumber/decBasic.c b/gcc-4.4.3/libdecnumber/decBasic.c deleted file mode 100644 index 66d98c773..000000000 --- a/gcc-4.4.3/libdecnumber/decBasic.c +++ /dev/null @@ -1,3768 +0,0 @@ -/* Common base code for the decNumber C Library. - 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 3, or (at your option) any later - version. - - GCC is distributed in the hope that it will be useful, but WITHOUT ANY - WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - -Under Section 7 of GPL version 3, you are granted additional -permissions described in the GCC Runtime Library Exception, version -3.1, as published by the Free Software Foundation. - -You should have received a copy of the GNU General Public License and -a copy of the GCC Runtime Library Exception along with this program; -see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -<http://www.gnu.org/licenses/>. */ - -/* ------------------------------------------------------------------ */ -/* 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("ient, &num, set); /* lay out the integer so far */ - DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */ - sign=DFWORD(dfl, 0); /* save sign of dfl */ - decFloatFMA(result, "ient, dfr, dfl, set); - if (!DFISZERO(result)) return result; - /* if the result is zero the sign shall be sign of dfl */ - DFWORD("ient, 0)=sign; /* construct decFloat of sign */ - return decFloatCopySign(result, result, "ient); - } /* 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 - ub = bufl; /* avoid type-pun violation */ - USHORTAT(ub)=0; - uc = bufr; /* avoid type-pun violation */ - USHORTAT(uc)=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 - ub=bufl; /* avoid type-pun violation */ - UINTAT(ub)=0; - uc=bufr; /* avoid type-pun violation */ - UINTAT(uc)=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 */ |