/* ------------------------------------------------------------------ */ /* decBasic.c -- common base code for Basic decimal types */ /* ------------------------------------------------------------------ */ /* Copyright (c) IBM Corporation, 2000, 2010. All rights reserved. */ /* */ /* This software is made available under the terms of the */ /* ICU License -- ICU 1.8.1 and later. */ /* */ /* The description and User's Guide ("The decNumber C Library") for */ /* this software is included in the package as decNumber.pdf. This */ /* document is also available in HTML, together with specifications, */ /* testcases, and Web links, on the General Decimal Arithmetic page. */ /* */ /* Please send comments, suggestions, and corrections to the author: */ /* mfc@uk.ibm.com */ /* Mike Cowlishaw, IBM Fellow */ /* IBM UK, PO Box 31, Birmingham Road, Warwick CV34 5JL, UK */ /* ------------------------------------------------------------------ */ /* 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<>(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 ((uInt)BILLION) // the base used for divide #define DIVOPLEN DECPMAX9 // operand length ('digits' base 10**9) #define DIVACCLEN (DIVOPLEN*3) // accumulator length (ditto) static decFloat * decDivide(decFloat *result, const decFloat *dfl, 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 uiwork; // for macros 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=lsua && *msua==0;) 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 (accunitsdiv: 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 && *msuadivunits // 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>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-=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=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==lsua && *msua==0) 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 UBFROMUI(ub, 0); // clear 9 BCD8s UBFROMUI(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) UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); } // BCD conversion loop ub--; // -> lsu // 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 (lsuamaxcount) { // 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 re-round digit uByte reround; // reround value *(num.msd-1)=0; // in case of left carry, or make 0 if (drop 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 (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); for (; *ub==9; ub--) *ub=0; // at most 3 more *ub+=1; if (ub9 #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; // .. uInt uiwork; // for macros #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; pl1 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=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 { // *pl1 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)=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) UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); } else { // *pa==0 UBFROMUI(ub, 0); // clear 9 BCD8s UBFROMUI(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 */ /* */ /* ------------------------------------------------------------------ */ #if QUAD // Table for testing MSDs for fastpath elimination; returns the MSD of // a decDouble or decQuad (top 6 bits tested) ignoring the sign. // Infinities return -32 and NaNs return -128 so that summing the two // MSDs also allows rapid tests for the Specials (see code below). const Int DECTESTMSD[64]={ 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128}; #else // The table for testing MSDs is shared between the modules extern const Int DECTESTMSD[64]; #endif decFloat * decFloatAdd(decFloat *result, const decFloat *dfl, const decFloat *dfr, decContext *set) { bcdnum num; // for final conversion Int bexpl, bexpr; // left and right biased exponents uByte *ub, *us, *ut; // work uInt uiwork; // for macros #if QUAD uShort uswork; // .. #endif uInt sourhil, sourhir; // top words from source decFloats // [valid only through end of // fastpath code -- before swap] uInt diffsign; // non-zero if signs differ uInt carry; // carry: 0 or 1 before add loop Int overlap; // coefficient overlap (if full) Int summ; // sum of the MSDs // the following buffers hold coefficients with various alignments // (see commentary and diagrams below) uByte acc[4+2+DECPMAX*3+8]; uByte 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). // The calculated exponents will be nonsense when the arguments are // Special, but are not used in that path sourhil=DFWORD(dfl, 0); // LHS top word summ=DECTESTMSD[sourhil>>26]; // get first MSD for testing bexpr=DECCOMBEXP[sourhil>>26]; // get exponent high bits (in place) bexpr+=GETECON(dfl); // .. + continuation sourhir=DFWORD(dfr, 0); // RHS top word summ+=DECTESTMSD[sourhir>>26]; // sum MSDs for testing bexpl=DECCOMBEXP[sourhir>>26]; bexpl+=GETECON(dfr); // here bexpr has biased exponent from lhs, and vice versa diffsign=(sourhil^sourhir)&DECFLOAT_Sign; // now determine whether to take a fast path or the full-function // slow path. The slow path must be taken when: // -- both numbers are finite, and: // the exponents are different, or // the signs are different, or // the sum of the MSDs is >8 (hence might overflow) // specialness and the sum of the MSDs can be tested at once using // the summ value just calculated, so the test for specials is no // longer on the worst-case path (as of 3.60) if (summ<=8) { // MSD+MSD is good, or there is a special if (summ<0) { // there is a special // Inf+Inf would give -64; Inf+finite is -32 or higher if (summ<-64) return decNaNs(result, dfl, dfr, set); // one or two NaNs // two infinities with different signs is invalid if (summ==-64 && diffsign) return decInvalid(result, set); if (DFISINF(dfl)) return decInfinity(result, dfl); // LHS is infinite return decInfinity(result, dfr); // RHS must be Inf } // Here when both arguments are finite; fast path is possible // (currently only for aligned and same-sign) if (bexpr==bexpl && !diffsign) { uInt tac[DECLETS+1]; // base-1000 coefficient uInt encode; // work // Get one coefficient as base-1000 and add the other GETCOEFFTHOU(dfl, tac); // least-significant goes to [0] ADDCOEFFTHOU(dfr, tac); // here the sum of the MSDs (plus any carry) will be <10 due to // the fastpath test earlier // construct the result; low word is the same for both formats encode =BIN2DPD[tac[0]]; encode|=BIN2DPD[tac[1]]<<10; encode|=BIN2DPD[tac[2]]<<20; encode|=BIN2DPD[tac[3]]<<30; DFWORD(result, (DECBYTES/4)-1)=encode; // collect next two declets (all that remains, for Double) encode =BIN2DPD[tac[3]]>>2; encode|=BIN2DPD[tac[4]]<<8; #if QUAD // complete and lay out middling words encode|=BIN2DPD[tac[5]]<<18; encode|=BIN2DPD[tac[6]]<<28; DFWORD(result, 2)=encode; encode =BIN2DPD[tac[6]]>>4; encode|=BIN2DPD[tac[7]]<<6; encode|=BIN2DPD[tac[8]]<<16; encode|=BIN2DPD[tac[9]]<<26; DFWORD(result, 1)=encode; // and final two declets encode =BIN2DPD[tac[9]]>>6; encode|=BIN2DPD[tac[10]]<<4; #endif // add exponent continuation and sign (from either argument) encode|=sourhil & (ECONMASK | DECFLOAT_Sign); // create lookup index = MSD + top two bits of biased exponent <<4 tac[DECLETS]|=(bexpl>>DECECONL)<<4; encode|=DECCOMBFROM[tac[DECLETS]]; // add constructed combination field DFWORD(result, 0)=encode; // complete // decFloatShow(result, ">"); return result; } // fast path OK // drop through to slow path } // low sum or Special(s) /* Slow path required -- arguments are finite and might overflow, */ /* or require alignment, or might have different signs */ // now swap either exponents or argument pointers if (bexpl<=bexpr) { // original left is bigger Int bexpswap=bexpl; bexpl=bexpr; bexpr=bexpswap; // printf("left bigger\n"); } else { const decFloat *dfswap=dfl; dfl=dfr; dfr=dfswap; // printf("right bigger\n"); } // [here dfl and bexpl refer to the datum with the larger exponent, // of if the exponents are equal then the original LHS argument] // if lhs is zero then result will be the rhs (now known to have // 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 UBFROMUS(acc+4, 0); // prefix 00 #define COFF 6 // offset into acc #endif GETCOEFF(dfl, acc+COFF); // decode from decFloat ulsd=acc+COFF+DECPMAX-1; umsd=acc+4; // [having this here avoids #if DECTRACE {bcdnum tum; tum.msd=umsd; tum.lsd=ulsd; tum.exponent=bexpl-DECBIAS; 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 UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4)); UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8)); UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12)); UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16)); #if QUAD UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20)); UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24)); UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28)); UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32)); UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36)); #endif } // diffsign // 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-(bexpl-bexpr); #if DECTRACE printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS)); 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) { bexpr+=gap-1; gap=DECPMAX; } ub=ulsd+gap+1; // where MSD will go // Fill the gap with 0s; note that there is no addition to do ut=acc+COFF+DECPMAX; // start of gap for (; ut 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) if (overlap==DECPMAX) { // aligned ub=buf+COFF; // where msd will go #if QUAD UBFROMUS(buf+4, 0); // clear quad's 00 #endif GETCOEFF(dfr, ub); // decode from decFloat } else { // unaligned ub=buf+COFF+DECPMAX-overlap; // where MSD will go // Fill the prefix gap with 0s; 8 will cover most common // unalignments, so start with direct assignments (a loop is // then used for any remaining -- the loop (and the one in a // moment) is not then on the critical path because the number // of additions is reduced by (at least) two in this case) UBFROMUI(buf+4, 0); // [clears decQuad 00 too] UBFROMUI(buf+8, 0); if (ub>buf+12) { ut=buf+12; // start any remaining for (; ut=acc+4; ut-=4, us-=4) { // big-endian add loop // bcd8 add carry+=UBTOUI(us); // rhs + carry if (carry==0) continue; // no-op carry+=UBTOUI(ut); // lhs // Big-endian BCD adjust (uses internal carry) carry+=0x76f6f6f6; // note top nibble not all bits // apply BCD adjust and save UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4)); carry>>=31; // true carry was at far left } // add loop #else for (; ut>=acc+4; ut-=4, us-=4) { // little-endian add loop // bcd8 add carry+=UBTOUI(us); // rhs + carry if (carry==0) continue; // no-op [common if unaligned] carry+=UBTOUI(ut); // lhs // Little-endian BCD adjust; inter-digit carry must be manual // because the lsb from the array will be in the most-significant // byte of carry carry+=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 UBFROMUI(ut, carry & 0x0f0f0f0f); // clear debris and save // here, final carry-out bit is at 0x00000080; move it ready // for next word-add (i.e., to 0x01000000) carry=(carry & 0x00000080)<<17; } // add loop #endif #if DECTRACE {bcdnum tum; printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign); 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=umsd+BNEXT) { // unaligned // eight will handle most unaligments for Double; 16 for Quad UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT)); UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4)); #if DOUBLE #define BNEXTY (BNEXT+8) #elif QUAD UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8)); UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12)); #define BNEXTY (BNEXT+16) #endif if (ulsd>=umsd+BNEXTY) { // very unaligned ut=umsd+BNEXTY; // -> continue for (;;ut+=4) { UBFROMUI(ut, 0x09090909-UBTOUI(ut)); // invert four digits if (ut>=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 (; umsd+3round==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=bexpr-DECBIAS; // set exponent to smaller, unbiassed #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 uInt uiwork; // for macros #if QUAD uShort uswork; // .. #endif if (DFISNAN(dfl) || DFISNAN(dfr)) { Int nanl, nanr; // work // morph NaNs to +/- 1 or 2, leave numbers as 0 nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; // quiet > signalling if (DFISSIGNED(dfl)) nanl=-nanl; nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2; if (DFISSIGNED(dfr)) nanr=-nanr; if (nanl>nanr) comp=+1; else if (nanl*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 (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2)) uByte acc[FMALEN]; // for multiplied coefficient in BCD // .. and for final result bcdnum mul; // for multiplication result bcdnum fin; // for final operand, expanded uByte coe[ROUNDUP4(DECPMAX)]; // dff coefficient in BCD bcdnum *hi, *lo; // bcdnum with higher/lower exponent uInt diffsign; // non-zero if signs differ uInt hipad; // pad digit for hi if needed Int padding; // excess exponent uInt carry; // +1 for ten's complement and during add uByte *ub, *uh, *ul; // work uInt uiwork; // for macros // handle all the special values [any special operand leads to a // special result] 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 points to the other [don't care, // if the same]. One coefficient will be in acc, the other in coe. if (mul.exponent>=fin.exponent) { hi=&mul; lo=&fin; } else { hi=&fin; lo=&mul; } // remove leading zeros on both operands; this will save time later // and make testing for zero trivial (tests are safe because acc // and coe are rounded up to uInts) for (; hi->msd+3lsd && UBTOUI(hi->msd)==0;) hi->msd+=4; for (; hi->msdlsd && *hi->msd==0;) hi->msd++; for (; lo->msd+3lsd && UBTOUI(lo->msd)==0;) lo->msd+=4; for (; lo->msdlsd && *lo->msd==0;) 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 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 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 lo is zero then padding with zeros may be needed, below) // if signs differ, take the ten's complement of hi (zeros to the // right do not matter because the complement of zero is zero); the // +1 is done later, as part of the addition, inserted at the // 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; uhlsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(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, and // will be aligned 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); // length of hi Int lolen=(Int)(lo->lsd-lo->msd+1); // and of lo if (hilen+padding-lolen > DECPMAX+2) { // can reduce lo to single // make sure it is virtually at least DECPMAX from hi->msd, at // least to right of hi->lsd (in case of destructive subtract), // and separated by at least two digits from either of those // (the tricky DOUBLE case is when hi is a 1 that will become a // 0.9999... by subtraction: // hi: 1 E+16 // lo: .................1000000000000000 E-16 // which for the addition pads to: // hi: 1000000000000000000 E-16 // lo: .................1000000000000000 E-16 Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3; // printf("FMA reduce: %ld\n", (LI)reduce); lo->lsd=lo->msd; // to single digit [maybe 0] lo->exponent=newexp; // new lowest exponent padding=hi->exponent-lo->exponent; // recalculate ul=lo->lsd; // .. and repoint } // padding is still > 0, but will fit in acc (less leading carry slot) #if DECCHECK if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding); if (hilen+padding+1>FMALEN) printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding); // printf("FMA padding: %ld\n", (LI)padding); #endif // padding digits can now be set in the result; one or more of // these will come from lo; others will be zeros in the gap for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) { UBFROMUI(ub-3, UBTOUI(ul-3)); // [cannot overlap] } for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul; for (;padding>0; padding--, ub--) *ub=0; // mind the gap } // addition now complete to the right of the rightmost digit of hi uh=hi->lsd; // dow do the add from hi->lsd to the left // [bytewise, because either operand can run out at any time] // carry was set up depending on ten's complement above // first assume both operands have some digits for (;; ub--) { if (uhmsd || ulmsd) break; *ub=(uByte)(carry+(*uh--)+(*ul--)); carry=0; if (*ub<10) continue; *ub-=10; carry=1; } // both loop if (ulmsd) { // to left of lo for (;; ub--) { if (uhmsd) break; *ub=(uByte)(carry+(*uh--)); // [+0] carry=0; if (*ub<10) continue; *ub-=10; carry=1; } // hi loop } else { // to left of hi for (;; ub--) { if (ulmsd) break; *ub=(uByte)(carry+hipad+(*ul--)); carry=0; if (*ub<10) continue; *ub-=10; carry=1; } // lo 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 hisign=hi->sign; // sign is lhs sign for (ul=lo->msd; ullsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul)); for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); // [leaves ul at lsd+1] // complete the ten's complement by adding 1 [cannot overrun] for (ul--; *ul==9; ul--) *ul=0; *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 (; lo->msd+3lsd && UBTOUI(lo->msd)==0;) lo->msd+=4; for (; lo->msdlsd && *lo->msd==0;) 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 #if DECCHECK // assert no left underrun if (lo->msdmsd)); } #endif 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 a uInt */ /* */ /* Many of these could be macros, but having them as real functions */ /* is a little 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 decFloatIsLogical(const decFloat *df) { return DFISUINT01(df) & DFISCC01(df); } uInt decFloatIsNaN(const decFloat *df) { return DFISNAN(df); } uInt decFloatIsNegative(const decFloat *df) { return DFISSIGNED(df) && !DFISZERO(df) && !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 decFloatIsPositive(const decFloat *df) { return !DFISSIGNED(df) && !DFISZERO(df) && !DFISNAN(df); } uInt decFloatIsSignaling(const decFloat *df) { return DFISSNAN(df); } uInt decFloatIsSignalling(const decFloat *df) { return DFISSNAN(df); } uInt decFloatIsSigned(const decFloat *df) { return DFISSIGNED(df)!=0; } uInt decFloatIsSubnormal(const decFloat *df) { if (DFISSPECIAL(df)) return 0; // is finite if (decFloatIsNormal(df)) return 0; // it is 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 754 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 754 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 754 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 754-1985 nextafter, as modified during revision (dropped */ /* from 754-2008); status may be set unless the result is a normal */ /* number. */ /* ------------------------------------------------------------------ */ decFloat * decFloatNextToward(decFloat *result, const decFloat *dfl, const decFloat *dfr, 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) { // lhsround; // 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 uByte *ub, *uc; // work Int drop; // .. uInt dpd; // .. uInt encode; // encoding accumulator uInt sourhil, sourhir; // top words from source decFloats uInt uiwork; // for macros #if QUAD uShort uswork; // .. #endif // the following buffer holds the coefficient for manipulation uByte buf[4+DECPMAX*3+2*QUAD]; // + 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 UBFROMUI(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 (dropstatus|=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 (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); // now at most 3 digits left to non-9 (usually just the one) for (; *ub==9; ub--) *ub=0; *ub+=1; // [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) { UBFROMUI(BUFOFF-8, 0); // must be at least 5 for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0); } } // need round (drop>0) else { // 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 // note that here zeros to the right are added by fours, so in // the Quad case this could write 36 zeros if the coefficient has // fewer than three significant digits (hence the +2*QUAD for buf) for (uc=BUFOFF+DECPMAX;; uc+=4) { UBFROMUI(uc, 0); if (UBTOUI(uc-DECPMAX)!=0) { // could be bad // if all four digits should be zero, definitely bad if (uc<=BUFOFF+DECPMAX+(-drop)-4) return decInvalid(result, set); // must be a 1- to 3-digit sequence; check more carefully if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0) return decInvalid(result, set); break; // no need for loop end test } if (uc>=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 754 */ /* */ /* result gets the result of the operation */ /* dfl is the first decFloat (lhs) */ /* dfr is the second decFloat (rhs), am integer (with q=0) */ /* 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; // .. uInt uiwork; // for macros if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); if (!DFISINT(dfr)) return decInvalid(result, set); digits=decFloatDigits(dfr); // calculate digits if (digits>2) return decInvalid(result, set); // definitely out of range shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; // is in bottom declet if (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 UBFROMUI(buf+DECPMAX, 0); // 8 will handle most cases UBFROMUI(buf+DECPMAX+4, 0); // .. if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); // all other cases num.msd+=shift; num.lsd=num.msd+DECPMAX-1; } savestat=set->status; // record 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 dfldfr */ /* */ /* No error is possible; status and mode are unchanged. */ /* ------------------------------------------------------------------ */ static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) { Int sigl, sigr; // LHS and RHS non-0 signums Int shift; // shift needed to align operands uByte *ub, *uc; // work uInt uiwork; // for macros // buffers +2 if Quad (36 digits), need double plus 4 for safe padding uByte bufl[DECPMAX*2+QUAD*2+4]; // for LHS coefficient + padding uByte bufr[DECPMAX*2+QUAD*2+4]; // for RHS coefficient + padding 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 UBFROMUI(bufl, 0); UBFROMUI(bufr, 0); #endif GETCOEFF(dfl, bufl+QUAD*2); // decode from decFloat GETCOEFF(dfr, bufr+QUAD*2); // .. if (shift==0) { // aligned; common and easy // all multiples of four, here for (ub=bufl, uc=bufr; ub*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 UBFROMUI(bufl+DECPMAX+QUAD*2, 0); // add eight zeros UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); // .. if (shift>8) { // more than eight; fill the rest, and also worth doing the // lead-in by fours uByte *up; // work uByte *upend=bufl+DECPMAX+QUAD*2+shift; for (up=bufl+DECPMAX+QUAD*2+8; upbufl+shift-4) break; } } // check remaining leading digits for (; 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 UBFROMUI(bufr+DECPMAX+QUAD*2, 0); // add eight zeros UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); // .. if (shift<-8) { // more than eight; fill the rest, and also worth doing the // lead-in by fours uByte *up; // work uByte *upend=bufr+DECPMAX+QUAD*2-shift; for (up=bufr+DECPMAX+QUAD*2+8; upbufr-shift-4) break; } } // check remaining leading digits for (; uc*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