163d1a8abSmrg /* Common base code for the decNumber C Library.
2*ec02198aSmrg Copyright (C) 2007-2020 Free Software Foundation, Inc.
363d1a8abSmrg Contributed by IBM Corporation. Author Mike Cowlishaw.
463d1a8abSmrg
563d1a8abSmrg This file is part of GCC.
663d1a8abSmrg
763d1a8abSmrg GCC is free software; you can redistribute it and/or modify it under
863d1a8abSmrg the terms of the GNU General Public License as published by the Free
963d1a8abSmrg Software Foundation; either version 3, or (at your option) any later
1063d1a8abSmrg version.
1163d1a8abSmrg
1263d1a8abSmrg GCC is distributed in the hope that it will be useful, but WITHOUT ANY
1363d1a8abSmrg WARRANTY; without even the implied warranty of MERCHANTABILITY or
1463d1a8abSmrg FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
1563d1a8abSmrg for more details.
1663d1a8abSmrg
1763d1a8abSmrg Under Section 7 of GPL version 3, you are granted additional
1863d1a8abSmrg permissions described in the GCC Runtime Library Exception, version
1963d1a8abSmrg 3.1, as published by the Free Software Foundation.
2063d1a8abSmrg
2163d1a8abSmrg You should have received a copy of the GNU General Public License and
2263d1a8abSmrg a copy of the GCC Runtime Library Exception along with this program;
2363d1a8abSmrg see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
2463d1a8abSmrg <http://www.gnu.org/licenses/>. */
2563d1a8abSmrg
2663d1a8abSmrg /* ------------------------------------------------------------------ */
2763d1a8abSmrg /* decBasic.c -- common base code for Basic decimal types */
2863d1a8abSmrg /* ------------------------------------------------------------------ */
2963d1a8abSmrg /* This module comprises code that is shared between decDouble and */
3063d1a8abSmrg /* decQuad (but not decSingle). The main arithmetic operations are */
3163d1a8abSmrg /* here (Add, Subtract, Multiply, FMA, and Division operators). */
3263d1a8abSmrg /* */
3363d1a8abSmrg /* Unlike decNumber, parameterization takes place at compile time */
3463d1a8abSmrg /* rather than at runtime. The parameters are set in the decDouble.c */
3563d1a8abSmrg /* (etc.) files, which then include this one to produce the compiled */
3663d1a8abSmrg /* code. The functions here, therefore, are code shared between */
3763d1a8abSmrg /* multiple formats. */
3863d1a8abSmrg /* */
3963d1a8abSmrg /* This must be included after decCommon.c. */
4063d1a8abSmrg /* ------------------------------------------------------------------ */
4163d1a8abSmrg /* Names here refer to decFloat rather than to decDouble, etc., and */
4263d1a8abSmrg /* the functions are in strict alphabetical order. */
4363d1a8abSmrg
4463d1a8abSmrg /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
4563d1a8abSmrg /* decCommon.c */
4663d1a8abSmrg #if !defined(QUAD)
4763d1a8abSmrg #error decBasic.c must be included after decCommon.c
4863d1a8abSmrg #endif
4963d1a8abSmrg #if SINGLE
5063d1a8abSmrg #error Routines in decBasic.c are for decDouble and decQuad only
5163d1a8abSmrg #endif
5263d1a8abSmrg
5363d1a8abSmrg /* Private constants */
5463d1a8abSmrg #define DIVIDE 0x80000000 /* Divide operations [as flags] */
5563d1a8abSmrg #define REMAINDER 0x40000000 /* .. */
5663d1a8abSmrg #define DIVIDEINT 0x20000000 /* .. */
5763d1a8abSmrg #define REMNEAR 0x10000000 /* .. */
5863d1a8abSmrg
5963d1a8abSmrg /* Private functions (local, used only by routines in this module) */
6063d1a8abSmrg static decFloat *decDivide(decFloat *, const decFloat *,
6163d1a8abSmrg const decFloat *, decContext *, uInt);
6263d1a8abSmrg static decFloat *decCanonical(decFloat *, const decFloat *);
6363d1a8abSmrg static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
6463d1a8abSmrg const decFloat *);
6563d1a8abSmrg static decFloat *decInfinity(decFloat *, const decFloat *);
6663d1a8abSmrg static decFloat *decInvalid(decFloat *, decContext *);
6763d1a8abSmrg static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
6863d1a8abSmrg decContext *);
6963d1a8abSmrg static Int decNumCompare(const decFloat *, const decFloat *, Flag);
7063d1a8abSmrg static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
7163d1a8abSmrg enum rounding, Flag);
7263d1a8abSmrg static uInt decToInt32(const decFloat *, decContext *, enum rounding,
7363d1a8abSmrg Flag, Flag);
7463d1a8abSmrg
7563d1a8abSmrg /* ------------------------------------------------------------------ */
7663d1a8abSmrg /* decCanonical -- copy a decFloat, making canonical */
7763d1a8abSmrg /* */
7863d1a8abSmrg /* result gets the canonicalized df */
7963d1a8abSmrg /* df is the decFloat to copy and make canonical */
8063d1a8abSmrg /* returns result */
8163d1a8abSmrg /* */
8263d1a8abSmrg /* This is exposed via decFloatCanonical for Double and Quad only. */
8363d1a8abSmrg /* This works on specials, too; no error or exception is possible. */
8463d1a8abSmrg /* ------------------------------------------------------------------ */
decCanonical(decFloat * result,const decFloat * df)8563d1a8abSmrg static decFloat * decCanonical(decFloat *result, const decFloat *df) {
8663d1a8abSmrg uInt encode, precode, dpd; /* work */
8763d1a8abSmrg uInt inword, uoff, canon; /* .. */
8863d1a8abSmrg Int n; /* counter (down) */
8963d1a8abSmrg if (df!=result) *result=*df; /* effect copy if needed */
9063d1a8abSmrg if (DFISSPECIAL(result)) {
9163d1a8abSmrg if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
9263d1a8abSmrg /* is a NaN */
9363d1a8abSmrg DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */
9463d1a8abSmrg if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */
9563d1a8abSmrg /* drop through to check payload */
9663d1a8abSmrg }
9763d1a8abSmrg /* return quickly if the coefficient continuation is canonical */
9863d1a8abSmrg { /* declare block */
9963d1a8abSmrg #if DOUBLE
10063d1a8abSmrg uInt sourhi=DFWORD(df, 0);
10163d1a8abSmrg uInt sourlo=DFWORD(df, 1);
10263d1a8abSmrg if (CANONDPDOFF(sourhi, 8)
10363d1a8abSmrg && CANONDPDTWO(sourhi, sourlo, 30)
10463d1a8abSmrg && CANONDPDOFF(sourlo, 20)
10563d1a8abSmrg && CANONDPDOFF(sourlo, 10)
10663d1a8abSmrg && CANONDPDOFF(sourlo, 0)) return result;
10763d1a8abSmrg #elif QUAD
10863d1a8abSmrg uInt sourhi=DFWORD(df, 0);
10963d1a8abSmrg uInt sourmh=DFWORD(df, 1);
11063d1a8abSmrg uInt sourml=DFWORD(df, 2);
11163d1a8abSmrg uInt sourlo=DFWORD(df, 3);
11263d1a8abSmrg if (CANONDPDOFF(sourhi, 4)
11363d1a8abSmrg && CANONDPDTWO(sourhi, sourmh, 26)
11463d1a8abSmrg && CANONDPDOFF(sourmh, 16)
11563d1a8abSmrg && CANONDPDOFF(sourmh, 6)
11663d1a8abSmrg && CANONDPDTWO(sourmh, sourml, 28)
11763d1a8abSmrg && CANONDPDOFF(sourml, 18)
11863d1a8abSmrg && CANONDPDOFF(sourml, 8)
11963d1a8abSmrg && CANONDPDTWO(sourml, sourlo, 30)
12063d1a8abSmrg && CANONDPDOFF(sourlo, 20)
12163d1a8abSmrg && CANONDPDOFF(sourlo, 10)
12263d1a8abSmrg && CANONDPDOFF(sourlo, 0)) return result;
12363d1a8abSmrg #endif
12463d1a8abSmrg } /* block */
12563d1a8abSmrg
12663d1a8abSmrg /* Loop to repair a non-canonical coefficent, as needed */
12763d1a8abSmrg inword=DECWORDS-1; /* current input word */
12863d1a8abSmrg uoff=0; /* bit offset of declet */
12963d1a8abSmrg encode=DFWORD(result, inword);
13063d1a8abSmrg for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */
13163d1a8abSmrg dpd=encode>>uoff;
13263d1a8abSmrg uoff+=10;
13363d1a8abSmrg if (uoff>32) { /* crossed uInt boundary */
13463d1a8abSmrg inword--;
13563d1a8abSmrg encode=DFWORD(result, inword);
13663d1a8abSmrg uoff-=32;
13763d1a8abSmrg dpd|=encode<<(10-uoff); /* get pending bits */
13863d1a8abSmrg }
13963d1a8abSmrg dpd&=0x3ff; /* clear uninteresting bits */
14063d1a8abSmrg if (dpd<0x16e) continue; /* must be canonical */
14163d1a8abSmrg canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */
14263d1a8abSmrg if (canon==dpd) continue; /* have canonical declet */
14363d1a8abSmrg /* need to replace declet */
14463d1a8abSmrg if (uoff>=10) { /* all within current word */
14563d1a8abSmrg encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
14663d1a8abSmrg encode|=canon<<(uoff-10); /* insert the canonical form */
14763d1a8abSmrg DFWORD(result, inword)=encode; /* .. and save */
14863d1a8abSmrg continue;
14963d1a8abSmrg }
15063d1a8abSmrg /* straddled words */
15163d1a8abSmrg precode=DFWORD(result, inword+1); /* get previous */
15263d1a8abSmrg precode&=0xffffffff>>(10-uoff); /* clear top bits */
15363d1a8abSmrg DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
15463d1a8abSmrg encode&=0xffffffff<<uoff; /* clear bottom bits */
15563d1a8abSmrg encode|=canon>>(10-uoff); /* insert canonical */
15663d1a8abSmrg DFWORD(result, inword)=encode; /* .. and save */
15763d1a8abSmrg } /* n */
15863d1a8abSmrg return result;
15963d1a8abSmrg } /* decCanonical */
16063d1a8abSmrg
16163d1a8abSmrg /* ------------------------------------------------------------------ */
16263d1a8abSmrg /* decDivide -- divide operations */
16363d1a8abSmrg /* */
16463d1a8abSmrg /* result gets the result of dividing dfl by dfr: */
16563d1a8abSmrg /* dfl is the first decFloat (lhs) */
16663d1a8abSmrg /* dfr is the second decFloat (rhs) */
16763d1a8abSmrg /* set is the context */
16863d1a8abSmrg /* op is the operation selector */
16963d1a8abSmrg /* returns result */
17063d1a8abSmrg /* */
17163d1a8abSmrg /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
17263d1a8abSmrg /* ------------------------------------------------------------------ */
17363d1a8abSmrg #define DIVCOUNT 0 /* 1 to instrument subtractions counter */
17463d1a8abSmrg #define DIVBASE ((uInt)BILLION) /* the base used for divide */
17563d1a8abSmrg #define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
17663d1a8abSmrg #define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */
decDivide(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set,uInt op)17763d1a8abSmrg static decFloat * decDivide(decFloat *result, const decFloat *dfl,
17863d1a8abSmrg const decFloat *dfr, decContext *set, uInt op) {
17963d1a8abSmrg decFloat quotient; /* for remainders */
18063d1a8abSmrg bcdnum num; /* for final conversion */
18163d1a8abSmrg uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
18263d1a8abSmrg uInt div[DIVOPLEN]; /* divisor in base-billion .. */
18363d1a8abSmrg uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */
18463d1a8abSmrg uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
18563d1a8abSmrg uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */
18663d1a8abSmrg Int divunits, accunits; /* lengths */
18763d1a8abSmrg Int quodigits; /* digits in quotient */
18863d1a8abSmrg uInt *lsua, *lsuq; /* -> current acc and quo lsus */
18963d1a8abSmrg Int length, multiplier; /* work */
19063d1a8abSmrg uInt carry, sign; /* .. */
19163d1a8abSmrg uInt *ua, *ud, *uq; /* .. */
19263d1a8abSmrg uByte *ub; /* .. */
19363d1a8abSmrg uInt uiwork; /* for macros */
19463d1a8abSmrg uInt divtop; /* top unit of div adjusted for estimating */
19563d1a8abSmrg #if DIVCOUNT
19663d1a8abSmrg static uInt maxcount=0; /* worst-seen subtractions count */
19763d1a8abSmrg uInt divcount=0; /* subtractions count [this divide] */
19863d1a8abSmrg #endif
19963d1a8abSmrg
20063d1a8abSmrg /* calculate sign */
20163d1a8abSmrg num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
20263d1a8abSmrg
20363d1a8abSmrg if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
20463d1a8abSmrg /* NaNs are handled as usual */
20563d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
20663d1a8abSmrg /* one or two infinities */
20763d1a8abSmrg if (DFISINF(dfl)) {
20863d1a8abSmrg if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
20963d1a8abSmrg if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
21063d1a8abSmrg /* Infinity/x is infinite and quiet, even if x=0 */
21163d1a8abSmrg DFWORD(result, 0)=num.sign;
21263d1a8abSmrg return decInfinity(result, result);
21363d1a8abSmrg }
21463d1a8abSmrg /* must be x/Infinity -- remainders are lhs */
21563d1a8abSmrg if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
21663d1a8abSmrg /* divides: return zero with correct sign and exponent depending */
21763d1a8abSmrg /* on op (Etiny for divide, 0 for divideInt) */
21863d1a8abSmrg decFloatZero(result);
21963d1a8abSmrg if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
22063d1a8abSmrg else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */
22163d1a8abSmrg return result;
22263d1a8abSmrg }
22363d1a8abSmrg /* next, handle zero operands (x/0 and 0/x) */
22463d1a8abSmrg if (DFISZERO(dfr)) { /* x/0 */
22563d1a8abSmrg if (DFISZERO(dfl)) { /* 0/0 is undefined */
22663d1a8abSmrg decFloatZero(result);
22763d1a8abSmrg DFWORD(result, 0)=DECFLOAT_qNaN;
22863d1a8abSmrg set->status|=DEC_Division_undefined;
22963d1a8abSmrg return result;
23063d1a8abSmrg }
23163d1a8abSmrg if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
23263d1a8abSmrg set->status|=DEC_Division_by_zero;
23363d1a8abSmrg DFWORD(result, 0)=num.sign;
23463d1a8abSmrg return decInfinity(result, result); /* x/0 -> signed Infinity */
23563d1a8abSmrg }
23663d1a8abSmrg num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
23763d1a8abSmrg if (DFISZERO(dfl)) { /* 0/x (x!=0) */
23863d1a8abSmrg /* if divide, result is 0 with ideal exponent; divideInt has */
23963d1a8abSmrg /* exponent=0, remainders give zero with lower exponent */
24063d1a8abSmrg if (op&DIVIDEINT) {
24163d1a8abSmrg decFloatZero(result);
24263d1a8abSmrg DFWORD(result, 0)|=num.sign; /* add sign */
24363d1a8abSmrg return result;
24463d1a8abSmrg }
24563d1a8abSmrg if (!(op&DIVIDE)) { /* a remainder */
24663d1a8abSmrg /* exponent is the minimum of the operands */
24763d1a8abSmrg num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
24863d1a8abSmrg /* if the result is zero the sign shall be sign of dfl */
24963d1a8abSmrg num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
25063d1a8abSmrg }
25163d1a8abSmrg bcdacc[0]=0;
25263d1a8abSmrg num.msd=bcdacc; /* -> 0 */
25363d1a8abSmrg num.lsd=bcdacc; /* .. */
25463d1a8abSmrg return decFinalize(result, &num, set); /* [divide may clamp exponent] */
25563d1a8abSmrg } /* 0/x */
25663d1a8abSmrg /* [here, both operands are known to be finite and non-zero] */
25763d1a8abSmrg
25863d1a8abSmrg /* extract the operand coefficents into 'units' which are */
25963d1a8abSmrg /* base-billion; the lhs is high-aligned in acc and the msu of both */
26063d1a8abSmrg /* acc and div is at the right-hand end of array (offset length-1); */
26163d1a8abSmrg /* the quotient can need one more unit than the operands as digits */
26263d1a8abSmrg /* in it are not necessarily aligned neatly; further, the quotient */
26363d1a8abSmrg /* may not start accumulating until after the end of the initial */
26463d1a8abSmrg /* operand in acc if that is small (e.g., 1) so the accumulator */
26563d1a8abSmrg /* must have at least that number of units extra (at the ls end) */
26663d1a8abSmrg GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
26763d1a8abSmrg GETCOEFFBILL(dfr, div);
26863d1a8abSmrg /* zero the low uInts of acc */
26963d1a8abSmrg acc[0]=0;
27063d1a8abSmrg acc[1]=0;
27163d1a8abSmrg acc[2]=0;
27263d1a8abSmrg acc[3]=0;
27363d1a8abSmrg #if DOUBLE
27463d1a8abSmrg #if DIVOPLEN!=2
27563d1a8abSmrg #error Unexpected Double DIVOPLEN
27663d1a8abSmrg #endif
27763d1a8abSmrg #elif QUAD
27863d1a8abSmrg acc[4]=0;
27963d1a8abSmrg acc[5]=0;
28063d1a8abSmrg acc[6]=0;
28163d1a8abSmrg acc[7]=0;
28263d1a8abSmrg #if DIVOPLEN!=4
28363d1a8abSmrg #error Unexpected Quad DIVOPLEN
28463d1a8abSmrg #endif
28563d1a8abSmrg #endif
28663d1a8abSmrg
28763d1a8abSmrg /* set msu and lsu pointers */
28863d1a8abSmrg msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
28963d1a8abSmrg msuq=quo+DIVOPLEN;
29063d1a8abSmrg /*[loop for div will terminate because operands are non-zero] */
29163d1a8abSmrg for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
29263d1a8abSmrg /* the initial least-significant unit of acc is set so acc appears */
29363d1a8abSmrg /* to have the same length as div. */
29463d1a8abSmrg /* This moves one position towards the least possible for each */
29563d1a8abSmrg /* iteration */
29663d1a8abSmrg divunits=(Int)(msud-div+1); /* precalculate */
29763d1a8abSmrg lsua=msua-divunits+1; /* initial working lsu of acc */
29863d1a8abSmrg lsuq=msuq; /* and of quo */
29963d1a8abSmrg
30063d1a8abSmrg /* set up the estimator for the multiplier; this is the msu of div, */
30163d1a8abSmrg /* plus two bits from the unit below (if any) rounded up by one if */
30263d1a8abSmrg /* there are any non-zero bits or units below that [the extra two */
30363d1a8abSmrg /* bits makes for a much better estimate when the top unit is small] */
30463d1a8abSmrg divtop=*msud<<2;
30563d1a8abSmrg if (divunits>1) {
30663d1a8abSmrg uInt *um=msud-1;
30763d1a8abSmrg uInt d=*um;
30863d1a8abSmrg if (d>=750000000) {divtop+=3; d-=750000000;}
30963d1a8abSmrg else if (d>=500000000) {divtop+=2; d-=500000000;}
31063d1a8abSmrg else if (d>=250000000) {divtop++; d-=250000000;}
31163d1a8abSmrg if (d) divtop++;
31263d1a8abSmrg else for (um--; um>=div; um--) if (*um) {
31363d1a8abSmrg divtop++;
31463d1a8abSmrg break;
31563d1a8abSmrg }
31663d1a8abSmrg } /* >1 unit */
31763d1a8abSmrg
31863d1a8abSmrg #if DECTRACE
31963d1a8abSmrg {Int i;
32063d1a8abSmrg printf("----- div=");
32163d1a8abSmrg for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
32263d1a8abSmrg printf("\n");}
32363d1a8abSmrg #endif
32463d1a8abSmrg
32563d1a8abSmrg /* now collect up to DECPMAX+1 digits in the quotient (this may */
32663d1a8abSmrg /* need OPLEN+1 uInts if unaligned) */
32763d1a8abSmrg quodigits=0; /* no digits yet */
32863d1a8abSmrg for (;; lsua--) { /* outer loop -- each input position */
32963d1a8abSmrg #if DECCHECK
33063d1a8abSmrg if (lsua<acc) {
33163d1a8abSmrg printf("Acc underrun...\n");
33263d1a8abSmrg break;
33363d1a8abSmrg }
33463d1a8abSmrg #endif
33563d1a8abSmrg #if DECTRACE
33663d1a8abSmrg printf("Outer: quodigits=%ld acc=", (LI)quodigits);
33763d1a8abSmrg for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
33863d1a8abSmrg printf("\n");
33963d1a8abSmrg #endif
34063d1a8abSmrg *lsuq=0; /* default unit result is 0 */
34163d1a8abSmrg for (;;) { /* inner loop -- calculate quotient unit */
34263d1a8abSmrg /* strip leading zero units from acc (either there initially or */
34363d1a8abSmrg /* from subtraction below); this may strip all if exactly 0 */
34463d1a8abSmrg for (; *msua==0 && msua>=lsua;) msua--;
34563d1a8abSmrg accunits=(Int)(msua-lsua+1); /* [maybe 0] */
34663d1a8abSmrg /* subtraction is only necessary and possible if there are as */
34763d1a8abSmrg /* least as many units remaining in acc for this iteration as */
34863d1a8abSmrg /* there are in div */
34963d1a8abSmrg if (accunits<divunits) {
35063d1a8abSmrg if (accunits==0) msua++; /* restore */
35163d1a8abSmrg break;
35263d1a8abSmrg }
35363d1a8abSmrg
35463d1a8abSmrg /* If acc is longer than div then subtraction is definitely */
35563d1a8abSmrg /* possible (as msu of both is non-zero), but if they are the */
35663d1a8abSmrg /* same length a comparison is needed. */
35763d1a8abSmrg /* If a subtraction is needed then a good estimate of the */
35863d1a8abSmrg /* multiplier for the subtraction is also needed in order to */
35963d1a8abSmrg /* minimise the iterations of this inner loop because the */
36063d1a8abSmrg /* subtractions needed dominate division performance. */
36163d1a8abSmrg if (accunits==divunits) {
36263d1a8abSmrg /* compare the high divunits of acc and div: */
36363d1a8abSmrg /* acc<div: this quotient unit is unchanged; subtraction */
36463d1a8abSmrg /* will be possible on the next iteration */
36563d1a8abSmrg /* acc==div: quotient gains 1, set acc=0 */
36663d1a8abSmrg /* acc>div: subtraction necessary at this position */
36763d1a8abSmrg for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
36863d1a8abSmrg /* [now at first mismatch or lsu] */
36963d1a8abSmrg if (*ud>*ua) break; /* next time... */
37063d1a8abSmrg if (*ud==*ua) { /* all compared equal */
37163d1a8abSmrg *lsuq+=1; /* increment result */
37263d1a8abSmrg msua=lsua; /* collapse acc units */
37363d1a8abSmrg *msua=0; /* .. to a zero */
37463d1a8abSmrg break;
37563d1a8abSmrg }
37663d1a8abSmrg
37763d1a8abSmrg /* subtraction necessary; estimate multiplier [see above] */
37863d1a8abSmrg /* if both *msud and *msua are small it is cost-effective to */
37963d1a8abSmrg /* bring in part of the following units (if any) to get a */
38063d1a8abSmrg /* better estimate (assume some other non-zero in div) */
38163d1a8abSmrg #define DIVLO 1000000U
38263d1a8abSmrg #define DIVHI (DIVBASE/DIVLO)
38363d1a8abSmrg #if DECUSE64
38463d1a8abSmrg if (divunits>1) {
38563d1a8abSmrg /* there cannot be a *(msud-2) for DECDOUBLE so next is */
38663d1a8abSmrg /* an exact calculation unless DECQUAD (which needs to */
38763d1a8abSmrg /* assume bits out there if divunits>2) */
38863d1a8abSmrg uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
38963d1a8abSmrg uLong div=(uLong)*msud * DIVBASE + *(msud-1);
39063d1a8abSmrg #if QUAD
39163d1a8abSmrg if (divunits>2) div++;
39263d1a8abSmrg #endif
39363d1a8abSmrg mul/=div;
39463d1a8abSmrg multiplier=(Int)mul;
39563d1a8abSmrg }
39663d1a8abSmrg else multiplier=*msua/(*msud);
39763d1a8abSmrg #else
39863d1a8abSmrg if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
39963d1a8abSmrg multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
40063d1a8abSmrg /(*msud*DIVHI + *(msud-1)/DIVLO +1);
40163d1a8abSmrg }
40263d1a8abSmrg else multiplier=(*msua<<2)/divtop;
40363d1a8abSmrg #endif
40463d1a8abSmrg }
40563d1a8abSmrg else { /* accunits>divunits */
40663d1a8abSmrg /* msud is one unit 'lower' than msua, so estimate differently */
40763d1a8abSmrg #if DECUSE64
40863d1a8abSmrg uLong mul;
40963d1a8abSmrg /* as before, bring in extra digits if possible */
41063d1a8abSmrg if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
41163d1a8abSmrg mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
41263d1a8abSmrg + *(msua-2)/DIVLO;
41363d1a8abSmrg mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
41463d1a8abSmrg }
41563d1a8abSmrg else if (divunits==1) {
41663d1a8abSmrg mul=(uLong)*msua * DIVBASE + *(msua-1);
41763d1a8abSmrg mul/=*msud; /* no more to the right */
41863d1a8abSmrg }
41963d1a8abSmrg else {
42063d1a8abSmrg mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
42163d1a8abSmrg + (*(msua-1)<<2);
42263d1a8abSmrg mul/=divtop; /* [divtop already allows for sticky bits] */
42363d1a8abSmrg }
42463d1a8abSmrg multiplier=(Int)mul;
42563d1a8abSmrg #else
42663d1a8abSmrg multiplier=*msua * ((DIVBASE<<2)/divtop);
42763d1a8abSmrg #endif
42863d1a8abSmrg }
42963d1a8abSmrg if (multiplier==0) multiplier=1; /* marginal case */
43063d1a8abSmrg *lsuq+=multiplier;
43163d1a8abSmrg
43263d1a8abSmrg #if DIVCOUNT
43363d1a8abSmrg /* printf("Multiplier: %ld\n", (LI)multiplier); */
43463d1a8abSmrg divcount++;
43563d1a8abSmrg #endif
43663d1a8abSmrg
43763d1a8abSmrg /* Carry out the subtraction acc-(div*multiplier); for each */
43863d1a8abSmrg /* unit in div, do the multiply, split to units (see */
43963d1a8abSmrg /* decFloatMultiply for the algorithm), and subtract from acc */
44063d1a8abSmrg #define DIVMAGIC 2305843009U /* 2**61/10**9 */
44163d1a8abSmrg #define DIVSHIFTA 29
44263d1a8abSmrg #define DIVSHIFTB 32
44363d1a8abSmrg carry=0;
44463d1a8abSmrg for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
44563d1a8abSmrg uInt lo, hop;
44663d1a8abSmrg #if DECUSE64
44763d1a8abSmrg uLong sub=(uLong)multiplier*(*ud)+carry;
44863d1a8abSmrg if (sub<DIVBASE) {
44963d1a8abSmrg carry=0;
45063d1a8abSmrg lo=(uInt)sub;
45163d1a8abSmrg }
45263d1a8abSmrg else {
45363d1a8abSmrg hop=(uInt)(sub>>DIVSHIFTA);
45463d1a8abSmrg carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
45563d1a8abSmrg /* the estimate is now in hi; now calculate sub-hi*10**9 */
45663d1a8abSmrg /* to get the remainder (which will be <DIVBASE)) */
45763d1a8abSmrg lo=(uInt)sub;
45863d1a8abSmrg lo-=carry*DIVBASE; /* low word of result */
45963d1a8abSmrg if (lo>=DIVBASE) {
46063d1a8abSmrg lo-=DIVBASE; /* correct by +1 */
46163d1a8abSmrg carry++;
46263d1a8abSmrg }
46363d1a8abSmrg }
46463d1a8abSmrg #else /* 32-bit */
46563d1a8abSmrg uInt hi;
46663d1a8abSmrg /* calculate multiplier*(*ud) into hi and lo */
46763d1a8abSmrg LONGMUL32HI(hi, *ud, multiplier); /* get the high word */
46863d1a8abSmrg lo=multiplier*(*ud); /* .. and the low */
46963d1a8abSmrg lo+=carry; /* add the old hi */
47063d1a8abSmrg carry=hi+(lo<carry); /* .. with any carry */
47163d1a8abSmrg if (carry || lo>=DIVBASE) { /* split is needed */
47263d1a8abSmrg hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */
47363d1a8abSmrg LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */
47463d1a8abSmrg /* [DIVSHIFTB is 32, so carry can be used directly] */
47563d1a8abSmrg /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
47663d1a8abSmrg /* happily the top word of the result is irrelevant because it */
47763d1a8abSmrg /* will always be zero so this needs only one multiplication */
47863d1a8abSmrg lo-=(carry*DIVBASE);
47963d1a8abSmrg /* the correction here will be at most +1; do it */
48063d1a8abSmrg if (lo>=DIVBASE) {
48163d1a8abSmrg lo-=DIVBASE;
48263d1a8abSmrg carry++;
48363d1a8abSmrg }
48463d1a8abSmrg }
48563d1a8abSmrg #endif
48663d1a8abSmrg if (lo>*ua) { /* borrow needed */
48763d1a8abSmrg *ua+=DIVBASE;
48863d1a8abSmrg carry++;
48963d1a8abSmrg }
49063d1a8abSmrg *ua-=lo;
49163d1a8abSmrg } /* ud loop */
49263d1a8abSmrg if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */
49363d1a8abSmrg } /* inner loop */
49463d1a8abSmrg
49563d1a8abSmrg /* the outer loop terminates when there is either an exact result */
49663d1a8abSmrg /* or enough digits; first update the quotient digit count and */
49763d1a8abSmrg /* pointer (if any significant digits) */
49863d1a8abSmrg #if DECTRACE
49963d1a8abSmrg if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
50063d1a8abSmrg #endif
50163d1a8abSmrg if (quodigits) {
50263d1a8abSmrg quodigits+=9; /* had leading unit earlier */
50363d1a8abSmrg lsuq--;
50463d1a8abSmrg if (quodigits>DECPMAX+1) break; /* have enough */
50563d1a8abSmrg }
50663d1a8abSmrg else if (*lsuq) { /* first quotient digits */
50763d1a8abSmrg const uInt *pow;
50863d1a8abSmrg for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
50963d1a8abSmrg lsuq--;
51063d1a8abSmrg /* [cannot have >DECPMAX+1 on first unit] */
51163d1a8abSmrg }
51263d1a8abSmrg
51363d1a8abSmrg if (*msua!=0) continue; /* not an exact result */
51463d1a8abSmrg /* acc is zero iff used all of original units and zero down to lsua */
51563d1a8abSmrg /* (must also continue to original lsu for correct quotient length) */
51663d1a8abSmrg if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
51763d1a8abSmrg for (; msua>lsua && *msua==0;) msua--;
51863d1a8abSmrg if (*msua==0 && msua==lsua) break;
51963d1a8abSmrg } /* outer loop */
52063d1a8abSmrg
52163d1a8abSmrg /* all of the original operand in acc has been covered at this point */
52263d1a8abSmrg /* quotient now has at least DECPMAX+2 digits */
52363d1a8abSmrg /* *msua is now non-0 if inexact and sticky bits */
52463d1a8abSmrg /* lsuq is one below the last uint of the quotient */
52563d1a8abSmrg lsuq++; /* set -> true lsu of quo */
52663d1a8abSmrg if (*msua) *lsuq|=1; /* apply sticky bit */
52763d1a8abSmrg
52863d1a8abSmrg /* quo now holds the (unrounded) quotient in base-billion; one */
52963d1a8abSmrg /* base-billion 'digit' per uInt. */
53063d1a8abSmrg #if DECTRACE
53163d1a8abSmrg printf("DivQuo:");
53263d1a8abSmrg for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
53363d1a8abSmrg printf("\n");
53463d1a8abSmrg #endif
53563d1a8abSmrg
53663d1a8abSmrg /* Now convert to BCD for rounding and cleanup, starting from the */
53763d1a8abSmrg /* most significant end [offset by one into bcdacc to leave room */
53863d1a8abSmrg /* for a possible carry digit if rounding for REMNEAR is needed] */
53963d1a8abSmrg for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
54063d1a8abSmrg uInt top, mid, rem; /* work */
54163d1a8abSmrg if (*uq==0) { /* no split needed */
54263d1a8abSmrg UBFROMUI(ub, 0); /* clear 9 BCD8s */
54363d1a8abSmrg UBFROMUI(ub+4, 0); /* .. */
54463d1a8abSmrg *(ub+8)=0; /* .. */
54563d1a8abSmrg continue;
54663d1a8abSmrg }
54763d1a8abSmrg /* *uq is non-zero -- split the base-billion digit into */
54863d1a8abSmrg /* hi, mid, and low three-digits */
54963d1a8abSmrg #define divsplit9 1000000 /* divisor */
55063d1a8abSmrg #define divsplit6 1000 /* divisor */
55163d1a8abSmrg /* The splitting is done by simple divides and remainders, */
55263d1a8abSmrg /* assuming the compiler will optimize these [GCC does] */
55363d1a8abSmrg top=*uq/divsplit9;
55463d1a8abSmrg rem=*uq%divsplit9;
55563d1a8abSmrg mid=rem/divsplit6;
55663d1a8abSmrg rem=rem%divsplit6;
55763d1a8abSmrg /* lay out the nine BCD digits (plus one unwanted byte) */
55863d1a8abSmrg UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
55963d1a8abSmrg UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
56063d1a8abSmrg UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
56163d1a8abSmrg } /* BCD conversion loop */
56263d1a8abSmrg ub--; /* -> lsu */
56363d1a8abSmrg
56463d1a8abSmrg /* complete the bcdnum; quodigits is correct, so the position of */
56563d1a8abSmrg /* the first non-zero is known */
56663d1a8abSmrg num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
56763d1a8abSmrg num.lsd=ub;
56863d1a8abSmrg
56963d1a8abSmrg /* make exponent adjustments, etc */
57063d1a8abSmrg if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */
57163d1a8abSmrg num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
57263d1a8abSmrg /* if the result was exact then there may be up to 8 extra */
57363d1a8abSmrg /* trailing zeros in the overflowed quotient final unit */
57463d1a8abSmrg if (*msua==0) {
57563d1a8abSmrg for (; *ub==0;) ub--; /* drop zeros */
57663d1a8abSmrg num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */
57763d1a8abSmrg num.lsd=ub;
57863d1a8abSmrg }
57963d1a8abSmrg } /* adjustment needed */
58063d1a8abSmrg
58163d1a8abSmrg #if DIVCOUNT
58263d1a8abSmrg if (divcount>maxcount) { /* new high-water nark */
58363d1a8abSmrg maxcount=divcount;
58463d1a8abSmrg printf("DivNewMaxCount: %ld\n", (LI)maxcount);
58563d1a8abSmrg }
58663d1a8abSmrg #endif
58763d1a8abSmrg
58863d1a8abSmrg if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
58963d1a8abSmrg
59063d1a8abSmrg /* Is DIVIDEINT or a remainder; there is more to do -- first form */
59163d1a8abSmrg /* the integer (this is done 'after the fact', unlike as in */
59263d1a8abSmrg /* decNumber, so as not to tax DIVIDE) */
59363d1a8abSmrg
59463d1a8abSmrg /* The first non-zero digit will be in the first 9 digits, known */
59563d1a8abSmrg /* from quodigits and num.msd, so there is always space for DECPMAX */
59663d1a8abSmrg /* digits */
59763d1a8abSmrg
59863d1a8abSmrg length=(Int)(num.lsd-num.msd+1);
59963d1a8abSmrg /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
60063d1a8abSmrg
60163d1a8abSmrg if (length+num.exponent>DECPMAX) { /* cannot fit */
60263d1a8abSmrg decFloatZero(result);
60363d1a8abSmrg DFWORD(result, 0)=DECFLOAT_qNaN;
60463d1a8abSmrg set->status|=DEC_Division_impossible;
60563d1a8abSmrg return result;
60663d1a8abSmrg }
60763d1a8abSmrg
60863d1a8abSmrg if (num.exponent>=0) { /* already an int, or need pad zeros */
60963d1a8abSmrg for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
61063d1a8abSmrg num.lsd+=num.exponent;
61163d1a8abSmrg }
61263d1a8abSmrg else { /* too long: round or truncate needed */
61363d1a8abSmrg Int drop=-num.exponent;
61463d1a8abSmrg if (!(op&REMNEAR)) { /* simple truncate */
61563d1a8abSmrg num.lsd-=drop;
61663d1a8abSmrg if (num.lsd<num.msd) { /* truncated all */
61763d1a8abSmrg num.lsd=num.msd; /* make 0 */
61863d1a8abSmrg *num.lsd=0; /* .. [sign still relevant] */
61963d1a8abSmrg }
62063d1a8abSmrg }
62163d1a8abSmrg else { /* round to nearest even [sigh] */
62263d1a8abSmrg /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
62363d1a8abSmrg /* (this is a special case of Quantize -- q.v. for commentary) */
62463d1a8abSmrg uByte *roundat; /* -> re-round digit */
62563d1a8abSmrg uByte reround; /* reround value */
62663d1a8abSmrg *(num.msd-1)=0; /* in case of left carry, or make 0 */
62763d1a8abSmrg if (drop<length) roundat=num.lsd-drop+1;
62863d1a8abSmrg else if (drop==length) roundat=num.msd;
62963d1a8abSmrg else roundat=num.msd-1; /* [-> 0] */
63063d1a8abSmrg reround=*roundat;
63163d1a8abSmrg for (ub=roundat+1; ub<=num.lsd; ub++) {
63263d1a8abSmrg if (*ub!=0) {
63363d1a8abSmrg reround=DECSTICKYTAB[reround];
63463d1a8abSmrg break;
63563d1a8abSmrg }
63663d1a8abSmrg } /* check stickies */
63763d1a8abSmrg if (roundat>num.msd) num.lsd=roundat-1;
63863d1a8abSmrg else {
63963d1a8abSmrg num.msd--; /* use the 0 .. */
64063d1a8abSmrg num.lsd=num.msd; /* .. at the new MSD place */
64163d1a8abSmrg }
64263d1a8abSmrg if (reround!=0) { /* discarding non-zero */
64363d1a8abSmrg uInt bump=0;
64463d1a8abSmrg /* rounding is DEC_ROUND_HALF_EVEN always */
64563d1a8abSmrg if (reround>5) bump=1; /* >0.5 goes up */
64663d1a8abSmrg else if (reround==5) /* exactly 0.5000 .. */
64763d1a8abSmrg bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */
64863d1a8abSmrg if (bump!=0) { /* need increment */
64963d1a8abSmrg /* increment the coefficient; this might end up with 1000... */
65063d1a8abSmrg ub=num.lsd;
65163d1a8abSmrg for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
65263d1a8abSmrg for (; *ub==9; ub--) *ub=0; /* at most 3 more */
65363d1a8abSmrg *ub+=1;
65463d1a8abSmrg if (ub<num.msd) num.msd--; /* carried */
65563d1a8abSmrg } /* bump needed */
65663d1a8abSmrg } /* reround!=0 */
65763d1a8abSmrg } /* remnear */
65863d1a8abSmrg } /* round or truncate needed */
65963d1a8abSmrg num.exponent=0; /* all paths */
66063d1a8abSmrg /*decShowNum(&num, "int"); */
66163d1a8abSmrg
66263d1a8abSmrg if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
66363d1a8abSmrg
66463d1a8abSmrg /* Have a remainder to calculate */
66563d1a8abSmrg decFinalize("ient, &num, set); /* lay out the integer so far */
66663d1a8abSmrg DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */
66763d1a8abSmrg sign=DFWORD(dfl, 0); /* save sign of dfl */
66863d1a8abSmrg decFloatFMA(result, "ient, dfr, dfl, set);
66963d1a8abSmrg if (!DFISZERO(result)) return result;
67063d1a8abSmrg /* if the result is zero the sign shall be sign of dfl */
67163d1a8abSmrg DFWORD("ient, 0)=sign; /* construct decFloat of sign */
67263d1a8abSmrg return decFloatCopySign(result, result, "ient);
67363d1a8abSmrg } /* decDivide */
67463d1a8abSmrg
67563d1a8abSmrg /* ------------------------------------------------------------------ */
67663d1a8abSmrg /* decFiniteMultiply -- multiply two finite decFloats */
67763d1a8abSmrg /* */
67863d1a8abSmrg /* num gets the result of multiplying dfl and dfr */
67963d1a8abSmrg /* bcdacc .. with the coefficient in this array */
68063d1a8abSmrg /* dfl is the first decFloat (lhs) */
68163d1a8abSmrg /* dfr is the second decFloat (rhs) */
68263d1a8abSmrg /* */
68363d1a8abSmrg /* This effects the multiplication of two decFloats, both known to be */
68463d1a8abSmrg /* finite, leaving the result in a bcdnum ready for decFinalize (for */
68563d1a8abSmrg /* use in Multiply) or in a following addition (FMA). */
68663d1a8abSmrg /* */
68763d1a8abSmrg /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
68863d1a8abSmrg /* No error is possible and no status is set. */
68963d1a8abSmrg /* ------------------------------------------------------------------ */
69063d1a8abSmrg /* This routine has two separate implementations of the core */
69163d1a8abSmrg /* multiplication; both using base-billion. One uses only 32-bit */
69263d1a8abSmrg /* variables (Ints and uInts) or smaller; the other uses uLongs (for */
69363d1a8abSmrg /* multiplication and addition only). Both implementations cover */
69463d1a8abSmrg /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
69563d1a8abSmrg /* comparisons. In any one compilation only one implementation for */
69663d1a8abSmrg /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
69763d1a8abSmrg /* version is forced. */
69863d1a8abSmrg /* */
69963d1a8abSmrg /* Historical note: an earlier version of this code also supported the */
70063d1a8abSmrg /* 256-bit format and has been preserved. That is somewhat trickier */
70163d1a8abSmrg /* during lazy carry splitting because the initial quotient estimate */
70263d1a8abSmrg /* (est) can exceed 32 bits. */
70363d1a8abSmrg
70463d1a8abSmrg #define MULTBASE ((uInt)BILLION) /* the base used for multiply */
70563d1a8abSmrg #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
70663d1a8abSmrg #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */
70763d1a8abSmrg #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
70863d1a8abSmrg
70963d1a8abSmrg /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
71063d1a8abSmrg #if DECEMAXD>9
71163d1a8abSmrg #error Exponent may overflow when doubled for Multiply
71263d1a8abSmrg #endif
71363d1a8abSmrg #if MULACCLEN!=(MULACCLEN/4)*4
71463d1a8abSmrg /* This assumption is used below only for initialization */
71563d1a8abSmrg #error MULACCLEN is not a multiple of 4
71663d1a8abSmrg #endif
71763d1a8abSmrg
decFiniteMultiply(bcdnum * num,uByte * bcdacc,const decFloat * dfl,const decFloat * dfr)71863d1a8abSmrg static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
71963d1a8abSmrg const decFloat *dfl, const decFloat *dfr) {
72063d1a8abSmrg uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
72163d1a8abSmrg uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
72263d1a8abSmrg uInt *ui, *uj; /* work */
72363d1a8abSmrg uByte *ub; /* .. */
72463d1a8abSmrg uInt uiwork; /* for macros */
72563d1a8abSmrg
72663d1a8abSmrg #if DECUSE64
72763d1a8abSmrg uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
72863d1a8abSmrg uLong *pl; /* work -> lazy accumulator */
72963d1a8abSmrg uInt acc[MULACCLEN]; /* coefficent in base-billion .. */
73063d1a8abSmrg #else
73163d1a8abSmrg uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */
73263d1a8abSmrg #endif
73363d1a8abSmrg uInt *pa; /* work -> accumulator */
73463d1a8abSmrg /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
73563d1a8abSmrg
73663d1a8abSmrg /* Calculate sign and exponent */
73763d1a8abSmrg num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
73863d1a8abSmrg num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
73963d1a8abSmrg
74063d1a8abSmrg /* Extract the coefficients and prepare the accumulator */
74163d1a8abSmrg /* the coefficients of the operands are decoded into base-billion */
74263d1a8abSmrg /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
74363d1a8abSmrg /* appropriate size. */
74463d1a8abSmrg GETCOEFFBILL(dfl, bufl);
74563d1a8abSmrg GETCOEFFBILL(dfr, bufr);
74663d1a8abSmrg #if DECTRACE && 0
74763d1a8abSmrg printf("CoeffbL:");
74863d1a8abSmrg for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
74963d1a8abSmrg printf("\n");
75063d1a8abSmrg printf("CoeffbR:");
75163d1a8abSmrg for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
75263d1a8abSmrg printf("\n");
75363d1a8abSmrg #endif
75463d1a8abSmrg
75563d1a8abSmrg /* start the 64-bit/32-bit differing paths... */
75663d1a8abSmrg #if DECUSE64
75763d1a8abSmrg
75863d1a8abSmrg /* zero the accumulator */
75963d1a8abSmrg #if MULACCLEN==4
76063d1a8abSmrg accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
76163d1a8abSmrg #else /* use a loop */
76263d1a8abSmrg /* MULACCLEN is a multiple of four, asserted above */
76363d1a8abSmrg for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
76463d1a8abSmrg *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
76563d1a8abSmrg } /* pl */
76663d1a8abSmrg #endif
76763d1a8abSmrg
76863d1a8abSmrg /* Effect the multiplication */
76963d1a8abSmrg /* The multiplcation proceeds using MFC's lazy-carry resolution */
77063d1a8abSmrg /* algorithm from decNumber. First, the multiplication is */
77163d1a8abSmrg /* effected, allowing accumulation of the partial products (which */
77263d1a8abSmrg /* are in base-billion at each column position) into 64 bits */
77363d1a8abSmrg /* without resolving back to base=billion after each addition. */
77463d1a8abSmrg /* These 64-bit numbers (which may contain up to 19 decimal digits) */
77563d1a8abSmrg /* are then split using the Clark & Cowlishaw algorithm (see below). */
77663d1a8abSmrg /* [Testing for 0 in the inner loop is not really a 'win'] */
77763d1a8abSmrg for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
77863d1a8abSmrg if (*ui==0) continue; /* product cannot affect result */
77963d1a8abSmrg pl=accl+(ui-bufr); /* where to add the lhs */
78063d1a8abSmrg for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
78163d1a8abSmrg /* if (*uj==0) continue; // product cannot affect result */
78263d1a8abSmrg *pl+=((uLong)*ui)*(*uj);
78363d1a8abSmrg } /* uj */
78463d1a8abSmrg } /* ui */
78563d1a8abSmrg
78663d1a8abSmrg /* The 64-bit carries must now be resolved; this means that a */
78763d1a8abSmrg /* quotient/remainder has to be calculated for base-billion (1E+9). */
78863d1a8abSmrg /* For this, Clark & Cowlishaw's quotient estimation approach (also */
78963d1a8abSmrg /* used in decNumber) is needed, because 64-bit divide is generally */
79063d1a8abSmrg /* extremely slow on 32-bit machines, and may be slower than this */
79163d1a8abSmrg /* approach even on 64-bit machines. This algorithm splits X */
79263d1a8abSmrg /* using: */
79363d1a8abSmrg /* */
79463d1a8abSmrg /* magic=2**(A+B)/1E+9; // 'magic number' */
79563d1a8abSmrg /* hop=X/2**A; // high order part of X (by shift) */
79663d1a8abSmrg /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
79763d1a8abSmrg /* */
79863d1a8abSmrg /* A and B are quite constrained; hop and magic must fit in 32 bits, */
79963d1a8abSmrg /* and 2**(A+B) must be as large as possible (which is 2**61 if */
80063d1a8abSmrg /* magic is to fit). Further, maxX increases with the length of */
80163d1a8abSmrg /* the operands (and hence the number of partial products */
80263d1a8abSmrg /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
80363d1a8abSmrg /* */
80463d1a8abSmrg /* It can be shown that when OPLEN is 2 then the maximum error in */
80563d1a8abSmrg /* the estimated quotient is <1, but for larger maximum x the */
80663d1a8abSmrg /* maximum error is above 1 so a correction that is >1 may be */
80763d1a8abSmrg /* needed. Values of A and B are chosen to satisfy the constraints */
80863d1a8abSmrg /* just mentioned while minimizing the maximum error (and hence the */
80963d1a8abSmrg /* maximum correction), as shown in the following table: */
81063d1a8abSmrg /* */
81163d1a8abSmrg /* Type OPLEN A B maxX maxError maxCorrection */
81263d1a8abSmrg /* --------------------------------------------------------- */
81363d1a8abSmrg /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
81463d1a8abSmrg /* QUAD 4 30 31 <4*10**18 1.17 2 */
81563d1a8abSmrg /* */
81663d1a8abSmrg /* In the OPLEN==2 case there is most choice, but the value for B */
81763d1a8abSmrg /* of 32 has a big advantage as then the calculation of the */
81863d1a8abSmrg /* estimate requires no shifting; the compiler can extract the high */
81963d1a8abSmrg /* word directly after multiplying magic*hop. */
82063d1a8abSmrg #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
82163d1a8abSmrg #if DOUBLE
82263d1a8abSmrg #define MULSHIFTA 29
82363d1a8abSmrg #define MULSHIFTB 32
82463d1a8abSmrg #elif QUAD
82563d1a8abSmrg #define MULSHIFTA 30
82663d1a8abSmrg #define MULSHIFTB 31
82763d1a8abSmrg #else
82863d1a8abSmrg #error Unexpected type
82963d1a8abSmrg #endif
83063d1a8abSmrg
83163d1a8abSmrg #if DECTRACE
83263d1a8abSmrg printf("MulAccl:");
83363d1a8abSmrg for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
83463d1a8abSmrg printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
83563d1a8abSmrg printf("\n");
83663d1a8abSmrg #endif
83763d1a8abSmrg
83863d1a8abSmrg for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
83963d1a8abSmrg uInt lo, hop; /* work */
84063d1a8abSmrg uInt est; /* cannot exceed 4E+9 */
84163d1a8abSmrg if (*pl>=MULTBASE) {
84263d1a8abSmrg /* *pl holds a binary number which needs to be split */
84363d1a8abSmrg hop=(uInt)(*pl>>MULSHIFTA);
84463d1a8abSmrg est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
84563d1a8abSmrg /* the estimate is now in est; now calculate hi:lo-est*10**9; */
84663d1a8abSmrg /* happily the top word of the result is irrelevant because it */
84763d1a8abSmrg /* will always be zero so this needs only one multiplication */
84863d1a8abSmrg lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */
84963d1a8abSmrg /* If QUAD, the correction here could be +2 */
85063d1a8abSmrg if (lo>=MULTBASE) {
85163d1a8abSmrg lo-=MULTBASE; /* correct by +1 */
85263d1a8abSmrg est++;
85363d1a8abSmrg #if QUAD
85463d1a8abSmrg /* may need to correct by +2 */
85563d1a8abSmrg if (lo>=MULTBASE) {
85663d1a8abSmrg lo-=MULTBASE;
85763d1a8abSmrg est++;
85863d1a8abSmrg }
85963d1a8abSmrg #endif
86063d1a8abSmrg }
86163d1a8abSmrg /* finally place lo as the new coefficient 'digit' and add est to */
86263d1a8abSmrg /* the next place up [this is safe because this path is never */
86363d1a8abSmrg /* taken on the final iteration as *pl will fit] */
86463d1a8abSmrg *pa=lo;
86563d1a8abSmrg *(pl+1)+=est;
86663d1a8abSmrg } /* *pl needed split */
86763d1a8abSmrg else { /* *pl<MULTBASE */
86863d1a8abSmrg *pa=(uInt)*pl; /* just copy across */
86963d1a8abSmrg }
87063d1a8abSmrg } /* pl loop */
87163d1a8abSmrg
87263d1a8abSmrg #else /* 32-bit */
87363d1a8abSmrg for (pa=acc;; pa+=4) { /* zero the accumulator */
87463d1a8abSmrg *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */
87563d1a8abSmrg if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */
87663d1a8abSmrg } /* pa */
87763d1a8abSmrg
87863d1a8abSmrg /* Effect the multiplication */
87963d1a8abSmrg /* uLongs are not available (and in particular, there is no uLong */
88063d1a8abSmrg /* divide) but it is still possible to use MFC's lazy-carry */
88163d1a8abSmrg /* resolution algorithm from decNumber. First, the multiplication */
88263d1a8abSmrg /* is effected, allowing accumulation of the partial products */
88363d1a8abSmrg /* (which are in base-billion at each column position) into 64 bits */
88463d1a8abSmrg /* [with the high-order 32 bits in each position being held at */
88563d1a8abSmrg /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
88663d1a8abSmrg /* These 64-bit numbers (which may contain up to 19 decimal digits) */
88763d1a8abSmrg /* are then split using the Clark & Cowlishaw algorithm (see */
88863d1a8abSmrg /* below). */
88963d1a8abSmrg for (ui=bufr;; ui++) { /* over each item in rhs */
89063d1a8abSmrg uInt hi, lo; /* words of exact multiply result */
89163d1a8abSmrg pa=acc+(ui-bufr); /* where to add the lhs */
89263d1a8abSmrg for (uj=bufl;; uj++, pa++) { /* over each item in lhs */
89363d1a8abSmrg LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */
89463d1a8abSmrg lo=(*ui)*(*uj); /* .. */
89563d1a8abSmrg *pa+=lo; /* accumulate low bits and .. */
89663d1a8abSmrg *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */
89763d1a8abSmrg if (uj==bufl+MULOPLEN-1) break;
89863d1a8abSmrg }
89963d1a8abSmrg if (ui==bufr+MULOPLEN-1) break;
90063d1a8abSmrg }
90163d1a8abSmrg
90263d1a8abSmrg /* The 64-bit carries must now be resolved; this means that a */
90363d1a8abSmrg /* quotient/remainder has to be calculated for base-billion (1E+9). */
90463d1a8abSmrg /* For this, Clark & Cowlishaw's quotient estimation approach (also */
90563d1a8abSmrg /* used in decNumber) is needed, because 64-bit divide is generally */
90663d1a8abSmrg /* extremely slow on 32-bit machines. This algorithm splits X */
90763d1a8abSmrg /* using: */
90863d1a8abSmrg /* */
90963d1a8abSmrg /* magic=2**(A+B)/1E+9; // 'magic number' */
91063d1a8abSmrg /* hop=X/2**A; // high order part of X (by shift) */
91163d1a8abSmrg /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
91263d1a8abSmrg /* */
91363d1a8abSmrg /* A and B are quite constrained; hop and magic must fit in 32 bits, */
91463d1a8abSmrg /* and 2**(A+B) must be as large as possible (which is 2**61 if */
91563d1a8abSmrg /* magic is to fit). Further, maxX increases with the length of */
91663d1a8abSmrg /* the operands (and hence the number of partial products */
91763d1a8abSmrg /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
91863d1a8abSmrg /* */
91963d1a8abSmrg /* It can be shown that when OPLEN is 2 then the maximum error in */
92063d1a8abSmrg /* the estimated quotient is <1, but for larger maximum x the */
92163d1a8abSmrg /* maximum error is above 1 so a correction that is >1 may be */
92263d1a8abSmrg /* needed. Values of A and B are chosen to satisfy the constraints */
92363d1a8abSmrg /* just mentioned while minimizing the maximum error (and hence the */
92463d1a8abSmrg /* maximum correction), as shown in the following table: */
92563d1a8abSmrg /* */
92663d1a8abSmrg /* Type OPLEN A B maxX maxError maxCorrection */
92763d1a8abSmrg /* --------------------------------------------------------- */
92863d1a8abSmrg /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
92963d1a8abSmrg /* QUAD 4 30 31 <4*10**18 1.17 2 */
93063d1a8abSmrg /* */
93163d1a8abSmrg /* In the OPLEN==2 case there is most choice, but the value for B */
93263d1a8abSmrg /* of 32 has a big advantage as then the calculation of the */
93363d1a8abSmrg /* estimate requires no shifting; the high word is simply */
93463d1a8abSmrg /* calculated from multiplying magic*hop. */
93563d1a8abSmrg #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
93663d1a8abSmrg #if DOUBLE
93763d1a8abSmrg #define MULSHIFTA 29
93863d1a8abSmrg #define MULSHIFTB 32
93963d1a8abSmrg #elif QUAD
94063d1a8abSmrg #define MULSHIFTA 30
94163d1a8abSmrg #define MULSHIFTB 31
94263d1a8abSmrg #else
94363d1a8abSmrg #error Unexpected type
94463d1a8abSmrg #endif
94563d1a8abSmrg
94663d1a8abSmrg #if DECTRACE
94763d1a8abSmrg printf("MulHiLo:");
94863d1a8abSmrg for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
94963d1a8abSmrg printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
95063d1a8abSmrg printf("\n");
95163d1a8abSmrg #endif
95263d1a8abSmrg
95363d1a8abSmrg for (pa=acc;; pa++) { /* each low uInt */
95463d1a8abSmrg uInt hi, lo; /* words of exact multiply result */
95563d1a8abSmrg uInt hop, estlo; /* work */
95663d1a8abSmrg #if QUAD
95763d1a8abSmrg uInt esthi; /* .. */
95863d1a8abSmrg #endif
95963d1a8abSmrg
96063d1a8abSmrg lo=*pa;
96163d1a8abSmrg hi=*(pa+MULACCLEN); /* top 32 bits */
96263d1a8abSmrg /* hi and lo now hold a binary number which needs to be split */
96363d1a8abSmrg
96463d1a8abSmrg #if DOUBLE
96563d1a8abSmrg hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */
96663d1a8abSmrg LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
96763d1a8abSmrg /* [MULSHIFTB is 32, so estlo can be used directly] */
96863d1a8abSmrg /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
96963d1a8abSmrg /* happily the top word of the result is irrelevant because it */
97063d1a8abSmrg /* will always be zero so this needs only one multiplication */
97163d1a8abSmrg lo-=(estlo*MULTBASE);
97263d1a8abSmrg /* esthi=0; // high word is ignored below */
97363d1a8abSmrg /* the correction here will be at most +1; do it */
97463d1a8abSmrg if (lo>=MULTBASE) {
97563d1a8abSmrg lo-=MULTBASE;
97663d1a8abSmrg estlo++;
97763d1a8abSmrg }
97863d1a8abSmrg #elif QUAD
97963d1a8abSmrg hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */
98063d1a8abSmrg LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
98163d1a8abSmrg estlo=hop*MULMAGIC; /* .. so low word needed */
98263d1a8abSmrg estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
98363d1a8abSmrg /* esthi=0; // high word is ignored below */
98463d1a8abSmrg lo-=(estlo*MULTBASE); /* as above */
98563d1a8abSmrg /* the correction here could be +1 or +2 */
98663d1a8abSmrg if (lo>=MULTBASE) {
98763d1a8abSmrg lo-=MULTBASE;
98863d1a8abSmrg estlo++;
98963d1a8abSmrg }
99063d1a8abSmrg if (lo>=MULTBASE) {
99163d1a8abSmrg lo-=MULTBASE;
99263d1a8abSmrg estlo++;
99363d1a8abSmrg }
99463d1a8abSmrg #else
99563d1a8abSmrg #error Unexpected type
99663d1a8abSmrg #endif
99763d1a8abSmrg
99863d1a8abSmrg /* finally place lo as the new accumulator digit and add est to */
99963d1a8abSmrg /* the next place up; this latter add could cause a carry of 1 */
100063d1a8abSmrg /* to the high word of the next place */
100163d1a8abSmrg *pa=lo;
100263d1a8abSmrg *(pa+1)+=estlo;
100363d1a8abSmrg /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
100463d1a8abSmrg /* *(pa+1+MULACCLEN)+=esthi; */
100563d1a8abSmrg if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
100663d1a8abSmrg if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */
100763d1a8abSmrg } /* pa loop */
100863d1a8abSmrg #endif
100963d1a8abSmrg
101063d1a8abSmrg /* At this point, whether using the 64-bit or the 32-bit paths, the */
101163d1a8abSmrg /* accumulator now holds the (unrounded) result in base-billion; */
101263d1a8abSmrg /* one base-billion 'digit' per uInt. */
101363d1a8abSmrg #if DECTRACE
101463d1a8abSmrg printf("MultAcc:");
101563d1a8abSmrg for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
101663d1a8abSmrg printf("\n");
101763d1a8abSmrg #endif
101863d1a8abSmrg
101963d1a8abSmrg /* Now convert to BCD for rounding and cleanup, starting from the */
102063d1a8abSmrg /* most significant end */
102163d1a8abSmrg pa=acc+MULACCLEN-1;
102263d1a8abSmrg if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
102363d1a8abSmrg else { /* >=1 word of leading zeros */
102463d1a8abSmrg num->msd=bcdacc; /* known leading zeros are gone */
102563d1a8abSmrg pa--; /* skip first word .. */
102663d1a8abSmrg for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
102763d1a8abSmrg }
102863d1a8abSmrg for (ub=bcdacc;; pa--, ub+=9) {
102963d1a8abSmrg if (*pa!=0) { /* split(s) needed */
103063d1a8abSmrg uInt top, mid, rem; /* work */
103163d1a8abSmrg /* *pa is non-zero -- split the base-billion acc digit into */
103263d1a8abSmrg /* hi, mid, and low three-digits */
103363d1a8abSmrg #define mulsplit9 1000000 /* divisor */
103463d1a8abSmrg #define mulsplit6 1000 /* divisor */
103563d1a8abSmrg /* The splitting is done by simple divides and remainders, */
103663d1a8abSmrg /* assuming the compiler will optimize these where useful */
103763d1a8abSmrg /* [GCC does] */
103863d1a8abSmrg top=*pa/mulsplit9;
103963d1a8abSmrg rem=*pa%mulsplit9;
104063d1a8abSmrg mid=rem/mulsplit6;
104163d1a8abSmrg rem=rem%mulsplit6;
104263d1a8abSmrg /* lay out the nine BCD digits (plus one unwanted byte) */
104363d1a8abSmrg UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
104463d1a8abSmrg UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
104563d1a8abSmrg UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
104663d1a8abSmrg }
104763d1a8abSmrg else { /* *pa==0 */
104863d1a8abSmrg UBFROMUI(ub, 0); /* clear 9 BCD8s */
104963d1a8abSmrg UBFROMUI(ub+4, 0); /* .. */
105063d1a8abSmrg *(ub+8)=0; /* .. */
105163d1a8abSmrg }
105263d1a8abSmrg if (pa==acc) break;
105363d1a8abSmrg } /* BCD conversion loop */
105463d1a8abSmrg
105563d1a8abSmrg num->lsd=ub+8; /* complete the bcdnum .. */
105663d1a8abSmrg
105763d1a8abSmrg #if DECTRACE
105863d1a8abSmrg decShowNum(num, "postmult");
105963d1a8abSmrg decFloatShow(dfl, "dfl");
106063d1a8abSmrg decFloatShow(dfr, "dfr");
106163d1a8abSmrg #endif
106263d1a8abSmrg return;
106363d1a8abSmrg } /* decFiniteMultiply */
106463d1a8abSmrg
106563d1a8abSmrg /* ------------------------------------------------------------------ */
106663d1a8abSmrg /* decFloatAbs -- absolute value, heeding NaNs, etc. */
106763d1a8abSmrg /* */
106863d1a8abSmrg /* result gets the canonicalized df with sign 0 */
106963d1a8abSmrg /* df is the decFloat to abs */
107063d1a8abSmrg /* set is the context */
107163d1a8abSmrg /* returns result */
107263d1a8abSmrg /* */
107363d1a8abSmrg /* This has the same effect as decFloatPlus unless df is negative, */
107463d1a8abSmrg /* in which case it has the same effect as decFloatMinus. The */
107563d1a8abSmrg /* effect is also the same as decFloatCopyAbs except that NaNs are */
107663d1a8abSmrg /* handled normally (the sign of a NaN is not affected, and an sNaN */
107763d1a8abSmrg /* will signal) and the result will be canonical. */
107863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatAbs(decFloat * result,const decFloat * df,decContext * set)107963d1a8abSmrg decFloat * decFloatAbs(decFloat *result, const decFloat *df,
108063d1a8abSmrg decContext *set) {
108163d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
108263d1a8abSmrg decCanonical(result, df); /* copy and check */
108363d1a8abSmrg DFBYTE(result, 0)&=~0x80; /* zero sign bit */
108463d1a8abSmrg return result;
108563d1a8abSmrg } /* decFloatAbs */
108663d1a8abSmrg
108763d1a8abSmrg /* ------------------------------------------------------------------ */
108863d1a8abSmrg /* decFloatAdd -- add two decFloats */
108963d1a8abSmrg /* */
109063d1a8abSmrg /* result gets the result of adding dfl and dfr: */
109163d1a8abSmrg /* dfl is the first decFloat (lhs) */
109263d1a8abSmrg /* dfr is the second decFloat (rhs) */
109363d1a8abSmrg /* set is the context */
109463d1a8abSmrg /* returns result */
109563d1a8abSmrg /* */
109663d1a8abSmrg /* ------------------------------------------------------------------ */
109763d1a8abSmrg #if QUAD
109863d1a8abSmrg /* Table for testing MSDs for fastpath elimination; returns the MSD of */
109963d1a8abSmrg /* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
110063d1a8abSmrg /* Infinities return -32 and NaNs return -128 so that summing the two */
110163d1a8abSmrg /* MSDs also allows rapid tests for the Specials (see code below). */
110263d1a8abSmrg const Int DECTESTMSD[64]={
110363d1a8abSmrg 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
110463d1a8abSmrg 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
110563d1a8abSmrg 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
110663d1a8abSmrg 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
110763d1a8abSmrg #else
110863d1a8abSmrg /* The table for testing MSDs is shared between the modules */
110963d1a8abSmrg extern const Int DECTESTMSD[64];
111063d1a8abSmrg #endif
111163d1a8abSmrg
decFloatAdd(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)111263d1a8abSmrg decFloat * decFloatAdd(decFloat *result,
111363d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
111463d1a8abSmrg decContext *set) {
111563d1a8abSmrg bcdnum num; /* for final conversion */
111663d1a8abSmrg Int bexpl, bexpr; /* left and right biased exponents */
111763d1a8abSmrg uByte *ub, *us, *ut; /* work */
111863d1a8abSmrg uInt uiwork; /* for macros */
111963d1a8abSmrg #if QUAD
112063d1a8abSmrg uShort uswork; /* .. */
112163d1a8abSmrg #endif
112263d1a8abSmrg
112363d1a8abSmrg uInt sourhil, sourhir; /* top words from source decFloats */
112463d1a8abSmrg /* [valid only through end of */
112563d1a8abSmrg /* fastpath code -- before swap] */
112663d1a8abSmrg uInt diffsign; /* non-zero if signs differ */
112763d1a8abSmrg uInt carry; /* carry: 0 or 1 before add loop */
112863d1a8abSmrg Int overlap; /* coefficient overlap (if full) */
112963d1a8abSmrg Int summ; /* sum of the MSDs */
113063d1a8abSmrg /* the following buffers hold coefficients with various alignments */
113163d1a8abSmrg /* (see commentary and diagrams below) */
113263d1a8abSmrg uByte acc[4+2+DECPMAX*3+8];
113363d1a8abSmrg uByte buf[4+2+DECPMAX*2];
113463d1a8abSmrg uByte *umsd, *ulsd; /* local MSD and LSD pointers */
113563d1a8abSmrg
113663d1a8abSmrg #if DECLITEND
113763d1a8abSmrg #define CARRYPAT 0x01000000 /* carry=1 pattern */
113863d1a8abSmrg #else
113963d1a8abSmrg #define CARRYPAT 0x00000001 /* carry=1 pattern */
114063d1a8abSmrg #endif
114163d1a8abSmrg
114263d1a8abSmrg /* Start decoding the arguments */
114363d1a8abSmrg /* The initial exponents are placed into the opposite Ints to */
114463d1a8abSmrg /* that which might be expected; there are two sets of data to */
114563d1a8abSmrg /* keep track of (each decFloat and the corresponding exponent), */
114663d1a8abSmrg /* and this scheme means that at the swap point (after comparing */
114763d1a8abSmrg /* exponents) only one pair of words needs to be swapped */
114863d1a8abSmrg /* whichever path is taken (thereby minimising worst-case path). */
114963d1a8abSmrg /* The calculated exponents will be nonsense when the arguments are */
115063d1a8abSmrg /* Special, but are not used in that path */
115163d1a8abSmrg sourhil=DFWORD(dfl, 0); /* LHS top word */
115263d1a8abSmrg summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */
115363d1a8abSmrg bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
115463d1a8abSmrg bexpr+=GETECON(dfl); /* .. + continuation */
115563d1a8abSmrg
115663d1a8abSmrg sourhir=DFWORD(dfr, 0); /* RHS top word */
115763d1a8abSmrg summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */
115863d1a8abSmrg bexpl=DECCOMBEXP[sourhir>>26];
115963d1a8abSmrg bexpl+=GETECON(dfr);
116063d1a8abSmrg
116163d1a8abSmrg /* here bexpr has biased exponent from lhs, and vice versa */
116263d1a8abSmrg
116363d1a8abSmrg diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
116463d1a8abSmrg
116563d1a8abSmrg /* now determine whether to take a fast path or the full-function */
116663d1a8abSmrg /* slow path. The slow path must be taken when: */
116763d1a8abSmrg /* -- both numbers are finite, and: */
116863d1a8abSmrg /* the exponents are different, or */
116963d1a8abSmrg /* the signs are different, or */
117063d1a8abSmrg /* the sum of the MSDs is >8 (hence might overflow) */
117163d1a8abSmrg /* specialness and the sum of the MSDs can be tested at once using */
117263d1a8abSmrg /* the summ value just calculated, so the test for specials is no */
117363d1a8abSmrg /* longer on the worst-case path (as of 3.60) */
117463d1a8abSmrg
117563d1a8abSmrg if (summ<=8) { /* MSD+MSD is good, or there is a special */
117663d1a8abSmrg if (summ<0) { /* there is a special */
117763d1a8abSmrg /* Inf+Inf would give -64; Inf+finite is -32 or higher */
117863d1a8abSmrg if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */
117963d1a8abSmrg /* two infinities with different signs is invalid */
118063d1a8abSmrg if (summ==-64 && diffsign) return decInvalid(result, set);
118163d1a8abSmrg if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
118263d1a8abSmrg return decInfinity(result, dfr); /* RHS must be Inf */
118363d1a8abSmrg }
118463d1a8abSmrg /* Here when both arguments are finite; fast path is possible */
118563d1a8abSmrg /* (currently only for aligned and same-sign) */
118663d1a8abSmrg if (bexpr==bexpl && !diffsign) {
118763d1a8abSmrg uInt tac[DECLETS+1]; /* base-1000 coefficient */
118863d1a8abSmrg uInt encode; /* work */
118963d1a8abSmrg
119063d1a8abSmrg /* Get one coefficient as base-1000 and add the other */
119163d1a8abSmrg GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */
119263d1a8abSmrg ADDCOEFFTHOU(dfr, tac);
119363d1a8abSmrg /* here the sum of the MSDs (plus any carry) will be <10 due to */
119463d1a8abSmrg /* the fastpath test earlier */
119563d1a8abSmrg
119663d1a8abSmrg /* construct the result; low word is the same for both formats */
119763d1a8abSmrg encode =BIN2DPD[tac[0]];
119863d1a8abSmrg encode|=BIN2DPD[tac[1]]<<10;
119963d1a8abSmrg encode|=BIN2DPD[tac[2]]<<20;
120063d1a8abSmrg encode|=BIN2DPD[tac[3]]<<30;
120163d1a8abSmrg DFWORD(result, (DECBYTES/4)-1)=encode;
120263d1a8abSmrg
120363d1a8abSmrg /* collect next two declets (all that remains, for Double) */
120463d1a8abSmrg encode =BIN2DPD[tac[3]]>>2;
120563d1a8abSmrg encode|=BIN2DPD[tac[4]]<<8;
120663d1a8abSmrg
120763d1a8abSmrg #if QUAD
120863d1a8abSmrg /* complete and lay out middling words */
120963d1a8abSmrg encode|=BIN2DPD[tac[5]]<<18;
121063d1a8abSmrg encode|=BIN2DPD[tac[6]]<<28;
121163d1a8abSmrg DFWORD(result, 2)=encode;
121263d1a8abSmrg
121363d1a8abSmrg encode =BIN2DPD[tac[6]]>>4;
121463d1a8abSmrg encode|=BIN2DPD[tac[7]]<<6;
121563d1a8abSmrg encode|=BIN2DPD[tac[8]]<<16;
121663d1a8abSmrg encode|=BIN2DPD[tac[9]]<<26;
121763d1a8abSmrg DFWORD(result, 1)=encode;
121863d1a8abSmrg
121963d1a8abSmrg /* and final two declets */
122063d1a8abSmrg encode =BIN2DPD[tac[9]]>>6;
122163d1a8abSmrg encode|=BIN2DPD[tac[10]]<<4;
122263d1a8abSmrg #endif
122363d1a8abSmrg
122463d1a8abSmrg /* add exponent continuation and sign (from either argument) */
122563d1a8abSmrg encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
122663d1a8abSmrg
122763d1a8abSmrg /* create lookup index = MSD + top two bits of biased exponent <<4 */
122863d1a8abSmrg tac[DECLETS]|=(bexpl>>DECECONL)<<4;
122963d1a8abSmrg encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
123063d1a8abSmrg DFWORD(result, 0)=encode; /* complete */
123163d1a8abSmrg
123263d1a8abSmrg /* decFloatShow(result, ">"); */
123363d1a8abSmrg return result;
123463d1a8abSmrg } /* fast path OK */
123563d1a8abSmrg /* drop through to slow path */
123663d1a8abSmrg } /* low sum or Special(s) */
123763d1a8abSmrg
123863d1a8abSmrg /* Slow path required -- arguments are finite and might overflow, */
123963d1a8abSmrg /* or require alignment, or might have different signs */
124063d1a8abSmrg
124163d1a8abSmrg /* now swap either exponents or argument pointers */
124263d1a8abSmrg if (bexpl<=bexpr) {
124363d1a8abSmrg /* original left is bigger */
124463d1a8abSmrg Int bexpswap=bexpl;
124563d1a8abSmrg bexpl=bexpr;
124663d1a8abSmrg bexpr=bexpswap;
124763d1a8abSmrg /* printf("left bigger\n"); */
124863d1a8abSmrg }
124963d1a8abSmrg else {
125063d1a8abSmrg const decFloat *dfswap=dfl;
125163d1a8abSmrg dfl=dfr;
125263d1a8abSmrg dfr=dfswap;
125363d1a8abSmrg /* printf("right bigger\n"); */
125463d1a8abSmrg }
125563d1a8abSmrg /* [here dfl and bexpl refer to the datum with the larger exponent, */
125663d1a8abSmrg /* of if the exponents are equal then the original LHS argument] */
125763d1a8abSmrg
125863d1a8abSmrg /* if lhs is zero then result will be the rhs (now known to have */
125963d1a8abSmrg /* the smaller exponent), which also may need to be tested for zero */
126063d1a8abSmrg /* for the weird IEEE 754 sign rules */
126163d1a8abSmrg if (DFISZERO(dfl)) {
126263d1a8abSmrg decCanonical(result, dfr); /* clean copy */
126363d1a8abSmrg /* "When the sum of two operands with opposite signs is */
126463d1a8abSmrg /* exactly zero, the sign of that sum shall be '+' in all */
126563d1a8abSmrg /* rounding modes except round toward -Infinity, in which */
126663d1a8abSmrg /* mode that sign shall be '-'." */
126763d1a8abSmrg if (diffsign && DFISZERO(result)) {
126863d1a8abSmrg DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */
126963d1a8abSmrg if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
127063d1a8abSmrg }
127163d1a8abSmrg return result;
127263d1a8abSmrg } /* numfl is zero */
127363d1a8abSmrg /* [here, LHS is non-zero; code below assumes that] */
127463d1a8abSmrg
127563d1a8abSmrg /* Coefficients layout during the calculations to follow: */
127663d1a8abSmrg /* */
127763d1a8abSmrg /* Overlap case: */
127863d1a8abSmrg /* +------------------------------------------------+ */
127963d1a8abSmrg /* acc: |0000| coeffa | tail B | | */
128063d1a8abSmrg /* +------------------------------------------------+ */
128163d1a8abSmrg /* buf: |0000| pad0s | coeffb | | */
128263d1a8abSmrg /* +------------------------------------------------+ */
128363d1a8abSmrg /* */
128463d1a8abSmrg /* Touching coefficients or gap: */
128563d1a8abSmrg /* +------------------------------------------------+ */
128663d1a8abSmrg /* acc: |0000| coeffa | gap | coeffb | */
128763d1a8abSmrg /* +------------------------------------------------+ */
128863d1a8abSmrg /* [buf not used or needed; gap clamped to Pmax] */
128963d1a8abSmrg
129063d1a8abSmrg /* lay out lhs coefficient into accumulator; this starts at acc+4 */
129163d1a8abSmrg /* for decDouble or acc+6 for decQuad so the LSD is word- */
129263d1a8abSmrg /* aligned; the top word gap is there only in case a carry digit */
129363d1a8abSmrg /* is prefixed after the add -- it does not need to be zeroed */
129463d1a8abSmrg #if DOUBLE
129563d1a8abSmrg #define COFF 4 /* offset into acc */
129663d1a8abSmrg #elif QUAD
129763d1a8abSmrg UBFROMUS(acc+4, 0); /* prefix 00 */
129863d1a8abSmrg #define COFF 6 /* offset into acc */
129963d1a8abSmrg #endif
130063d1a8abSmrg
130163d1a8abSmrg GETCOEFF(dfl, acc+COFF); /* decode from decFloat */
130263d1a8abSmrg ulsd=acc+COFF+DECPMAX-1;
130363d1a8abSmrg umsd=acc+4; /* [having this here avoids */
130463d1a8abSmrg
130563d1a8abSmrg #if DECTRACE
130663d1a8abSmrg {bcdnum tum;
130763d1a8abSmrg tum.msd=umsd;
130863d1a8abSmrg tum.lsd=ulsd;
130963d1a8abSmrg tum.exponent=bexpl-DECBIAS;
131063d1a8abSmrg tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
131163d1a8abSmrg decShowNum(&tum, "dflx");}
131263d1a8abSmrg #endif
131363d1a8abSmrg
131463d1a8abSmrg /* if signs differ, take ten's complement of lhs (here the */
131563d1a8abSmrg /* coefficient is subtracted from all-nines; the 1 is added during */
131663d1a8abSmrg /* the later add cycle -- zeros to the right do not matter because */
131763d1a8abSmrg /* the complement of zero is zero); these are fixed-length inverts */
131863d1a8abSmrg /* where the lsd is known to be at a 4-byte boundary (so no borrow */
131963d1a8abSmrg /* possible) */
132063d1a8abSmrg carry=0; /* assume no carry */
132163d1a8abSmrg if (diffsign) {
132263d1a8abSmrg carry=CARRYPAT; /* for +1 during add */
132363d1a8abSmrg UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
132463d1a8abSmrg UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
132563d1a8abSmrg UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
132663d1a8abSmrg UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
132763d1a8abSmrg #if QUAD
132863d1a8abSmrg UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
132963d1a8abSmrg UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
133063d1a8abSmrg UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
133163d1a8abSmrg UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
133263d1a8abSmrg UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
133363d1a8abSmrg #endif
133463d1a8abSmrg } /* diffsign */
133563d1a8abSmrg
133663d1a8abSmrg /* now process the rhs coefficient; if it cannot overlap lhs then */
133763d1a8abSmrg /* it can be put straight into acc (with an appropriate gap, if */
133863d1a8abSmrg /* needed) because no actual addition will be needed (except */
133963d1a8abSmrg /* possibly to complete ten's complement) */
134063d1a8abSmrg overlap=DECPMAX-(bexpl-bexpr);
134163d1a8abSmrg #if DECTRACE
134263d1a8abSmrg printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
134363d1a8abSmrg printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
134463d1a8abSmrg #endif
134563d1a8abSmrg
134663d1a8abSmrg if (overlap<=0) { /* no overlap possible */
134763d1a8abSmrg uInt gap; /* local work */
134863d1a8abSmrg /* since a full addition is not needed, a ten's complement */
134963d1a8abSmrg /* calculation started above may need to be completed */
135063d1a8abSmrg if (carry) {
135163d1a8abSmrg for (ub=ulsd; *ub==9; ub--) *ub=0;
135263d1a8abSmrg *ub+=1;
135363d1a8abSmrg carry=0; /* taken care of */
135463d1a8abSmrg }
135563d1a8abSmrg /* up to DECPMAX-1 digits of the final result can extend down */
135663d1a8abSmrg /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
135763d1a8abSmrg /* rhs will be simply sticky bits. In this case the gap is */
135863d1a8abSmrg /* clamped to DECPMAX and the exponent adjusted to suit [this is */
135963d1a8abSmrg /* safe because the lhs is non-zero]. */
136063d1a8abSmrg gap=-overlap;
136163d1a8abSmrg if (gap>DECPMAX) {
136263d1a8abSmrg bexpr+=gap-1;
136363d1a8abSmrg gap=DECPMAX;
136463d1a8abSmrg }
136563d1a8abSmrg ub=ulsd+gap+1; /* where MSD will go */
136663d1a8abSmrg /* Fill the gap with 0s; note that there is no addition to do */
136763d1a8abSmrg ut=acc+COFF+DECPMAX; /* start of gap */
136863d1a8abSmrg for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */
136963d1a8abSmrg if (overlap<-DECPMAX) { /* gap was > DECPMAX */
137063d1a8abSmrg *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
137163d1a8abSmrg }
137263d1a8abSmrg else { /* need full coefficient */
137363d1a8abSmrg GETCOEFF(dfr, ub); /* decode from decFloat */
137463d1a8abSmrg ub+=DECPMAX-1; /* new LSD... */
137563d1a8abSmrg }
137663d1a8abSmrg ulsd=ub; /* save new LSD */
137763d1a8abSmrg } /* no overlap possible */
137863d1a8abSmrg
137963d1a8abSmrg else { /* overlap>0 */
138063d1a8abSmrg /* coefficients overlap (perhaps completely, although also */
138163d1a8abSmrg /* perhaps only where zeros) */
138263d1a8abSmrg if (overlap==DECPMAX) { /* aligned */
138363d1a8abSmrg ub=buf+COFF; /* where msd will go */
138463d1a8abSmrg #if QUAD
138563d1a8abSmrg UBFROMUS(buf+4, 0); /* clear quad's 00 */
138663d1a8abSmrg #endif
138763d1a8abSmrg GETCOEFF(dfr, ub); /* decode from decFloat */
138863d1a8abSmrg }
138963d1a8abSmrg else { /* unaligned */
139063d1a8abSmrg ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */
139163d1a8abSmrg /* Fill the prefix gap with 0s; 8 will cover most common */
139263d1a8abSmrg /* unalignments, so start with direct assignments (a loop is */
139363d1a8abSmrg /* then used for any remaining -- the loop (and the one in a */
139463d1a8abSmrg /* moment) is not then on the critical path because the number */
139563d1a8abSmrg /* of additions is reduced by (at least) two in this case) */
139663d1a8abSmrg UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */
139763d1a8abSmrg UBFROMUI(buf+8, 0);
139863d1a8abSmrg if (ub>buf+12) {
139963d1a8abSmrg ut=buf+12; /* start any remaining */
140063d1a8abSmrg for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
140163d1a8abSmrg }
140263d1a8abSmrg GETCOEFF(dfr, ub); /* decode from decFloat */
140363d1a8abSmrg
140463d1a8abSmrg /* now move tail of rhs across to main acc; again use direct */
140563d1a8abSmrg /* copies for 8 digits-worth */
140663d1a8abSmrg UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX));
140763d1a8abSmrg UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
140863d1a8abSmrg if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
140963d1a8abSmrg us=buf+COFF+DECPMAX+8; /* source */
141063d1a8abSmrg ut=acc+COFF+DECPMAX+8; /* target */
141163d1a8abSmrg for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
141263d1a8abSmrg }
141363d1a8abSmrg } /* unaligned */
141463d1a8abSmrg
141563d1a8abSmrg ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
141663d1a8abSmrg
141763d1a8abSmrg /* Now do the add of the non-tail; this is all nicely aligned, */
141863d1a8abSmrg /* and is over a multiple of four digits (because for Quad two */
141963d1a8abSmrg /* zero digits were added on the left); words in both acc and */
142063d1a8abSmrg /* buf (buf especially) will often be zero */
142163d1a8abSmrg /* [byte-by-byte add, here, is about 15% slower total effect than */
142263d1a8abSmrg /* the by-fours] */
142363d1a8abSmrg
142463d1a8abSmrg /* Now effect the add; this is harder on a little-endian */
142563d1a8abSmrg /* machine as the inter-digit carry cannot use the usual BCD */
142663d1a8abSmrg /* addition trick because the bytes are loaded in the wrong order */
142763d1a8abSmrg /* [this loop could be unrolled, but probably scarcely worth it] */
142863d1a8abSmrg
142963d1a8abSmrg ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */
143063d1a8abSmrg us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */
143163d1a8abSmrg
143263d1a8abSmrg #if !DECLITEND
143363d1a8abSmrg for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */
143463d1a8abSmrg /* bcd8 add */
143563d1a8abSmrg carry+=UBTOUI(us); /* rhs + carry */
143663d1a8abSmrg if (carry==0) continue; /* no-op */
143763d1a8abSmrg carry+=UBTOUI(ut); /* lhs */
143863d1a8abSmrg /* Big-endian BCD adjust (uses internal carry) */
143963d1a8abSmrg carry+=0x76f6f6f6; /* note top nibble not all bits */
144063d1a8abSmrg /* apply BCD adjust and save */
144163d1a8abSmrg UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
144263d1a8abSmrg carry>>=31; /* true carry was at far left */
144363d1a8abSmrg } /* add loop */
144463d1a8abSmrg #else
144563d1a8abSmrg for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */
144663d1a8abSmrg /* bcd8 add */
144763d1a8abSmrg carry+=UBTOUI(us); /* rhs + carry */
144863d1a8abSmrg if (carry==0) continue; /* no-op [common if unaligned] */
144963d1a8abSmrg carry+=UBTOUI(ut); /* lhs */
145063d1a8abSmrg /* Little-endian BCD adjust; inter-digit carry must be manual */
145163d1a8abSmrg /* because the lsb from the array will be in the most-significant */
145263d1a8abSmrg /* byte of carry */
145363d1a8abSmrg carry+=0x76767676; /* note no inter-byte carries */
145463d1a8abSmrg carry+=(carry & 0x80000000)>>15;
145563d1a8abSmrg carry+=(carry & 0x00800000)>>15;
145663d1a8abSmrg carry+=(carry & 0x00008000)>>15;
145763d1a8abSmrg carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
145863d1a8abSmrg UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */
145963d1a8abSmrg /* here, final carry-out bit is at 0x00000080; move it ready */
146063d1a8abSmrg /* for next word-add (i.e., to 0x01000000) */
146163d1a8abSmrg carry=(carry & 0x00000080)<<17;
146263d1a8abSmrg } /* add loop */
146363d1a8abSmrg #endif
146463d1a8abSmrg
146563d1a8abSmrg #if DECTRACE
146663d1a8abSmrg {bcdnum tum;
146763d1a8abSmrg printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
146863d1a8abSmrg tum.msd=umsd; /* acc+4; */
146963d1a8abSmrg tum.lsd=ulsd;
147063d1a8abSmrg tum.exponent=0;
147163d1a8abSmrg tum.sign=0;
147263d1a8abSmrg decShowNum(&tum, "dfadd");}
147363d1a8abSmrg #endif
147463d1a8abSmrg } /* overlap possible */
147563d1a8abSmrg
147663d1a8abSmrg /* ordering here is a little strange in order to have slowest path */
147763d1a8abSmrg /* first in GCC asm listing */
147863d1a8abSmrg if (diffsign) { /* subtraction */
147963d1a8abSmrg if (!carry) { /* no carry out means RHS<LHS */
148063d1a8abSmrg /* borrowed -- take ten's complement */
148163d1a8abSmrg /* sign is lhs sign */
148263d1a8abSmrg num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
148363d1a8abSmrg
148463d1a8abSmrg /* invert the coefficient first by fours, then add one; space */
148563d1a8abSmrg /* at the end of the buffer ensures the by-fours is always */
148663d1a8abSmrg /* safe, but lsd+1 must be cleared to prevent a borrow */
148763d1a8abSmrg /* if big-endian */
148863d1a8abSmrg #if !DECLITEND
148963d1a8abSmrg *(ulsd+1)=0;
149063d1a8abSmrg #endif
149163d1a8abSmrg /* there are always at least four coefficient words */
149263d1a8abSmrg UBFROMUI(umsd, 0x09090909-UBTOUI(umsd));
149363d1a8abSmrg UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4));
149463d1a8abSmrg UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8));
149563d1a8abSmrg UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
149663d1a8abSmrg #if DOUBLE
149763d1a8abSmrg #define BNEXT 16
149863d1a8abSmrg #elif QUAD
149963d1a8abSmrg UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
150063d1a8abSmrg UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
150163d1a8abSmrg UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
150263d1a8abSmrg UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
150363d1a8abSmrg UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32));
150463d1a8abSmrg #define BNEXT 36
150563d1a8abSmrg #endif
150663d1a8abSmrg if (ulsd>=umsd+BNEXT) { /* unaligned */
150763d1a8abSmrg /* eight will handle most unaligments for Double; 16 for Quad */
150863d1a8abSmrg UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT));
150963d1a8abSmrg UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
151063d1a8abSmrg #if DOUBLE
151163d1a8abSmrg #define BNEXTY (BNEXT+8)
151263d1a8abSmrg #elif QUAD
151363d1a8abSmrg UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8));
151463d1a8abSmrg UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
151563d1a8abSmrg #define BNEXTY (BNEXT+16)
151663d1a8abSmrg #endif
151763d1a8abSmrg if (ulsd>=umsd+BNEXTY) { /* very unaligned */
151863d1a8abSmrg ut=umsd+BNEXTY; /* -> continue */
151963d1a8abSmrg for (;;ut+=4) {
152063d1a8abSmrg UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
152163d1a8abSmrg if (ut>=ulsd-3) break; /* all done */
152263d1a8abSmrg }
152363d1a8abSmrg }
152463d1a8abSmrg }
152563d1a8abSmrg /* complete the ten's complement by adding 1 */
152663d1a8abSmrg for (ub=ulsd; *ub==9; ub--) *ub=0;
152763d1a8abSmrg *ub+=1;
152863d1a8abSmrg } /* borrowed */
152963d1a8abSmrg
153063d1a8abSmrg else { /* carry out means RHS>=LHS */
153163d1a8abSmrg num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
153263d1a8abSmrg /* all done except for the special IEEE 754 exact-zero-result */
153363d1a8abSmrg /* rule (see above); while testing for zero, strip leading */
153463d1a8abSmrg /* zeros (which will save decFinalize doing it) (this is in */
153563d1a8abSmrg /* diffsign path, so carry impossible and true umsd is */
153663d1a8abSmrg /* acc+COFF) */
153763d1a8abSmrg
153863d1a8abSmrg /* Check the initial coefficient area using the fast macro; */
153963d1a8abSmrg /* this will often be all that needs to be done (as on the */
154063d1a8abSmrg /* worst-case path when the subtraction was aligned and */
154163d1a8abSmrg /* full-length) */
154263d1a8abSmrg if (ISCOEFFZERO(acc+COFF)) {
154363d1a8abSmrg umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
154463d1a8abSmrg if (ulsd>umsd) { /* more to check */
154563d1a8abSmrg umsd++; /* to align after checked area */
154663d1a8abSmrg for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
154763d1a8abSmrg for (; *umsd==0 && umsd<ulsd;) umsd++;
154863d1a8abSmrg }
154963d1a8abSmrg if (*umsd==0) { /* must be true zero (and diffsign) */
155063d1a8abSmrg num.sign=0; /* assume + */
155163d1a8abSmrg if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
155263d1a8abSmrg }
155363d1a8abSmrg }
155463d1a8abSmrg /* [else was not zero, might still have leading zeros] */
155563d1a8abSmrg } /* subtraction gave positive result */
155663d1a8abSmrg } /* diffsign */
155763d1a8abSmrg
155863d1a8abSmrg else { /* same-sign addition */
155963d1a8abSmrg num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
156063d1a8abSmrg #if DOUBLE
156163d1a8abSmrg if (carry) { /* only possible with decDouble */
156263d1a8abSmrg *(acc+3)=1; /* [Quad has leading 00] */
156363d1a8abSmrg umsd=acc+3;
156463d1a8abSmrg }
156563d1a8abSmrg #endif
156663d1a8abSmrg } /* same sign */
156763d1a8abSmrg
156863d1a8abSmrg num.msd=umsd; /* set MSD .. */
156963d1a8abSmrg num.lsd=ulsd; /* .. and LSD */
157063d1a8abSmrg num.exponent=bexpr-DECBIAS; /* set exponent to smaller, unbiassed */
157163d1a8abSmrg
157263d1a8abSmrg #if DECTRACE
157363d1a8abSmrg decFloatShow(dfl, "dfl");
157463d1a8abSmrg decFloatShow(dfr, "dfr");
157563d1a8abSmrg decShowNum(&num, "postadd");
157663d1a8abSmrg #endif
157763d1a8abSmrg return decFinalize(result, &num, set); /* round, check, and lay out */
157863d1a8abSmrg } /* decFloatAdd */
157963d1a8abSmrg
158063d1a8abSmrg /* ------------------------------------------------------------------ */
158163d1a8abSmrg /* decFloatAnd -- logical digitwise AND of two decFloats */
158263d1a8abSmrg /* */
158363d1a8abSmrg /* result gets the result of ANDing dfl and dfr */
158463d1a8abSmrg /* dfl is the first decFloat (lhs) */
158563d1a8abSmrg /* dfr is the second decFloat (rhs) */
158663d1a8abSmrg /* set is the context */
158763d1a8abSmrg /* returns result, which will be canonical with sign=0 */
158863d1a8abSmrg /* */
158963d1a8abSmrg /* The operands must be positive, finite with exponent q=0, and */
159063d1a8abSmrg /* comprise just zeros and ones; if not, Invalid operation results. */
159163d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatAnd(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)159263d1a8abSmrg decFloat * decFloatAnd(decFloat *result,
159363d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
159463d1a8abSmrg decContext *set) {
159563d1a8abSmrg if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
159663d1a8abSmrg || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
159763d1a8abSmrg /* the operands are positive finite integers (q=0) with just 0s and 1s */
159863d1a8abSmrg #if DOUBLE
159963d1a8abSmrg DFWORD(result, 0)=ZEROWORD
160063d1a8abSmrg |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
160163d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
160263d1a8abSmrg #elif QUAD
160363d1a8abSmrg DFWORD(result, 0)=ZEROWORD
160463d1a8abSmrg |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
160563d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
160663d1a8abSmrg DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
160763d1a8abSmrg DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
160863d1a8abSmrg #endif
160963d1a8abSmrg return result;
161063d1a8abSmrg } /* decFloatAnd */
161163d1a8abSmrg
161263d1a8abSmrg /* ------------------------------------------------------------------ */
161363d1a8abSmrg /* decFloatCanonical -- copy a decFloat, making canonical */
161463d1a8abSmrg /* */
161563d1a8abSmrg /* result gets the canonicalized df */
161663d1a8abSmrg /* df is the decFloat to copy and make canonical */
161763d1a8abSmrg /* returns result */
161863d1a8abSmrg /* */
161963d1a8abSmrg /* This works on specials, too; no error or exception is possible. */
162063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCanonical(decFloat * result,const decFloat * df)162163d1a8abSmrg decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
162263d1a8abSmrg return decCanonical(result, df);
162363d1a8abSmrg } /* decFloatCanonical */
162463d1a8abSmrg
162563d1a8abSmrg /* ------------------------------------------------------------------ */
162663d1a8abSmrg /* decFloatClass -- return the class of a decFloat */
162763d1a8abSmrg /* */
162863d1a8abSmrg /* df is the decFloat to test */
162963d1a8abSmrg /* returns the decClass that df falls into */
163063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatClass(const decFloat * df)163163d1a8abSmrg enum decClass decFloatClass(const decFloat *df) {
163263d1a8abSmrg Int exp; /* exponent */
163363d1a8abSmrg if (DFISSPECIAL(df)) {
163463d1a8abSmrg if (DFISQNAN(df)) return DEC_CLASS_QNAN;
163563d1a8abSmrg if (DFISSNAN(df)) return DEC_CLASS_SNAN;
163663d1a8abSmrg /* must be an infinity */
163763d1a8abSmrg if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
163863d1a8abSmrg return DEC_CLASS_POS_INF;
163963d1a8abSmrg }
164063d1a8abSmrg if (DFISZERO(df)) { /* quite common */
164163d1a8abSmrg if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
164263d1a8abSmrg return DEC_CLASS_POS_ZERO;
164363d1a8abSmrg }
164463d1a8abSmrg /* is finite and non-zero; similar code to decFloatIsNormal, here */
164563d1a8abSmrg /* [this could be speeded up slightly by in-lining decFloatDigits] */
164663d1a8abSmrg exp=GETEXPUN(df) /* get unbiased exponent .. */
164763d1a8abSmrg +decFloatDigits(df)-1; /* .. and make adjusted exponent */
164863d1a8abSmrg if (exp>=DECEMIN) { /* is normal */
164963d1a8abSmrg if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
165063d1a8abSmrg return DEC_CLASS_POS_NORMAL;
165163d1a8abSmrg }
165263d1a8abSmrg /* is subnormal */
165363d1a8abSmrg if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
165463d1a8abSmrg return DEC_CLASS_POS_SUBNORMAL;
165563d1a8abSmrg } /* decFloatClass */
165663d1a8abSmrg
165763d1a8abSmrg /* ------------------------------------------------------------------ */
165863d1a8abSmrg /* decFloatClassString -- return the class of a decFloat as a string */
165963d1a8abSmrg /* */
166063d1a8abSmrg /* df is the decFloat to test */
166163d1a8abSmrg /* returns a constant string describing the class df falls into */
166263d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatClassString(const decFloat * df)166363d1a8abSmrg const char *decFloatClassString(const decFloat *df) {
166463d1a8abSmrg enum decClass eclass=decFloatClass(df);
166563d1a8abSmrg if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
166663d1a8abSmrg if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
166763d1a8abSmrg if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
166863d1a8abSmrg if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
166963d1a8abSmrg if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
167063d1a8abSmrg if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
167163d1a8abSmrg if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
167263d1a8abSmrg if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
167363d1a8abSmrg if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
167463d1a8abSmrg if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
167563d1a8abSmrg return DEC_ClassString_UN; /* Unknown */
167663d1a8abSmrg } /* decFloatClassString */
167763d1a8abSmrg
167863d1a8abSmrg /* ------------------------------------------------------------------ */
167963d1a8abSmrg /* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
168063d1a8abSmrg /* */
168163d1a8abSmrg /* result gets the result of comparing dfl and dfr */
168263d1a8abSmrg /* dfl is the first decFloat (lhs) */
168363d1a8abSmrg /* dfr is the second decFloat (rhs) */
168463d1a8abSmrg /* set is the context */
168563d1a8abSmrg /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
168663d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCompare(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)168763d1a8abSmrg decFloat * decFloatCompare(decFloat *result,
168863d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
168963d1a8abSmrg decContext *set) {
169063d1a8abSmrg Int comp; /* work */
169163d1a8abSmrg /* NaNs are handled as usual */
169263d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
169363d1a8abSmrg /* numeric comparison needed */
169463d1a8abSmrg comp=decNumCompare(dfl, dfr, 0);
169563d1a8abSmrg decFloatZero(result);
169663d1a8abSmrg if (comp==0) return result;
169763d1a8abSmrg DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
169863d1a8abSmrg if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
169963d1a8abSmrg return result;
170063d1a8abSmrg } /* decFloatCompare */
170163d1a8abSmrg
170263d1a8abSmrg /* ------------------------------------------------------------------ */
170363d1a8abSmrg /* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
170463d1a8abSmrg /* */
170563d1a8abSmrg /* result gets the result of comparing dfl and dfr */
170663d1a8abSmrg /* dfl is the first decFloat (lhs) */
170763d1a8abSmrg /* dfr is the second decFloat (rhs) */
170863d1a8abSmrg /* set is the context */
170963d1a8abSmrg /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
171063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCompareSignal(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)171163d1a8abSmrg decFloat * decFloatCompareSignal(decFloat *result,
171263d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
171363d1a8abSmrg decContext *set) {
171463d1a8abSmrg Int comp; /* work */
171563d1a8abSmrg /* NaNs are handled as usual, except that all NaNs signal */
171663d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) {
171763d1a8abSmrg set->status|=DEC_Invalid_operation;
171863d1a8abSmrg return decNaNs(result, dfl, dfr, set);
171963d1a8abSmrg }
172063d1a8abSmrg /* numeric comparison needed */
172163d1a8abSmrg comp=decNumCompare(dfl, dfr, 0);
172263d1a8abSmrg decFloatZero(result);
172363d1a8abSmrg if (comp==0) return result;
172463d1a8abSmrg DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
172563d1a8abSmrg if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
172663d1a8abSmrg return result;
172763d1a8abSmrg } /* decFloatCompareSignal */
172863d1a8abSmrg
172963d1a8abSmrg /* ------------------------------------------------------------------ */
173063d1a8abSmrg /* decFloatCompareTotal -- compare two decFloats with total ordering */
173163d1a8abSmrg /* */
173263d1a8abSmrg /* result gets the result of comparing dfl and dfr */
173363d1a8abSmrg /* dfl is the first decFloat (lhs) */
173463d1a8abSmrg /* dfr is the second decFloat (rhs) */
173563d1a8abSmrg /* returns result, which may be -1, 0, or 1 */
173663d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCompareTotal(decFloat * result,const decFloat * dfl,const decFloat * dfr)173763d1a8abSmrg decFloat * decFloatCompareTotal(decFloat *result,
173863d1a8abSmrg const decFloat *dfl, const decFloat *dfr) {
173963d1a8abSmrg Int comp; /* work */
174063d1a8abSmrg uInt uiwork; /* for macros */
174163d1a8abSmrg #if QUAD
174263d1a8abSmrg uShort uswork; /* .. */
174363d1a8abSmrg #endif
174463d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) {
174563d1a8abSmrg Int nanl, nanr; /* work */
174663d1a8abSmrg /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
174763d1a8abSmrg nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
174863d1a8abSmrg if (DFISSIGNED(dfl)) nanl=-nanl;
174963d1a8abSmrg nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
175063d1a8abSmrg if (DFISSIGNED(dfr)) nanr=-nanr;
175163d1a8abSmrg if (nanl>nanr) comp=+1;
175263d1a8abSmrg else if (nanl<nanr) comp=-1;
175363d1a8abSmrg else { /* NaNs are the same type and sign .. must compare payload */
175463d1a8abSmrg /* buffers need +2 for QUAD */
175563d1a8abSmrg uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */
175663d1a8abSmrg uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */
175763d1a8abSmrg uByte *ub, *uc; /* work */
175863d1a8abSmrg Int sigl; /* signum of LHS */
175963d1a8abSmrg sigl=(DFISSIGNED(dfl) ? -1 : +1);
176063d1a8abSmrg
176163d1a8abSmrg /* decode the coefficients */
176263d1a8abSmrg /* (shift both right two if Quad to make a multiple of four) */
176363d1a8abSmrg #if QUAD
176463d1a8abSmrg UBFROMUS(bufl, 0);
176563d1a8abSmrg UBFROMUS(bufr, 0);
176663d1a8abSmrg #endif
176763d1a8abSmrg GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
176863d1a8abSmrg GETCOEFF(dfr, bufr+QUAD*2); /* .. */
176963d1a8abSmrg /* all multiples of four, here */
177063d1a8abSmrg comp=0; /* assume equal */
177163d1a8abSmrg for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
177263d1a8abSmrg uInt ui=UBTOUI(ub);
177363d1a8abSmrg if (ui==UBTOUI(uc)) continue; /* so far so same */
177463d1a8abSmrg /* about to find a winner; go by bytes in case little-endian */
177563d1a8abSmrg for (;; ub++, uc++) {
177663d1a8abSmrg if (*ub==*uc) continue;
177763d1a8abSmrg if (*ub>*uc) comp=sigl; /* difference found */
177863d1a8abSmrg else comp=-sigl; /* .. */
177963d1a8abSmrg break;
178063d1a8abSmrg }
178163d1a8abSmrg }
178263d1a8abSmrg } /* same NaN type and sign */
178363d1a8abSmrg }
178463d1a8abSmrg else {
178563d1a8abSmrg /* numeric comparison needed */
178663d1a8abSmrg comp=decNumCompare(dfl, dfr, 1); /* total ordering */
178763d1a8abSmrg }
178863d1a8abSmrg decFloatZero(result);
178963d1a8abSmrg if (comp==0) return result;
179063d1a8abSmrg DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
179163d1a8abSmrg if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
179263d1a8abSmrg return result;
179363d1a8abSmrg } /* decFloatCompareTotal */
179463d1a8abSmrg
179563d1a8abSmrg /* ------------------------------------------------------------------ */
179663d1a8abSmrg /* decFloatCompareTotalMag -- compare magnitudes with total ordering */
179763d1a8abSmrg /* */
179863d1a8abSmrg /* result gets the result of comparing abs(dfl) and abs(dfr) */
179963d1a8abSmrg /* dfl is the first decFloat (lhs) */
180063d1a8abSmrg /* dfr is the second decFloat (rhs) */
180163d1a8abSmrg /* returns result, which may be -1, 0, or 1 */
180263d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCompareTotalMag(decFloat * result,const decFloat * dfl,const decFloat * dfr)180363d1a8abSmrg decFloat * decFloatCompareTotalMag(decFloat *result,
180463d1a8abSmrg const decFloat *dfl, const decFloat *dfr) {
180563d1a8abSmrg decFloat a, b; /* for copy if needed */
180663d1a8abSmrg /* copy and redirect signed operand(s) */
180763d1a8abSmrg if (DFISSIGNED(dfl)) {
180863d1a8abSmrg decFloatCopyAbs(&a, dfl);
180963d1a8abSmrg dfl=&a;
181063d1a8abSmrg }
181163d1a8abSmrg if (DFISSIGNED(dfr)) {
181263d1a8abSmrg decFloatCopyAbs(&b, dfr);
181363d1a8abSmrg dfr=&b;
181463d1a8abSmrg }
181563d1a8abSmrg return decFloatCompareTotal(result, dfl, dfr);
181663d1a8abSmrg } /* decFloatCompareTotalMag */
181763d1a8abSmrg
181863d1a8abSmrg /* ------------------------------------------------------------------ */
181963d1a8abSmrg /* decFloatCopy -- copy a decFloat as-is */
182063d1a8abSmrg /* */
182163d1a8abSmrg /* result gets the copy of dfl */
182263d1a8abSmrg /* dfl is the decFloat to copy */
182363d1a8abSmrg /* returns result */
182463d1a8abSmrg /* */
182563d1a8abSmrg /* This is a bitwise operation; no errors or exceptions are possible. */
182663d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCopy(decFloat * result,const decFloat * dfl)182763d1a8abSmrg decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
182863d1a8abSmrg if (dfl!=result) *result=*dfl; /* copy needed */
182963d1a8abSmrg return result;
183063d1a8abSmrg } /* decFloatCopy */
183163d1a8abSmrg
183263d1a8abSmrg /* ------------------------------------------------------------------ */
183363d1a8abSmrg /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
183463d1a8abSmrg /* */
183563d1a8abSmrg /* result gets the copy of dfl with sign bit 0 */
183663d1a8abSmrg /* dfl is the decFloat to copy */
183763d1a8abSmrg /* returns result */
183863d1a8abSmrg /* */
183963d1a8abSmrg /* This is a bitwise operation; no errors or exceptions are possible. */
184063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCopyAbs(decFloat * result,const decFloat * dfl)184163d1a8abSmrg decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
184263d1a8abSmrg if (dfl!=result) *result=*dfl; /* copy needed */
184363d1a8abSmrg DFBYTE(result, 0)&=~0x80; /* zero sign bit */
184463d1a8abSmrg return result;
184563d1a8abSmrg } /* decFloatCopyAbs */
184663d1a8abSmrg
184763d1a8abSmrg /* ------------------------------------------------------------------ */
184863d1a8abSmrg /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
184963d1a8abSmrg /* */
185063d1a8abSmrg /* result gets the copy of dfl with sign bit inverted */
185163d1a8abSmrg /* dfl is the decFloat to copy */
185263d1a8abSmrg /* returns result */
185363d1a8abSmrg /* */
185463d1a8abSmrg /* This is a bitwise operation; no errors or exceptions are possible. */
185563d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCopyNegate(decFloat * result,const decFloat * dfl)185663d1a8abSmrg decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
185763d1a8abSmrg if (dfl!=result) *result=*dfl; /* copy needed */
185863d1a8abSmrg DFBYTE(result, 0)^=0x80; /* invert sign bit */
185963d1a8abSmrg return result;
186063d1a8abSmrg } /* decFloatCopyNegate */
186163d1a8abSmrg
186263d1a8abSmrg /* ------------------------------------------------------------------ */
186363d1a8abSmrg /* decFloatCopySign -- copy a decFloat with the sign of another */
186463d1a8abSmrg /* */
186563d1a8abSmrg /* result gets the result of copying dfl with the sign of dfr */
186663d1a8abSmrg /* dfl is the first decFloat (lhs) */
186763d1a8abSmrg /* dfr is the second decFloat (rhs) */
186863d1a8abSmrg /* returns result */
186963d1a8abSmrg /* */
187063d1a8abSmrg /* This is a bitwise operation; no errors or exceptions are possible. */
187163d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatCopySign(decFloat * result,const decFloat * dfl,const decFloat * dfr)187263d1a8abSmrg decFloat * decFloatCopySign(decFloat *result,
187363d1a8abSmrg const decFloat *dfl, const decFloat *dfr) {
187463d1a8abSmrg uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
187563d1a8abSmrg if (dfl!=result) *result=*dfl; /* copy needed */
187663d1a8abSmrg DFBYTE(result, 0)&=~0x80; /* clear sign .. */
187763d1a8abSmrg DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
187863d1a8abSmrg return result;
187963d1a8abSmrg } /* decFloatCopySign */
188063d1a8abSmrg
188163d1a8abSmrg /* ------------------------------------------------------------------ */
188263d1a8abSmrg /* decFloatDigits -- return the number of digits in a decFloat */
188363d1a8abSmrg /* */
188463d1a8abSmrg /* df is the decFloat to investigate */
188563d1a8abSmrg /* returns the number of significant digits in the decFloat; a */
188663d1a8abSmrg /* zero coefficient returns 1 as does an infinity (a NaN returns */
188763d1a8abSmrg /* the number of digits in the payload) */
188863d1a8abSmrg /* ------------------------------------------------------------------ */
188963d1a8abSmrg /* private macro to extract a declet according to provided formula */
189063d1a8abSmrg /* (form), and if it is non-zero then return the calculated digits */
189163d1a8abSmrg /* depending on the declet number (n), where n=0 for the most */
189263d1a8abSmrg /* significant declet; uses uInt dpd for work */
189363d1a8abSmrg #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
189463d1a8abSmrg if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
189563d1a8abSmrg /* next one is used when it is known that the declet must be */
189663d1a8abSmrg /* non-zero, or is the final zero declet */
189763d1a8abSmrg #define dpdlendun(n, form) {dpd=(form)&0x3ff; \
189863d1a8abSmrg if (dpd==0) return 1; \
189963d1a8abSmrg return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
190063d1a8abSmrg
decFloatDigits(const decFloat * df)190163d1a8abSmrg uInt decFloatDigits(const decFloat *df) {
190263d1a8abSmrg uInt dpd; /* work */
190363d1a8abSmrg uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
190463d1a8abSmrg #if QUAD
190563d1a8abSmrg uInt sourmh, sourml;
190663d1a8abSmrg #endif
190763d1a8abSmrg uInt sourlo;
190863d1a8abSmrg
190963d1a8abSmrg if (DFISINF(df)) return 1;
191063d1a8abSmrg /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
191163d1a8abSmrg /* then the coefficient is full-length */
191263d1a8abSmrg if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
191363d1a8abSmrg
191463d1a8abSmrg #if DOUBLE
191563d1a8abSmrg if (sourhi&0x0003ffff) { /* ends in first */
191663d1a8abSmrg dpdlenchk(0, sourhi>>8);
191763d1a8abSmrg sourlo=DFWORD(df, 1);
191863d1a8abSmrg dpdlendun(1, (sourhi<<2) | (sourlo>>30));
191963d1a8abSmrg } /* [cannot drop through] */
192063d1a8abSmrg sourlo=DFWORD(df, 1); /* sourhi not involved now */
192163d1a8abSmrg if (sourlo&0xfff00000) { /* in one of first two */
192263d1a8abSmrg dpdlenchk(1, sourlo>>30); /* very rare */
192363d1a8abSmrg dpdlendun(2, sourlo>>20);
192463d1a8abSmrg } /* [cannot drop through] */
192563d1a8abSmrg dpdlenchk(3, sourlo>>10);
192663d1a8abSmrg dpdlendun(4, sourlo);
192763d1a8abSmrg /* [cannot drop through] */
192863d1a8abSmrg
192963d1a8abSmrg #elif QUAD
193063d1a8abSmrg if (sourhi&0x00003fff) { /* ends in first */
193163d1a8abSmrg dpdlenchk(0, sourhi>>4);
193263d1a8abSmrg sourmh=DFWORD(df, 1);
193363d1a8abSmrg dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
193463d1a8abSmrg } /* [cannot drop through] */
193563d1a8abSmrg sourmh=DFWORD(df, 1);
193663d1a8abSmrg if (sourmh) {
193763d1a8abSmrg dpdlenchk(1, sourmh>>26);
193863d1a8abSmrg dpdlenchk(2, sourmh>>16);
193963d1a8abSmrg dpdlenchk(3, sourmh>>6);
194063d1a8abSmrg sourml=DFWORD(df, 2);
194163d1a8abSmrg dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
194263d1a8abSmrg } /* [cannot drop through] */
194363d1a8abSmrg sourml=DFWORD(df, 2);
194463d1a8abSmrg if (sourml) {
194563d1a8abSmrg dpdlenchk(4, sourml>>28);
194663d1a8abSmrg dpdlenchk(5, sourml>>18);
194763d1a8abSmrg dpdlenchk(6, sourml>>8);
194863d1a8abSmrg sourlo=DFWORD(df, 3);
194963d1a8abSmrg dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
195063d1a8abSmrg } /* [cannot drop through] */
195163d1a8abSmrg sourlo=DFWORD(df, 3);
195263d1a8abSmrg if (sourlo&0xfff00000) { /* in one of first two */
195363d1a8abSmrg dpdlenchk(7, sourlo>>30); /* very rare */
195463d1a8abSmrg dpdlendun(8, sourlo>>20);
195563d1a8abSmrg } /* [cannot drop through] */
195663d1a8abSmrg dpdlenchk(9, sourlo>>10);
195763d1a8abSmrg dpdlendun(10, sourlo);
195863d1a8abSmrg /* [cannot drop through] */
195963d1a8abSmrg #endif
196063d1a8abSmrg } /* decFloatDigits */
196163d1a8abSmrg
196263d1a8abSmrg /* ------------------------------------------------------------------ */
196363d1a8abSmrg /* decFloatDivide -- divide a decFloat by another */
196463d1a8abSmrg /* */
196563d1a8abSmrg /* result gets the result of dividing dfl by dfr: */
196663d1a8abSmrg /* dfl is the first decFloat (lhs) */
196763d1a8abSmrg /* dfr is the second decFloat (rhs) */
196863d1a8abSmrg /* set is the context */
196963d1a8abSmrg /* returns result */
197063d1a8abSmrg /* */
197163d1a8abSmrg /* ------------------------------------------------------------------ */
197263d1a8abSmrg /* This is just a wrapper. */
decFloatDivide(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)197363d1a8abSmrg decFloat * decFloatDivide(decFloat *result,
197463d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
197563d1a8abSmrg decContext *set) {
197663d1a8abSmrg return decDivide(result, dfl, dfr, set, DIVIDE);
197763d1a8abSmrg } /* decFloatDivide */
197863d1a8abSmrg
197963d1a8abSmrg /* ------------------------------------------------------------------ */
198063d1a8abSmrg /* decFloatDivideInteger -- integer divide a decFloat by another */
198163d1a8abSmrg /* */
198263d1a8abSmrg /* result gets the result of dividing dfl by dfr: */
198363d1a8abSmrg /* dfl is the first decFloat (lhs) */
198463d1a8abSmrg /* dfr is the second decFloat (rhs) */
198563d1a8abSmrg /* set is the context */
198663d1a8abSmrg /* returns result */
198763d1a8abSmrg /* */
198863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatDivideInteger(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)198963d1a8abSmrg decFloat * decFloatDivideInteger(decFloat *result,
199063d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
199163d1a8abSmrg decContext *set) {
199263d1a8abSmrg return decDivide(result, dfl, dfr, set, DIVIDEINT);
199363d1a8abSmrg } /* decFloatDivideInteger */
199463d1a8abSmrg
199563d1a8abSmrg /* ------------------------------------------------------------------ */
199663d1a8abSmrg /* decFloatFMA -- multiply and add three decFloats, fused */
199763d1a8abSmrg /* */
199863d1a8abSmrg /* result gets the result of (dfl*dfr)+dff with a single rounding */
199963d1a8abSmrg /* dfl is the first decFloat (lhs) */
200063d1a8abSmrg /* dfr is the second decFloat (rhs) */
200163d1a8abSmrg /* dff is the final decFloat (fhs) */
200263d1a8abSmrg /* set is the context */
200363d1a8abSmrg /* returns result */
200463d1a8abSmrg /* */
200563d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatFMA(decFloat * result,const decFloat * dfl,const decFloat * dfr,const decFloat * dff,decContext * set)200663d1a8abSmrg decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
200763d1a8abSmrg const decFloat *dfr, const decFloat *dff,
200863d1a8abSmrg decContext *set) {
200963d1a8abSmrg
201063d1a8abSmrg /* The accumulator has the bytes needed for FiniteMultiply, plus */
201163d1a8abSmrg /* one byte to the left in case of carry, plus DECPMAX+2 to the */
201263d1a8abSmrg /* right for the final addition (up to full fhs + round & sticky) */
201363d1a8abSmrg #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
201463d1a8abSmrg uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
201563d1a8abSmrg /* .. and for final result */
201663d1a8abSmrg bcdnum mul; /* for multiplication result */
201763d1a8abSmrg bcdnum fin; /* for final operand, expanded */
201863d1a8abSmrg uByte coe[ROUNDUP4(DECPMAX)]; /* dff coefficient in BCD */
201963d1a8abSmrg bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
202063d1a8abSmrg uInt diffsign; /* non-zero if signs differ */
202163d1a8abSmrg uInt hipad; /* pad digit for hi if needed */
202263d1a8abSmrg Int padding; /* excess exponent */
202363d1a8abSmrg uInt carry; /* +1 for ten's complement and during add */
202463d1a8abSmrg uByte *ub, *uh, *ul; /* work */
202563d1a8abSmrg uInt uiwork; /* for macros */
202663d1a8abSmrg
202763d1a8abSmrg /* handle all the special values [any special operand leads to a */
202863d1a8abSmrg /* special result] */
202963d1a8abSmrg if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
203063d1a8abSmrg decFloat proxy; /* multiplication result proxy */
203163d1a8abSmrg /* NaNs are handled as usual, giving priority to sNaNs */
203263d1a8abSmrg if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
203363d1a8abSmrg if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
203463d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
203563d1a8abSmrg if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
203663d1a8abSmrg /* One or more of the three is infinite */
203763d1a8abSmrg /* infinity times zero is bad */
203863d1a8abSmrg decFloatZero(&proxy);
203963d1a8abSmrg if (DFISINF(dfl)) {
204063d1a8abSmrg if (DFISZERO(dfr)) return decInvalid(result, set);
204163d1a8abSmrg decInfinity(&proxy, &proxy);
204263d1a8abSmrg }
204363d1a8abSmrg else if (DFISINF(dfr)) {
204463d1a8abSmrg if (DFISZERO(dfl)) return decInvalid(result, set);
204563d1a8abSmrg decInfinity(&proxy, &proxy);
204663d1a8abSmrg }
204763d1a8abSmrg /* compute sign of multiplication and place in proxy */
204863d1a8abSmrg DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
204963d1a8abSmrg if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
205063d1a8abSmrg /* dff is Infinite */
205163d1a8abSmrg if (!DFISINF(&proxy)) return decInfinity(result, dff);
205263d1a8abSmrg /* both sides of addition are infinite; different sign is bad */
205363d1a8abSmrg if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
205463d1a8abSmrg return decInvalid(result, set);
205563d1a8abSmrg return decFloatCopy(result, &proxy);
205663d1a8abSmrg }
205763d1a8abSmrg
205863d1a8abSmrg /* Here when all operands are finite */
205963d1a8abSmrg
206063d1a8abSmrg /* First multiply dfl*dfr */
206163d1a8abSmrg decFiniteMultiply(&mul, acc+1, dfl, dfr);
206263d1a8abSmrg /* The multiply is complete, exact and unbounded, and described in */
206363d1a8abSmrg /* mul with the coefficient held in acc[1...] */
206463d1a8abSmrg
206563d1a8abSmrg /* now add in dff; the algorithm is essentially the same as */
206663d1a8abSmrg /* decFloatAdd, but the code is different because the code there */
206763d1a8abSmrg /* is highly optimized for adding two numbers of the same size */
206863d1a8abSmrg fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
206963d1a8abSmrg fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
207063d1a8abSmrg diffsign=mul.sign^fin.sign; /* note if signs differ */
207163d1a8abSmrg fin.msd=coe;
207263d1a8abSmrg fin.lsd=coe+DECPMAX-1;
207363d1a8abSmrg GETCOEFF(dff, coe); /* extract the coefficient */
207463d1a8abSmrg
207563d1a8abSmrg /* now set hi and lo so that hi points to whichever of mul and fin */
207663d1a8abSmrg /* has the higher exponent and lo points to the other [don't care, */
207763d1a8abSmrg /* if the same]. One coefficient will be in acc, the other in coe. */
207863d1a8abSmrg if (mul.exponent>=fin.exponent) {
207963d1a8abSmrg hi=&mul;
208063d1a8abSmrg lo=&fin;
208163d1a8abSmrg }
208263d1a8abSmrg else {
208363d1a8abSmrg hi=&fin;
208463d1a8abSmrg lo=&mul;
208563d1a8abSmrg }
208663d1a8abSmrg
208763d1a8abSmrg /* remove leading zeros on both operands; this will save time later */
208863d1a8abSmrg /* and make testing for zero trivial (tests are safe because acc */
208963d1a8abSmrg /* and coe are rounded up to uInts) */
209063d1a8abSmrg for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
209163d1a8abSmrg for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
209263d1a8abSmrg for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
209363d1a8abSmrg for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
209463d1a8abSmrg
209563d1a8abSmrg /* if hi is zero then result will be lo (which has the smaller */
209663d1a8abSmrg /* exponent), which also may need to be tested for zero for the */
209763d1a8abSmrg /* weird IEEE 754 sign rules */
209863d1a8abSmrg if (*hi->msd==0) { /* hi is zero */
209963d1a8abSmrg /* "When the sum of two operands with opposite signs is */
210063d1a8abSmrg /* exactly zero, the sign of that sum shall be '+' in all */
210163d1a8abSmrg /* rounding modes except round toward -Infinity, in which */
210263d1a8abSmrg /* mode that sign shall be '-'." */
210363d1a8abSmrg if (diffsign) {
210463d1a8abSmrg if (*lo->msd==0) { /* lo is zero */
210563d1a8abSmrg lo->sign=0;
210663d1a8abSmrg if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
210763d1a8abSmrg } /* diffsign && lo=0 */
210863d1a8abSmrg } /* diffsign */
210963d1a8abSmrg return decFinalize(result, lo, set); /* may need clamping */
211063d1a8abSmrg } /* numfl is zero */
211163d1a8abSmrg /* [here, both are minimal length and hi is non-zero] */
211263d1a8abSmrg /* (if lo is zero then padding with zeros may be needed, below) */
211363d1a8abSmrg
211463d1a8abSmrg /* if signs differ, take the ten's complement of hi (zeros to the */
211563d1a8abSmrg /* right do not matter because the complement of zero is zero); the */
211663d1a8abSmrg /* +1 is done later, as part of the addition, inserted at the */
211763d1a8abSmrg /* correct digit */
211863d1a8abSmrg hipad=0;
211963d1a8abSmrg carry=0;
212063d1a8abSmrg if (diffsign) {
212163d1a8abSmrg hipad=9;
212263d1a8abSmrg carry=1;
212363d1a8abSmrg /* exactly the correct number of digits must be inverted */
212463d1a8abSmrg for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
212563d1a8abSmrg for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
212663d1a8abSmrg }
212763d1a8abSmrg
212863d1a8abSmrg /* ready to add; note that hi has no leading zeros so gap */
212963d1a8abSmrg /* calculation does not have to be as pessimistic as in decFloatAdd */
213063d1a8abSmrg /* (this is much more like the arbitrary-precision algorithm in */
213163d1a8abSmrg /* Rexx and decNumber) */
213263d1a8abSmrg
213363d1a8abSmrg /* padding is the number of zeros that would need to be added to hi */
213463d1a8abSmrg /* for its lsd to be aligned with the lsd of lo */
213563d1a8abSmrg padding=hi->exponent-lo->exponent;
213663d1a8abSmrg /* printf("FMA pad %ld\n", (LI)padding); */
213763d1a8abSmrg
213863d1a8abSmrg /* the result of the addition will be built into the accumulator, */
213963d1a8abSmrg /* starting from the far right; this could be either hi or lo, and */
214063d1a8abSmrg /* will be aligned */
214163d1a8abSmrg ub=acc+FMALEN-1; /* where lsd of result will go */
214263d1a8abSmrg ul=lo->lsd; /* lsd of rhs */
214363d1a8abSmrg
214463d1a8abSmrg if (padding!=0) { /* unaligned */
214563d1a8abSmrg /* if the msd of lo is more than DECPMAX+2 digits to the right of */
214663d1a8abSmrg /* the original msd of hi then it can be reduced to a single */
214763d1a8abSmrg /* digit at the right place, as it stays clear of hi digits */
214863d1a8abSmrg /* [it must be DECPMAX+2 because during a subtraction the msd */
214963d1a8abSmrg /* could become 0 after a borrow from 1.000 to 0.9999...] */
215063d1a8abSmrg
215163d1a8abSmrg Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
215263d1a8abSmrg Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
215363d1a8abSmrg
215463d1a8abSmrg if (hilen+padding-lolen > DECPMAX+2) { /* can reduce lo to single */
215563d1a8abSmrg /* make sure it is virtually at least DECPMAX from hi->msd, at */
215663d1a8abSmrg /* least to right of hi->lsd (in case of destructive subtract), */
215763d1a8abSmrg /* and separated by at least two digits from either of those */
215863d1a8abSmrg /* (the tricky DOUBLE case is when hi is a 1 that will become a */
215963d1a8abSmrg /* 0.9999... by subtraction: */
216063d1a8abSmrg /* hi: 1 E+16 */
216163d1a8abSmrg /* lo: .................1000000000000000 E-16 */
216263d1a8abSmrg /* which for the addition pads to: */
216363d1a8abSmrg /* hi: 1000000000000000000 E-16 */
216463d1a8abSmrg /* lo: .................1000000000000000 E-16 */
216563d1a8abSmrg Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
216663d1a8abSmrg
216763d1a8abSmrg /* printf("FMA reduce: %ld\n", (LI)reduce); */
216863d1a8abSmrg lo->lsd=lo->msd; /* to single digit [maybe 0] */
216963d1a8abSmrg lo->exponent=newexp; /* new lowest exponent */
217063d1a8abSmrg padding=hi->exponent-lo->exponent; /* recalculate */
217163d1a8abSmrg ul=lo->lsd; /* .. and repoint */
217263d1a8abSmrg }
217363d1a8abSmrg
217463d1a8abSmrg /* padding is still > 0, but will fit in acc (less leading carry slot) */
217563d1a8abSmrg #if DECCHECK
217663d1a8abSmrg if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
217763d1a8abSmrg if (hilen+padding+1>FMALEN)
217863d1a8abSmrg printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding);
217963d1a8abSmrg /* printf("FMA padding: %ld\n", (LI)padding); */
218063d1a8abSmrg #endif
218163d1a8abSmrg
218263d1a8abSmrg /* padding digits can now be set in the result; one or more of */
218363d1a8abSmrg /* these will come from lo; others will be zeros in the gap */
218463d1a8abSmrg for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
218563d1a8abSmrg UBFROMUI(ub-3, UBTOUI(ul-3)); /* [cannot overlap] */
218663d1a8abSmrg }
218763d1a8abSmrg for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
218863d1a8abSmrg for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
218963d1a8abSmrg }
219063d1a8abSmrg
219163d1a8abSmrg /* addition now complete to the right of the rightmost digit of hi */
219263d1a8abSmrg uh=hi->lsd;
219363d1a8abSmrg
219463d1a8abSmrg /* dow do the add from hi->lsd to the left */
219563d1a8abSmrg /* [bytewise, because either operand can run out at any time] */
219663d1a8abSmrg /* carry was set up depending on ten's complement above */
219763d1a8abSmrg /* first assume both operands have some digits */
219863d1a8abSmrg for (;; ub--) {
219963d1a8abSmrg if (uh<hi->msd || ul<lo->msd) break;
220063d1a8abSmrg *ub=(uByte)(carry+(*uh--)+(*ul--));
220163d1a8abSmrg carry=0;
220263d1a8abSmrg if (*ub<10) continue;
220363d1a8abSmrg *ub-=10;
220463d1a8abSmrg carry=1;
220563d1a8abSmrg } /* both loop */
220663d1a8abSmrg
220763d1a8abSmrg if (ul<lo->msd) { /* to left of lo */
220863d1a8abSmrg for (;; ub--) {
220963d1a8abSmrg if (uh<hi->msd) break;
221063d1a8abSmrg *ub=(uByte)(carry+(*uh--)); /* [+0] */
221163d1a8abSmrg carry=0;
221263d1a8abSmrg if (*ub<10) continue;
221363d1a8abSmrg *ub-=10;
221463d1a8abSmrg carry=1;
221563d1a8abSmrg } /* hi loop */
221663d1a8abSmrg }
221763d1a8abSmrg else { /* to left of hi */
221863d1a8abSmrg for (;; ub--) {
221963d1a8abSmrg if (ul<lo->msd) break;
222063d1a8abSmrg *ub=(uByte)(carry+hipad+(*ul--));
222163d1a8abSmrg carry=0;
222263d1a8abSmrg if (*ub<10) continue;
222363d1a8abSmrg *ub-=10;
222463d1a8abSmrg carry=1;
222563d1a8abSmrg } /* lo loop */
222663d1a8abSmrg }
222763d1a8abSmrg
222863d1a8abSmrg /* addition complete -- now handle carry, borrow, etc. */
222963d1a8abSmrg /* use lo to set up the num (its exponent is already correct, and */
223063d1a8abSmrg /* sign usually is) */
223163d1a8abSmrg lo->msd=ub+1;
223263d1a8abSmrg lo->lsd=acc+FMALEN-1;
223363d1a8abSmrg /* decShowNum(lo, "lo"); */
223463d1a8abSmrg if (!diffsign) { /* same-sign addition */
223563d1a8abSmrg if (carry) { /* carry out */
223663d1a8abSmrg *ub=1; /* place the 1 .. */
223763d1a8abSmrg lo->msd--; /* .. and update */
223863d1a8abSmrg }
223963d1a8abSmrg } /* same sign */
224063d1a8abSmrg else { /* signs differed (subtraction) */
224163d1a8abSmrg if (!carry) { /* no carry out means hi<lo */
224263d1a8abSmrg /* borrowed -- take ten's complement of the right digits */
224363d1a8abSmrg lo->sign=hi->sign; /* sign is lhs sign */
224463d1a8abSmrg for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul));
224563d1a8abSmrg for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
224663d1a8abSmrg /* complete the ten's complement by adding 1 [cannot overrun] */
224763d1a8abSmrg for (ul--; *ul==9; ul--) *ul=0;
224863d1a8abSmrg *ul+=1;
224963d1a8abSmrg } /* borrowed */
225063d1a8abSmrg else { /* carry out means hi>=lo */
225163d1a8abSmrg /* sign to use is lo->sign */
225263d1a8abSmrg /* all done except for the special IEEE 754 exact-zero-result */
225363d1a8abSmrg /* rule (see above); while testing for zero, strip leading */
225463d1a8abSmrg /* zeros (which will save decFinalize doing it) */
225563d1a8abSmrg for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
225663d1a8abSmrg for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
225763d1a8abSmrg if (*lo->msd==0) { /* must be true zero (and diffsign) */
225863d1a8abSmrg lo->sign=0; /* assume + */
225963d1a8abSmrg if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
226063d1a8abSmrg }
226163d1a8abSmrg /* [else was not zero, might still have leading zeros] */
226263d1a8abSmrg } /* subtraction gave positive result */
226363d1a8abSmrg } /* diffsign */
226463d1a8abSmrg
226563d1a8abSmrg #if DECCHECK
226663d1a8abSmrg /* assert no left underrun */
226763d1a8abSmrg if (lo->msd<acc) {
226863d1a8abSmrg printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
226963d1a8abSmrg }
227063d1a8abSmrg #endif
227163d1a8abSmrg
227263d1a8abSmrg return decFinalize(result, lo, set); /* round, check, and lay out */
227363d1a8abSmrg } /* decFloatFMA */
227463d1a8abSmrg
227563d1a8abSmrg /* ------------------------------------------------------------------ */
227663d1a8abSmrg /* decFloatFromInt -- initialise a decFloat from an Int */
227763d1a8abSmrg /* */
227863d1a8abSmrg /* result gets the converted Int */
227963d1a8abSmrg /* n is the Int to convert */
228063d1a8abSmrg /* returns result */
228163d1a8abSmrg /* */
228263d1a8abSmrg /* The result is Exact; no errors or exceptions are possible. */
228363d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatFromInt32(decFloat * result,Int n)228463d1a8abSmrg decFloat * decFloatFromInt32(decFloat *result, Int n) {
228563d1a8abSmrg uInt u=(uInt)n; /* copy as bits */
228663d1a8abSmrg uInt encode; /* work */
228763d1a8abSmrg DFWORD(result, 0)=ZEROWORD; /* always */
228863d1a8abSmrg #if QUAD
228963d1a8abSmrg DFWORD(result, 1)=0;
229063d1a8abSmrg DFWORD(result, 2)=0;
229163d1a8abSmrg #endif
229263d1a8abSmrg if (n<0) { /* handle -n with care */
229363d1a8abSmrg /* [This can be done without the test, but is then slightly slower] */
229463d1a8abSmrg u=(~u)+1;
229563d1a8abSmrg DFWORD(result, 0)|=DECFLOAT_Sign;
229663d1a8abSmrg }
229763d1a8abSmrg /* Since the maximum value of u now is 2**31, only the low word of */
229863d1a8abSmrg /* result is affected */
229963d1a8abSmrg encode=BIN2DPD[u%1000];
230063d1a8abSmrg u/=1000;
230163d1a8abSmrg encode|=BIN2DPD[u%1000]<<10;
230263d1a8abSmrg u/=1000;
230363d1a8abSmrg encode|=BIN2DPD[u%1000]<<20;
230463d1a8abSmrg u/=1000; /* now 0, 1, or 2 */
230563d1a8abSmrg encode|=u<<30;
230663d1a8abSmrg DFWORD(result, DECWORDS-1)=encode;
230763d1a8abSmrg return result;
230863d1a8abSmrg } /* decFloatFromInt32 */
230963d1a8abSmrg
231063d1a8abSmrg /* ------------------------------------------------------------------ */
231163d1a8abSmrg /* decFloatFromUInt -- initialise a decFloat from a uInt */
231263d1a8abSmrg /* */
231363d1a8abSmrg /* result gets the converted uInt */
231463d1a8abSmrg /* n is the uInt to convert */
231563d1a8abSmrg /* returns result */
231663d1a8abSmrg /* */
231763d1a8abSmrg /* The result is Exact; no errors or exceptions are possible. */
231863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatFromUInt32(decFloat * result,uInt u)231963d1a8abSmrg decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
232063d1a8abSmrg uInt encode; /* work */
232163d1a8abSmrg DFWORD(result, 0)=ZEROWORD; /* always */
232263d1a8abSmrg #if QUAD
232363d1a8abSmrg DFWORD(result, 1)=0;
232463d1a8abSmrg DFWORD(result, 2)=0;
232563d1a8abSmrg #endif
232663d1a8abSmrg encode=BIN2DPD[u%1000];
232763d1a8abSmrg u/=1000;
232863d1a8abSmrg encode|=BIN2DPD[u%1000]<<10;
232963d1a8abSmrg u/=1000;
233063d1a8abSmrg encode|=BIN2DPD[u%1000]<<20;
233163d1a8abSmrg u/=1000; /* now 0 -> 4 */
233263d1a8abSmrg encode|=u<<30;
233363d1a8abSmrg DFWORD(result, DECWORDS-1)=encode;
233463d1a8abSmrg DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
233563d1a8abSmrg return result;
233663d1a8abSmrg } /* decFloatFromUInt32 */
233763d1a8abSmrg
233863d1a8abSmrg /* ------------------------------------------------------------------ */
233963d1a8abSmrg /* decFloatInvert -- logical digitwise INVERT of a decFloat */
234063d1a8abSmrg /* */
234163d1a8abSmrg /* result gets the result of INVERTing df */
234263d1a8abSmrg /* df is the decFloat to invert */
234363d1a8abSmrg /* set is the context */
234463d1a8abSmrg /* returns result, which will be canonical with sign=0 */
234563d1a8abSmrg /* */
234663d1a8abSmrg /* The operand must be positive, finite with exponent q=0, and */
234763d1a8abSmrg /* comprise just zeros and ones; if not, Invalid operation results. */
234863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatInvert(decFloat * result,const decFloat * df,decContext * set)234963d1a8abSmrg decFloat * decFloatInvert(decFloat *result, const decFloat *df,
235063d1a8abSmrg decContext *set) {
235163d1a8abSmrg uInt sourhi=DFWORD(df, 0); /* top word of dfs */
235263d1a8abSmrg
235363d1a8abSmrg if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
235463d1a8abSmrg /* the operand is a finite integer (q=0) */
235563d1a8abSmrg #if DOUBLE
235663d1a8abSmrg DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
235763d1a8abSmrg DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
235863d1a8abSmrg #elif QUAD
235963d1a8abSmrg DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
236063d1a8abSmrg DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
236163d1a8abSmrg DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
236263d1a8abSmrg DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
236363d1a8abSmrg #endif
236463d1a8abSmrg return result;
236563d1a8abSmrg } /* decFloatInvert */
236663d1a8abSmrg
236763d1a8abSmrg /* ------------------------------------------------------------------ */
236863d1a8abSmrg /* decFloatIs -- decFloat tests (IsSigned, etc.) */
236963d1a8abSmrg /* */
237063d1a8abSmrg /* df is the decFloat to test */
237163d1a8abSmrg /* returns 0 or 1 in a uInt */
237263d1a8abSmrg /* */
237363d1a8abSmrg /* Many of these could be macros, but having them as real functions */
237463d1a8abSmrg /* is a little cleaner (and they can be referred to here by the */
237563d1a8abSmrg /* generic names) */
237663d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatIsCanonical(const decFloat * df)237763d1a8abSmrg uInt decFloatIsCanonical(const decFloat *df) {
237863d1a8abSmrg if (DFISSPECIAL(df)) {
237963d1a8abSmrg if (DFISINF(df)) {
238063d1a8abSmrg if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
238163d1a8abSmrg if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
238263d1a8abSmrg return 1;
238363d1a8abSmrg }
238463d1a8abSmrg /* is a NaN */
238563d1a8abSmrg if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
238663d1a8abSmrg if (DFISCCZERO(df)) return 1; /* coefficient continuation */
238763d1a8abSmrg /* drop through to check payload */
238863d1a8abSmrg }
238963d1a8abSmrg { /* declare block */
239063d1a8abSmrg #if DOUBLE
239163d1a8abSmrg uInt sourhi=DFWORD(df, 0);
239263d1a8abSmrg uInt sourlo=DFWORD(df, 1);
239363d1a8abSmrg if (CANONDPDOFF(sourhi, 8)
239463d1a8abSmrg && CANONDPDTWO(sourhi, sourlo, 30)
239563d1a8abSmrg && CANONDPDOFF(sourlo, 20)
239663d1a8abSmrg && CANONDPDOFF(sourlo, 10)
239763d1a8abSmrg && CANONDPDOFF(sourlo, 0)) return 1;
239863d1a8abSmrg #elif QUAD
239963d1a8abSmrg uInt sourhi=DFWORD(df, 0);
240063d1a8abSmrg uInt sourmh=DFWORD(df, 1);
240163d1a8abSmrg uInt sourml=DFWORD(df, 2);
240263d1a8abSmrg uInt sourlo=DFWORD(df, 3);
240363d1a8abSmrg if (CANONDPDOFF(sourhi, 4)
240463d1a8abSmrg && CANONDPDTWO(sourhi, sourmh, 26)
240563d1a8abSmrg && CANONDPDOFF(sourmh, 16)
240663d1a8abSmrg && CANONDPDOFF(sourmh, 6)
240763d1a8abSmrg && CANONDPDTWO(sourmh, sourml, 28)
240863d1a8abSmrg && CANONDPDOFF(sourml, 18)
240963d1a8abSmrg && CANONDPDOFF(sourml, 8)
241063d1a8abSmrg && CANONDPDTWO(sourml, sourlo, 30)
241163d1a8abSmrg && CANONDPDOFF(sourlo, 20)
241263d1a8abSmrg && CANONDPDOFF(sourlo, 10)
241363d1a8abSmrg && CANONDPDOFF(sourlo, 0)) return 1;
241463d1a8abSmrg #endif
241563d1a8abSmrg } /* block */
241663d1a8abSmrg return 0; /* a declet is non-canonical */
241763d1a8abSmrg }
241863d1a8abSmrg
decFloatIsFinite(const decFloat * df)241963d1a8abSmrg uInt decFloatIsFinite(const decFloat *df) {
242063d1a8abSmrg return !DFISSPECIAL(df);
242163d1a8abSmrg }
decFloatIsInfinite(const decFloat * df)242263d1a8abSmrg uInt decFloatIsInfinite(const decFloat *df) {
242363d1a8abSmrg return DFISINF(df);
242463d1a8abSmrg }
decFloatIsInteger(const decFloat * df)242563d1a8abSmrg uInt decFloatIsInteger(const decFloat *df) {
242663d1a8abSmrg return DFISINT(df);
242763d1a8abSmrg }
decFloatIsNaN(const decFloat * df)242863d1a8abSmrg uInt decFloatIsNaN(const decFloat *df) {
242963d1a8abSmrg return DFISNAN(df);
243063d1a8abSmrg }
decFloatIsNormal(const decFloat * df)243163d1a8abSmrg uInt decFloatIsNormal(const decFloat *df) {
243263d1a8abSmrg Int exp; /* exponent */
243363d1a8abSmrg if (DFISSPECIAL(df)) return 0;
243463d1a8abSmrg if (DFISZERO(df)) return 0;
243563d1a8abSmrg /* is finite and non-zero */
243663d1a8abSmrg exp=GETEXPUN(df) /* get unbiased exponent .. */
243763d1a8abSmrg +decFloatDigits(df)-1; /* .. and make adjusted exponent */
243863d1a8abSmrg return (exp>=DECEMIN); /* < DECEMIN is subnormal */
243963d1a8abSmrg }
decFloatIsSignaling(const decFloat * df)244063d1a8abSmrg uInt decFloatIsSignaling(const decFloat *df) {
244163d1a8abSmrg return DFISSNAN(df);
244263d1a8abSmrg }
decFloatIsSignalling(const decFloat * df)244363d1a8abSmrg uInt decFloatIsSignalling(const decFloat *df) {
244463d1a8abSmrg return DFISSNAN(df);
244563d1a8abSmrg }
decFloatIsSigned(const decFloat * df)244663d1a8abSmrg uInt decFloatIsSigned(const decFloat *df) {
244763d1a8abSmrg return DFISSIGNED(df);
244863d1a8abSmrg }
decFloatIsSubnormal(const decFloat * df)244963d1a8abSmrg uInt decFloatIsSubnormal(const decFloat *df) {
245063d1a8abSmrg if (DFISSPECIAL(df)) return 0;
245163d1a8abSmrg /* is finite */
245263d1a8abSmrg if (decFloatIsNormal(df)) return 0;
245363d1a8abSmrg /* it is <Nmin, but could be zero */
245463d1a8abSmrg if (DFISZERO(df)) return 0;
245563d1a8abSmrg return 1; /* is subnormal */
245663d1a8abSmrg }
decFloatIsZero(const decFloat * df)245763d1a8abSmrg uInt decFloatIsZero(const decFloat *df) {
245863d1a8abSmrg return DFISZERO(df);
245963d1a8abSmrg } /* decFloatIs... */
246063d1a8abSmrg
246163d1a8abSmrg /* ------------------------------------------------------------------ */
246263d1a8abSmrg /* decFloatLogB -- return adjusted exponent, by 754 rules */
246363d1a8abSmrg /* */
246463d1a8abSmrg /* result gets the adjusted exponent as an integer, or a NaN etc. */
246563d1a8abSmrg /* df is the decFloat to be examined */
246663d1a8abSmrg /* set is the context */
246763d1a8abSmrg /* returns result */
246863d1a8abSmrg /* */
246963d1a8abSmrg /* Notable cases: */
247063d1a8abSmrg /* A<0 -> Use |A| */
247163d1a8abSmrg /* A=0 -> -Infinity (Division by zero) */
247263d1a8abSmrg /* A=Infinite -> +Infinity (Exact) */
247363d1a8abSmrg /* A=1 exactly -> 0 (Exact) */
247463d1a8abSmrg /* NaNs are propagated as usual */
247563d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatLogB(decFloat * result,const decFloat * df,decContext * set)247663d1a8abSmrg decFloat * decFloatLogB(decFloat *result, const decFloat *df,
247763d1a8abSmrg decContext *set) {
247863d1a8abSmrg Int ae; /* adjusted exponent */
247963d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
248063d1a8abSmrg if (DFISINF(df)) {
248163d1a8abSmrg DFWORD(result, 0)=0; /* need +ve */
248263d1a8abSmrg return decInfinity(result, result); /* canonical +Infinity */
248363d1a8abSmrg }
248463d1a8abSmrg if (DFISZERO(df)) {
248563d1a8abSmrg set->status|=DEC_Division_by_zero; /* as per 754 */
248663d1a8abSmrg DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
248763d1a8abSmrg return decInfinity(result, result); /* canonical -Infinity */
248863d1a8abSmrg }
248963d1a8abSmrg ae=GETEXPUN(df) /* get unbiased exponent .. */
249063d1a8abSmrg +decFloatDigits(df)-1; /* .. and make adjusted exponent */
249163d1a8abSmrg /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
249263d1a8abSmrg /* it is worth using a special case of decFloatFromInt32 */
249363d1a8abSmrg DFWORD(result, 0)=ZEROWORD; /* always */
249463d1a8abSmrg if (ae<0) {
249563d1a8abSmrg DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
249663d1a8abSmrg ae=-ae;
249763d1a8abSmrg }
249863d1a8abSmrg #if DOUBLE
249963d1a8abSmrg DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
250063d1a8abSmrg #elif QUAD
250163d1a8abSmrg DFWORD(result, 1)=0;
250263d1a8abSmrg DFWORD(result, 2)=0;
250363d1a8abSmrg DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
250463d1a8abSmrg DFWORD(result, 3)|=BIN2DPD[ae%1000];
250563d1a8abSmrg #endif
250663d1a8abSmrg return result;
250763d1a8abSmrg } /* decFloatLogB */
250863d1a8abSmrg
250963d1a8abSmrg /* ------------------------------------------------------------------ */
251063d1a8abSmrg /* decFloatMax -- return maxnum of two operands */
251163d1a8abSmrg /* */
251263d1a8abSmrg /* result gets the chosen decFloat */
251363d1a8abSmrg /* dfl is the first decFloat (lhs) */
251463d1a8abSmrg /* dfr is the second decFloat (rhs) */
251563d1a8abSmrg /* set is the context */
251663d1a8abSmrg /* returns result */
251763d1a8abSmrg /* */
251863d1a8abSmrg /* If just one operand is a quiet NaN it is ignored. */
251963d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMax(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)252063d1a8abSmrg decFloat * decFloatMax(decFloat *result,
252163d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
252263d1a8abSmrg decContext *set) {
252363d1a8abSmrg Int comp;
252463d1a8abSmrg if (DFISNAN(dfl)) {
252563d1a8abSmrg /* sNaN or both NaNs leads to normal NaN processing */
252663d1a8abSmrg if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
252763d1a8abSmrg return decCanonical(result, dfr); /* RHS is numeric */
252863d1a8abSmrg }
252963d1a8abSmrg if (DFISNAN(dfr)) {
253063d1a8abSmrg /* sNaN leads to normal NaN processing (both NaN handled above) */
253163d1a8abSmrg if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
253263d1a8abSmrg return decCanonical(result, dfl); /* LHS is numeric */
253363d1a8abSmrg }
253463d1a8abSmrg /* Both operands are numeric; numeric comparison needed -- use */
253563d1a8abSmrg /* total order for a well-defined choice (and +0 > -0) */
253663d1a8abSmrg comp=decNumCompare(dfl, dfr, 1);
253763d1a8abSmrg if (comp>=0) return decCanonical(result, dfl);
253863d1a8abSmrg return decCanonical(result, dfr);
253963d1a8abSmrg } /* decFloatMax */
254063d1a8abSmrg
254163d1a8abSmrg /* ------------------------------------------------------------------ */
254263d1a8abSmrg /* decFloatMaxMag -- return maxnummag of two operands */
254363d1a8abSmrg /* */
254463d1a8abSmrg /* result gets the chosen decFloat */
254563d1a8abSmrg /* dfl is the first decFloat (lhs) */
254663d1a8abSmrg /* dfr is the second decFloat (rhs) */
254763d1a8abSmrg /* set is the context */
254863d1a8abSmrg /* returns result */
254963d1a8abSmrg /* */
255063d1a8abSmrg /* Returns according to the magnitude comparisons if both numeric and */
255163d1a8abSmrg /* unequal, otherwise returns maxnum */
255263d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMaxMag(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)255363d1a8abSmrg decFloat * decFloatMaxMag(decFloat *result,
255463d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
255563d1a8abSmrg decContext *set) {
255663d1a8abSmrg Int comp;
255763d1a8abSmrg decFloat absl, absr;
255863d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
255963d1a8abSmrg
256063d1a8abSmrg decFloatCopyAbs(&absl, dfl);
256163d1a8abSmrg decFloatCopyAbs(&absr, dfr);
256263d1a8abSmrg comp=decNumCompare(&absl, &absr, 0);
256363d1a8abSmrg if (comp>0) return decCanonical(result, dfl);
256463d1a8abSmrg if (comp<0) return decCanonical(result, dfr);
256563d1a8abSmrg return decFloatMax(result, dfl, dfr, set);
256663d1a8abSmrg } /* decFloatMaxMag */
256763d1a8abSmrg
256863d1a8abSmrg /* ------------------------------------------------------------------ */
256963d1a8abSmrg /* decFloatMin -- return minnum of two operands */
257063d1a8abSmrg /* */
257163d1a8abSmrg /* result gets the chosen decFloat */
257263d1a8abSmrg /* dfl is the first decFloat (lhs) */
257363d1a8abSmrg /* dfr is the second decFloat (rhs) */
257463d1a8abSmrg /* set is the context */
257563d1a8abSmrg /* returns result */
257663d1a8abSmrg /* */
257763d1a8abSmrg /* If just one operand is a quiet NaN it is ignored. */
257863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMin(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)257963d1a8abSmrg decFloat * decFloatMin(decFloat *result,
258063d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
258163d1a8abSmrg decContext *set) {
258263d1a8abSmrg Int comp;
258363d1a8abSmrg if (DFISNAN(dfl)) {
258463d1a8abSmrg /* sNaN or both NaNs leads to normal NaN processing */
258563d1a8abSmrg if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
258663d1a8abSmrg return decCanonical(result, dfr); /* RHS is numeric */
258763d1a8abSmrg }
258863d1a8abSmrg if (DFISNAN(dfr)) {
258963d1a8abSmrg /* sNaN leads to normal NaN processing (both NaN handled above) */
259063d1a8abSmrg if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
259163d1a8abSmrg return decCanonical(result, dfl); /* LHS is numeric */
259263d1a8abSmrg }
259363d1a8abSmrg /* Both operands are numeric; numeric comparison needed -- use */
259463d1a8abSmrg /* total order for a well-defined choice (and +0 > -0) */
259563d1a8abSmrg comp=decNumCompare(dfl, dfr, 1);
259663d1a8abSmrg if (comp<=0) return decCanonical(result, dfl);
259763d1a8abSmrg return decCanonical(result, dfr);
259863d1a8abSmrg } /* decFloatMin */
259963d1a8abSmrg
260063d1a8abSmrg /* ------------------------------------------------------------------ */
260163d1a8abSmrg /* decFloatMinMag -- return minnummag of two operands */
260263d1a8abSmrg /* */
260363d1a8abSmrg /* result gets the chosen decFloat */
260463d1a8abSmrg /* dfl is the first decFloat (lhs) */
260563d1a8abSmrg /* dfr is the second decFloat (rhs) */
260663d1a8abSmrg /* set is the context */
260763d1a8abSmrg /* returns result */
260863d1a8abSmrg /* */
260963d1a8abSmrg /* Returns according to the magnitude comparisons if both numeric and */
261063d1a8abSmrg /* unequal, otherwise returns minnum */
261163d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMinMag(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)261263d1a8abSmrg decFloat * decFloatMinMag(decFloat *result,
261363d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
261463d1a8abSmrg decContext *set) {
261563d1a8abSmrg Int comp;
261663d1a8abSmrg decFloat absl, absr;
261763d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
261863d1a8abSmrg
261963d1a8abSmrg decFloatCopyAbs(&absl, dfl);
262063d1a8abSmrg decFloatCopyAbs(&absr, dfr);
262163d1a8abSmrg comp=decNumCompare(&absl, &absr, 0);
262263d1a8abSmrg if (comp<0) return decCanonical(result, dfl);
262363d1a8abSmrg if (comp>0) return decCanonical(result, dfr);
262463d1a8abSmrg return decFloatMin(result, dfl, dfr, set);
262563d1a8abSmrg } /* decFloatMinMag */
262663d1a8abSmrg
262763d1a8abSmrg /* ------------------------------------------------------------------ */
262863d1a8abSmrg /* decFloatMinus -- negate value, heeding NaNs, etc. */
262963d1a8abSmrg /* */
263063d1a8abSmrg /* result gets the canonicalized 0-df */
263163d1a8abSmrg /* df is the decFloat to minus */
263263d1a8abSmrg /* set is the context */
263363d1a8abSmrg /* returns result */
263463d1a8abSmrg /* */
263563d1a8abSmrg /* This has the same effect as 0-df where the exponent of the zero is */
263663d1a8abSmrg /* the same as that of df (if df is finite). */
263763d1a8abSmrg /* The effect is also the same as decFloatCopyNegate except that NaNs */
263863d1a8abSmrg /* are handled normally (the sign of a NaN is not affected, and an */
263963d1a8abSmrg /* sNaN will signal), the result is canonical, and zero gets sign 0. */
264063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMinus(decFloat * result,const decFloat * df,decContext * set)264163d1a8abSmrg decFloat * decFloatMinus(decFloat *result, const decFloat *df,
264263d1a8abSmrg decContext *set) {
264363d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
264463d1a8abSmrg decCanonical(result, df); /* copy and check */
264563d1a8abSmrg if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
264663d1a8abSmrg else DFBYTE(result, 0)^=0x80; /* flip sign bit */
264763d1a8abSmrg return result;
264863d1a8abSmrg } /* decFloatMinus */
264963d1a8abSmrg
265063d1a8abSmrg /* ------------------------------------------------------------------ */
265163d1a8abSmrg /* decFloatMultiply -- multiply two decFloats */
265263d1a8abSmrg /* */
265363d1a8abSmrg /* result gets the result of multiplying dfl and dfr: */
265463d1a8abSmrg /* dfl is the first decFloat (lhs) */
265563d1a8abSmrg /* dfr is the second decFloat (rhs) */
265663d1a8abSmrg /* set is the context */
265763d1a8abSmrg /* returns result */
265863d1a8abSmrg /* */
265963d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatMultiply(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)266063d1a8abSmrg decFloat * decFloatMultiply(decFloat *result,
266163d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
266263d1a8abSmrg decContext *set) {
266363d1a8abSmrg bcdnum num; /* for final conversion */
266463d1a8abSmrg uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
266563d1a8abSmrg
266663d1a8abSmrg if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
266763d1a8abSmrg /* NaNs are handled as usual */
266863d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
266963d1a8abSmrg /* infinity times zero is bad */
267063d1a8abSmrg if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
267163d1a8abSmrg if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
267263d1a8abSmrg /* both infinite; return canonical infinity with computed sign */
267363d1a8abSmrg DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
267463d1a8abSmrg return decInfinity(result, result);
267563d1a8abSmrg }
267663d1a8abSmrg
267763d1a8abSmrg /* Here when both operands are finite */
267863d1a8abSmrg decFiniteMultiply(&num, bcdacc, dfl, dfr);
267963d1a8abSmrg return decFinalize(result, &num, set); /* round, check, and lay out */
268063d1a8abSmrg } /* decFloatMultiply */
268163d1a8abSmrg
268263d1a8abSmrg /* ------------------------------------------------------------------ */
268363d1a8abSmrg /* decFloatNextMinus -- next towards -Infinity */
268463d1a8abSmrg /* */
268563d1a8abSmrg /* result gets the next lesser decFloat */
268663d1a8abSmrg /* dfl is the decFloat to start with */
268763d1a8abSmrg /* set is the context */
268863d1a8abSmrg /* returns result */
268963d1a8abSmrg /* */
269063d1a8abSmrg /* This is 754 nextdown; Invalid is the only status possible (from */
269163d1a8abSmrg /* an sNaN). */
269263d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatNextMinus(decFloat * result,const decFloat * dfl,decContext * set)269363d1a8abSmrg decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
269463d1a8abSmrg decContext *set) {
269563d1a8abSmrg decFloat delta; /* tiny increment */
269663d1a8abSmrg uInt savestat; /* saves status */
269763d1a8abSmrg enum rounding saveround; /* .. and mode */
269863d1a8abSmrg
269963d1a8abSmrg /* +Infinity is the special case */
270063d1a8abSmrg if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
270163d1a8abSmrg DFSETNMAX(result);
270263d1a8abSmrg return result; /* [no status to set] */
270363d1a8abSmrg }
270463d1a8abSmrg /* other cases are effected by sutracting a tiny delta -- this */
270563d1a8abSmrg /* should be done in a wider format as the delta is unrepresentable */
270663d1a8abSmrg /* here (but can be done with normal add if the sign of zero is */
270763d1a8abSmrg /* treated carefully, because no Inexactitude is interesting); */
270863d1a8abSmrg /* rounding to -Infinity then pushes the result to next below */
270963d1a8abSmrg decFloatZero(&delta); /* set up tiny delta */
271063d1a8abSmrg DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
271163d1a8abSmrg DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
271263d1a8abSmrg /* set up for the directional round */
271363d1a8abSmrg saveround=set->round; /* save mode */
271463d1a8abSmrg set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
271563d1a8abSmrg savestat=set->status; /* save status */
271663d1a8abSmrg decFloatAdd(result, dfl, &delta, set);
271763d1a8abSmrg /* Add rules mess up the sign when going from +Ntiny to 0 */
271863d1a8abSmrg if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
271963d1a8abSmrg set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
272063d1a8abSmrg set->status|=savestat; /* restore pending flags */
272163d1a8abSmrg set->round=saveround; /* .. and mode */
272263d1a8abSmrg return result;
272363d1a8abSmrg } /* decFloatNextMinus */
272463d1a8abSmrg
272563d1a8abSmrg /* ------------------------------------------------------------------ */
272663d1a8abSmrg /* decFloatNextPlus -- next towards +Infinity */
272763d1a8abSmrg /* */
272863d1a8abSmrg /* result gets the next larger decFloat */
272963d1a8abSmrg /* dfl is the decFloat to start with */
273063d1a8abSmrg /* set is the context */
273163d1a8abSmrg /* returns result */
273263d1a8abSmrg /* */
273363d1a8abSmrg /* This is 754 nextup; Invalid is the only status possible (from */
273463d1a8abSmrg /* an sNaN). */
273563d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatNextPlus(decFloat * result,const decFloat * dfl,decContext * set)273663d1a8abSmrg decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
273763d1a8abSmrg decContext *set) {
273863d1a8abSmrg uInt savestat; /* saves status */
273963d1a8abSmrg enum rounding saveround; /* .. and mode */
274063d1a8abSmrg decFloat delta; /* tiny increment */
274163d1a8abSmrg
274263d1a8abSmrg /* -Infinity is the special case */
274363d1a8abSmrg if (DFISINF(dfl) && DFISSIGNED(dfl)) {
274463d1a8abSmrg DFSETNMAX(result);
274563d1a8abSmrg DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
274663d1a8abSmrg return result; /* [no status to set] */
274763d1a8abSmrg }
274863d1a8abSmrg /* other cases are effected by sutracting a tiny delta -- this */
274963d1a8abSmrg /* should be done in a wider format as the delta is unrepresentable */
275063d1a8abSmrg /* here (but can be done with normal add if the sign of zero is */
275163d1a8abSmrg /* treated carefully, because no Inexactitude is interesting); */
275263d1a8abSmrg /* rounding to +Infinity then pushes the result to next above */
275363d1a8abSmrg decFloatZero(&delta); /* set up tiny delta */
275463d1a8abSmrg DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
275563d1a8abSmrg DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
275663d1a8abSmrg /* set up for the directional round */
275763d1a8abSmrg saveround=set->round; /* save mode */
275863d1a8abSmrg set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
275963d1a8abSmrg savestat=set->status; /* save status */
276063d1a8abSmrg decFloatAdd(result, dfl, &delta, set);
276163d1a8abSmrg /* Add rules mess up the sign when going from -Ntiny to -0 */
276263d1a8abSmrg if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
276363d1a8abSmrg set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
276463d1a8abSmrg set->status|=savestat; /* restore pending flags */
276563d1a8abSmrg set->round=saveround; /* .. and mode */
276663d1a8abSmrg return result;
276763d1a8abSmrg } /* decFloatNextPlus */
276863d1a8abSmrg
276963d1a8abSmrg /* ------------------------------------------------------------------ */
277063d1a8abSmrg /* decFloatNextToward -- next towards a decFloat */
277163d1a8abSmrg /* */
277263d1a8abSmrg /* result gets the next decFloat */
277363d1a8abSmrg /* dfl is the decFloat to start with */
277463d1a8abSmrg /* dfr is the decFloat to move toward */
277563d1a8abSmrg /* set is the context */
277663d1a8abSmrg /* returns result */
277763d1a8abSmrg /* */
277863d1a8abSmrg /* This is 754-1985 nextafter, as modified during revision (dropped */
277963d1a8abSmrg /* from 754-2008); status may be set unless the result is a normal */
278063d1a8abSmrg /* number. */
278163d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatNextToward(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)278263d1a8abSmrg decFloat * decFloatNextToward(decFloat *result,
278363d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
278463d1a8abSmrg decContext *set) {
278563d1a8abSmrg decFloat delta; /* tiny increment or decrement */
278663d1a8abSmrg decFloat pointone; /* 1e-1 */
278763d1a8abSmrg uInt savestat; /* saves status */
278863d1a8abSmrg enum rounding saveround; /* .. and mode */
278963d1a8abSmrg uInt deltatop; /* top word for delta */
279063d1a8abSmrg Int comp; /* work */
279163d1a8abSmrg
279263d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
279363d1a8abSmrg /* Both are numeric, so Invalid no longer a possibility */
279463d1a8abSmrg comp=decNumCompare(dfl, dfr, 0);
279563d1a8abSmrg if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
279663d1a8abSmrg /* unequal; do NextPlus or NextMinus but with different status rules */
279763d1a8abSmrg
279863d1a8abSmrg if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
279963d1a8abSmrg if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
280063d1a8abSmrg DFSETNMAX(result);
280163d1a8abSmrg DFWORD(result, 0)|=DECFLOAT_Sign;
280263d1a8abSmrg return result;
280363d1a8abSmrg }
280463d1a8abSmrg saveround=set->round; /* save mode */
280563d1a8abSmrg set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
280663d1a8abSmrg deltatop=0; /* positive delta */
280763d1a8abSmrg }
280863d1a8abSmrg else { /* lhs>rhs, do NextMinus, see above for commentary */
280963d1a8abSmrg if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
281063d1a8abSmrg DFSETNMAX(result);
281163d1a8abSmrg return result;
281263d1a8abSmrg }
281363d1a8abSmrg saveround=set->round; /* save mode */
281463d1a8abSmrg set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
281563d1a8abSmrg deltatop=DECFLOAT_Sign; /* negative delta */
281663d1a8abSmrg }
281763d1a8abSmrg savestat=set->status; /* save status */
281863d1a8abSmrg /* Here, Inexact is needed where appropriate (and hence Underflow, */
281963d1a8abSmrg /* etc.). Therefore the tiny delta which is otherwise */
282063d1a8abSmrg /* unrepresentable (see NextPlus and NextMinus) is constructed */
282163d1a8abSmrg /* using the multiplication of FMA. */
282263d1a8abSmrg decFloatZero(&delta); /* set up tiny delta */
282363d1a8abSmrg DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
282463d1a8abSmrg DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
282563d1a8abSmrg decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
282663d1a8abSmrg decFloatFMA(result, &delta, &pointone, dfl, set);
282763d1a8abSmrg /* [Delta is truly tiny, so no need to correct sign of zero] */
282863d1a8abSmrg /* use new status unless the result is normal */
282963d1a8abSmrg if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
283063d1a8abSmrg set->round=saveround; /* restore mode */
283163d1a8abSmrg return result;
283263d1a8abSmrg } /* decFloatNextToward */
283363d1a8abSmrg
283463d1a8abSmrg /* ------------------------------------------------------------------ */
283563d1a8abSmrg /* decFloatOr -- logical digitwise OR of two decFloats */
283663d1a8abSmrg /* */
283763d1a8abSmrg /* result gets the result of ORing dfl and dfr */
283863d1a8abSmrg /* dfl is the first decFloat (lhs) */
283963d1a8abSmrg /* dfr is the second decFloat (rhs) */
284063d1a8abSmrg /* set is the context */
284163d1a8abSmrg /* returns result, which will be canonical with sign=0 */
284263d1a8abSmrg /* */
284363d1a8abSmrg /* The operands must be positive, finite with exponent q=0, and */
284463d1a8abSmrg /* comprise just zeros and ones; if not, Invalid operation results. */
284563d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatOr(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)284663d1a8abSmrg decFloat * decFloatOr(decFloat *result,
284763d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
284863d1a8abSmrg decContext *set) {
284963d1a8abSmrg if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
285063d1a8abSmrg || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
285163d1a8abSmrg /* the operands are positive finite integers (q=0) with just 0s and 1s */
285263d1a8abSmrg #if DOUBLE
285363d1a8abSmrg DFWORD(result, 0)=ZEROWORD
285463d1a8abSmrg |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
285563d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
285663d1a8abSmrg #elif QUAD
285763d1a8abSmrg DFWORD(result, 0)=ZEROWORD
285863d1a8abSmrg |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
285963d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
286063d1a8abSmrg DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
286163d1a8abSmrg DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
286263d1a8abSmrg #endif
286363d1a8abSmrg return result;
286463d1a8abSmrg } /* decFloatOr */
286563d1a8abSmrg
286663d1a8abSmrg /* ------------------------------------------------------------------ */
286763d1a8abSmrg /* decFloatPlus -- add value to 0, heeding NaNs, etc. */
286863d1a8abSmrg /* */
286963d1a8abSmrg /* result gets the canonicalized 0+df */
287063d1a8abSmrg /* df is the decFloat to plus */
287163d1a8abSmrg /* set is the context */
287263d1a8abSmrg /* returns result */
287363d1a8abSmrg /* */
287463d1a8abSmrg /* This has the same effect as 0+df where the exponent of the zero is */
287563d1a8abSmrg /* the same as that of df (if df is finite). */
287663d1a8abSmrg /* The effect is also the same as decFloatCopy except that NaNs */
287763d1a8abSmrg /* are handled normally (the sign of a NaN is not affected, and an */
287863d1a8abSmrg /* sNaN will signal), the result is canonical, and zero gets sign 0. */
287963d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatPlus(decFloat * result,const decFloat * df,decContext * set)288063d1a8abSmrg decFloat * decFloatPlus(decFloat *result, const decFloat *df,
288163d1a8abSmrg decContext *set) {
288263d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
288363d1a8abSmrg decCanonical(result, df); /* copy and check */
288463d1a8abSmrg if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
288563d1a8abSmrg return result;
288663d1a8abSmrg } /* decFloatPlus */
288763d1a8abSmrg
288863d1a8abSmrg /* ------------------------------------------------------------------ */
288963d1a8abSmrg /* decFloatQuantize -- quantize a decFloat */
289063d1a8abSmrg /* */
289163d1a8abSmrg /* result gets the result of quantizing dfl to match dfr */
289263d1a8abSmrg /* dfl is the first decFloat (lhs) */
289363d1a8abSmrg /* dfr is the second decFloat (rhs), which sets the exponent */
289463d1a8abSmrg /* set is the context */
289563d1a8abSmrg /* returns result */
289663d1a8abSmrg /* */
289763d1a8abSmrg /* Unless there is an error or the result is infinite, the exponent */
289863d1a8abSmrg /* of result is guaranteed to be the same as that of dfr. */
289963d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatQuantize(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)290063d1a8abSmrg decFloat * decFloatQuantize(decFloat *result,
290163d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
290263d1a8abSmrg decContext *set) {
290363d1a8abSmrg Int explb, exprb; /* left and right biased exponents */
290463d1a8abSmrg uByte *ulsd; /* local LSD pointer */
290563d1a8abSmrg uByte *ub, *uc; /* work */
290663d1a8abSmrg Int drop; /* .. */
290763d1a8abSmrg uInt dpd; /* .. */
290863d1a8abSmrg uInt encode; /* encoding accumulator */
290963d1a8abSmrg uInt sourhil, sourhir; /* top words from source decFloats */
291063d1a8abSmrg uInt uiwork; /* for macros */
291163d1a8abSmrg #if QUAD
291263d1a8abSmrg uShort uswork; /* .. */
291363d1a8abSmrg #endif
291463d1a8abSmrg /* the following buffer holds the coefficient for manipulation */
291563d1a8abSmrg uByte buf[4+DECPMAX*3+2*QUAD]; /* + space for zeros to left or right */
291663d1a8abSmrg #if DECTRACE
291763d1a8abSmrg bcdnum num; /* for trace displays */
291863d1a8abSmrg #endif
291963d1a8abSmrg
292063d1a8abSmrg /* Start decoding the arguments */
292163d1a8abSmrg sourhil=DFWORD(dfl, 0); /* LHS top word */
292263d1a8abSmrg explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
292363d1a8abSmrg sourhir=DFWORD(dfr, 0); /* RHS top word */
292463d1a8abSmrg exprb=DECCOMBEXP[sourhir>>26];
292563d1a8abSmrg
292663d1a8abSmrg if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
292763d1a8abSmrg /* NaNs are handled as usual */
292863d1a8abSmrg if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
292963d1a8abSmrg /* one infinity but not both is bad */
293063d1a8abSmrg if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
293163d1a8abSmrg /* both infinite; return canonical infinity with sign of LHS */
293263d1a8abSmrg return decInfinity(result, dfl);
293363d1a8abSmrg }
293463d1a8abSmrg
293563d1a8abSmrg /* Here when both arguments are finite */
293663d1a8abSmrg /* complete extraction of the exponents [no need to unbias] */
293763d1a8abSmrg explb+=GETECON(dfl); /* + continuation */
293863d1a8abSmrg exprb+=GETECON(dfr); /* .. */
293963d1a8abSmrg
294063d1a8abSmrg /* calculate the number of digits to drop from the coefficient */
294163d1a8abSmrg drop=exprb-explb; /* 0 if nothing to do */
294263d1a8abSmrg if (drop==0) return decCanonical(result, dfl); /* return canonical */
294363d1a8abSmrg
294463d1a8abSmrg /* the coefficient is needed; lay it out into buf, offset so zeros */
294563d1a8abSmrg /* can be added before or after as needed -- an extra heading is */
294663d1a8abSmrg /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
294763d1a8abSmrg /* fours */
294863d1a8abSmrg #define BUFOFF (buf+4+DECPMAX)
294963d1a8abSmrg GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
295063d1a8abSmrg /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
295163d1a8abSmrg
295263d1a8abSmrg #if DECTRACE
295363d1a8abSmrg num.msd=BUFOFF;
295463d1a8abSmrg num.lsd=BUFOFF+DECPMAX-1;
295563d1a8abSmrg num.exponent=explb-DECBIAS;
295663d1a8abSmrg num.sign=sourhil & DECFLOAT_Sign;
295763d1a8abSmrg decShowNum(&num, "dfl");
295863d1a8abSmrg #endif
295963d1a8abSmrg
296063d1a8abSmrg if (drop>0) { /* [most common case] */
296163d1a8abSmrg /* (this code is very similar to that in decFloatFinalize, but */
296263d1a8abSmrg /* has many differences so is duplicated here -- so any changes */
296363d1a8abSmrg /* may need to be made there, too) */
296463d1a8abSmrg uByte *roundat; /* -> re-round digit */
296563d1a8abSmrg uByte reround; /* reround value */
296663d1a8abSmrg /* printf("Rounding; drop=%ld\n", (LI)drop); */
296763d1a8abSmrg
296863d1a8abSmrg /* there is at least one zero needed to the left, in all but one */
296963d1a8abSmrg /* exceptional (all-nines) case, so place four zeros now; this is */
297063d1a8abSmrg /* needed almost always and makes rounding all-nines by fours safe */
297163d1a8abSmrg UBFROMUI(BUFOFF-4, 0);
297263d1a8abSmrg
297363d1a8abSmrg /* Three cases here: */
297463d1a8abSmrg /* 1. new LSD is in coefficient (almost always) */
297563d1a8abSmrg /* 2. new LSD is digit to left of coefficient (so MSD is */
297663d1a8abSmrg /* round-for-reround digit) */
297763d1a8abSmrg /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
297863d1a8abSmrg /* Note that leading zeros can safely be treated as useful digits */
297963d1a8abSmrg
298063d1a8abSmrg /* [duplicate check-stickies code to save a test] */
298163d1a8abSmrg /* [by-digit check for stickies as runs of zeros are rare] */
298263d1a8abSmrg if (drop<DECPMAX) { /* NB lengths not addresses */
298363d1a8abSmrg roundat=BUFOFF+DECPMAX-drop;
298463d1a8abSmrg reround=*roundat;
298563d1a8abSmrg for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
298663d1a8abSmrg if (*ub!=0) { /* non-zero to be discarded */
298763d1a8abSmrg reround=DECSTICKYTAB[reround]; /* apply sticky bit */
298863d1a8abSmrg break; /* [remainder don't-care] */
298963d1a8abSmrg }
299063d1a8abSmrg } /* check stickies */
299163d1a8abSmrg ulsd=roundat-1; /* set LSD */
299263d1a8abSmrg }
299363d1a8abSmrg else { /* edge case */
299463d1a8abSmrg if (drop==DECPMAX) {
299563d1a8abSmrg roundat=BUFOFF;
299663d1a8abSmrg reround=*roundat;
299763d1a8abSmrg }
299863d1a8abSmrg else {
299963d1a8abSmrg roundat=BUFOFF-1;
300063d1a8abSmrg reround=0;
300163d1a8abSmrg }
300263d1a8abSmrg for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
300363d1a8abSmrg if (*ub!=0) { /* non-zero to be discarded */
300463d1a8abSmrg reround=DECSTICKYTAB[reround]; /* apply sticky bit */
300563d1a8abSmrg break; /* [remainder don't-care] */
300663d1a8abSmrg }
300763d1a8abSmrg } /* check stickies */
300863d1a8abSmrg *BUFOFF=0; /* make a coefficient of 0 */
300963d1a8abSmrg ulsd=BUFOFF; /* .. at the MSD place */
301063d1a8abSmrg }
301163d1a8abSmrg
301263d1a8abSmrg if (reround!=0) { /* discarding non-zero */
301363d1a8abSmrg uInt bump=0;
301463d1a8abSmrg set->status|=DEC_Inexact;
301563d1a8abSmrg
301663d1a8abSmrg /* next decide whether to increment the coefficient */
301763d1a8abSmrg if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
301863d1a8abSmrg if (reround>5) bump=1; /* >0.5 goes up */
301963d1a8abSmrg else if (reround==5) /* exactly 0.5000 .. */
302063d1a8abSmrg bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */
302163d1a8abSmrg } /* r-h-e */
302263d1a8abSmrg else switch (set->round) {
302363d1a8abSmrg case DEC_ROUND_DOWN: {
302463d1a8abSmrg /* no change */
302563d1a8abSmrg break;} /* r-d */
302663d1a8abSmrg case DEC_ROUND_HALF_DOWN: {
302763d1a8abSmrg if (reround>5) bump=1;
302863d1a8abSmrg break;} /* r-h-d */
302963d1a8abSmrg case DEC_ROUND_HALF_UP: {
303063d1a8abSmrg if (reround>=5) bump=1;
303163d1a8abSmrg break;} /* r-h-u */
303263d1a8abSmrg case DEC_ROUND_UP: {
303363d1a8abSmrg if (reround>0) bump=1;
303463d1a8abSmrg break;} /* r-u */
303563d1a8abSmrg case DEC_ROUND_CEILING: {
303663d1a8abSmrg /* same as _UP for positive numbers, and as _DOWN for negatives */
303763d1a8abSmrg if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
303863d1a8abSmrg break;} /* r-c */
303963d1a8abSmrg case DEC_ROUND_FLOOR: {
304063d1a8abSmrg /* same as _UP for negative numbers, and as _DOWN for positive */
304163d1a8abSmrg /* [negative reround cannot occur on 0] */
304263d1a8abSmrg if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
304363d1a8abSmrg break;} /* r-f */
304463d1a8abSmrg case DEC_ROUND_05UP: {
304563d1a8abSmrg if (reround>0) { /* anything out there is 'sticky' */
304663d1a8abSmrg /* bump iff lsd=0 or 5; this cannot carry so it could be */
304763d1a8abSmrg /* effected immediately with no bump -- but the code */
304863d1a8abSmrg /* is clearer if this is done the same way as the others */
304963d1a8abSmrg if (*ulsd==0 || *ulsd==5) bump=1;
305063d1a8abSmrg }
305163d1a8abSmrg break;} /* r-r */
305263d1a8abSmrg default: { /* e.g., DEC_ROUND_MAX */
305363d1a8abSmrg set->status|=DEC_Invalid_context;
305463d1a8abSmrg #if DECCHECK
305563d1a8abSmrg printf("Unknown rounding mode: %ld\n", (LI)set->round);
305663d1a8abSmrg #endif
305763d1a8abSmrg break;}
305863d1a8abSmrg } /* switch (not r-h-e) */
305963d1a8abSmrg /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
306063d1a8abSmrg
306163d1a8abSmrg if (bump!=0) { /* need increment */
306263d1a8abSmrg /* increment the coefficient; this could give 1000... (after */
306363d1a8abSmrg /* the all nines case) */
306463d1a8abSmrg ub=ulsd;
306563d1a8abSmrg for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
306663d1a8abSmrg /* now at most 3 digits left to non-9 (usually just the one) */
306763d1a8abSmrg for (; *ub==9; ub--) *ub=0;
306863d1a8abSmrg *ub+=1;
306963d1a8abSmrg /* [the all-nines case will have carried one digit to the */
307063d1a8abSmrg /* left of the original MSD -- just where it is needed] */
307163d1a8abSmrg } /* bump needed */
307263d1a8abSmrg } /* inexact rounding */
307363d1a8abSmrg
307463d1a8abSmrg /* now clear zeros to the left so exactly DECPMAX digits will be */
307563d1a8abSmrg /* available in the coefficent -- the first word to the left was */
307663d1a8abSmrg /* cleared earlier for safe carry; now add any more needed */
307763d1a8abSmrg if (drop>4) {
307863d1a8abSmrg UBFROMUI(BUFOFF-8, 0); /* must be at least 5 */
307963d1a8abSmrg for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0);
308063d1a8abSmrg }
308163d1a8abSmrg } /* need round (drop>0) */
308263d1a8abSmrg
308363d1a8abSmrg else { /* drop<0; padding with -drop digits is needed */
308463d1a8abSmrg /* This is the case where an error can occur if the padded */
308563d1a8abSmrg /* coefficient will not fit; checking for this can be done in the */
308663d1a8abSmrg /* same loop as padding for zeros if the no-hope and zero cases */
308763d1a8abSmrg /* are checked first */
308863d1a8abSmrg if (-drop>DECPMAX-1) { /* cannot fit unless 0 */
308963d1a8abSmrg if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
309063d1a8abSmrg /* a zero can have any exponent; just drop through and use it */
309163d1a8abSmrg ulsd=BUFOFF+DECPMAX-1;
309263d1a8abSmrg }
309363d1a8abSmrg else { /* padding will fit (but may still be too long) */
309463d1a8abSmrg /* final-word mask depends on endianess */
309563d1a8abSmrg #if DECLITEND
309663d1a8abSmrg static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
309763d1a8abSmrg #else
309863d1a8abSmrg static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
309963d1a8abSmrg #endif
310063d1a8abSmrg /* note that here zeros to the right are added by fours, so in */
310163d1a8abSmrg /* the Quad case this could write 36 zeros if the coefficient has */
310263d1a8abSmrg /* fewer than three significant digits (hence the +2*QUAD for buf) */
310363d1a8abSmrg for (uc=BUFOFF+DECPMAX;; uc+=4) {
310463d1a8abSmrg UBFROMUI(uc, 0);
310563d1a8abSmrg if (UBTOUI(uc-DECPMAX)!=0) { /* could be bad */
310663d1a8abSmrg /* if all four digits should be zero, definitely bad */
310763d1a8abSmrg if (uc<=BUFOFF+DECPMAX+(-drop)-4)
310863d1a8abSmrg return decInvalid(result, set);
310963d1a8abSmrg /* must be a 1- to 3-digit sequence; check more carefully */
311063d1a8abSmrg if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0)
311163d1a8abSmrg return decInvalid(result, set);
311263d1a8abSmrg break; /* no need for loop end test */
311363d1a8abSmrg }
311463d1a8abSmrg if (uc>=BUFOFF+DECPMAX+(-drop)-4) break; /* done */
311563d1a8abSmrg }
311663d1a8abSmrg ulsd=BUFOFF+DECPMAX+(-drop)-1;
311763d1a8abSmrg } /* pad and check leading zeros */
311863d1a8abSmrg } /* drop<0 */
311963d1a8abSmrg
312063d1a8abSmrg #if DECTRACE
312163d1a8abSmrg num.msd=ulsd-DECPMAX+1;
312263d1a8abSmrg num.lsd=ulsd;
312363d1a8abSmrg num.exponent=explb-DECBIAS;
312463d1a8abSmrg num.sign=sourhil & DECFLOAT_Sign;
312563d1a8abSmrg decShowNum(&num, "res");
312663d1a8abSmrg #endif
312763d1a8abSmrg
312863d1a8abSmrg /*------------------------------------------------------------------*/
312963d1a8abSmrg /* At this point the result is DECPMAX digits, ending at ulsd, so */
313063d1a8abSmrg /* fits the encoding exactly; there is no possibility of error */
313163d1a8abSmrg /*------------------------------------------------------------------*/
313263d1a8abSmrg encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
313363d1a8abSmrg encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */
313463d1a8abSmrg /* the exponent continuation can be extracted from the original RHS */
313563d1a8abSmrg encode|=sourhir & ECONMASK;
313663d1a8abSmrg encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */
313763d1a8abSmrg
313863d1a8abSmrg /* finally encode the coefficient */
313963d1a8abSmrg /* private macro to encode a declet; this version can be used */
314063d1a8abSmrg /* because all coefficient digits exist */
314163d1a8abSmrg #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
314263d1a8abSmrg dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
314363d1a8abSmrg
314463d1a8abSmrg #if DOUBLE
314563d1a8abSmrg getDPD3q(dpd, 4); encode|=dpd<<8;
314663d1a8abSmrg getDPD3q(dpd, 3); encode|=dpd>>2;
314763d1a8abSmrg DFWORD(result, 0)=encode;
314863d1a8abSmrg encode=dpd<<30;
314963d1a8abSmrg getDPD3q(dpd, 2); encode|=dpd<<20;
315063d1a8abSmrg getDPD3q(dpd, 1); encode|=dpd<<10;
315163d1a8abSmrg getDPD3q(dpd, 0); encode|=dpd;
315263d1a8abSmrg DFWORD(result, 1)=encode;
315363d1a8abSmrg
315463d1a8abSmrg #elif QUAD
315563d1a8abSmrg getDPD3q(dpd,10); encode|=dpd<<4;
315663d1a8abSmrg getDPD3q(dpd, 9); encode|=dpd>>6;
315763d1a8abSmrg DFWORD(result, 0)=encode;
315863d1a8abSmrg encode=dpd<<26;
315963d1a8abSmrg getDPD3q(dpd, 8); encode|=dpd<<16;
316063d1a8abSmrg getDPD3q(dpd, 7); encode|=dpd<<6;
316163d1a8abSmrg getDPD3q(dpd, 6); encode|=dpd>>4;
316263d1a8abSmrg DFWORD(result, 1)=encode;
316363d1a8abSmrg encode=dpd<<28;
316463d1a8abSmrg getDPD3q(dpd, 5); encode|=dpd<<18;
316563d1a8abSmrg getDPD3q(dpd, 4); encode|=dpd<<8;
316663d1a8abSmrg getDPD3q(dpd, 3); encode|=dpd>>2;
316763d1a8abSmrg DFWORD(result, 2)=encode;
316863d1a8abSmrg encode=dpd<<30;
316963d1a8abSmrg getDPD3q(dpd, 2); encode|=dpd<<20;
317063d1a8abSmrg getDPD3q(dpd, 1); encode|=dpd<<10;
317163d1a8abSmrg getDPD3q(dpd, 0); encode|=dpd;
317263d1a8abSmrg DFWORD(result, 3)=encode;
317363d1a8abSmrg #endif
317463d1a8abSmrg return result;
317563d1a8abSmrg } /* decFloatQuantize */
317663d1a8abSmrg
317763d1a8abSmrg /* ------------------------------------------------------------------ */
317863d1a8abSmrg /* decFloatReduce -- reduce finite coefficient to minimum length */
317963d1a8abSmrg /* */
318063d1a8abSmrg /* result gets the reduced decFloat */
318163d1a8abSmrg /* df is the source decFloat */
318263d1a8abSmrg /* set is the context */
318363d1a8abSmrg /* returns result, which will be canonical */
318463d1a8abSmrg /* */
318563d1a8abSmrg /* This removes all possible trailing zeros from the coefficient; */
318663d1a8abSmrg /* some may remain when the number is very close to Nmax. */
318763d1a8abSmrg /* Special values are unchanged and no status is set unless df=sNaN. */
318863d1a8abSmrg /* Reduced zero has an exponent q=0. */
318963d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatReduce(decFloat * result,const decFloat * df,decContext * set)319063d1a8abSmrg decFloat * decFloatReduce(decFloat *result, const decFloat *df,
319163d1a8abSmrg decContext *set) {
319263d1a8abSmrg bcdnum num; /* work */
319363d1a8abSmrg uByte buf[DECPMAX], *ub; /* coefficient and pointer */
319463d1a8abSmrg if (df!=result) *result=*df; /* copy, if needed */
319563d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */
319663d1a8abSmrg /* zeros and infinites propagate too */
319763d1a8abSmrg if (DFISINF(df)) return decInfinity(result, df); /* canonical */
319863d1a8abSmrg if (DFISZERO(df)) {
319963d1a8abSmrg uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
320063d1a8abSmrg decFloatZero(result);
320163d1a8abSmrg DFWORD(result, 0)|=sign;
320263d1a8abSmrg return result; /* exponent dropped, sign OK */
320363d1a8abSmrg }
320463d1a8abSmrg /* non-zero finite */
320563d1a8abSmrg GETCOEFF(df, buf);
320663d1a8abSmrg ub=buf+DECPMAX-1; /* -> lsd */
320763d1a8abSmrg if (*ub) return result; /* no trailing zeros */
320863d1a8abSmrg for (ub--; *ub==0;) ub--; /* terminates because non-zero */
320963d1a8abSmrg /* *ub is the first non-zero from the right */
321063d1a8abSmrg num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
321163d1a8abSmrg num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
321263d1a8abSmrg num.msd=buf;
321363d1a8abSmrg num.lsd=ub;
321463d1a8abSmrg return decFinalize(result, &num, set);
321563d1a8abSmrg } /* decFloatReduce */
321663d1a8abSmrg
321763d1a8abSmrg /* ------------------------------------------------------------------ */
321863d1a8abSmrg /* decFloatRemainder -- integer divide and return remainder */
321963d1a8abSmrg /* */
322063d1a8abSmrg /* result gets the remainder of dividing dfl by dfr: */
322163d1a8abSmrg /* dfl is the first decFloat (lhs) */
322263d1a8abSmrg /* dfr is the second decFloat (rhs) */
322363d1a8abSmrg /* set is the context */
322463d1a8abSmrg /* returns result */
322563d1a8abSmrg /* */
322663d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatRemainder(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)322763d1a8abSmrg decFloat * decFloatRemainder(decFloat *result,
322863d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
322963d1a8abSmrg decContext *set) {
323063d1a8abSmrg return decDivide(result, dfl, dfr, set, REMAINDER);
323163d1a8abSmrg } /* decFloatRemainder */
323263d1a8abSmrg
323363d1a8abSmrg /* ------------------------------------------------------------------ */
323463d1a8abSmrg /* decFloatRemainderNear -- integer divide to nearest and remainder */
323563d1a8abSmrg /* */
323663d1a8abSmrg /* result gets the remainder of dividing dfl by dfr: */
323763d1a8abSmrg /* dfl is the first decFloat (lhs) */
323863d1a8abSmrg /* dfr is the second decFloat (rhs) */
323963d1a8abSmrg /* set is the context */
324063d1a8abSmrg /* returns result */
324163d1a8abSmrg /* */
324263d1a8abSmrg /* This is the IEEE remainder, where the nearest integer is used. */
324363d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatRemainderNear(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)324463d1a8abSmrg decFloat * decFloatRemainderNear(decFloat *result,
324563d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
324663d1a8abSmrg decContext *set) {
324763d1a8abSmrg return decDivide(result, dfl, dfr, set, REMNEAR);
324863d1a8abSmrg } /* decFloatRemainderNear */
324963d1a8abSmrg
325063d1a8abSmrg /* ------------------------------------------------------------------ */
325163d1a8abSmrg /* decFloatRotate -- rotate the coefficient of a decFloat left/right */
325263d1a8abSmrg /* */
325363d1a8abSmrg /* result gets the result of rotating dfl */
325463d1a8abSmrg /* dfl is the source decFloat to rotate */
325563d1a8abSmrg /* dfr is the count of digits to rotate, an integer (with q=0) */
325663d1a8abSmrg /* set is the context */
325763d1a8abSmrg /* returns result */
325863d1a8abSmrg /* */
325963d1a8abSmrg /* The digits of the coefficient of dfl are rotated to the left (if */
326063d1a8abSmrg /* dfr is positive) or to the right (if dfr is negative) without */
326163d1a8abSmrg /* adjusting the exponent or the sign of dfl. */
326263d1a8abSmrg /* */
326363d1a8abSmrg /* dfr must be in the range -DECPMAX through +DECPMAX. */
326463d1a8abSmrg /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
326563d1a8abSmrg /* dfr must be valid). No status is set unless dfr is invalid or an */
326663d1a8abSmrg /* operand is an sNaN. The result is canonical. */
326763d1a8abSmrg /* ------------------------------------------------------------------ */
326863d1a8abSmrg #define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
decFloatRotate(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)326963d1a8abSmrg decFloat * decFloatRotate(decFloat *result,
327063d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
327163d1a8abSmrg decContext *set) {
327263d1a8abSmrg Int rotate; /* dfr as an Int */
327363d1a8abSmrg uByte buf[DECPMAX+PHALF]; /* coefficient + half */
327463d1a8abSmrg uInt digits, savestat; /* work */
327563d1a8abSmrg bcdnum num; /* .. */
327663d1a8abSmrg uByte *ub; /* .. */
327763d1a8abSmrg
327863d1a8abSmrg if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
327963d1a8abSmrg if (!DFISINT(dfr)) return decInvalid(result, set);
328063d1a8abSmrg digits=decFloatDigits(dfr); /* calculate digits */
328163d1a8abSmrg if (digits>2) return decInvalid(result, set); /* definitely out of range */
328263d1a8abSmrg rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
328363d1a8abSmrg if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
328463d1a8abSmrg /* [from here on no error or status change is possible] */
328563d1a8abSmrg if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
328663d1a8abSmrg /* handle no-rotate cases */
328763d1a8abSmrg if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
328863d1a8abSmrg /* a real rotate is needed: 0 < rotate < DECPMAX */
328963d1a8abSmrg /* reduce the rotation to no more than half to reduce copying later */
329063d1a8abSmrg /* (for QUAD in fact half + 2 digits) */
329163d1a8abSmrg if (DFISSIGNED(dfr)) rotate=-rotate;
329263d1a8abSmrg if (abs(rotate)>PHALF) {
329363d1a8abSmrg if (rotate<0) rotate=DECPMAX+rotate;
329463d1a8abSmrg else rotate=rotate-DECPMAX;
329563d1a8abSmrg }
329663d1a8abSmrg /* now lay out the coefficient, leaving room to the right or the */
329763d1a8abSmrg /* left depending on the direction of rotation */
329863d1a8abSmrg ub=buf;
329963d1a8abSmrg if (rotate<0) ub+=PHALF; /* rotate right, so space to left */
330063d1a8abSmrg GETCOEFF(dfl, ub);
330163d1a8abSmrg /* copy half the digits to left or right, and set num.msd */
330263d1a8abSmrg if (rotate<0) {
330363d1a8abSmrg memcpy(buf, buf+DECPMAX, PHALF);
330463d1a8abSmrg num.msd=buf+PHALF+rotate;
330563d1a8abSmrg }
330663d1a8abSmrg else {
330763d1a8abSmrg memcpy(buf+DECPMAX, buf, PHALF);
330863d1a8abSmrg num.msd=buf+rotate;
330963d1a8abSmrg }
331063d1a8abSmrg /* fill in rest of num */
331163d1a8abSmrg num.lsd=num.msd+DECPMAX-1;
331263d1a8abSmrg num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
331363d1a8abSmrg num.exponent=GETEXPUN(dfl);
331463d1a8abSmrg savestat=set->status; /* record */
331563d1a8abSmrg decFinalize(result, &num, set);
331663d1a8abSmrg set->status=savestat; /* restore */
331763d1a8abSmrg return result;
331863d1a8abSmrg } /* decFloatRotate */
331963d1a8abSmrg
332063d1a8abSmrg /* ------------------------------------------------------------------ */
332163d1a8abSmrg /* decFloatSameQuantum -- test decFloats for same quantum */
332263d1a8abSmrg /* */
332363d1a8abSmrg /* dfl is the first decFloat (lhs) */
332463d1a8abSmrg /* dfr is the second decFloat (rhs) */
332563d1a8abSmrg /* returns 1 if the operands have the same quantum, 0 otherwise */
332663d1a8abSmrg /* */
332763d1a8abSmrg /* No error is possible and no status results. */
332863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatSameQuantum(const decFloat * dfl,const decFloat * dfr)332963d1a8abSmrg uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
333063d1a8abSmrg if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
333163d1a8abSmrg if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
333263d1a8abSmrg if (DFISINF(dfl) && DFISINF(dfr)) return 1;
333363d1a8abSmrg return 0; /* any other special mixture gives false */
333463d1a8abSmrg }
333563d1a8abSmrg if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
333663d1a8abSmrg return 0;
333763d1a8abSmrg } /* decFloatSameQuantum */
333863d1a8abSmrg
333963d1a8abSmrg /* ------------------------------------------------------------------ */
334063d1a8abSmrg /* decFloatScaleB -- multiply by a power of 10, as per 754 */
334163d1a8abSmrg /* */
334263d1a8abSmrg /* result gets the result of the operation */
334363d1a8abSmrg /* dfl is the first decFloat (lhs) */
334463d1a8abSmrg /* dfr is the second decFloat (rhs), am integer (with q=0) */
334563d1a8abSmrg /* set is the context */
334663d1a8abSmrg /* returns result */
334763d1a8abSmrg /* */
334863d1a8abSmrg /* This computes result=dfl x 10**dfr where dfr is an integer in the */
334963d1a8abSmrg /* range +/-2*(emax+pmax), typically resulting from LogB. */
335063d1a8abSmrg /* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
335163d1a8abSmrg /* as usual. */
335263d1a8abSmrg /* ------------------------------------------------------------------ */
335363d1a8abSmrg #define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
decFloatScaleB(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)335463d1a8abSmrg decFloat * decFloatScaleB(decFloat *result,
335563d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
335663d1a8abSmrg decContext *set) {
335763d1a8abSmrg uInt digits; /* work */
335863d1a8abSmrg Int expr; /* dfr as an Int */
335963d1a8abSmrg
336063d1a8abSmrg if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
336163d1a8abSmrg if (!DFISINT(dfr)) return decInvalid(result, set);
336263d1a8abSmrg digits=decFloatDigits(dfr); /* calculate digits */
336363d1a8abSmrg
336463d1a8abSmrg #if DOUBLE
336563d1a8abSmrg if (digits>3) return decInvalid(result, set); /* definitely out of range */
336663d1a8abSmrg expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */
336763d1a8abSmrg #elif QUAD
336863d1a8abSmrg if (digits>5) return decInvalid(result, set); /* definitely out of range */
336963d1a8abSmrg expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */
337063d1a8abSmrg +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */
337163d1a8abSmrg #endif
337263d1a8abSmrg if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */
337363d1a8abSmrg /* [from now on no error possible] */
337463d1a8abSmrg if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
337563d1a8abSmrg if (DFISSIGNED(dfr)) expr=-expr;
337663d1a8abSmrg /* dfl is finite and expr is valid */
337763d1a8abSmrg *result=*dfl; /* copy to target */
337863d1a8abSmrg return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
337963d1a8abSmrg } /* decFloatScaleB */
338063d1a8abSmrg
338163d1a8abSmrg /* ------------------------------------------------------------------ */
338263d1a8abSmrg /* decFloatShift -- shift the coefficient of a decFloat left or right */
338363d1a8abSmrg /* */
338463d1a8abSmrg /* result gets the result of shifting dfl */
338563d1a8abSmrg /* dfl is the source decFloat to shift */
338663d1a8abSmrg /* dfr is the count of digits to shift, an integer (with q=0) */
338763d1a8abSmrg /* set is the context */
338863d1a8abSmrg /* returns result */
338963d1a8abSmrg /* */
339063d1a8abSmrg /* The digits of the coefficient of dfl are shifted to the left (if */
339163d1a8abSmrg /* dfr is positive) or to the right (if dfr is negative) without */
339263d1a8abSmrg /* adjusting the exponent or the sign of dfl. */
339363d1a8abSmrg /* */
339463d1a8abSmrg /* dfr must be in the range -DECPMAX through +DECPMAX. */
339563d1a8abSmrg /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
339663d1a8abSmrg /* dfr must be valid). No status is set unless dfr is invalid or an */
339763d1a8abSmrg /* operand is an sNaN. The result is canonical. */
339863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatShift(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)339963d1a8abSmrg decFloat * decFloatShift(decFloat *result,
340063d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
340163d1a8abSmrg decContext *set) {
340263d1a8abSmrg Int shift; /* dfr as an Int */
340363d1a8abSmrg uByte buf[DECPMAX*2]; /* coefficient + padding */
340463d1a8abSmrg uInt digits, savestat; /* work */
340563d1a8abSmrg bcdnum num; /* .. */
340663d1a8abSmrg uInt uiwork; /* for macros */
340763d1a8abSmrg
340863d1a8abSmrg if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
340963d1a8abSmrg if (!DFISINT(dfr)) return decInvalid(result, set);
341063d1a8abSmrg digits=decFloatDigits(dfr); /* calculate digits */
341163d1a8abSmrg if (digits>2) return decInvalid(result, set); /* definitely out of range */
341263d1a8abSmrg shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
341363d1a8abSmrg if (shift>DECPMAX) return decInvalid(result, set); /* too big */
341463d1a8abSmrg /* [from here on no error or status change is possible] */
341563d1a8abSmrg
341663d1a8abSmrg if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
341763d1a8abSmrg /* handle no-shift and all-shift (clear to zero) cases */
341863d1a8abSmrg if (shift==0) return decCanonical(result, dfl);
341963d1a8abSmrg if (shift==DECPMAX) { /* zero with sign */
342063d1a8abSmrg uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
342163d1a8abSmrg decFloatZero(result); /* make +0 */
342263d1a8abSmrg DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
342363d1a8abSmrg /* [cannot safely use CopySign] */
342463d1a8abSmrg return result;
342563d1a8abSmrg }
342663d1a8abSmrg /* a real shift is needed: 0 < shift < DECPMAX */
342763d1a8abSmrg num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
342863d1a8abSmrg num.exponent=GETEXPUN(dfl);
342963d1a8abSmrg num.msd=buf;
343063d1a8abSmrg GETCOEFF(dfl, buf);
343163d1a8abSmrg if (DFISSIGNED(dfr)) { /* shift right */
343263d1a8abSmrg /* edge cases are taken care of, so this is easy */
343363d1a8abSmrg num.lsd=buf+DECPMAX-shift-1;
343463d1a8abSmrg }
343563d1a8abSmrg else { /* shift left -- zero padding needed to right */
343663d1a8abSmrg UBFROMUI(buf+DECPMAX, 0); /* 8 will handle most cases */
343763d1a8abSmrg UBFROMUI(buf+DECPMAX+4, 0); /* .. */
343863d1a8abSmrg if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
343963d1a8abSmrg num.msd+=shift;
344063d1a8abSmrg num.lsd=num.msd+DECPMAX-1;
344163d1a8abSmrg }
344263d1a8abSmrg savestat=set->status; /* record */
344363d1a8abSmrg decFinalize(result, &num, set);
344463d1a8abSmrg set->status=savestat; /* restore */
344563d1a8abSmrg return result;
344663d1a8abSmrg } /* decFloatShift */
344763d1a8abSmrg
344863d1a8abSmrg /* ------------------------------------------------------------------ */
344963d1a8abSmrg /* decFloatSubtract -- subtract a decFloat from another */
345063d1a8abSmrg /* */
345163d1a8abSmrg /* result gets the result of subtracting dfr from dfl: */
345263d1a8abSmrg /* dfl is the first decFloat (lhs) */
345363d1a8abSmrg /* dfr is the second decFloat (rhs) */
345463d1a8abSmrg /* set is the context */
345563d1a8abSmrg /* returns result */
345663d1a8abSmrg /* */
345763d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatSubtract(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)345863d1a8abSmrg decFloat * decFloatSubtract(decFloat *result,
345963d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
346063d1a8abSmrg decContext *set) {
346163d1a8abSmrg decFloat temp;
346263d1a8abSmrg /* NaNs must propagate without sign change */
346363d1a8abSmrg if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
346463d1a8abSmrg temp=*dfr; /* make a copy */
346563d1a8abSmrg DFBYTE(&temp, 0)^=0x80; /* flip sign */
346663d1a8abSmrg return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
346763d1a8abSmrg } /* decFloatSubtract */
346863d1a8abSmrg
346963d1a8abSmrg /* ------------------------------------------------------------------ */
347063d1a8abSmrg /* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
347163d1a8abSmrg /* */
347263d1a8abSmrg /* df is the decFloat to round */
347363d1a8abSmrg /* set is the context */
347463d1a8abSmrg /* round is the rounding mode to use */
347563d1a8abSmrg /* returns a uInt or an Int, rounded according to the name */
347663d1a8abSmrg /* */
347763d1a8abSmrg /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
347863d1a8abSmrg /* outside the range of the target; Inexact will not be signaled for */
347963d1a8abSmrg /* simple rounding unless 'Exact' appears in the name. */
348063d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatToUInt32(const decFloat * df,decContext * set,enum rounding round)348163d1a8abSmrg uInt decFloatToUInt32(const decFloat *df, decContext *set,
348263d1a8abSmrg enum rounding round) {
348363d1a8abSmrg return decToInt32(df, set, round, 0, 1);}
348463d1a8abSmrg
decFloatToUInt32Exact(const decFloat * df,decContext * set,enum rounding round)348563d1a8abSmrg uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
348663d1a8abSmrg enum rounding round) {
348763d1a8abSmrg return decToInt32(df, set, round, 1, 1);}
348863d1a8abSmrg
decFloatToInt32(const decFloat * df,decContext * set,enum rounding round)348963d1a8abSmrg Int decFloatToInt32(const decFloat *df, decContext *set,
349063d1a8abSmrg enum rounding round) {
349163d1a8abSmrg return (Int)decToInt32(df, set, round, 0, 0);}
349263d1a8abSmrg
decFloatToInt32Exact(const decFloat * df,decContext * set,enum rounding round)349363d1a8abSmrg Int decFloatToInt32Exact(const decFloat *df, decContext *set,
349463d1a8abSmrg enum rounding round) {
349563d1a8abSmrg return (Int)decToInt32(df, set, round, 1, 0);}
349663d1a8abSmrg
349763d1a8abSmrg /* ------------------------------------------------------------------ */
349863d1a8abSmrg /* decFloatToIntegral -- round to integral value (two flavours) */
349963d1a8abSmrg /* */
350063d1a8abSmrg /* result gets the result */
350163d1a8abSmrg /* df is the decFloat to round */
350263d1a8abSmrg /* set is the context */
350363d1a8abSmrg /* round is the rounding mode to use */
350463d1a8abSmrg /* returns result */
350563d1a8abSmrg /* */
350663d1a8abSmrg /* No exceptions, even Inexact, are raised except for sNaN input, or */
350763d1a8abSmrg /* if 'Exact' appears in the name. */
350863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatToIntegralValue(decFloat * result,const decFloat * df,decContext * set,enum rounding round)350963d1a8abSmrg decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
351063d1a8abSmrg decContext *set, enum rounding round) {
351163d1a8abSmrg return decToIntegral(result, df, set, round, 0);}
351263d1a8abSmrg
decFloatToIntegralExact(decFloat * result,const decFloat * df,decContext * set)351363d1a8abSmrg decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
351463d1a8abSmrg decContext *set) {
351563d1a8abSmrg return decToIntegral(result, df, set, set->round, 1);}
351663d1a8abSmrg
351763d1a8abSmrg /* ------------------------------------------------------------------ */
351863d1a8abSmrg /* decFloatXor -- logical digitwise XOR of two decFloats */
351963d1a8abSmrg /* */
352063d1a8abSmrg /* result gets the result of XORing dfl and dfr */
352163d1a8abSmrg /* dfl is the first decFloat (lhs) */
352263d1a8abSmrg /* dfr is the second decFloat (rhs) */
352363d1a8abSmrg /* set is the context */
352463d1a8abSmrg /* returns result, which will be canonical with sign=0 */
352563d1a8abSmrg /* */
352663d1a8abSmrg /* The operands must be positive, finite with exponent q=0, and */
352763d1a8abSmrg /* comprise just zeros and ones; if not, Invalid operation results. */
352863d1a8abSmrg /* ------------------------------------------------------------------ */
decFloatXor(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)352963d1a8abSmrg decFloat * decFloatXor(decFloat *result,
353063d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
353163d1a8abSmrg decContext *set) {
353263d1a8abSmrg if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
353363d1a8abSmrg || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
353463d1a8abSmrg /* the operands are positive finite integers (q=0) with just 0s and 1s */
353563d1a8abSmrg #if DOUBLE
353663d1a8abSmrg DFWORD(result, 0)=ZEROWORD
353763d1a8abSmrg |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
353863d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
353963d1a8abSmrg #elif QUAD
354063d1a8abSmrg DFWORD(result, 0)=ZEROWORD
354163d1a8abSmrg |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
354263d1a8abSmrg DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
354363d1a8abSmrg DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
354463d1a8abSmrg DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
354563d1a8abSmrg #endif
354663d1a8abSmrg return result;
354763d1a8abSmrg } /* decFloatXor */
354863d1a8abSmrg
354963d1a8abSmrg /* ------------------------------------------------------------------ */
355063d1a8abSmrg /* decInvalid -- set Invalid_operation result */
355163d1a8abSmrg /* */
355263d1a8abSmrg /* result gets a canonical NaN */
355363d1a8abSmrg /* set is the context */
355463d1a8abSmrg /* returns result */
355563d1a8abSmrg /* */
355663d1a8abSmrg /* status has Invalid_operation added */
355763d1a8abSmrg /* ------------------------------------------------------------------ */
decInvalid(decFloat * result,decContext * set)355863d1a8abSmrg static decFloat *decInvalid(decFloat *result, decContext *set) {
355963d1a8abSmrg decFloatZero(result);
356063d1a8abSmrg DFWORD(result, 0)=DECFLOAT_qNaN;
356163d1a8abSmrg set->status|=DEC_Invalid_operation;
356263d1a8abSmrg return result;
356363d1a8abSmrg } /* decInvalid */
356463d1a8abSmrg
356563d1a8abSmrg /* ------------------------------------------------------------------ */
356663d1a8abSmrg /* decInfinity -- set canonical Infinity with sign from a decFloat */
356763d1a8abSmrg /* */
356863d1a8abSmrg /* result gets a canonical Infinity */
356963d1a8abSmrg /* df is source decFloat (only the sign is used) */
357063d1a8abSmrg /* returns result */
357163d1a8abSmrg /* */
357263d1a8abSmrg /* df may be the same as result */
357363d1a8abSmrg /* ------------------------------------------------------------------ */
decInfinity(decFloat * result,const decFloat * df)357463d1a8abSmrg static decFloat *decInfinity(decFloat *result, const decFloat *df) {
357563d1a8abSmrg uInt sign=DFWORD(df, 0); /* save source signword */
357663d1a8abSmrg decFloatZero(result); /* clear everything */
357763d1a8abSmrg DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
357863d1a8abSmrg return result;
357963d1a8abSmrg } /* decInfinity */
358063d1a8abSmrg
358163d1a8abSmrg /* ------------------------------------------------------------------ */
358263d1a8abSmrg /* decNaNs -- handle NaN argument(s) */
358363d1a8abSmrg /* */
358463d1a8abSmrg /* result gets the result of handling dfl and dfr, one or both of */
358563d1a8abSmrg /* which is a NaN */
358663d1a8abSmrg /* dfl is the first decFloat (lhs) */
358763d1a8abSmrg /* dfr is the second decFloat (rhs) -- may be NULL for a single- */
358863d1a8abSmrg /* operand operation */
358963d1a8abSmrg /* set is the context */
359063d1a8abSmrg /* returns result */
359163d1a8abSmrg /* */
359263d1a8abSmrg /* Called when one or both operands is a NaN, and propagates the */
359363d1a8abSmrg /* appropriate result to res. When an sNaN is found, it is changed */
359463d1a8abSmrg /* to a qNaN and Invalid operation is set. */
359563d1a8abSmrg /* ------------------------------------------------------------------ */
decNaNs(decFloat * result,const decFloat * dfl,const decFloat * dfr,decContext * set)359663d1a8abSmrg static decFloat *decNaNs(decFloat *result,
359763d1a8abSmrg const decFloat *dfl, const decFloat *dfr,
359863d1a8abSmrg decContext *set) {
359963d1a8abSmrg /* handle sNaNs first */
360063d1a8abSmrg if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
360163d1a8abSmrg if (DFISSNAN(dfl)) {
360263d1a8abSmrg decCanonical(result, dfl); /* propagate canonical sNaN */
360363d1a8abSmrg DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
360463d1a8abSmrg set->status|=DEC_Invalid_operation;
360563d1a8abSmrg return result;
360663d1a8abSmrg }
360763d1a8abSmrg /* one or both is a quiet NaN */
360863d1a8abSmrg if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */
360963d1a8abSmrg return decCanonical(result, dfl); /* propagate canonical qNaN */
361063d1a8abSmrg } /* decNaNs */
361163d1a8abSmrg
361263d1a8abSmrg /* ------------------------------------------------------------------ */
361363d1a8abSmrg /* decNumCompare -- numeric comparison of two decFloats */
361463d1a8abSmrg /* */
361563d1a8abSmrg /* dfl is the left-hand decFloat, which is not a NaN */
361663d1a8abSmrg /* dfr is the right-hand decFloat, which is not a NaN */
361763d1a8abSmrg /* tot is 1 for total order compare, 0 for simple numeric */
361863d1a8abSmrg /* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
361963d1a8abSmrg /* */
362063d1a8abSmrg /* No error is possible; status and mode are unchanged. */
362163d1a8abSmrg /* ------------------------------------------------------------------ */
decNumCompare(const decFloat * dfl,const decFloat * dfr,Flag tot)362263d1a8abSmrg static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
362363d1a8abSmrg Int sigl, sigr; /* LHS and RHS non-0 signums */
362463d1a8abSmrg Int shift; /* shift needed to align operands */
362563d1a8abSmrg uByte *ub, *uc; /* work */
362663d1a8abSmrg uInt uiwork; /* for macros */
362763d1a8abSmrg /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
362863d1a8abSmrg uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */
362963d1a8abSmrg uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */
363063d1a8abSmrg
363163d1a8abSmrg sigl=1;
363263d1a8abSmrg if (DFISSIGNED(dfl)) {
363363d1a8abSmrg if (!DFISSIGNED(dfr)) { /* -LHS +RHS */
363463d1a8abSmrg if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
363563d1a8abSmrg return -1; /* RHS wins */
363663d1a8abSmrg }
363763d1a8abSmrg sigl=-1;
363863d1a8abSmrg }
363963d1a8abSmrg if (DFISSIGNED(dfr)) {
364063d1a8abSmrg if (!DFISSIGNED(dfl)) { /* +LHS -RHS */
364163d1a8abSmrg if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
364263d1a8abSmrg return +1; /* LHS wins */
364363d1a8abSmrg }
364463d1a8abSmrg }
364563d1a8abSmrg
364663d1a8abSmrg /* signs are the same; operand(s) could be zero */
364763d1a8abSmrg sigr=-sigl; /* sign to return if abs(RHS) wins */
364863d1a8abSmrg
364963d1a8abSmrg if (DFISINF(dfl)) {
365063d1a8abSmrg if (DFISINF(dfr)) return 0; /* both infinite & same sign */
365163d1a8abSmrg return sigl; /* inf > n */
365263d1a8abSmrg }
365363d1a8abSmrg if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
365463d1a8abSmrg
365563d1a8abSmrg /* here, both are same sign and finite; calculate their offset */
365663d1a8abSmrg shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */
365763d1a8abSmrg /* [bias can be ignored -- the absolute exponent is not relevant] */
365863d1a8abSmrg
365963d1a8abSmrg if (DFISZERO(dfl)) {
366063d1a8abSmrg if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */
366163d1a8abSmrg /* both are zero, return 0 if both same exponent or numeric compare */
366263d1a8abSmrg if (shift==0 || !tot) return 0;
366363d1a8abSmrg if (shift>0) return sigl;
366463d1a8abSmrg return sigr; /* [shift<0] */
366563d1a8abSmrg }
366663d1a8abSmrg else { /* LHS!=0 */
366763d1a8abSmrg if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */
366863d1a8abSmrg }
366963d1a8abSmrg /* both are known to be non-zero at this point */
367063d1a8abSmrg
367163d1a8abSmrg /* if the exponents are so different that the coefficients do not */
367263d1a8abSmrg /* overlap (by even one digit) then a full comparison is not needed */
367363d1a8abSmrg if (abs(shift)>=DECPMAX) { /* no overlap */
367463d1a8abSmrg /* coefficients are known to be non-zero */
367563d1a8abSmrg if (shift>0) return sigl;
367663d1a8abSmrg return sigr; /* [shift<0] */
367763d1a8abSmrg }
367863d1a8abSmrg
367963d1a8abSmrg /* decode the coefficients */
368063d1a8abSmrg /* (shift both right two if Quad to make a multiple of four) */
368163d1a8abSmrg #if QUAD
368263d1a8abSmrg UBFROMUI(bufl, 0);
368363d1a8abSmrg UBFROMUI(bufr, 0);
368463d1a8abSmrg #endif
368563d1a8abSmrg GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
368663d1a8abSmrg GETCOEFF(dfr, bufr+QUAD*2); /* .. */
368763d1a8abSmrg if (shift==0) { /* aligned; common and easy */
368863d1a8abSmrg /* all multiples of four, here */
368963d1a8abSmrg for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
369063d1a8abSmrg uInt ui=UBTOUI(ub);
369163d1a8abSmrg if (ui==UBTOUI(uc)) continue; /* so far so same */
369263d1a8abSmrg /* about to find a winner; go by bytes in case little-endian */
369363d1a8abSmrg for (;; ub++, uc++) {
369463d1a8abSmrg if (*ub>*uc) return sigl; /* difference found */
369563d1a8abSmrg if (*ub<*uc) return sigr; /* .. */
369663d1a8abSmrg }
369763d1a8abSmrg }
369863d1a8abSmrg } /* aligned */
369963d1a8abSmrg else if (shift>0) { /* lhs to left */
370063d1a8abSmrg ub=bufl; /* RHS pointer */
370163d1a8abSmrg /* pad bufl so right-aligned; most shifts will fit in 8 */
370263d1a8abSmrg UBFROMUI(bufl+DECPMAX+QUAD*2, 0); /* add eight zeros */
370363d1a8abSmrg UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); /* .. */
370463d1a8abSmrg if (shift>8) {
370563d1a8abSmrg /* more than eight; fill the rest, and also worth doing the */
370663d1a8abSmrg /* lead-in by fours */
370763d1a8abSmrg uByte *up; /* work */
370863d1a8abSmrg uByte *upend=bufl+DECPMAX+QUAD*2+shift;
370963d1a8abSmrg for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
371063d1a8abSmrg /* [pads up to 36 in all for Quad] */
371163d1a8abSmrg for (;; ub+=4) {
371263d1a8abSmrg if (UBTOUI(ub)!=0) return sigl;
371363d1a8abSmrg if (ub+4>bufl+shift-4) break;
371463d1a8abSmrg }
371563d1a8abSmrg }
371663d1a8abSmrg /* check remaining leading digits */
371763d1a8abSmrg for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
371863d1a8abSmrg /* now start the overlapped part; bufl has been padded, so the */
371963d1a8abSmrg /* comparison can go for the full length of bufr, which is a */
372063d1a8abSmrg /* multiple of 4 bytes */
372163d1a8abSmrg for (uc=bufr; ; uc+=4, ub+=4) {
372263d1a8abSmrg uInt ui=UBTOUI(ub);
372363d1a8abSmrg if (ui!=UBTOUI(uc)) { /* mismatch found */
372463d1a8abSmrg for (;; uc++, ub++) { /* check from left [little-endian?] */
372563d1a8abSmrg if (*ub>*uc) return sigl; /* difference found */
372663d1a8abSmrg if (*ub<*uc) return sigr; /* .. */
372763d1a8abSmrg }
372863d1a8abSmrg } /* mismatch */
372963d1a8abSmrg if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
373063d1a8abSmrg }
373163d1a8abSmrg } /* shift>0 */
373263d1a8abSmrg
373363d1a8abSmrg else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
373463d1a8abSmrg uc=bufr; /* RHS pointer */
373563d1a8abSmrg /* pad bufr so right-aligned; most shifts will fit in 8 */
373663d1a8abSmrg UBFROMUI(bufr+DECPMAX+QUAD*2, 0); /* add eight zeros */
373763d1a8abSmrg UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); /* .. */
373863d1a8abSmrg if (shift<-8) {
373963d1a8abSmrg /* more than eight; fill the rest, and also worth doing the */
374063d1a8abSmrg /* lead-in by fours */
374163d1a8abSmrg uByte *up; /* work */
374263d1a8abSmrg uByte *upend=bufr+DECPMAX+QUAD*2-shift;
374363d1a8abSmrg for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
374463d1a8abSmrg /* [pads up to 36 in all for Quad] */
374563d1a8abSmrg for (;; uc+=4) {
374663d1a8abSmrg if (UBTOUI(uc)!=0) return sigr;
374763d1a8abSmrg if (uc+4>bufr-shift-4) break;
374863d1a8abSmrg }
374963d1a8abSmrg }
375063d1a8abSmrg /* check remaining leading digits */
375163d1a8abSmrg for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
375263d1a8abSmrg /* now start the overlapped part; bufr has been padded, so the */
375363d1a8abSmrg /* comparison can go for the full length of bufl, which is a */
375463d1a8abSmrg /* multiple of 4 bytes */
375563d1a8abSmrg for (ub=bufl; ; ub+=4, uc+=4) {
375663d1a8abSmrg uInt ui=UBTOUI(ub);
375763d1a8abSmrg if (ui!=UBTOUI(uc)) { /* mismatch found */
375863d1a8abSmrg for (;; ub++, uc++) { /* check from left [little-endian?] */
375963d1a8abSmrg if (*ub>*uc) return sigl; /* difference found */
376063d1a8abSmrg if (*ub<*uc) return sigr; /* .. */
376163d1a8abSmrg }
376263d1a8abSmrg } /* mismatch */
376363d1a8abSmrg if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
376463d1a8abSmrg }
376563d1a8abSmrg } /* shift<0 */
376663d1a8abSmrg
376763d1a8abSmrg /* Here when compare equal */
376863d1a8abSmrg if (!tot) return 0; /* numerically equal */
376963d1a8abSmrg /* total ordering .. exponent matters */
377063d1a8abSmrg if (shift>0) return sigl; /* total order by exponent */
377163d1a8abSmrg if (shift<0) return sigr; /* .. */
377263d1a8abSmrg return 0;
377363d1a8abSmrg } /* decNumCompare */
377463d1a8abSmrg
377563d1a8abSmrg /* ------------------------------------------------------------------ */
377663d1a8abSmrg /* decToInt32 -- local routine to effect ToInteger conversions */
377763d1a8abSmrg /* */
377863d1a8abSmrg /* df is the decFloat to convert */
377963d1a8abSmrg /* set is the context */
378063d1a8abSmrg /* rmode is the rounding mode to use */
378163d1a8abSmrg /* exact is 1 if Inexact should be signalled */
378263d1a8abSmrg /* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
378363d1a8abSmrg /* returns 32-bit result as a uInt */
378463d1a8abSmrg /* */
378563d1a8abSmrg /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
378663d1a8abSmrg /* these cases 0 is returned. */
378763d1a8abSmrg /* ------------------------------------------------------------------ */
decToInt32(const decFloat * df,decContext * set,enum rounding rmode,Flag exact,Flag unsign)378863d1a8abSmrg static uInt decToInt32(const decFloat *df, decContext *set,
378963d1a8abSmrg enum rounding rmode, Flag exact, Flag unsign) {
379063d1a8abSmrg Int exp; /* exponent */
379163d1a8abSmrg uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */
379263d1a8abSmrg uInt hi, lo; /* .. penultimate, least, etc. */
379363d1a8abSmrg decFloat zero, result; /* work */
379463d1a8abSmrg Int i; /* .. */
379563d1a8abSmrg
379663d1a8abSmrg /* Start decoding the argument */
379763d1a8abSmrg sourhi=DFWORD(df, 0); /* top word */
379863d1a8abSmrg exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
379963d1a8abSmrg if (EXPISSPECIAL(exp)) { /* is special? */
380063d1a8abSmrg set->status|=DEC_Invalid_operation; /* signal */
380163d1a8abSmrg return 0;
380263d1a8abSmrg }
380363d1a8abSmrg
380463d1a8abSmrg /* Here when the argument is finite */
380563d1a8abSmrg if (GETEXPUN(df)==0) result=*df; /* already a true integer */
380663d1a8abSmrg else { /* need to round to integer */
380763d1a8abSmrg enum rounding saveround; /* saver */
380863d1a8abSmrg uInt savestatus; /* .. */
380963d1a8abSmrg saveround=set->round; /* save rounding mode .. */
381063d1a8abSmrg savestatus=set->status; /* .. and status */
381163d1a8abSmrg set->round=rmode; /* set mode */
381263d1a8abSmrg decFloatZero(&zero); /* make 0E+0 */
381363d1a8abSmrg set->status=0; /* clear */
381463d1a8abSmrg decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
381563d1a8abSmrg set->round=saveround; /* restore rounding mode .. */
381663d1a8abSmrg if (exact) set->status|=savestatus; /* include Inexact */
381763d1a8abSmrg else set->status=savestatus; /* .. or just original status */
381863d1a8abSmrg }
381963d1a8abSmrg
382063d1a8abSmrg /* only the last four declets of the coefficient can contain */
382163d1a8abSmrg /* non-zero; check for others (and also NaN or Infinity from the */
382263d1a8abSmrg /* Quantize) first (see DFISZERO for explanation): */
382363d1a8abSmrg /* decFloatShow(&result, "sofar"); */
382463d1a8abSmrg #if DOUBLE
382563d1a8abSmrg if ((DFWORD(&result, 0)&0x1c03ff00)!=0
382663d1a8abSmrg || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
382763d1a8abSmrg #elif QUAD
382863d1a8abSmrg if ((DFWORD(&result, 2)&0xffffff00)!=0
382963d1a8abSmrg || DFWORD(&result, 1)!=0
383063d1a8abSmrg || (DFWORD(&result, 0)&0x1c003fff)!=0
383163d1a8abSmrg || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
383263d1a8abSmrg #endif
383363d1a8abSmrg set->status|=DEC_Invalid_operation; /* Invalid or out of range */
383463d1a8abSmrg return 0;
383563d1a8abSmrg }
383663d1a8abSmrg /* get last twelve digits of the coefficent into hi & ho, base */
383763d1a8abSmrg /* 10**9 (see GETCOEFFBILL): */
383863d1a8abSmrg sourlo=DFWORD(&result, DECWORDS-1);
383963d1a8abSmrg lo=DPD2BIN0[sourlo&0x3ff]
384063d1a8abSmrg +DPD2BINK[(sourlo>>10)&0x3ff]
384163d1a8abSmrg +DPD2BINM[(sourlo>>20)&0x3ff];
384263d1a8abSmrg sourpen=DFWORD(&result, DECWORDS-2);
384363d1a8abSmrg hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
384463d1a8abSmrg
384563d1a8abSmrg /* according to request, check range carefully */
384663d1a8abSmrg if (unsign) {
384763d1a8abSmrg if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
384863d1a8abSmrg set->status|=DEC_Invalid_operation; /* out of range */
384963d1a8abSmrg return 0;
385063d1a8abSmrg }
385163d1a8abSmrg return hi*BILLION+lo;
385263d1a8abSmrg }
385363d1a8abSmrg /* signed */
385463d1a8abSmrg if (hi>2 || (hi==2 && lo>147483647)) {
385563d1a8abSmrg /* handle the usual edge case */
385663d1a8abSmrg if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
385763d1a8abSmrg set->status|=DEC_Invalid_operation; /* truly out of range */
385863d1a8abSmrg return 0;
385963d1a8abSmrg }
386063d1a8abSmrg i=hi*BILLION+lo;
386163d1a8abSmrg if (DFISSIGNED(&result)) i=-i;
386263d1a8abSmrg return (uInt)i;
386363d1a8abSmrg } /* decToInt32 */
386463d1a8abSmrg
386563d1a8abSmrg /* ------------------------------------------------------------------ */
386663d1a8abSmrg /* decToIntegral -- local routine to effect ToIntegral value */
386763d1a8abSmrg /* */
386863d1a8abSmrg /* result gets the result */
386963d1a8abSmrg /* df is the decFloat to round */
387063d1a8abSmrg /* set is the context */
387163d1a8abSmrg /* rmode is the rounding mode to use */
387263d1a8abSmrg /* exact is 1 if Inexact should be signalled */
387363d1a8abSmrg /* returns result */
387463d1a8abSmrg /* ------------------------------------------------------------------ */
387563d1a8abSmrg static decFloat * decToIntegral(decFloat *result, const decFloat *df,
387663d1a8abSmrg decContext *set, enum rounding rmode,
387763d1a8abSmrg Flag exact) {
387863d1a8abSmrg Int exp; /* exponent */
387963d1a8abSmrg uInt sourhi; /* top word from source decFloat */
388063d1a8abSmrg enum rounding saveround; /* saver */
388163d1a8abSmrg uInt savestatus; /* .. */
388263d1a8abSmrg decFloat zero; /* work */
388363d1a8abSmrg
388463d1a8abSmrg /* Start decoding the argument */
388563d1a8abSmrg sourhi=DFWORD(df, 0); /* top word */
388663d1a8abSmrg exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
388763d1a8abSmrg
388863d1a8abSmrg if (EXPISSPECIAL(exp)) { /* is special? */
388963d1a8abSmrg /* NaNs are handled as usual */
389063d1a8abSmrg if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
389163d1a8abSmrg /* must be infinite; return canonical infinity with sign of df */
389263d1a8abSmrg return decInfinity(result, df);
389363d1a8abSmrg }
389463d1a8abSmrg
389563d1a8abSmrg /* Here when the argument is finite */
389663d1a8abSmrg /* complete extraction of the exponent */
389763d1a8abSmrg exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */
389863d1a8abSmrg
389963d1a8abSmrg if (exp>=0) return decCanonical(result, df); /* already integral */
390063d1a8abSmrg
390163d1a8abSmrg saveround=set->round; /* save rounding mode .. */
390263d1a8abSmrg savestatus=set->status; /* .. and status */
390363d1a8abSmrg set->round=rmode; /* set mode */
390463d1a8abSmrg decFloatZero(&zero); /* make 0E+0 */
390563d1a8abSmrg decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
390663d1a8abSmrg set->round=saveround; /* restore rounding mode .. */
390763d1a8abSmrg if (!exact) set->status=savestatus; /* .. and status, unless exact */
390863d1a8abSmrg return result;
390963d1a8abSmrg } /* decToIntegral */
3910