Google Git
Sign in
gnu / gcc / a8404c07e7fca388c02c39077865f7d5fa928430 / . / libdecnumber / decBasic.c
blob: 27f98263b4e3a53864360c29b23f278030d27ad0 [file] [log] [blame]
/* Common base code for the decNumber C Library.
Copyright (C) 2007-2022 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(&quotient, &num, set); /* lay out the integer so far */
DFWORD(&quotient, 0)^=DECFLOAT_Sign; /* negate it */
sign=DFWORD(dfl, 0); /* save sign of dfl */
decFloatFMA(result, &quotient, dfr, dfl, set);
if (!DFISZERO(result)) return result;
/* if the result is zero the sign shall be sign of dfl */
DFWORD(&quotient, 0)=sign; /* construct decFloat of sign */
return decFloatCopySign(result, result, &quotient);
} /* decDivide */
/* ------------------------------------------------------------------ */
/* decFiniteMultiply -- multiply two finite decFloats */
/* */
/* num gets the result of multiplying dfl and dfr */
/* bcdacc .. with the coefficient in this array */
/* dfl is the first decFloat (lhs) */
/* dfr is the second decFloat (rhs) */
/* */
/* This effects the multiplication of two decFloats, both known to be */
/* finite, leaving the result in a bcdnum ready for decFinalize (for */
/* use in Multiply) or in a following addition (FMA). */
/* */
/* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
/* No error is possible and no status is set. */
/* ------------------------------------------------------------------ */
/* This routine has two separate implementations of the core */
/* multiplication; both using base-billion. One uses only 32-bit */
/* variables (Ints and uInts) or smaller; the other uses uLongs (for */
/* multiplication and addition only). Both implementations cover */
/* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
/* comparisons. In any one compilation only one implementation for */
/* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
/* version is forced. */
/* */
/* Historical note: an earlier version of this code also supported the */
/* 256-bit format and has been preserved. That is somewhat trickier */
/* during lazy carry splitting because the initial quotient estimate */
/* (est) can exceed 32 bits. */
#define MULTBASE ((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;
}