| /* Common base code for the decNumber C Library. |
| Copyright (C) 2007-2014 Free Software Foundation, Inc. |
| Contributed by IBM Corporation. Author Mike Cowlishaw. |
| |
| This file is part of GCC. |
| |
| GCC is free software; you can redistribute it and/or modify it under |
| the terms of the GNU General Public License as published by the Free |
| Software Foundation; either version 3, or (at your option) any later |
| version. |
| |
| GCC is distributed in the hope that it will be useful, but WITHOUT ANY |
| WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| for more details. |
| |
| Under Section 7 of GPL version 3, you are granted additional |
| permissions described in the GCC Runtime Library Exception, version |
| 3.1, as published by the Free Software Foundation. |
| |
| You should have received a copy of the GNU General Public License and |
| a copy of the GCC Runtime Library Exception along with this program; |
| see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| <http://www.gnu.org/licenses/>. */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decBasic.c -- common base code for Basic decimal types */ |
| /* ------------------------------------------------------------------ */ |
| /* This module comprises code that is shared between decDouble and */ |
| /* decQuad (but not decSingle). The main arithmetic operations are */ |
| /* here (Add, Subtract, Multiply, FMA, and Division operators). */ |
| /* */ |
| /* Unlike decNumber, parameterization takes place at compile time */ |
| /* rather than at runtime. The parameters are set in the decDouble.c */ |
| /* (etc.) files, which then include this one to produce the compiled */ |
| /* code. The functions here, therefore, are code shared between */ |
| /* multiple formats. */ |
| /* */ |
| /* This must be included after decCommon.c. */ |
| /* ------------------------------------------------------------------ */ |
| /* Names here refer to decFloat rather than to decDouble, etc., and */ |
| /* the functions are in strict alphabetical order. */ |
| |
| /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */ |
| /* decCommon.c */ |
| #if !defined(QUAD) |
| #error decBasic.c must be included after decCommon.c |
| #endif |
| #if SINGLE |
| #error Routines in decBasic.c are for decDouble and decQuad only |
| #endif |
| |
| /* Private constants */ |
| #define DIVIDE 0x80000000 /* Divide operations [as flags] */ |
| #define REMAINDER 0x40000000 /* .. */ |
| #define DIVIDEINT 0x20000000 /* .. */ |
| #define REMNEAR 0x10000000 /* .. */ |
| |
| /* Private functions (local, used only by routines in this module) */ |
| static decFloat *decDivide(decFloat *, const decFloat *, |
| const decFloat *, decContext *, uInt); |
| static decFloat *decCanonical(decFloat *, const decFloat *); |
| static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *, |
| const decFloat *); |
| static decFloat *decInfinity(decFloat *, const decFloat *); |
| static decFloat *decInvalid(decFloat *, decContext *); |
| static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *, |
| decContext *); |
| static Int decNumCompare(const decFloat *, const decFloat *, Flag); |
| static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *, |
| enum rounding, Flag); |
| static uInt decToInt32(const decFloat *, decContext *, enum rounding, |
| Flag, Flag); |
| |
| /* ------------------------------------------------------------------ */ |
| /* decCanonical -- copy a decFloat, making canonical */ |
| /* */ |
| /* result gets the canonicalized df */ |
| /* df is the decFloat to copy and make canonical */ |
| /* returns result */ |
| /* */ |
| /* This is exposed via decFloatCanonical for Double and Quad only. */ |
| /* This works on specials, too; no error or exception is possible. */ |
| /* ------------------------------------------------------------------ */ |
| static decFloat * decCanonical(decFloat *result, const decFloat *df) { |
| uInt encode, precode, dpd; /* work */ |
| uInt inword, uoff, canon; /* .. */ |
| Int n; /* counter (down) */ |
| if (df!=result) *result=*df; /* effect copy if needed */ |
| if (DFISSPECIAL(result)) { |
| if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */ |
| /* is a NaN */ |
| DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */ |
| if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */ |
| /* drop through to check payload */ |
| } |
| /* return quickly if the coefficient continuation is canonical */ |
| { /* declare block */ |
| #if DOUBLE |
| uInt sourhi=DFWORD(df, 0); |
| uInt sourlo=DFWORD(df, 1); |
| if (CANONDPDOFF(sourhi, 8) |
| && CANONDPDTWO(sourhi, sourlo, 30) |
| && CANONDPDOFF(sourlo, 20) |
| && CANONDPDOFF(sourlo, 10) |
| && CANONDPDOFF(sourlo, 0)) return result; |
| #elif QUAD |
| uInt sourhi=DFWORD(df, 0); |
| uInt sourmh=DFWORD(df, 1); |
| uInt sourml=DFWORD(df, 2); |
| uInt sourlo=DFWORD(df, 3); |
| if (CANONDPDOFF(sourhi, 4) |
| && CANONDPDTWO(sourhi, sourmh, 26) |
| && CANONDPDOFF(sourmh, 16) |
| && CANONDPDOFF(sourmh, 6) |
| && CANONDPDTWO(sourmh, sourml, 28) |
| && CANONDPDOFF(sourml, 18) |
| && CANONDPDOFF(sourml, 8) |
| && CANONDPDTWO(sourml, sourlo, 30) |
| && CANONDPDOFF(sourlo, 20) |
| && CANONDPDOFF(sourlo, 10) |
| && CANONDPDOFF(sourlo, 0)) return result; |
| #endif |
| } /* block */ |
| |
| /* Loop to repair a non-canonical coefficent, as needed */ |
| inword=DECWORDS-1; /* current input word */ |
| uoff=0; /* bit offset of declet */ |
| encode=DFWORD(result, inword); |
| for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */ |
| dpd=encode>>uoff; |
| uoff+=10; |
| if (uoff>32) { /* crossed uInt boundary */ |
| inword--; |
| encode=DFWORD(result, inword); |
| uoff-=32; |
| dpd|=encode<<(10-uoff); /* get pending bits */ |
| } |
| dpd&=0x3ff; /* clear uninteresting bits */ |
| if (dpd<0x16e) continue; /* must be canonical */ |
| canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */ |
| if (canon==dpd) continue; /* have canonical declet */ |
| /* need to replace declet */ |
| if (uoff>=10) { /* all within current word */ |
| encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */ |
| encode|=canon<<(uoff-10); /* insert the canonical form */ |
| DFWORD(result, inword)=encode; /* .. and save */ |
| continue; |
| } |
| /* straddled words */ |
| precode=DFWORD(result, inword+1); /* get previous */ |
| precode&=0xffffffff>>(10-uoff); /* clear top bits */ |
| DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff))); |
| encode&=0xffffffff<<uoff; /* clear bottom bits */ |
| encode|=canon>>(10-uoff); /* insert canonical */ |
| DFWORD(result, inword)=encode; /* .. and save */ |
| } /* n */ |
| return result; |
| } /* decCanonical */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decDivide -- divide operations */ |
| /* */ |
| /* result gets the result of dividing dfl by dfr: */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* set is the context */ |
| /* op is the operation selector */ |
| /* returns result */ |
| /* */ |
| /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */ |
| /* ------------------------------------------------------------------ */ |
| #define DIVCOUNT 0 /* 1 to instrument subtractions counter */ |
| #define DIVBASE ((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<acc) { |
| printf("Acc underrun...\n"); |
| break; |
| } |
| #endif |
| #if DECTRACE |
| printf("Outer: quodigits=%ld acc=", (LI)quodigits); |
| for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua); |
| printf("\n"); |
| #endif |
| *lsuq=0; /* default unit result is 0 */ |
| for (;;) { /* inner loop -- calculate quotient unit */ |
| /* strip leading zero units from acc (either there initially or */ |
| /* from subtraction below); this may strip all if exactly 0 */ |
| for (; *msua==0 && msua>=lsua;) msua--; |
| accunits=(Int)(msua-lsua+1); /* [maybe 0] */ |
| /* subtraction is only necessary and possible if there are as */ |
| /* least as many units remaining in acc for this iteration as */ |
| /* there are in div */ |
| if (accunits<divunits) { |
| if (accunits==0) msua++; /* restore */ |
| break; |
| } |
| |
| /* If acc is longer than div then subtraction is definitely */ |
| /* possible (as msu of both is non-zero), but if they are the */ |
| /* same length a comparison is needed. */ |
| /* If a subtraction is needed then a good estimate of the */ |
| /* multiplier for the subtraction is also needed in order to */ |
| /* minimise the iterations of this inner loop because the */ |
| /* subtractions needed dominate division performance. */ |
| if (accunits==divunits) { |
| /* compare the high divunits of acc and div: */ |
| /* acc<div: this quotient unit is unchanged; subtraction */ |
| /* will be possible on the next iteration */ |
| /* acc==div: quotient gains 1, set acc=0 */ |
| /* acc>div: subtraction necessary at this position */ |
| for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break; |
| /* [now at first mismatch or lsu] */ |
| if (*ud>*ua) break; /* next time... */ |
| if (*ud==*ua) { /* all compared equal */ |
| *lsuq+=1; /* increment result */ |
| msua=lsua; /* collapse acc units */ |
| *msua=0; /* .. to a zero */ |
| break; |
| } |
| |
| /* subtraction necessary; estimate multiplier [see above] */ |
| /* if both *msud and *msua are small it is cost-effective to */ |
| /* bring in part of the following units (if any) to get a */ |
| /* better estimate (assume some other non-zero in div) */ |
| #define DIVLO 1000000U |
| #define DIVHI (DIVBASE/DIVLO) |
| #if DECUSE64 |
| if (divunits>1) { |
| /* there cannot be a *(msud-2) for DECDOUBLE so next is */ |
| /* an exact calculation unless DECQUAD (which needs to */ |
| /* assume bits out there if divunits>2) */ |
| uLong mul=(uLong)*msua * DIVBASE + *(msua-1); |
| uLong div=(uLong)*msud * DIVBASE + *(msud-1); |
| #if QUAD |
| if (divunits>2) div++; |
| #endif |
| mul/=div; |
| multiplier=(Int)mul; |
| } |
| else multiplier=*msua/(*msud); |
| #else |
| if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { |
| multiplier=(*msua*DIVHI + *(msua-1)/DIVLO) |
| /(*msud*DIVHI + *(msud-1)/DIVLO +1); |
| } |
| else multiplier=(*msua<<2)/divtop; |
| #endif |
| } |
| else { /* accunits>divunits */ |
| /* msud is one unit 'lower' than msua, so estimate differently */ |
| #if DECUSE64 |
| uLong mul; |
| /* as before, bring in extra digits if possible */ |
| if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { |
| mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI |
| + *(msua-2)/DIVLO; |
| mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1); |
| } |
| else if (divunits==1) { |
| mul=(uLong)*msua * DIVBASE + *(msua-1); |
| mul/=*msud; /* no more to the right */ |
| } |
| else { |
| mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) |
| + (*(msua-1)<<2); |
| mul/=divtop; /* [divtop already allows for sticky bits] */ |
| } |
| multiplier=(Int)mul; |
| #else |
| multiplier=*msua * ((DIVBASE<<2)/divtop); |
| #endif |
| } |
| if (multiplier==0) multiplier=1; /* marginal case */ |
| *lsuq+=multiplier; |
| |
| #if DIVCOUNT |
| /* printf("Multiplier: %ld\n", (LI)multiplier); */ |
| divcount++; |
| #endif |
| |
| /* Carry out the subtraction acc-(div*multiplier); for each */ |
| /* unit in div, do the multiply, split to units (see */ |
| /* decFloatMultiply for the algorithm), and subtract from acc */ |
| #define DIVMAGIC 2305843009U /* 2**61/10**9 */ |
| #define DIVSHIFTA 29 |
| #define DIVSHIFTB 32 |
| carry=0; |
| for (ud=div, ua=lsua; ud<=msud; ud++, ua++) { |
| uInt lo, hop; |
| #if DECUSE64 |
| uLong sub=(uLong)multiplier*(*ud)+carry; |
| if (sub<DIVBASE) { |
| carry=0; |
| lo=(uInt)sub; |
| } |
| else { |
| hop=(uInt)(sub>>DIVSHIFTA); |
| carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB); |
| /* the estimate is now in hi; now calculate sub-hi*10**9 */ |
| /* to get the remainder (which will be <DIVBASE)) */ |
| lo=(uInt)sub; |
| lo-=carry*DIVBASE; /* low word of result */ |
| if (lo>=DIVBASE) { |
| lo-=DIVBASE; /* correct by +1 */ |
| carry++; |
| } |
| } |
| #else /* 32-bit */ |
| uInt hi; |
| /* calculate multiplier*(*ud) into hi and lo */ |
| LONGMUL32HI(hi, *ud, multiplier); /* get the high word */ |
| lo=multiplier*(*ud); /* .. and the low */ |
| lo+=carry; /* add the old hi */ |
| carry=hi+(lo<carry); /* .. with any carry */ |
| if (carry || lo>=DIVBASE) { /* split is needed */ |
| hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */ |
| LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */ |
| /* [DIVSHIFTB is 32, so carry can be used directly] */ |
| /* the estimate is now in carry; now calculate hi:lo-est*10**9; */ |
| /* happily the top word of the result is irrelevant because it */ |
| /* will always be zero so this needs only one multiplication */ |
| lo-=(carry*DIVBASE); |
| /* the correction here will be at most +1; do it */ |
| if (lo>=DIVBASE) { |
| lo-=DIVBASE; |
| carry++; |
| } |
| } |
| #endif |
| if (lo>*ua) { /* borrow needed */ |
| *ua+=DIVBASE; |
| carry++; |
| } |
| *ua-=lo; |
| } /* ud loop */ |
| if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */ |
| } /* inner loop */ |
| |
| /* the outer loop terminates when there is either an exact result */ |
| /* or enough digits; first update the quotient digit count and */ |
| /* pointer (if any significant digits) */ |
| #if DECTRACE |
| if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq); |
| #endif |
| if (quodigits) { |
| quodigits+=9; /* had leading unit earlier */ |
| lsuq--; |
| if (quodigits>DECPMAX+1) break; /* have enough */ |
| } |
| else if (*lsuq) { /* first quotient digits */ |
| const uInt *pow; |
| for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++; |
| lsuq--; |
| /* [cannot have >DECPMAX+1 on first unit] */ |
| } |
| |
| if (*msua!=0) continue; /* not an exact result */ |
| /* acc is zero iff used all of original units and zero down to lsua */ |
| /* (must also continue to original lsu for correct quotient length) */ |
| if (lsua>acc+DIVACCLEN-DIVOPLEN) continue; |
| for (; msua>lsua && *msua==0;) msua--; |
| if (*msua==0 && msua==lsua) break; |
| } /* outer loop */ |
| |
| /* all of the original operand in acc has been covered at this point */ |
| /* quotient now has at least DECPMAX+2 digits */ |
| /* *msua is now non-0 if inexact and sticky bits */ |
| /* lsuq is one below the last uint of the quotient */ |
| lsuq++; /* set -> true lsu of quo */ |
| if (*msua) *lsuq|=1; /* apply sticky bit */ |
| |
| /* quo now holds the (unrounded) quotient in base-billion; one */ |
| /* base-billion 'digit' per uInt. */ |
| #if DECTRACE |
| printf("DivQuo:"); |
| for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq); |
| printf("\n"); |
| #endif |
| |
| /* Now convert to BCD for rounding and cleanup, starting from the */ |
| /* most significant end [offset by one into bcdacc to leave room */ |
| /* for a possible carry digit if rounding for REMNEAR is needed] */ |
| for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) { |
| uInt top, mid, rem; /* work */ |
| if (*uq==0) { /* no split needed */ |
| 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 (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */ |
| num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9); |
| /* if the result was exact then there may be up to 8 extra */ |
| /* trailing zeros in the overflowed quotient final unit */ |
| if (*msua==0) { |
| for (; *ub==0;) ub--; /* drop zeros */ |
| num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */ |
| num.lsd=ub; |
| } |
| } /* adjustment needed */ |
| |
| #if DIVCOUNT |
| if (divcount>maxcount) { /* new high-water nark */ |
| maxcount=divcount; |
| printf("DivNewMaxCount: %ld\n", (LI)maxcount); |
| } |
| #endif |
| |
| if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */ |
| |
| /* Is DIVIDEINT or a remainder; there is more to do -- first form */ |
| /* the integer (this is done 'after the fact', unlike as in */ |
| /* decNumber, so as not to tax DIVIDE) */ |
| |
| /* The first non-zero digit will be in the first 9 digits, known */ |
| /* from quodigits and num.msd, so there is always space for DECPMAX */ |
| /* digits */ |
| |
| length=(Int)(num.lsd-num.msd+1); |
| /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */ |
| |
| if (length+num.exponent>DECPMAX) { /* cannot fit */ |
| decFloatZero(result); |
| DFWORD(result, 0)=DECFLOAT_qNaN; |
| set->status|=DEC_Division_impossible; |
| return result; |
| } |
| |
| if (num.exponent>=0) { /* already an int, or need pad zeros */ |
| for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0; |
| num.lsd+=num.exponent; |
| } |
| else { /* too long: round or truncate needed */ |
| Int drop=-num.exponent; |
| if (!(op&REMNEAR)) { /* simple truncate */ |
| num.lsd-=drop; |
| if (num.lsd<num.msd) { /* truncated all */ |
| num.lsd=num.msd; /* make 0 */ |
| *num.lsd=0; /* .. [sign still relevant] */ |
| } |
| } |
| else { /* round to nearest even [sigh] */ |
| /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */ |
| /* (this is a special case of Quantize -- q.v. for commentary) */ |
| uByte *roundat; /* -> re-round digit */ |
| uByte reround; /* reround value */ |
| *(num.msd-1)=0; /* in case of left carry, or make 0 */ |
| if (drop<length) roundat=num.lsd-drop+1; |
| else if (drop==length) roundat=num.msd; |
| else roundat=num.msd-1; /* [-> 0] */ |
| reround=*roundat; |
| for (ub=roundat+1; ub<=num.lsd; ub++) { |
| if (*ub!=0) { |
| reround=DECSTICKYTAB[reround]; |
| break; |
| } |
| } /* check stickies */ |
| if (roundat>num.msd) num.lsd=roundat-1; |
| else { |
| num.msd--; /* use the 0 .. */ |
| num.lsd=num.msd; /* .. at the new MSD place */ |
| } |
| if (reround!=0) { /* discarding non-zero */ |
| uInt bump=0; |
| /* rounding is DEC_ROUND_HALF_EVEN always */ |
| if (reround>5) bump=1; /* >0.5 goes up */ |
| else if (reround==5) /* exactly 0.5000 .. */ |
| bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */ |
| if (bump!=0) { /* need increment */ |
| /* increment the coefficient; this might end up with 1000... */ |
| ub=num.lsd; |
| for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); |
| for (; *ub==9; ub--) *ub=0; /* at most 3 more */ |
| *ub+=1; |
| if (ub<num.msd) num.msd--; /* carried */ |
| } /* bump needed */ |
| } /* reround!=0 */ |
| } /* remnear */ |
| } /* round or truncate needed */ |
| num.exponent=0; /* all paths */ |
| /*decShowNum(&num, "int"); */ |
| |
| if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */ |
| |
| /* Have a remainder to calculate */ |
| decFinalize("ient, &num, set); /* lay out the integer so far */ |
| DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */ |
| sign=DFWORD(dfl, 0); /* save sign of dfl */ |
| decFloatFMA(result, "ient, dfr, dfl, set); |
| if (!DFISZERO(result)) return result; |
| /* if the result is zero the sign shall be sign of dfl */ |
| DFWORD("ient, 0)=sign; /* construct decFloat of sign */ |
| return decFloatCopySign(result, result, "ient); |
| } /* decDivide */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFiniteMultiply -- multiply two finite decFloats */ |
| /* */ |
| /* num gets the result of multiplying dfl and dfr */ |
| /* bcdacc .. with the coefficient in this array */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* */ |
| /* This effects the multiplication of two decFloats, both known to be */ |
| /* finite, leaving the result in a bcdnum ready for decFinalize (for */ |
| /* use in Multiply) or in a following addition (FMA). */ |
| /* */ |
| /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */ |
| /* No error is possible and no status is set. */ |
| /* ------------------------------------------------------------------ */ |
| /* This routine has two separate implementations of the core */ |
| /* multiplication; both using base-billion. One uses only 32-bit */ |
| /* variables (Ints and uInts) or smaller; the other uses uLongs (for */ |
| /* multiplication and addition only). Both implementations cover */ |
| /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */ |
| /* comparisons. In any one compilation only one implementation for */ |
| /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */ |
| /* version is forced. */ |
| /* */ |
| /* Historical note: an earlier version of this code also supported the */ |
| /* 256-bit format and has been preserved. That is somewhat trickier */ |
| /* during lazy carry splitting because the initial quotient estimate */ |
| /* (est) can exceed 32 bits. */ |
| |
| #define MULTBASE ((uInt)BILLION) /* the base used for multiply */ |
| #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */ |
| #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */ |
| #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */ |
| |
| /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */ |
| #if DECEMAXD>9 |
| #error Exponent may overflow when doubled for Multiply |
| #endif |
| #if MULACCLEN!=(MULACCLEN/4)*4 |
| /* This assumption is used below only for initialization */ |
| #error MULACCLEN is not a multiple of 4 |
| #endif |
| |
| static void decFiniteMultiply(bcdnum *num, uByte *bcdacc, |
| const decFloat *dfl, const decFloat *dfr) { |
| uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */ |
| uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */ |
| uInt *ui, *uj; /* work */ |
| uByte *ub; /* .. */ |
| 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; pl<accl+MULACCLEN; pl+=4) { |
| *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */ |
| } /* pl */ |
| #endif |
| |
| /* Effect the multiplication */ |
| /* The multiplcation proceeds using MFC's lazy-carry resolution */ |
| /* algorithm from decNumber. First, the multiplication is */ |
| /* effected, allowing accumulation of the partial products (which */ |
| /* are in base-billion at each column position) into 64 bits */ |
| /* without resolving back to base=billion after each addition. */ |
| /* These 64-bit numbers (which may contain up to 19 decimal digits) */ |
| /* are then split using the Clark & Cowlishaw algorithm (see below). */ |
| /* [Testing for 0 in the inner loop is not really a 'win'] */ |
| for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */ |
| if (*ui==0) continue; /* product cannot affect result */ |
| pl=accl+(ui-bufr); /* where to add the lhs */ |
| for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */ |
| /* if (*uj==0) continue; // product cannot affect result */ |
| *pl+=((uLong)*ui)*(*uj); |
| } /* uj */ |
| } /* ui */ |
| |
| /* The 64-bit carries must now be resolved; this means that a */ |
| /* quotient/remainder has to be calculated for base-billion (1E+9). */ |
| /* For this, Clark & Cowlishaw's quotient estimation approach (also */ |
| /* used in decNumber) is needed, because 64-bit divide is generally */ |
| /* extremely slow on 32-bit machines, and may be slower than this */ |
| /* approach even on 64-bit machines. This algorithm splits X */ |
| /* using: */ |
| /* */ |
| /* magic=2**(A+B)/1E+9; // 'magic number' */ |
| /* hop=X/2**A; // high order part of X (by shift) */ |
| /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ |
| /* */ |
| /* A and B are quite constrained; hop and magic must fit in 32 bits, */ |
| /* and 2**(A+B) must be as large as possible (which is 2**61 if */ |
| /* magic is to fit). Further, maxX increases with the length of */ |
| /* the operands (and hence the number of partial products */ |
| /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ |
| /* */ |
| /* It can be shown that when OPLEN is 2 then the maximum error in */ |
| /* the estimated quotient is <1, but for larger maximum x the */ |
| /* maximum error is above 1 so a correction that is >1 may be */ |
| /* needed. Values of A and B are chosen to satisfy the constraints */ |
| /* just mentioned while minimizing the maximum error (and hence the */ |
| /* maximum correction), as shown in the following table: */ |
| /* */ |
| /* Type OPLEN A B maxX maxError maxCorrection */ |
| /* --------------------------------------------------------- */ |
| /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ |
| /* QUAD 4 30 31 <4*10**18 1.17 2 */ |
| /* */ |
| /* In the OPLEN==2 case there is most choice, but the value for B */ |
| /* of 32 has a big advantage as then the calculation of the */ |
| /* estimate requires no shifting; the compiler can extract the high */ |
| /* word directly after multiplying magic*hop. */ |
| #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ |
| #if DOUBLE |
| #define MULSHIFTA 29 |
| #define MULSHIFTB 32 |
| #elif QUAD |
| #define MULSHIFTA 30 |
| #define MULSHIFTB 31 |
| #else |
| #error Unexpected type |
| #endif |
| |
| #if DECTRACE |
| printf("MulAccl:"); |
| for (pl=accl+MULACCLEN-1; pl>=accl; pl--) |
| printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff)); |
| printf("\n"); |
| #endif |
| |
| for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */ |
| uInt lo, hop; /* work */ |
| uInt est; /* cannot exceed 4E+9 */ |
| if (*pl>=MULTBASE) { |
| /* *pl holds a binary number which needs to be split */ |
| hop=(uInt)(*pl>>MULSHIFTA); |
| est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB); |
| /* the estimate is now in est; now calculate hi:lo-est*10**9; */ |
| /* happily the top word of the result is irrelevant because it */ |
| /* will always be zero so this needs only one multiplication */ |
| lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */ |
| /* If QUAD, the correction here could be +2 */ |
| if (lo>=MULTBASE) { |
| lo-=MULTBASE; /* correct by +1 */ |
| est++; |
| #if QUAD |
| /* may need to correct by +2 */ |
| if (lo>=MULTBASE) { |
| lo-=MULTBASE; |
| est++; |
| } |
| #endif |
| } |
| /* finally place lo as the new coefficient 'digit' and add est to */ |
| /* the next place up [this is safe because this path is never */ |
| /* taken on the final iteration as *pl will fit] */ |
| *pa=lo; |
| *(pl+1)+=est; |
| } /* *pl needed split */ |
| else { /* *pl<MULTBASE */ |
| *pa=(uInt)*pl; /* just copy across */ |
| } |
| } /* pl loop */ |
| |
| #else /* 32-bit */ |
| for (pa=acc;; pa+=4) { /* zero the accumulator */ |
| *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */ |
| if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */ |
| } /* pa */ |
| |
| /* Effect the multiplication */ |
| /* uLongs are not available (and in particular, there is no uLong */ |
| /* divide) but it is still possible to use MFC's lazy-carry */ |
| /* resolution algorithm from decNumber. First, the multiplication */ |
| /* is effected, allowing accumulation of the partial products */ |
| /* (which are in base-billion at each column position) into 64 bits */ |
| /* [with the high-order 32 bits in each position being held at */ |
| /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */ |
| /* These 64-bit numbers (which may contain up to 19 decimal digits) */ |
| /* are then split using the Clark & Cowlishaw algorithm (see */ |
| /* below). */ |
| for (ui=bufr;; ui++) { /* over each item in rhs */ |
| uInt hi, lo; /* words of exact multiply result */ |
| pa=acc+(ui-bufr); /* where to add the lhs */ |
| for (uj=bufl;; uj++, pa++) { /* over each item in lhs */ |
| LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */ |
| lo=(*ui)*(*uj); /* .. */ |
| *pa+=lo; /* accumulate low bits and .. */ |
| *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */ |
| if (uj==bufl+MULOPLEN-1) break; |
| } |
| if (ui==bufr+MULOPLEN-1) break; |
| } |
| |
| /* The 64-bit carries must now be resolved; this means that a */ |
| /* quotient/remainder has to be calculated for base-billion (1E+9). */ |
| /* For this, Clark & Cowlishaw's quotient estimation approach (also */ |
| /* used in decNumber) is needed, because 64-bit divide is generally */ |
| /* extremely slow on 32-bit machines. This algorithm splits X */ |
| /* using: */ |
| /* */ |
| /* magic=2**(A+B)/1E+9; // 'magic number' */ |
| /* hop=X/2**A; // high order part of X (by shift) */ |
| /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ |
| /* */ |
| /* A and B are quite constrained; hop and magic must fit in 32 bits, */ |
| /* and 2**(A+B) must be as large as possible (which is 2**61 if */ |
| /* magic is to fit). Further, maxX increases with the length of */ |
| /* the operands (and hence the number of partial products */ |
| /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ |
| /* */ |
| /* It can be shown that when OPLEN is 2 then the maximum error in */ |
| /* the estimated quotient is <1, but for larger maximum x the */ |
| /* maximum error is above 1 so a correction that is >1 may be */ |
| /* needed. Values of A and B are chosen to satisfy the constraints */ |
| /* just mentioned while minimizing the maximum error (and hence the */ |
| /* maximum correction), as shown in the following table: */ |
| /* */ |
| /* Type OPLEN A B maxX maxError maxCorrection */ |
| /* --------------------------------------------------------- */ |
| /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ |
| /* QUAD 4 30 31 <4*10**18 1.17 2 */ |
| /* */ |
| /* In the OPLEN==2 case there is most choice, but the value for B */ |
| /* of 32 has a big advantage as then the calculation of the */ |
| /* estimate requires no shifting; the high word is simply */ |
| /* calculated from multiplying magic*hop. */ |
| #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ |
| #if DOUBLE |
| #define MULSHIFTA 29 |
| #define MULSHIFTB 32 |
| #elif QUAD |
| #define MULSHIFTA 30 |
| #define MULSHIFTB 31 |
| #else |
| #error Unexpected type |
| #endif |
| |
| #if DECTRACE |
| printf("MulHiLo:"); |
| for (pa=acc+MULACCLEN-1; pa>=acc; pa--) |
| printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa); |
| printf("\n"); |
| #endif |
| |
| for (pa=acc;; pa++) { /* each low uInt */ |
| uInt hi, lo; /* words of exact multiply result */ |
| uInt hop, estlo; /* work */ |
| #if QUAD |
| uInt esthi; /* .. */ |
| #endif |
| |
| lo=*pa; |
| hi=*(pa+MULACCLEN); /* top 32 bits */ |
| /* hi and lo now hold a binary number which needs to be split */ |
| |
| #if DOUBLE |
| hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */ |
| LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */ |
| /* [MULSHIFTB is 32, so estlo can be used directly] */ |
| /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */ |
| /* happily the top word of the result is irrelevant because it */ |
| /* will always be zero so this needs only one multiplication */ |
| lo-=(estlo*MULTBASE); |
| /* esthi=0; // high word is ignored below */ |
| /* the correction here will be at most +1; do it */ |
| if (lo>=MULTBASE) { |
| lo-=MULTBASE; |
| estlo++; |
| } |
| #elif QUAD |
| hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */ |
| LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */ |
| estlo=hop*MULMAGIC; /* .. so low word needed */ |
| estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */ |
| /* esthi=0; // high word is ignored below */ |
| lo-=(estlo*MULTBASE); /* as above */ |
| /* the correction here could be +1 or +2 */ |
| if (lo>=MULTBASE) { |
| lo-=MULTBASE; |
| estlo++; |
| } |
| if (lo>=MULTBASE) { |
| lo-=MULTBASE; |
| estlo++; |
| } |
| #else |
| #error Unexpected type |
| #endif |
| |
| /* finally place lo as the new accumulator digit and add est to */ |
| /* the next place up; this latter add could cause a carry of 1 */ |
| /* to the high word of the next place */ |
| *pa=lo; |
| *(pa+1)+=estlo; |
| /* esthi is always 0 for DOUBLE and QUAD so this is skipped */ |
| /* *(pa+1+MULACCLEN)+=esthi; */ |
| if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */ |
| if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */ |
| } /* pa loop */ |
| #endif |
| |
| /* At this point, whether using the 64-bit or the 32-bit paths, the */ |
| /* accumulator now holds the (unrounded) result in base-billion; */ |
| /* one base-billion 'digit' per uInt. */ |
| #if DECTRACE |
| printf("MultAcc:"); |
| for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa); |
| printf("\n"); |
| #endif |
| |
| /* Now convert to BCD for rounding and cleanup, starting from the */ |
| /* most significant end */ |
| pa=acc+MULACCLEN-1; |
| if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */ |
| else { /* >=1 word of leading zeros */ |
| num->msd=bcdacc; /* known leading zeros are gone */ |
| pa--; /* skip first word .. */ |
| for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */ |
| } |
| for (ub=bcdacc;; pa--, ub+=9) { |
| if (*pa!=0) { /* split(s) needed */ |
| uInt top, mid, rem; /* work */ |
| /* *pa is non-zero -- split the base-billion acc digit into */ |
| /* hi, mid, and low three-digits */ |
| #define mulsplit9 1000000 /* divisor */ |
| #define mulsplit6 1000 /* divisor */ |
| /* The splitting is done by simple divides and remainders, */ |
| /* assuming the compiler will optimize these where useful */ |
| /* [GCC does] */ |
| top=*pa/mulsplit9; |
| rem=*pa%mulsplit9; |
| mid=rem/mulsplit6; |
| rem=rem%mulsplit6; |
| /* lay out the nine BCD digits (plus one unwanted byte) */ |
| 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<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */ |
| if (overlap<-DECPMAX) { /* gap was > DECPMAX */ |
| *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */ |
| } |
| else { /* need full coefficient */ |
| GETCOEFF(dfr, ub); /* decode from decFloat */ |
| ub+=DECPMAX-1; /* new LSD... */ |
| } |
| ulsd=ub; /* save new LSD */ |
| } /* no overlap possible */ |
| |
| else { /* overlap>0 */ |
| /* coefficients overlap (perhaps completely, although also */ |
| /* perhaps only where zeros) */ |
| if (overlap==DECPMAX) { /* aligned */ |
| ub=buf+COFF; /* where msd will go */ |
| #if QUAD |
| UBFROMUS(buf+4, 0); /* clear quad's 00 */ |
| #endif |
| GETCOEFF(dfr, ub); /* decode from decFloat */ |
| } |
| else { /* unaligned */ |
| ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */ |
| /* Fill the prefix gap with 0s; 8 will cover most common */ |
| /* unalignments, so start with direct assignments (a loop is */ |
| /* then used for any remaining -- the loop (and the one in a */ |
| /* moment) is not then on the critical path because the number */ |
| /* of additions is reduced by (at least) two in this case) */ |
| UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */ |
| UBFROMUI(buf+8, 0); |
| if (ub>buf+12) { |
| ut=buf+12; /* start any remaining */ |
| for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */ |
| } |
| GETCOEFF(dfr, ub); /* decode from decFloat */ |
| |
| /* now move tail of rhs across to main acc; again use direct */ |
| /* copies for 8 digits-worth */ |
| UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX)); |
| UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4)); |
| if (buf+COFF+DECPMAX+8<ub+DECPMAX) { |
| us=buf+COFF+DECPMAX+8; /* source */ |
| ut=acc+COFF+DECPMAX+8; /* target */ |
| for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us)); |
| } |
| } /* unaligned */ |
| |
| ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */ |
| |
| /* Now do the add of the non-tail; this is all nicely aligned, */ |
| /* and is over a multiple of four digits (because for Quad two */ |
| /* zero digits were added on the left); words in both acc and */ |
| /* buf (buf especially) will often be zero */ |
| /* [byte-by-byte add, here, is about 15% slower total effect than */ |
| /* the by-fours] */ |
| |
| /* Now effect the add; this is harder on a little-endian */ |
| /* machine as the inter-digit carry cannot use the usual BCD */ |
| /* addition trick because the bytes are loaded in the wrong order */ |
| /* [this loop could be unrolled, but probably scarcely worth it] */ |
| |
| ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */ |
| us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */ |
| |
| #if !DECLITEND |
| 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<LHS */ |
| /* borrowed -- take ten's complement */ |
| /* sign is lhs sign */ |
| num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign; |
| |
| /* invert the coefficient first by fours, then add one; space */ |
| /* at the end of the buffer ensures the by-fours is always */ |
| /* safe, but lsd+1 must be cleared to prevent a borrow */ |
| /* if big-endian */ |
| #if !DECLITEND |
| *(ulsd+1)=0; |
| #endif |
| /* there are always at least four coefficient words */ |
| UBFROMUI(umsd, 0x09090909-UBTOUI(umsd)); |
| UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4)); |
| UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8)); |
| UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12)); |
| #if DOUBLE |
| #define BNEXT 16 |
| #elif QUAD |
| UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16)); |
| UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20)); |
| UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24)); |
| UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28)); |
| UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32)); |
| #define BNEXT 36 |
| #endif |
| if (ulsd>=umsd+BNEXT) { /* unaligned */ |
| /* eight will handle most unaligments for Double; 16 for Quad */ |
| 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 (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4; |
| for (; *umsd==0 && umsd<ulsd;) umsd++; |
| } |
| if (*umsd==0) { /* must be true zero (and diffsign) */ |
| num.sign=0; /* assume + */ |
| if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign; |
| } |
| } |
| /* [else was not zero, might still have leading zeros] */ |
| } /* subtraction gave positive result */ |
| } /* diffsign */ |
| |
| else { /* same-sign addition */ |
| num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; |
| #if DOUBLE |
| if (carry) { /* only possible with decDouble */ |
| *(acc+3)=1; /* [Quad has leading 00] */ |
| umsd=acc+3; |
| } |
| #endif |
| } /* same sign */ |
| |
| num.msd=umsd; /* set MSD .. */ |
| num.lsd=ulsd; /* .. and LSD */ |
| num.exponent=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<nanr) comp=-1; |
| else { /* NaNs are the same type and sign .. must compare payload */ |
| /* buffers need +2 for QUAD */ |
| uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */ |
| uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */ |
| uByte *ub, *uc; /* work */ |
| Int sigl; /* signum of LHS */ |
| sigl=(DFISSIGNED(dfl) ? -1 : +1); |
| |
| /* decode the coefficients */ |
| /* (shift both right two if Quad to make a multiple of four) */ |
| #if QUAD |
| UBFROMUS(bufl, 0); |
| UBFROMUS(bufr, 0); |
| #endif |
| GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */ |
| GETCOEFF(dfr, bufr+QUAD*2); /* .. */ |
| /* all multiples of four, here */ |
| comp=0; /* assume equal */ |
| for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) { |
| uInt ui=UBTOUI(ub); |
| if (ui==UBTOUI(uc)) continue; /* so far so same */ |
| /* about to find a winner; go by bytes in case little-endian */ |
| for (;; ub++, uc++) { |
| if (*ub==*uc) continue; |
| if (*ub>*uc) comp=sigl; /* difference found */ |
| else comp=-sigl; /* .. */ |
| break; |
| } |
| } |
| } /* same NaN type and sign */ |
| } |
| else { |
| /* numeric comparison needed */ |
| comp=decNumCompare(dfl, dfr, 1); /* total ordering */ |
| } |
| decFloatZero(result); |
| if (comp==0) return result; |
| DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */ |
| if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */ |
| return result; |
| } /* decFloatCompareTotal */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatCompareTotalMag -- compare magnitudes with total ordering */ |
| /* */ |
| /* result gets the result of comparing abs(dfl) and abs(dfr) */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* returns result, which may be -1, 0, or 1 */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatCompareTotalMag(decFloat *result, |
| const decFloat *dfl, const decFloat *dfr) { |
| decFloat a, b; /* for copy if needed */ |
| /* copy and redirect signed operand(s) */ |
| if (DFISSIGNED(dfl)) { |
| decFloatCopyAbs(&a, dfl); |
| dfl=&a; |
| } |
| if (DFISSIGNED(dfr)) { |
| decFloatCopyAbs(&b, dfr); |
| dfr=&b; |
| } |
| return decFloatCompareTotal(result, dfl, dfr); |
| } /* decFloatCompareTotalMag */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatCopy -- copy a decFloat as-is */ |
| /* */ |
| /* result gets the copy of dfl */ |
| /* dfl is the decFloat to copy */ |
| /* returns result */ |
| /* */ |
| /* This is a bitwise operation; no errors or exceptions are possible. */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) { |
| if (dfl!=result) *result=*dfl; /* copy needed */ |
| return result; |
| } /* decFloatCopy */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */ |
| /* */ |
| /* result gets the copy of dfl with sign bit 0 */ |
| /* dfl is the decFloat to copy */ |
| /* returns result */ |
| /* */ |
| /* This is a bitwise operation; no errors or exceptions are possible. */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) { |
| if (dfl!=result) *result=*dfl; /* copy needed */ |
| DFBYTE(result, 0)&=~0x80; /* zero sign bit */ |
| return result; |
| } /* decFloatCopyAbs */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */ |
| /* */ |
| /* result gets the copy of dfl with sign bit inverted */ |
| /* dfl is the decFloat to copy */ |
| /* returns result */ |
| /* */ |
| /* This is a bitwise operation; no errors or exceptions are possible. */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) { |
| if (dfl!=result) *result=*dfl; /* copy needed */ |
| DFBYTE(result, 0)^=0x80; /* invert sign bit */ |
| return result; |
| } /* decFloatCopyNegate */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatCopySign -- copy a decFloat with the sign of another */ |
| /* */ |
| /* result gets the result of copying dfl with the sign of dfr */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* returns result */ |
| /* */ |
| /* This is a bitwise operation; no errors or exceptions are possible. */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatCopySign(decFloat *result, |
| const decFloat *dfl, const decFloat *dfr) { |
| uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */ |
| if (dfl!=result) *result=*dfl; /* copy needed */ |
| DFBYTE(result, 0)&=~0x80; /* clear sign .. */ |
| DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */ |
| return result; |
| } /* decFloatCopySign */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatDigits -- return the number of digits in a decFloat */ |
| /* */ |
| /* df is the decFloat to investigate */ |
| /* returns the number of significant digits in the decFloat; a */ |
| /* zero coefficient returns 1 as does an infinity (a NaN returns */ |
| /* the number of digits in the payload) */ |
| /* ------------------------------------------------------------------ */ |
| /* private macro to extract a declet according to provided formula */ |
| /* (form), and if it is non-zero then return the calculated digits */ |
| /* depending on the declet number (n), where n=0 for the most */ |
| /* significant declet; uses uInt dpd for work */ |
| #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \ |
| if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);} |
| /* next one is used when it is known that the declet must be */ |
| /* non-zero, or is the final zero declet */ |
| #define dpdlendun(n, form) {dpd=(form)&0x3ff; \ |
| if (dpd==0) return 1; \ |
| return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);} |
| |
| uInt decFloatDigits(const decFloat *df) { |
| uInt dpd; /* work */ |
| uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */ |
| #if QUAD |
| uInt sourmh, sourml; |
| #endif |
| uInt sourlo; |
| |
| if (DFISINF(df)) return 1; |
| /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */ |
| /* then the coefficient is full-length */ |
| if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX; |
| |
| #if DOUBLE |
| if (sourhi&0x0003ffff) { /* ends in first */ |
| dpdlenchk(0, sourhi>>8); |
| sourlo=DFWORD(df, 1); |
| dpdlendun(1, (sourhi<<2) | (sourlo>>30)); |
| } /* [cannot drop through] */ |
| sourlo=DFWORD(df, 1); /* sourhi not involved now */ |
| if (sourlo&0xfff00000) { /* in one of first two */ |
| dpdlenchk(1, sourlo>>30); /* very rare */ |
| dpdlendun(2, sourlo>>20); |
| } /* [cannot drop through] */ |
| dpdlenchk(3, sourlo>>10); |
| dpdlendun(4, sourlo); |
| /* [cannot drop through] */ |
| |
| #elif QUAD |
| if (sourhi&0x00003fff) { /* ends in first */ |
| dpdlenchk(0, sourhi>>4); |
| sourmh=DFWORD(df, 1); |
| dpdlendun(1, ((sourhi)<<6) | (sourmh>>26)); |
| } /* [cannot drop through] */ |
| sourmh=DFWORD(df, 1); |
| if (sourmh) { |
| dpdlenchk(1, sourmh>>26); |
| dpdlenchk(2, sourmh>>16); |
| dpdlenchk(3, sourmh>>6); |
| sourml=DFWORD(df, 2); |
| dpdlendun(4, ((sourmh)<<4) | (sourml>>28)); |
| } /* [cannot drop through] */ |
| sourml=DFWORD(df, 2); |
| if (sourml) { |
| dpdlenchk(4, sourml>>28); |
| dpdlenchk(5, sourml>>18); |
| dpdlenchk(6, sourml>>8); |
| sourlo=DFWORD(df, 3); |
| dpdlendun(7, ((sourml)<<2) | (sourlo>>30)); |
| } /* [cannot drop through] */ |
| sourlo=DFWORD(df, 3); |
| if (sourlo&0xfff00000) { /* in one of first two */ |
| dpdlenchk(7, sourlo>>30); /* very rare */ |
| dpdlendun(8, sourlo>>20); |
| } /* [cannot drop through] */ |
| dpdlenchk(9, sourlo>>10); |
| dpdlendun(10, sourlo); |
| /* [cannot drop through] */ |
| #endif |
| } /* decFloatDigits */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatDivide -- divide a decFloat by another */ |
| /* */ |
| /* result gets the result of dividing dfl by dfr: */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* set is the context */ |
| /* returns result */ |
| /* */ |
| /* ------------------------------------------------------------------ */ |
| /* This is just a wrapper. */ |
| decFloat * decFloatDivide(decFloat *result, |
| const decFloat *dfl, const decFloat *dfr, |
| decContext *set) { |
| return decDivide(result, dfl, dfr, set, DIVIDE); |
| } /* decFloatDivide */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatDivideInteger -- integer divide a decFloat by another */ |
| /* */ |
| /* result gets the result of dividing dfl by dfr: */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* set is the context */ |
| /* returns result */ |
| /* */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatDivideInteger(decFloat *result, |
| const decFloat *dfl, const decFloat *dfr, |
| decContext *set) { |
| return decDivide(result, dfl, dfr, set, DIVIDEINT); |
| } /* decFloatDivideInteger */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatFMA -- multiply and add three decFloats, fused */ |
| /* */ |
| /* result gets the result of (dfl*dfr)+dff with a single rounding */ |
| /* dfl is the first decFloat (lhs) */ |
| /* dfr is the second decFloat (rhs) */ |
| /* dff is the final decFloat (fhs) */ |
| /* set is the context */ |
| /* returns result */ |
| /* */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatFMA(decFloat *result, const decFloat *dfl, |
| const decFloat *dfr, const decFloat *dff, |
| decContext *set) { |
| |
| /* The accumulator has the bytes needed for FiniteMultiply, plus */ |
| /* one byte to the left in case of carry, plus DECPMAX+2 to the */ |
| /* right for the final addition (up to full fhs + round & sticky) */ |
| #define FMALEN (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 (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4; |
| for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++; |
| for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4; |
| for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++; |
| |
| /* if hi is zero then result will be lo (which has the smaller */ |
| /* exponent), which also may need to be tested for zero for the */ |
| /* weird IEEE 754 sign rules */ |
| if (*hi->msd==0) { /* hi 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; uh<hi->lsd-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 (uh<hi->msd || ul<lo->msd) break; |
| *ub=(uByte)(carry+(*uh--)+(*ul--)); |
| carry=0; |
| if (*ub<10) continue; |
| *ub-=10; |
| carry=1; |
| } /* both loop */ |
| |
| if (ul<lo->msd) { /* to left of lo */ |
| for (;; ub--) { |
| if (uh<hi->msd) break; |
| *ub=(uByte)(carry+(*uh--)); /* [+0] */ |
| carry=0; |
| if (*ub<10) continue; |
| *ub-=10; |
| carry=1; |
| } /* hi loop */ |
| } |
| else { /* to left of hi */ |
| for (;; ub--) { |
| if (ul<lo->msd) break; |
| *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 hi<lo */ |
| /* borrowed -- take ten's complement of the right digits */ |
| lo->sign=hi->sign; /* sign is lhs sign */ |
| for (ul=lo->msd; ul<lo->lsd-3; ul+=4) 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 (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4; |
| for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++; |
| if (*lo->msd==0) { /* must be true zero (and diffsign) */ |
| lo->sign=0; /* assume + */ |
| 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->msd<acc) { |
| printf("FMA underrun by %ld \n", (LI)(acc-lo->msd)); |
| } |
| #endif |
| |
| return decFinalize(result, lo, set); /* round, check, and lay out */ |
| } /* decFloatFMA */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatFromInt -- initialise a decFloat from an Int */ |
| /* */ |
| /* 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 decFloatIsNaN(const decFloat *df) { |
| return DFISNAN(df); |
| } |
| uInt decFloatIsNormal(const decFloat *df) { |
| Int exp; /* exponent */ |
| if (DFISSPECIAL(df)) return 0; |
| if (DFISZERO(df)) return 0; |
| /* is finite and non-zero */ |
| exp=GETEXPUN(df) /* get unbiased exponent .. */ |
| +decFloatDigits(df)-1; /* .. and make adjusted exponent */ |
| return (exp>=DECEMIN); /* < DECEMIN is subnormal */ |
| } |
| uInt decFloatIsSignaling(const decFloat *df) { |
| return DFISSNAN(df); |
| } |
| uInt decFloatIsSignalling(const decFloat *df) { |
| return DFISSNAN(df); |
| } |
| uInt decFloatIsSigned(const decFloat *df) { |
| return DFISSIGNED(df); |
| } |
| uInt decFloatIsSubnormal(const decFloat *df) { |
| if (DFISSPECIAL(df)) return 0; |
| /* is finite */ |
| if (decFloatIsNormal(df)) return 0; |
| /* it is <Nmin, but could be zero */ |
| if (DFISZERO(df)) return 0; |
| return 1; /* is subnormal */ |
| } |
| uInt decFloatIsZero(const decFloat *df) { |
| return DFISZERO(df); |
| } /* decFloatIs... */ |
| |
| /* ------------------------------------------------------------------ */ |
| /* decFloatLogB -- return adjusted exponent, by 754 rules */ |
| /* */ |
| /* result gets the adjusted exponent as an integer, or a NaN etc. */ |
| /* df is the decFloat to be examined */ |
| /* set is the context */ |
| /* returns result */ |
| /* */ |
| /* Notable cases: */ |
| /* A<0 -> Use |A| */ |
| /* A=0 -> -Infinity (Division by zero) */ |
| /* A=Infinite -> +Infinity (Exact) */ |
| /* A=1 exactly -> 0 (Exact) */ |
| /* NaNs are propagated as usual */ |
| /* ------------------------------------------------------------------ */ |
| decFloat * decFloatLogB(decFloat *result, const decFloat *df, |
| decContext *set) { |
| Int ae; /* adjusted exponent */ |
| if (DFISNAN(df)) return decNaNs(result, df, NULL, set); |
| if (DFISINF(df)) { |
| DFWORD(result, 0)=0; /* need +ve */ |
| return decInfinity(result, result); /* canonical +Infinity */ |
| } |
| if (DFISZERO(df)) { |
| set->status|=DEC_Division_by_zero; /* as per 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) { /* lhs<rhs, do NextPlus, see above for commentary */ |
| if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */ |
| DFSETNMAX(result); |
| DFWORD(result, 0)|=DECFLOAT_Sign; |
| return result; |
|