1 /* Common base code for the decNumber C Library.
2 Copyright (C) 2007, 2009 Free Software Foundation, Inc.
3 Contributed by IBM Corporation. Author Mike Cowlishaw.
5 This file is part of GCC.
7 GCC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU General Public License as published by the Free
9 Software Foundation; either version 3, or (at your option) any later
12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 Under Section 7 of GPL version 3, you are granted additional
18 permissions described in the GCC Runtime Library Exception, version
19 3.1, as published by the Free Software Foundation.
21 You should have received a copy of the GNU General Public License and
22 a copy of the GCC Runtime Library Exception along with this program;
23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 <http://www.gnu.org/licenses/>. */
26 /* ------------------------------------------------------------------ */
27 /* decBasic.c -- common base code for Basic decimal types */
28 /* ------------------------------------------------------------------ */
29 /* This module comprises code that is shared between decDouble and */
30 /* decQuad (but not decSingle). The main arithmetic operations are */
31 /* here (Add, Subtract, Multiply, FMA, and Division operators). */
33 /* Unlike decNumber, parameterization takes place at compile time */
34 /* rather than at runtime. The parameters are set in the decDouble.c */
35 /* (etc.) files, which then include this one to produce the compiled */
36 /* code. The functions here, therefore, are code shared between */
37 /* multiple formats. */
39 /* This must be included after decCommon.c. */
40 /* ------------------------------------------------------------------ */
41 /* Names here refer to decFloat rather than to decDouble, etc., and */
42 /* the functions are in strict alphabetical order. */
44 /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
47 #error decBasic.c must be included after decCommon.c
50 #error Routines in decBasic.c are for decDouble and decQuad only
53 /* Private constants */
54 #define DIVIDE 0x80000000 /* Divide operations [as flags] */
55 #define REMAINDER 0x40000000 /* .. */
56 #define DIVIDEINT 0x20000000 /* .. */
57 #define REMNEAR 0x10000000 /* .. */
59 /* Private functions (local, used only by routines in this module) */
60 static decFloat *decDivide(decFloat *, const decFloat *,
61 const decFloat *, decContext *, uInt);
62 static decFloat *decCanonical(decFloat *, const decFloat *);
63 static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
65 static decFloat *decInfinity(decFloat *, const decFloat *);
66 static decFloat *decInvalid(decFloat *, decContext *);
67 static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
69 static Int decNumCompare(const decFloat *, const decFloat *, Flag);
70 static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
72 static uInt decToInt32(const decFloat *, decContext *, enum rounding,
75 /* ------------------------------------------------------------------ */
76 /* decCanonical -- copy a decFloat, making canonical */
78 /* result gets the canonicalized df */
79 /* df is the decFloat to copy and make canonical */
82 /* This is exposed via decFloatCanonical for Double and Quad only. */
83 /* This works on specials, too; no error or exception is possible. */
84 /* ------------------------------------------------------------------ */
85 static decFloat * decCanonical(decFloat *result, const decFloat *df) {
86 uInt encode, precode, dpd; /* work */
87 uInt inword, uoff, canon; /* .. */
88 Int n; /* counter (down) */
89 if (df!=result) *result=*df; /* effect copy if needed */
90 if (DFISSPECIAL(result)) {
91 if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
93 DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */
94 if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */
95 /* drop through to check payload */
97 /* return quickly if the coefficient continuation is canonical */
100 uInt sourhi=DFWORD(df, 0);
101 uInt sourlo=DFWORD(df, 1);
102 if (CANONDPDOFF(sourhi, 8)
103 && CANONDPDTWO(sourhi, sourlo, 30)
104 && CANONDPDOFF(sourlo, 20)
105 && CANONDPDOFF(sourlo, 10)
106 && CANONDPDOFF(sourlo, 0)) return result;
108 uInt sourhi=DFWORD(df, 0);
109 uInt sourmh=DFWORD(df, 1);
110 uInt sourml=DFWORD(df, 2);
111 uInt sourlo=DFWORD(df, 3);
112 if (CANONDPDOFF(sourhi, 4)
113 && CANONDPDTWO(sourhi, sourmh, 26)
114 && CANONDPDOFF(sourmh, 16)
115 && CANONDPDOFF(sourmh, 6)
116 && CANONDPDTWO(sourmh, sourml, 28)
117 && CANONDPDOFF(sourml, 18)
118 && CANONDPDOFF(sourml, 8)
119 && CANONDPDTWO(sourml, sourlo, 30)
120 && CANONDPDOFF(sourlo, 20)
121 && CANONDPDOFF(sourlo, 10)
122 && CANONDPDOFF(sourlo, 0)) return result;
126 /* Loop to repair a non-canonical coefficent, as needed */
127 inword=DECWORDS-1; /* current input word */
128 uoff=0; /* bit offset of declet */
129 encode=DFWORD(result, inword);
130 for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */
133 if (uoff>32) { /* crossed uInt boundary */
135 encode=DFWORD(result, inword);
137 dpd|=encode<<(10-uoff); /* get pending bits */
139 dpd&=0x3ff; /* clear uninteresting bits */
140 if (dpd<0x16e) continue; /* must be canonical */
141 canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */
142 if (canon==dpd) continue; /* have canonical declet */
143 /* need to replace declet */
144 if (uoff>=10) { /* all within current word */
145 encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
146 encode|=canon<<(uoff-10); /* insert the canonical form */
147 DFWORD(result, inword)=encode; /* .. and save */
150 /* straddled words */
151 precode=DFWORD(result, inword+1); /* get previous */
152 precode&=0xffffffff>>(10-uoff); /* clear top bits */
153 DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
154 encode&=0xffffffff<<uoff; /* clear bottom bits */
155 encode|=canon>>(10-uoff); /* insert canonical */
156 DFWORD(result, inword)=encode; /* .. and save */
161 /* ------------------------------------------------------------------ */
162 /* decDivide -- divide operations */
164 /* result gets the result of dividing dfl by dfr: */
165 /* dfl is the first decFloat (lhs) */
166 /* dfr is the second decFloat (rhs) */
167 /* set is the context */
168 /* op is the operation selector */
171 /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
172 /* ------------------------------------------------------------------ */
173 #define DIVCOUNT 0 /* 1 to instrument subtractions counter */
174 #define DIVBASE ((uInt)BILLION) /* the base used for divide */
175 #define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
176 #define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */
177 static decFloat * decDivide(decFloat *result, const decFloat *dfl,
178 const decFloat *dfr, decContext *set, uInt op) {
179 decFloat quotient; /* for remainders */
180 bcdnum num; /* for final conversion */
181 uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
182 uInt div[DIVOPLEN]; /* divisor in base-billion .. */
183 uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */
184 uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
185 uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */
186 Int divunits, accunits; /* lengths */
187 Int quodigits; /* digits in quotient */
188 uInt *lsua, *lsuq; /* -> current acc and quo lsus */
189 Int length, multiplier; /* work */
190 uInt carry, sign; /* .. */
191 uInt *ua, *ud, *uq; /* .. */
193 uInt uiwork; /* for macros */
194 uInt divtop; /* top unit of div adjusted for estimating */
196 static uInt maxcount=0; /* worst-seen subtractions count */
197 uInt divcount=0; /* subtractions count [this divide] */
201 num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
203 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
204 /* NaNs are handled as usual */
205 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
206 /* one or two infinities */
208 if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
209 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
210 /* Infinity/x is infinite and quiet, even if x=0 */
211 DFWORD(result, 0)=num.sign;
212 return decInfinity(result, result);
214 /* must be x/Infinity -- remainders are lhs */
215 if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
216 /* divides: return zero with correct sign and exponent depending */
217 /* on op (Etiny for divide, 0 for divideInt) */
218 decFloatZero(result);
219 if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
220 else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */
223 /* next, handle zero operands (x/0 and 0/x) */
224 if (DFISZERO(dfr)) { /* x/0 */
225 if (DFISZERO(dfl)) { /* 0/0 is undefined */
226 decFloatZero(result);
227 DFWORD(result, 0)=DECFLOAT_qNaN;
228 set->status|=DEC_Division_undefined;
231 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
232 set->status|=DEC_Division_by_zero;
233 DFWORD(result, 0)=num.sign;
234 return decInfinity(result, result); /* x/0 -> signed Infinity */
236 num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
237 if (DFISZERO(dfl)) { /* 0/x (x!=0) */
238 /* if divide, result is 0 with ideal exponent; divideInt has */
239 /* exponent=0, remainders give zero with lower exponent */
241 decFloatZero(result);
242 DFWORD(result, 0)|=num.sign; /* add sign */
245 if (!(op&DIVIDE)) { /* a remainder */
246 /* exponent is the minimum of the operands */
247 num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
248 /* if the result is zero the sign shall be sign of dfl */
249 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
252 num.msd=bcdacc; /* -> 0 */
253 num.lsd=bcdacc; /* .. */
254 return decFinalize(result, &num, set); /* [divide may clamp exponent] */
256 /* [here, both operands are known to be finite and non-zero] */
258 /* extract the operand coefficents into 'units' which are */
259 /* base-billion; the lhs is high-aligned in acc and the msu of both */
260 /* acc and div is at the right-hand end of array (offset length-1); */
261 /* the quotient can need one more unit than the operands as digits */
262 /* in it are not necessarily aligned neatly; further, the quotient */
263 /* may not start accumulating until after the end of the initial */
264 /* operand in acc if that is small (e.g., 1) so the accumulator */
265 /* must have at least that number of units extra (at the ls end) */
266 GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
267 GETCOEFFBILL(dfr, div);
268 /* zero the low uInts of acc */
275 #error Unexpected Double DIVOPLEN
283 #error Unexpected Quad DIVOPLEN
287 /* set msu and lsu pointers */
288 msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
290 /*[loop for div will terminate because operands are non-zero] */
291 for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
292 /* the initial least-significant unit of acc is set so acc appears */
293 /* to have the same length as div. */
294 /* This moves one position towards the least possible for each */
296 divunits=(Int)(msud-div+1); /* precalculate */
297 lsua=msua-divunits+1; /* initial working lsu of acc */
298 lsuq=msuq; /* and of quo */
300 /* set up the estimator for the multiplier; this is the msu of div, */
301 /* plus two bits from the unit below (if any) rounded up by one if */
302 /* there are any non-zero bits or units below that [the extra two */
303 /* bits makes for a much better estimate when the top unit is small] */
308 if (d>=750000000) {divtop+=3; d-=750000000;}
309 else if (d>=500000000) {divtop+=2; d-=500000000;}
310 else if (d>=250000000) {divtop++; d-=250000000;}
312 else for (um--; um>=div; um--) if (*um) {
320 printf("----- div=");
321 for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
325 /* now collect up to DECPMAX+1 digits in the quotient (this may */
326 /* need OPLEN+1 uInts if unaligned) */
327 quodigits=0; /* no digits yet */
328 for (;; lsua--) { /* outer loop -- each input position */
331 printf("Acc underrun...\n");
336 printf("Outer: quodigits=%ld acc=", (LI)quodigits);
337 for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
340 *lsuq=0; /* default unit result is 0 */
341 for (;;) { /* inner loop -- calculate quotient unit */
342 /* strip leading zero units from acc (either there initially or */
343 /* from subtraction below); this may strip all if exactly 0 */
344 for (; *msua==0 && msua>=lsua;) msua--;
345 accunits=(Int)(msua-lsua+1); /* [maybe 0] */
346 /* subtraction is only necessary and possible if there are as */
347 /* least as many units remaining in acc for this iteration as */
348 /* there are in div */
349 if (accunits<divunits) {
350 if (accunits==0) msua++; /* restore */
354 /* If acc is longer than div then subtraction is definitely */
355 /* possible (as msu of both is non-zero), but if they are the */
356 /* same length a comparison is needed. */
357 /* If a subtraction is needed then a good estimate of the */
358 /* multiplier for the subtraction is also needed in order to */
359 /* minimise the iterations of this inner loop because the */
360 /* subtractions needed dominate division performance. */
361 if (accunits==divunits) {
362 /* compare the high divunits of acc and div: */
363 /* acc<div: this quotient unit is unchanged; subtraction */
364 /* will be possible on the next iteration */
365 /* acc==div: quotient gains 1, set acc=0 */
366 /* acc>div: subtraction necessary at this position */
367 for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
368 /* [now at first mismatch or lsu] */
369 if (*ud>*ua) break; /* next time... */
370 if (*ud==*ua) { /* all compared equal */
371 *lsuq+=1; /* increment result */
372 msua=lsua; /* collapse acc units */
373 *msua=0; /* .. to a zero */
377 /* subtraction necessary; estimate multiplier [see above] */
378 /* if both *msud and *msua are small it is cost-effective to */
379 /* bring in part of the following units (if any) to get a */
380 /* better estimate (assume some other non-zero in div) */
381 #define DIVLO 1000000U
382 #define DIVHI (DIVBASE/DIVLO)
385 /* there cannot be a *(msud-2) for DECDOUBLE so next is */
386 /* an exact calculation unless DECQUAD (which needs to */
387 /* assume bits out there if divunits>2) */
388 uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
389 uLong div=(uLong)*msud * DIVBASE + *(msud-1);
391 if (divunits>2) div++;
396 else multiplier=*msua/(*msud);
398 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
399 multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
400 /(*msud*DIVHI + *(msud-1)/DIVLO +1);
402 else multiplier=(*msua<<2)/divtop;
405 else { /* accunits>divunits */
406 /* msud is one unit 'lower' than msua, so estimate differently */
409 /* as before, bring in extra digits if possible */
410 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
411 mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
413 mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
415 else if (divunits==1) {
416 mul=(uLong)*msua * DIVBASE + *(msua-1);
417 mul/=*msud; /* no more to the right */
420 mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
422 mul/=divtop; /* [divtop already allows for sticky bits] */
426 multiplier=*msua * ((DIVBASE<<2)/divtop);
429 if (multiplier==0) multiplier=1; /* marginal case */
433 /* printf("Multiplier: %ld\n", (LI)multiplier); */
437 /* Carry out the subtraction acc-(div*multiplier); for each */
438 /* unit in div, do the multiply, split to units (see */
439 /* decFloatMultiply for the algorithm), and subtract from acc */
440 #define DIVMAGIC 2305843009U /* 2**61/10**9 */
444 for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
447 uLong sub=(uLong)multiplier*(*ud)+carry;
453 hop=(uInt)(sub>>DIVSHIFTA);
454 carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
455 /* the estimate is now in hi; now calculate sub-hi*10**9 */
456 /* to get the remainder (which will be <DIVBASE)) */
458 lo-=carry*DIVBASE; /* low word of result */
460 lo-=DIVBASE; /* correct by +1 */
466 /* calculate multiplier*(*ud) into hi and lo */
467 LONGMUL32HI(hi, *ud, multiplier); /* get the high word */
468 lo=multiplier*(*ud); /* .. and the low */
469 lo+=carry; /* add the old hi */
470 carry=hi+(lo<carry); /* .. with any carry */
471 if (carry || lo>=DIVBASE) { /* split is needed */
472 hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */
473 LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */
474 /* [DIVSHIFTB is 32, so carry can be used directly] */
475 /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
476 /* happily the top word of the result is irrelevant because it */
477 /* will always be zero so this needs only one multiplication */
479 /* the correction here will be at most +1; do it */
486 if (lo>*ua) { /* borrow needed */
492 if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */
495 /* the outer loop terminates when there is either an exact result */
496 /* or enough digits; first update the quotient digit count and */
497 /* pointer (if any significant digits) */
499 if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
502 quodigits+=9; /* had leading unit earlier */
504 if (quodigits>DECPMAX+1) break; /* have enough */
506 else if (*lsuq) { /* first quotient digits */
508 for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
510 /* [cannot have >DECPMAX+1 on first unit] */
513 if (*msua!=0) continue; /* not an exact result */
514 /* acc is zero iff used all of original units and zero down to lsua */
515 /* (must also continue to original lsu for correct quotient length) */
516 if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
517 for (; msua>lsua && *msua==0;) msua--;
518 if (*msua==0 && msua==lsua) break;
521 /* all of the original operand in acc has been covered at this point */
522 /* quotient now has at least DECPMAX+2 digits */
523 /* *msua is now non-0 if inexact and sticky bits */
524 /* lsuq is one below the last uint of the quotient */
525 lsuq++; /* set -> true lsu of quo */
526 if (*msua) *lsuq|=1; /* apply sticky bit */
528 /* quo now holds the (unrounded) quotient in base-billion; one */
529 /* base-billion 'digit' per uInt. */
532 for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
536 /* Now convert to BCD for rounding and cleanup, starting from the */
537 /* most significant end [offset by one into bcdacc to leave room */
538 /* for a possible carry digit if rounding for REMNEAR is needed] */
539 for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
540 uInt top, mid, rem; /* work */
541 if (*uq==0) { /* no split needed */
542 UBFROMUI(ub, 0); /* clear 9 BCD8s */
543 UBFROMUI(ub+4, 0); /* .. */
547 /* *uq is non-zero -- split the base-billion digit into */
548 /* hi, mid, and low three-digits */
549 #define divsplit9 1000000 /* divisor */
550 #define divsplit6 1000 /* divisor */
551 /* The splitting is done by simple divides and remainders, */
552 /* assuming the compiler will optimize these [GCC does] */
557 /* lay out the nine BCD digits (plus one unwanted byte) */
558 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
559 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
560 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
561 } /* BCD conversion loop */
564 /* complete the bcdnum; quodigits is correct, so the position of */
565 /* the first non-zero is known */
566 num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
569 /* make exponent adjustments, etc */
570 if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */
571 num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
572 /* if the result was exact then there may be up to 8 extra */
573 /* trailing zeros in the overflowed quotient final unit */
575 for (; *ub==0;) ub--; /* drop zeros */
576 num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */
579 } /* adjustment needed */
582 if (divcount>maxcount) { /* new high-water nark */
584 printf("DivNewMaxCount: %ld\n", (LI)maxcount);
588 if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
590 /* Is DIVIDEINT or a remainder; there is more to do -- first form */
591 /* the integer (this is done 'after the fact', unlike as in */
592 /* decNumber, so as not to tax DIVIDE) */
594 /* The first non-zero digit will be in the first 9 digits, known */
595 /* from quodigits and num.msd, so there is always space for DECPMAX */
598 length=(Int)(num.lsd-num.msd+1);
599 /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
601 if (length+num.exponent>DECPMAX) { /* cannot fit */
602 decFloatZero(result);
603 DFWORD(result, 0)=DECFLOAT_qNaN;
604 set->status|=DEC_Division_impossible;
608 if (num.exponent>=0) { /* already an int, or need pad zeros */
609 for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
610 num.lsd+=num.exponent;
612 else { /* too long: round or truncate needed */
613 Int drop=-num.exponent;
614 if (!(op&REMNEAR)) { /* simple truncate */
616 if (num.lsd<num.msd) { /* truncated all */
617 num.lsd=num.msd; /* make 0 */
618 *num.lsd=0; /* .. [sign still relevant] */
621 else { /* round to nearest even [sigh] */
622 /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
623 /* (this is a special case of Quantize -- q.v. for commentary) */
624 uByte *roundat; /* -> re-round digit */
625 uByte reround; /* reround value */
626 *(num.msd-1)=0; /* in case of left carry, or make 0 */
627 if (drop<length) roundat=num.lsd-drop+1;
628 else if (drop==length) roundat=num.msd;
629 else roundat=num.msd-1; /* [-> 0] */
631 for (ub=roundat+1; ub<=num.lsd; ub++) {
633 reround=DECSTICKYTAB[reround];
636 } /* check stickies */
637 if (roundat>num.msd) num.lsd=roundat-1;
639 num.msd--; /* use the 0 .. */
640 num.lsd=num.msd; /* .. at the new MSD place */
642 if (reround!=0) { /* discarding non-zero */
644 /* rounding is DEC_ROUND_HALF_EVEN always */
645 if (reround>5) bump=1; /* >0.5 goes up */
646 else if (reround==5) /* exactly 0.5000 .. */
647 bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */
648 if (bump!=0) { /* need increment */
649 /* increment the coefficient; this might end up with 1000... */
651 for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
652 for (; *ub==9; ub--) *ub=0; /* at most 3 more */
654 if (ub<num.msd) num.msd--; /* carried */
658 } /* round or truncate needed */
659 num.exponent=0; /* all paths */
660 /*decShowNum(&num, "int"); */
662 if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
664 /* Have a remainder to calculate */
665 decFinalize("ient, &num, set); /* lay out the integer so far */
666 DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */
667 sign=DFWORD(dfl, 0); /* save sign of dfl */
668 decFloatFMA(result, "ient, dfr, dfl, set);
669 if (!DFISZERO(result)) return result;
670 /* if the result is zero the sign shall be sign of dfl */
671 DFWORD("ient, 0)=sign; /* construct decFloat of sign */
672 return decFloatCopySign(result, result, "ient);
675 /* ------------------------------------------------------------------ */
676 /* decFiniteMultiply -- multiply two finite decFloats */
678 /* num gets the result of multiplying dfl and dfr */
679 /* bcdacc .. with the coefficient in this array */
680 /* dfl is the first decFloat (lhs) */
681 /* dfr is the second decFloat (rhs) */
683 /* This effects the multiplication of two decFloats, both known to be */
684 /* finite, leaving the result in a bcdnum ready for decFinalize (for */
685 /* use in Multiply) or in a following addition (FMA). */
687 /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
688 /* No error is possible and no status is set. */
689 /* ------------------------------------------------------------------ */
690 /* This routine has two separate implementations of the core */
691 /* multiplication; both using base-billion. One uses only 32-bit */
692 /* variables (Ints and uInts) or smaller; the other uses uLongs (for */
693 /* multiplication and addition only). Both implementations cover */
694 /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
695 /* comparisons. In any one compilation only one implementation for */
696 /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
697 /* version is forced. */
699 /* Historical note: an earlier version of this code also supported the */
700 /* 256-bit format and has been preserved. That is somewhat trickier */
701 /* during lazy carry splitting because the initial quotient estimate */
702 /* (est) can exceed 32 bits. */
704 #define MULTBASE ((uInt)BILLION) /* the base used for multiply */
705 #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
706 #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */
707 #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
709 /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
711 #error Exponent may overflow when doubled for Multiply
713 #if MULACCLEN!=(MULACCLEN/4)*4
714 /* This assumption is used below only for initialization */
715 #error MULACCLEN is not a multiple of 4
718 static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
719 const decFloat *dfl, const decFloat *dfr) {
720 uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
721 uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
722 uInt *ui, *uj; /* work */
724 uInt uiwork; /* for macros */
727 uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
728 uLong *pl; /* work -> lazy accumulator */
729 uInt acc[MULACCLEN]; /* coefficent in base-billion .. */
731 uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */
733 uInt *pa; /* work -> accumulator */
734 /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
736 /* Calculate sign and exponent */
737 num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
738 num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
740 /* Extract the coefficients and prepare the accumulator */
741 /* the coefficients of the operands are decoded into base-billion */
742 /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
743 /* appropriate size. */
744 GETCOEFFBILL(dfl, bufl);
745 GETCOEFFBILL(dfr, bufr);
748 for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
751 for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
755 /* start the 64-bit/32-bit differing paths... */
758 /* zero the accumulator */
760 accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
761 #else /* use a loop */
762 /* MULACCLEN is a multiple of four, asserted above */
763 for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
764 *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
768 /* Effect the multiplication */
769 /* The multiplcation proceeds using MFC's lazy-carry resolution */
770 /* algorithm from decNumber. First, the multiplication is */
771 /* effected, allowing accumulation of the partial products (which */
772 /* are in base-billion at each column position) into 64 bits */
773 /* without resolving back to base=billion after each addition. */
774 /* These 64-bit numbers (which may contain up to 19 decimal digits) */
775 /* are then split using the Clark & Cowlishaw algorithm (see below). */
776 /* [Testing for 0 in the inner loop is not really a 'win'] */
777 for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
778 if (*ui==0) continue; /* product cannot affect result */
779 pl=accl+(ui-bufr); /* where to add the lhs */
780 for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
781 /* if (*uj==0) continue; // product cannot affect result */
782 *pl+=((uLong)*ui)*(*uj);
786 /* The 64-bit carries must now be resolved; this means that a */
787 /* quotient/remainder has to be calculated for base-billion (1E+9). */
788 /* For this, Clark & Cowlishaw's quotient estimation approach (also */
789 /* used in decNumber) is needed, because 64-bit divide is generally */
790 /* extremely slow on 32-bit machines, and may be slower than this */
791 /* approach even on 64-bit machines. This algorithm splits X */
794 /* magic=2**(A+B)/1E+9; // 'magic number' */
795 /* hop=X/2**A; // high order part of X (by shift) */
796 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
798 /* A and B are quite constrained; hop and magic must fit in 32 bits, */
799 /* and 2**(A+B) must be as large as possible (which is 2**61 if */
800 /* magic is to fit). Further, maxX increases with the length of */
801 /* the operands (and hence the number of partial products */
802 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
804 /* It can be shown that when OPLEN is 2 then the maximum error in */
805 /* the estimated quotient is <1, but for larger maximum x the */
806 /* maximum error is above 1 so a correction that is >1 may be */
807 /* needed. Values of A and B are chosen to satisfy the constraints */
808 /* just mentioned while minimizing the maximum error (and hence the */
809 /* maximum correction), as shown in the following table: */
811 /* Type OPLEN A B maxX maxError maxCorrection */
812 /* --------------------------------------------------------- */
813 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
814 /* QUAD 4 30 31 <4*10**18 1.17 2 */
816 /* In the OPLEN==2 case there is most choice, but the value for B */
817 /* of 32 has a big advantage as then the calculation of the */
818 /* estimate requires no shifting; the compiler can extract the high */
819 /* word directly after multiplying magic*hop. */
820 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
828 #error Unexpected type
833 for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
834 printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
838 for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
839 uInt lo, hop; /* work */
840 uInt est; /* cannot exceed 4E+9 */
842 /* *pl holds a binary number which needs to be split */
843 hop=(uInt)(*pl>>MULSHIFTA);
844 est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
845 /* the estimate is now in est; now calculate hi:lo-est*10**9; */
846 /* happily the top word of the result is irrelevant because it */
847 /* will always be zero so this needs only one multiplication */
848 lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */
849 /* If QUAD, the correction here could be +2 */
851 lo-=MULTBASE; /* correct by +1 */
854 /* may need to correct by +2 */
861 /* finally place lo as the new coefficient 'digit' and add est to */
862 /* the next place up [this is safe because this path is never */
863 /* taken on the final iteration as *pl will fit] */
866 } /* *pl needed split */
867 else { /* *pl<MULTBASE */
868 *pa=(uInt)*pl; /* just copy across */
873 for (pa=acc;; pa+=4) { /* zero the accumulator */
874 *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */
875 if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */
878 /* Effect the multiplication */
879 /* uLongs are not available (and in particular, there is no uLong */
880 /* divide) but it is still possible to use MFC's lazy-carry */
881 /* resolution algorithm from decNumber. First, the multiplication */
882 /* is effected, allowing accumulation of the partial products */
883 /* (which are in base-billion at each column position) into 64 bits */
884 /* [with the high-order 32 bits in each position being held at */
885 /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
886 /* These 64-bit numbers (which may contain up to 19 decimal digits) */
887 /* are then split using the Clark & Cowlishaw algorithm (see */
889 for (ui=bufr;; ui++) { /* over each item in rhs */
890 uInt hi, lo; /* words of exact multiply result */
891 pa=acc+(ui-bufr); /* where to add the lhs */
892 for (uj=bufl;; uj++, pa++) { /* over each item in lhs */
893 LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */
894 lo=(*ui)*(*uj); /* .. */
895 *pa+=lo; /* accumulate low bits and .. */
896 *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */
897 if (uj==bufl+MULOPLEN-1) break;
899 if (ui==bufr+MULOPLEN-1) break;
902 /* The 64-bit carries must now be resolved; this means that a */
903 /* quotient/remainder has to be calculated for base-billion (1E+9). */
904 /* For this, Clark & Cowlishaw's quotient estimation approach (also */
905 /* used in decNumber) is needed, because 64-bit divide is generally */
906 /* extremely slow on 32-bit machines. This algorithm splits X */
909 /* magic=2**(A+B)/1E+9; // 'magic number' */
910 /* hop=X/2**A; // high order part of X (by shift) */
911 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
913 /* A and B are quite constrained; hop and magic must fit in 32 bits, */
914 /* and 2**(A+B) must be as large as possible (which is 2**61 if */
915 /* magic is to fit). Further, maxX increases with the length of */
916 /* the operands (and hence the number of partial products */
917 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
919 /* It can be shown that when OPLEN is 2 then the maximum error in */
920 /* the estimated quotient is <1, but for larger maximum x the */
921 /* maximum error is above 1 so a correction that is >1 may be */
922 /* needed. Values of A and B are chosen to satisfy the constraints */
923 /* just mentioned while minimizing the maximum error (and hence the */
924 /* maximum correction), as shown in the following table: */
926 /* Type OPLEN A B maxX maxError maxCorrection */
927 /* --------------------------------------------------------- */
928 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
929 /* QUAD 4 30 31 <4*10**18 1.17 2 */
931 /* In the OPLEN==2 case there is most choice, but the value for B */
932 /* of 32 has a big advantage as then the calculation of the */
933 /* estimate requires no shifting; the high word is simply */
934 /* calculated from multiplying magic*hop. */
935 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
943 #error Unexpected type
948 for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
949 printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
953 for (pa=acc;; pa++) { /* each low uInt */
954 uInt hi, lo; /* words of exact multiply result */
955 uInt hop, estlo; /* work */
961 hi=*(pa+MULACCLEN); /* top 32 bits */
962 /* hi and lo now hold a binary number which needs to be split */
965 hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */
966 LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
967 /* [MULSHIFTB is 32, so estlo can be used directly] */
968 /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
969 /* happily the top word of the result is irrelevant because it */
970 /* will always be zero so this needs only one multiplication */
971 lo-=(estlo*MULTBASE);
972 /* esthi=0; // high word is ignored below */
973 /* the correction here will be at most +1; do it */
979 hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */
980 LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
981 estlo=hop*MULMAGIC; /* .. so low word needed */
982 estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
983 /* esthi=0; // high word is ignored below */
984 lo-=(estlo*MULTBASE); /* as above */
985 /* the correction here could be +1 or +2 */
995 #error Unexpected type
998 /* finally place lo as the new accumulator digit and add est to */
999 /* the next place up; this latter add could cause a carry of 1 */
1000 /* to the high word of the next place */
1003 /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
1004 /* *(pa+1+MULACCLEN)+=esthi; */
1005 if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
1006 if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */
1010 /* At this point, whether using the 64-bit or the 32-bit paths, the */
1011 /* accumulator now holds the (unrounded) result in base-billion; */
1012 /* one base-billion 'digit' per uInt. */
1015 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
1019 /* Now convert to BCD for rounding and cleanup, starting from the */
1020 /* most significant end */
1022 if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
1023 else { /* >=1 word of leading zeros */
1024 num->msd=bcdacc; /* known leading zeros are gone */
1025 pa--; /* skip first word .. */
1026 for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
1028 for (ub=bcdacc;; pa--, ub+=9) {
1029 if (*pa!=0) { /* split(s) needed */
1030 uInt top, mid, rem; /* work */
1031 /* *pa is non-zero -- split the base-billion acc digit into */
1032 /* hi, mid, and low three-digits */
1033 #define mulsplit9 1000000 /* divisor */
1034 #define mulsplit6 1000 /* divisor */
1035 /* The splitting is done by simple divides and remainders, */
1036 /* assuming the compiler will optimize these where useful */
1042 /* lay out the nine BCD digits (plus one unwanted byte) */
1043 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
1044 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
1045 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
1048 UBFROMUI(ub, 0); /* clear 9 BCD8s */
1049 UBFROMUI(ub+4, 0); /* .. */
1053 } /* BCD conversion loop */
1055 num->lsd=ub+8; /* complete the bcdnum .. */
1058 decShowNum(num, "postmult");
1059 decFloatShow(dfl, "dfl");
1060 decFloatShow(dfr, "dfr");
1063 } /* decFiniteMultiply */
1065 /* ------------------------------------------------------------------ */
1066 /* decFloatAbs -- absolute value, heeding NaNs, etc. */
1068 /* result gets the canonicalized df with sign 0 */
1069 /* df is the decFloat to abs */
1070 /* set is the context */
1071 /* returns result */
1073 /* This has the same effect as decFloatPlus unless df is negative, */
1074 /* in which case it has the same effect as decFloatMinus. The */
1075 /* effect is also the same as decFloatCopyAbs except that NaNs are */
1076 /* handled normally (the sign of a NaN is not affected, and an sNaN */
1077 /* will signal) and the result will be canonical. */
1078 /* ------------------------------------------------------------------ */
1079 decFloat * decFloatAbs(decFloat *result, const decFloat *df,
1081 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
1082 decCanonical(result, df); /* copy and check */
1083 DFBYTE(result, 0)&=~0x80; /* zero sign bit */
1087 /* ------------------------------------------------------------------ */
1088 /* decFloatAdd -- add two decFloats */
1090 /* result gets the result of adding dfl and dfr: */
1091 /* dfl is the first decFloat (lhs) */
1092 /* dfr is the second decFloat (rhs) */
1093 /* set is the context */
1094 /* returns result */
1096 /* ------------------------------------------------------------------ */
1098 /* Table for testing MSDs for fastpath elimination; returns the MSD of */
1099 /* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
1100 /* Infinities return -32 and NaNs return -128 so that summing the two */
1101 /* MSDs also allows rapid tests for the Specials (see code below). */
1102 const Int DECTESTMSD[64]={
1103 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
1104 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
1105 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
1106 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
1108 /* The table for testing MSDs is shared between the modules */
1109 extern const Int DECTESTMSD[64];
1112 decFloat * decFloatAdd(decFloat *result,
1113 const decFloat *dfl, const decFloat *dfr,
1115 bcdnum num; /* for final conversion */
1116 Int bexpl, bexpr; /* left and right biased exponents */
1117 uByte *ub, *us, *ut; /* work */
1118 uInt uiwork; /* for macros */
1120 uShort uswork; /* .. */
1123 uInt sourhil, sourhir; /* top words from source decFloats */
1124 /* [valid only through end of */
1125 /* fastpath code -- before swap] */
1126 uInt diffsign; /* non-zero if signs differ */
1127 uInt carry; /* carry: 0 or 1 before add loop */
1128 Int overlap; /* coefficient overlap (if full) */
1129 Int summ; /* sum of the MSDs */
1130 /* the following buffers hold coefficients with various alignments */
1131 /* (see commentary and diagrams below) */
1132 uByte acc[4+2+DECPMAX*3+8];
1133 uByte buf[4+2+DECPMAX*2];
1134 uByte *umsd, *ulsd; /* local MSD and LSD pointers */
1137 #define CARRYPAT 0x01000000 /* carry=1 pattern */
1139 #define CARRYPAT 0x00000001 /* carry=1 pattern */
1142 /* Start decoding the arguments */
1143 /* The initial exponents are placed into the opposite Ints to */
1144 /* that which might be expected; there are two sets of data to */
1145 /* keep track of (each decFloat and the corresponding exponent), */
1146 /* and this scheme means that at the swap point (after comparing */
1147 /* exponents) only one pair of words needs to be swapped */
1148 /* whichever path is taken (thereby minimising worst-case path). */
1149 /* The calculated exponents will be nonsense when the arguments are */
1150 /* Special, but are not used in that path */
1151 sourhil=DFWORD(dfl, 0); /* LHS top word */
1152 summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */
1153 bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
1154 bexpr+=GETECON(dfl); /* .. + continuation */
1156 sourhir=DFWORD(dfr, 0); /* RHS top word */
1157 summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */
1158 bexpl=DECCOMBEXP[sourhir>>26];
1159 bexpl+=GETECON(dfr);
1161 /* here bexpr has biased exponent from lhs, and vice versa */
1163 diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
1165 /* now determine whether to take a fast path or the full-function */
1166 /* slow path. The slow path must be taken when: */
1167 /* -- both numbers are finite, and: */
1168 /* the exponents are different, or */
1169 /* the signs are different, or */
1170 /* the sum of the MSDs is >8 (hence might overflow) */
1171 /* specialness and the sum of the MSDs can be tested at once using */
1172 /* the summ value just calculated, so the test for specials is no */
1173 /* longer on the worst-case path (as of 3.60) */
1175 if (summ<=8) { /* MSD+MSD is good, or there is a special */
1176 if (summ<0) { /* there is a special */
1177 /* Inf+Inf would give -64; Inf+finite is -32 or higher */
1178 if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */
1179 /* two infinities with different signs is invalid */
1180 if (summ==-64 && diffsign) return decInvalid(result, set);
1181 if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
1182 return decInfinity(result, dfr); /* RHS must be Inf */
1184 /* Here when both arguments are finite; fast path is possible */
1185 /* (currently only for aligned and same-sign) */
1186 if (bexpr==bexpl && !diffsign) {
1187 uInt tac[DECLETS+1]; /* base-1000 coefficient */
1188 uInt encode; /* work */
1190 /* Get one coefficient as base-1000 and add the other */
1191 GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */
1192 ADDCOEFFTHOU(dfr, tac);
1193 /* here the sum of the MSDs (plus any carry) will be <10 due to */
1194 /* the fastpath test earlier */
1196 /* construct the result; low word is the same for both formats */
1197 encode =BIN2DPD[tac[0]];
1198 encode|=BIN2DPD[tac[1]]<<10;
1199 encode|=BIN2DPD[tac[2]]<<20;
1200 encode|=BIN2DPD[tac[3]]<<30;
1201 DFWORD(result, (DECBYTES/4)-1)=encode;
1203 /* collect next two declets (all that remains, for Double) */
1204 encode =BIN2DPD[tac[3]]>>2;
1205 encode|=BIN2DPD[tac[4]]<<8;
1208 /* complete and lay out middling words */
1209 encode|=BIN2DPD[tac[5]]<<18;
1210 encode|=BIN2DPD[tac[6]]<<28;
1211 DFWORD(result, 2)=encode;
1213 encode =BIN2DPD[tac[6]]>>4;
1214 encode|=BIN2DPD[tac[7]]<<6;
1215 encode|=BIN2DPD[tac[8]]<<16;
1216 encode|=BIN2DPD[tac[9]]<<26;
1217 DFWORD(result, 1)=encode;
1219 /* and final two declets */
1220 encode =BIN2DPD[tac[9]]>>6;
1221 encode|=BIN2DPD[tac[10]]<<4;
1224 /* add exponent continuation and sign (from either argument) */
1225 encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
1227 /* create lookup index = MSD + top two bits of biased exponent <<4 */
1228 tac[DECLETS]|=(bexpl>>DECECONL)<<4;
1229 encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
1230 DFWORD(result, 0)=encode; /* complete */
1232 /* decFloatShow(result, ">"); */
1234 } /* fast path OK */
1235 /* drop through to slow path */
1236 } /* low sum or Special(s) */
1238 /* Slow path required -- arguments are finite and might overflow, */
1239 /* or require alignment, or might have different signs */
1241 /* now swap either exponents or argument pointers */
1243 /* original left is bigger */
1247 /* printf("left bigger\n"); */
1250 const decFloat *dfswap=dfl;
1253 /* printf("right bigger\n"); */
1255 /* [here dfl and bexpl refer to the datum with the larger exponent, */
1256 /* of if the exponents are equal then the original LHS argument] */
1258 /* if lhs is zero then result will be the rhs (now known to have */
1259 /* the smaller exponent), which also may need to be tested for zero */
1260 /* for the weird IEEE 754 sign rules */
1261 if (DFISZERO(dfl)) {
1262 decCanonical(result, dfr); /* clean copy */
1263 /* "When the sum of two operands with opposite signs is */
1264 /* exactly zero, the sign of that sum shall be '+' in all */
1265 /* rounding modes except round toward -Infinity, in which */
1266 /* mode that sign shall be '-'." */
1267 if (diffsign && DFISZERO(result)) {
1268 DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */
1269 if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
1272 } /* numfl is zero */
1273 /* [here, LHS is non-zero; code below assumes that] */
1275 /* Coefficients layout during the calculations to follow: */
1278 /* +------------------------------------------------+ */
1279 /* acc: |0000| coeffa | tail B | | */
1280 /* +------------------------------------------------+ */
1281 /* buf: |0000| pad0s | coeffb | | */
1282 /* +------------------------------------------------+ */
1284 /* Touching coefficients or gap: */
1285 /* +------------------------------------------------+ */
1286 /* acc: |0000| coeffa | gap | coeffb | */
1287 /* +------------------------------------------------+ */
1288 /* [buf not used or needed; gap clamped to Pmax] */
1290 /* lay out lhs coefficient into accumulator; this starts at acc+4 */
1291 /* for decDouble or acc+6 for decQuad so the LSD is word- */
1292 /* aligned; the top word gap is there only in case a carry digit */
1293 /* is prefixed after the add -- it does not need to be zeroed */
1295 #define COFF 4 /* offset into acc */
1297 UBFROMUS(acc+4, 0); /* prefix 00 */
1298 #define COFF 6 /* offset into acc */
1301 GETCOEFF(dfl, acc+COFF); /* decode from decFloat */
1302 ulsd=acc+COFF+DECPMAX-1;
1303 umsd=acc+4; /* [having this here avoids */
1309 tum.exponent=bexpl-DECBIAS;
1310 tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1311 decShowNum(&tum, "dflx");}
1314 /* if signs differ, take ten's complement of lhs (here the */
1315 /* coefficient is subtracted from all-nines; the 1 is added during */
1316 /* the later add cycle -- zeros to the right do not matter because */
1317 /* the complement of zero is zero); these are fixed-length inverts */
1318 /* where the lsd is known to be at a 4-byte boundary (so no borrow */
1320 carry=0; /* assume no carry */
1322 carry=CARRYPAT; /* for +1 during add */
1323 UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
1324 UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
1325 UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
1326 UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
1328 UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
1329 UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
1330 UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
1331 UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
1332 UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
1336 /* now process the rhs coefficient; if it cannot overlap lhs then */
1337 /* it can be put straight into acc (with an appropriate gap, if */
1338 /* needed) because no actual addition will be needed (except */
1339 /* possibly to complete ten's complement) */
1340 overlap=DECPMAX-(bexpl-bexpr);
1342 printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
1343 printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
1346 if (overlap<=0) { /* no overlap possible */
1347 uInt gap; /* local work */
1348 /* since a full addition is not needed, a ten's complement */
1349 /* calculation started above may need to be completed */
1351 for (ub=ulsd; *ub==9; ub--) *ub=0;
1353 carry=0; /* taken care of */
1355 /* up to DECPMAX-1 digits of the final result can extend down */
1356 /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
1357 /* rhs will be simply sticky bits. In this case the gap is */
1358 /* clamped to DECPMAX and the exponent adjusted to suit [this is */
1359 /* safe because the lhs is non-zero]. */
1365 ub=ulsd+gap+1; /* where MSD will go */
1366 /* Fill the gap with 0s; note that there is no addition to do */
1367 ut=acc+COFF+DECPMAX; /* start of gap */
1368 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */
1369 if (overlap<-DECPMAX) { /* gap was > DECPMAX */
1370 *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
1372 else { /* need full coefficient */
1373 GETCOEFF(dfr, ub); /* decode from decFloat */
1374 ub+=DECPMAX-1; /* new LSD... */
1376 ulsd=ub; /* save new LSD */
1377 } /* no overlap possible */
1379 else { /* overlap>0 */
1380 /* coefficients overlap (perhaps completely, although also */
1381 /* perhaps only where zeros) */
1382 if (overlap==DECPMAX) { /* aligned */
1383 ub=buf+COFF; /* where msd will go */
1385 UBFROMUS(buf+4, 0); /* clear quad's 00 */
1387 GETCOEFF(dfr, ub); /* decode from decFloat */
1389 else { /* unaligned */
1390 ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */
1391 /* Fill the prefix gap with 0s; 8 will cover most common */
1392 /* unalignments, so start with direct assignments (a loop is */
1393 /* then used for any remaining -- the loop (and the one in a */
1394 /* moment) is not then on the critical path because the number */
1395 /* of additions is reduced by (at least) two in this case) */
1396 UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */
1399 ut=buf+12; /* start any remaining */
1400 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
1402 GETCOEFF(dfr, ub); /* decode from decFloat */
1404 /* now move tail of rhs across to main acc; again use direct */
1405 /* copies for 8 digits-worth */
1406 UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX));
1407 UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
1408 if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
1409 us=buf+COFF+DECPMAX+8; /* source */
1410 ut=acc+COFF+DECPMAX+8; /* target */
1411 for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
1415 ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
1417 /* Now do the add of the non-tail; this is all nicely aligned, */
1418 /* and is over a multiple of four digits (because for Quad two */
1419 /* zero digits were added on the left); words in both acc and */
1420 /* buf (buf especially) will often be zero */
1421 /* [byte-by-byte add, here, is about 15% slower total effect than */
1424 /* Now effect the add; this is harder on a little-endian */
1425 /* machine as the inter-digit carry cannot use the usual BCD */
1426 /* addition trick because the bytes are loaded in the wrong order */
1427 /* [this loop could be unrolled, but probably scarcely worth it] */
1429 ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */
1430 us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */
1433 for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */
1435 carry+=UBTOUI(us); /* rhs + carry */
1436 if (carry==0) continue; /* no-op */
1437 carry+=UBTOUI(ut); /* lhs */
1438 /* Big-endian BCD adjust (uses internal carry) */
1439 carry+=0x76f6f6f6; /* note top nibble not all bits */
1440 /* apply BCD adjust and save */
1441 UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
1442 carry>>=31; /* true carry was at far left */
1445 for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */
1447 carry+=UBTOUI(us); /* rhs + carry */
1448 if (carry==0) continue; /* no-op [common if unaligned] */
1449 carry+=UBTOUI(ut); /* lhs */
1450 /* Little-endian BCD adjust; inter-digit carry must be manual */
1451 /* because the lsb from the array will be in the most-significant */
1453 carry+=0x76767676; /* note no inter-byte carries */
1454 carry+=(carry & 0x80000000)>>15;
1455 carry+=(carry & 0x00800000)>>15;
1456 carry+=(carry & 0x00008000)>>15;
1457 carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
1458 UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */
1459 /* here, final carry-out bit is at 0x00000080; move it ready */
1460 /* for next word-add (i.e., to 0x01000000) */
1461 carry=(carry & 0x00000080)<<17;
1467 printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
1468 tum.msd=umsd; /* acc+4; */
1472 decShowNum(&tum, "dfadd");}
1474 } /* overlap possible */
1476 /* ordering here is a little strange in order to have slowest path */
1477 /* first in GCC asm listing */
1478 if (diffsign) { /* subtraction */
1479 if (!carry) { /* no carry out means RHS<LHS */
1480 /* borrowed -- take ten's complement */
1481 /* sign is lhs sign */
1482 num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1484 /* invert the coefficient first by fours, then add one; space */
1485 /* at the end of the buffer ensures the by-fours is always */
1486 /* safe, but lsd+1 must be cleared to prevent a borrow */
1491 /* there are always at least four coefficient words */
1492 UBFROMUI(umsd, 0x09090909-UBTOUI(umsd));
1493 UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4));
1494 UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8));
1495 UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
1499 UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
1500 UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
1501 UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
1502 UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
1503 UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32));
1506 if (ulsd>=umsd+BNEXT) { /* unaligned */
1507 /* eight will handle most unaligments for Double; 16 for Quad */
1508 UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT));
1509 UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
1511 #define BNEXTY (BNEXT+8)
1513 UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8));
1514 UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
1515 #define BNEXTY (BNEXT+16)
1517 if (ulsd>=umsd+BNEXTY) { /* very unaligned */
1518 ut=umsd+BNEXTY; /* -> continue */
1520 UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
1521 if (ut>=ulsd-3) break; /* all done */
1525 /* complete the ten's complement by adding 1 */
1526 for (ub=ulsd; *ub==9; ub--) *ub=0;
1530 else { /* carry out means RHS>=LHS */
1531 num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
1532 /* all done except for the special IEEE 754 exact-zero-result */
1533 /* rule (see above); while testing for zero, strip leading */
1534 /* zeros (which will save decFinalize doing it) (this is in */
1535 /* diffsign path, so carry impossible and true umsd is */
1538 /* Check the initial coefficient area using the fast macro; */
1539 /* this will often be all that needs to be done (as on the */
1540 /* worst-case path when the subtraction was aligned and */
1542 if (ISCOEFFZERO(acc+COFF)) {
1543 umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
1544 if (ulsd>umsd) { /* more to check */
1545 umsd++; /* to align after checked area */
1546 for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
1547 for (; *umsd==0 && umsd<ulsd;) umsd++;
1549 if (*umsd==0) { /* must be true zero (and diffsign) */
1550 num.sign=0; /* assume + */
1551 if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
1554 /* [else was not zero, might still have leading zeros] */
1555 } /* subtraction gave positive result */
1558 else { /* same-sign addition */
1559 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
1561 if (carry) { /* only possible with decDouble */
1562 *(acc+3)=1; /* [Quad has leading 00] */
1568 num.msd=umsd; /* set MSD .. */
1569 num.lsd=ulsd; /* .. and LSD */
1570 num.exponent=bexpr-DECBIAS; /* set exponent to smaller, unbiassed */
1573 decFloatShow(dfl, "dfl");
1574 decFloatShow(dfr, "dfr");
1575 decShowNum(&num, "postadd");
1577 return decFinalize(result, &num, set); /* round, check, and lay out */
1580 /* ------------------------------------------------------------------ */
1581 /* decFloatAnd -- logical digitwise AND of two decFloats */
1583 /* result gets the result of ANDing dfl and dfr */
1584 /* dfl is the first decFloat (lhs) */
1585 /* dfr is the second decFloat (rhs) */
1586 /* set is the context */
1587 /* returns result, which will be canonical with sign=0 */
1589 /* The operands must be positive, finite with exponent q=0, and */
1590 /* comprise just zeros and ones; if not, Invalid operation results. */
1591 /* ------------------------------------------------------------------ */
1592 decFloat * decFloatAnd(decFloat *result,
1593 const decFloat *dfl, const decFloat *dfr,
1595 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
1596 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
1597 /* the operands are positive finite integers (q=0) with just 0s and 1s */
1599 DFWORD(result, 0)=ZEROWORD
1600 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
1601 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
1603 DFWORD(result, 0)=ZEROWORD
1604 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
1605 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
1606 DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
1607 DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
1612 /* ------------------------------------------------------------------ */
1613 /* decFloatCanonical -- copy a decFloat, making canonical */
1615 /* result gets the canonicalized df */
1616 /* df is the decFloat to copy and make canonical */
1617 /* returns result */
1619 /* This works on specials, too; no error or exception is possible. */
1620 /* ------------------------------------------------------------------ */
1621 decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
1622 return decCanonical(result, df);
1623 } /* decFloatCanonical */
1625 /* ------------------------------------------------------------------ */
1626 /* decFloatClass -- return the class of a decFloat */
1628 /* df is the decFloat to test */
1629 /* returns the decClass that df falls into */
1630 /* ------------------------------------------------------------------ */
1631 enum decClass decFloatClass(const decFloat *df) {
1632 Int exp; /* exponent */
1633 if (DFISSPECIAL(df)) {
1634 if (DFISQNAN(df)) return DEC_CLASS_QNAN;
1635 if (DFISSNAN(df)) return DEC_CLASS_SNAN;
1636 /* must be an infinity */
1637 if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
1638 return DEC_CLASS_POS_INF;
1640 if (DFISZERO(df)) { /* quite common */
1641 if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
1642 return DEC_CLASS_POS_ZERO;
1644 /* is finite and non-zero; similar code to decFloatIsNormal, here */
1645 /* [this could be speeded up slightly by in-lining decFloatDigits] */
1646 exp=GETEXPUN(df) /* get unbiased exponent .. */
1647 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
1648 if (exp>=DECEMIN) { /* is normal */
1649 if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
1650 return DEC_CLASS_POS_NORMAL;
1653 if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
1654 return DEC_CLASS_POS_SUBNORMAL;
1655 } /* decFloatClass */
1657 /* ------------------------------------------------------------------ */
1658 /* decFloatClassString -- return the class of a decFloat as a string */
1660 /* df is the decFloat to test */
1661 /* returns a constant string describing the class df falls into */
1662 /* ------------------------------------------------------------------ */
1663 const char *decFloatClassString(const decFloat *df) {
1664 enum decClass eclass=decFloatClass(df);
1665 if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
1666 if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
1667 if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
1668 if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
1669 if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
1670 if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
1671 if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
1672 if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
1673 if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
1674 if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
1675 return DEC_ClassString_UN; /* Unknown */
1676 } /* decFloatClassString */
1678 /* ------------------------------------------------------------------ */
1679 /* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
1681 /* result gets the result of comparing dfl and dfr */
1682 /* dfl is the first decFloat (lhs) */
1683 /* dfr is the second decFloat (rhs) */
1684 /* set is the context */
1685 /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
1686 /* ------------------------------------------------------------------ */
1687 decFloat * decFloatCompare(decFloat *result,
1688 const decFloat *dfl, const decFloat *dfr,
1690 Int comp; /* work */
1691 /* NaNs are handled as usual */
1692 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1693 /* numeric comparison needed */
1694 comp=decNumCompare(dfl, dfr, 0);
1695 decFloatZero(result);
1696 if (comp==0) return result;
1697 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1698 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1700 } /* decFloatCompare */
1702 /* ------------------------------------------------------------------ */
1703 /* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
1705 /* result gets the result of comparing dfl and dfr */
1706 /* dfl is the first decFloat (lhs) */
1707 /* dfr is the second decFloat (rhs) */
1708 /* set is the context */
1709 /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
1710 /* ------------------------------------------------------------------ */
1711 decFloat * decFloatCompareSignal(decFloat *result,
1712 const decFloat *dfl, const decFloat *dfr,
1714 Int comp; /* work */
1715 /* NaNs are handled as usual, except that all NaNs signal */
1716 if (DFISNAN(dfl) || DFISNAN(dfr)) {
1717 set->status|=DEC_Invalid_operation;
1718 return decNaNs(result, dfl, dfr, set);
1720 /* numeric comparison needed */
1721 comp=decNumCompare(dfl, dfr, 0);
1722 decFloatZero(result);
1723 if (comp==0) return result;
1724 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1725 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1727 } /* decFloatCompareSignal */
1729 /* ------------------------------------------------------------------ */
1730 /* decFloatCompareTotal -- compare two decFloats with total ordering */
1732 /* result gets the result of comparing dfl and dfr */
1733 /* dfl is the first decFloat (lhs) */
1734 /* dfr is the second decFloat (rhs) */
1735 /* returns result, which may be -1, 0, or 1 */
1736 /* ------------------------------------------------------------------ */
1737 decFloat * decFloatCompareTotal(decFloat *result,
1738 const decFloat *dfl, const decFloat *dfr) {
1739 Int comp; /* work */
1740 uInt uiwork; /* for macros */
1742 uShort uswork; /* .. */
1744 if (DFISNAN(dfl) || DFISNAN(dfr)) {
1745 Int nanl, nanr; /* work */
1746 /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
1747 nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
1748 if (DFISSIGNED(dfl)) nanl=-nanl;
1749 nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
1750 if (DFISSIGNED(dfr)) nanr=-nanr;
1751 if (nanl>nanr) comp=+1;
1752 else if (nanl<nanr) comp=-1;
1753 else { /* NaNs are the same type and sign .. must compare payload */
1754 /* buffers need +2 for QUAD */
1755 uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */
1756 uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */
1757 uByte *ub, *uc; /* work */
1758 Int sigl; /* signum of LHS */
1759 sigl=(DFISSIGNED(dfl) ? -1 : +1);
1761 /* decode the coefficients */
1762 /* (shift both right two if Quad to make a multiple of four) */
1767 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
1768 GETCOEFF(dfr, bufr+QUAD*2); /* .. */
1769 /* all multiples of four, here */
1770 comp=0; /* assume equal */
1771 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
1773 if (ui==UBTOUI(uc)) continue; /* so far so same */
1774 /* about to find a winner; go by bytes in case little-endian */
1775 for (;; ub++, uc++) {
1776 if (*ub==*uc) continue;
1777 if (*ub>*uc) comp=sigl; /* difference found */
1778 else comp=-sigl; /* .. */
1782 } /* same NaN type and sign */
1785 /* numeric comparison needed */
1786 comp=decNumCompare(dfl, dfr, 1); /* total ordering */
1788 decFloatZero(result);
1789 if (comp==0) return result;
1790 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1791 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1793 } /* decFloatCompareTotal */
1795 /* ------------------------------------------------------------------ */
1796 /* decFloatCompareTotalMag -- compare magnitudes with total ordering */
1798 /* result gets the result of comparing abs(dfl) and abs(dfr) */
1799 /* dfl is the first decFloat (lhs) */
1800 /* dfr is the second decFloat (rhs) */
1801 /* returns result, which may be -1, 0, or 1 */
1802 /* ------------------------------------------------------------------ */
1803 decFloat * decFloatCompareTotalMag(decFloat *result,
1804 const decFloat *dfl, const decFloat *dfr) {
1805 decFloat a, b; /* for copy if needed */
1806 /* copy and redirect signed operand(s) */
1807 if (DFISSIGNED(dfl)) {
1808 decFloatCopyAbs(&a, dfl);
1811 if (DFISSIGNED(dfr)) {
1812 decFloatCopyAbs(&b, dfr);
1815 return decFloatCompareTotal(result, dfl, dfr);
1816 } /* decFloatCompareTotalMag */
1818 /* ------------------------------------------------------------------ */
1819 /* decFloatCopy -- copy a decFloat as-is */
1821 /* result gets the copy of dfl */
1822 /* dfl is the decFloat to copy */
1823 /* returns result */
1825 /* This is a bitwise operation; no errors or exceptions are possible. */
1826 /* ------------------------------------------------------------------ */
1827 decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
1828 if (dfl!=result) *result=*dfl; /* copy needed */
1830 } /* decFloatCopy */
1832 /* ------------------------------------------------------------------ */
1833 /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
1835 /* result gets the copy of dfl with sign bit 0 */
1836 /* dfl is the decFloat to copy */
1837 /* returns result */
1839 /* This is a bitwise operation; no errors or exceptions are possible. */
1840 /* ------------------------------------------------------------------ */
1841 decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
1842 if (dfl!=result) *result=*dfl; /* copy needed */
1843 DFBYTE(result, 0)&=~0x80; /* zero sign bit */
1845 } /* decFloatCopyAbs */
1847 /* ------------------------------------------------------------------ */
1848 /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
1850 /* result gets the copy of dfl with sign bit inverted */
1851 /* dfl is the decFloat to copy */
1852 /* returns result */
1854 /* This is a bitwise operation; no errors or exceptions are possible. */
1855 /* ------------------------------------------------------------------ */
1856 decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
1857 if (dfl!=result) *result=*dfl; /* copy needed */
1858 DFBYTE(result, 0)^=0x80; /* invert sign bit */
1860 } /* decFloatCopyNegate */
1862 /* ------------------------------------------------------------------ */
1863 /* decFloatCopySign -- copy a decFloat with the sign of another */
1865 /* result gets the result of copying dfl with the sign of dfr */
1866 /* dfl is the first decFloat (lhs) */
1867 /* dfr is the second decFloat (rhs) */
1868 /* returns result */
1870 /* This is a bitwise operation; no errors or exceptions are possible. */
1871 /* ------------------------------------------------------------------ */
1872 decFloat * decFloatCopySign(decFloat *result,
1873 const decFloat *dfl, const decFloat *dfr) {
1874 uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
1875 if (dfl!=result) *result=*dfl; /* copy needed */
1876 DFBYTE(result, 0)&=~0x80; /* clear sign .. */
1877 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
1879 } /* decFloatCopySign */
1881 /* ------------------------------------------------------------------ */
1882 /* decFloatDigits -- return the number of digits in a decFloat */
1884 /* df is the decFloat to investigate */
1885 /* returns the number of significant digits in the decFloat; a */
1886 /* zero coefficient returns 1 as does an infinity (a NaN returns */
1887 /* the number of digits in the payload) */
1888 /* ------------------------------------------------------------------ */
1889 /* private macro to extract a declet according to provided formula */
1890 /* (form), and if it is non-zero then return the calculated digits */
1891 /* depending on the declet number (n), where n=0 for the most */
1892 /* significant declet; uses uInt dpd for work */
1893 #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
1894 if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1895 /* next one is used when it is known that the declet must be */
1896 /* non-zero, or is the final zero declet */
1897 #define dpdlendun(n, form) {dpd=(form)&0x3ff; \
1898 if (dpd==0) return 1; \
1899 return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1901 uInt decFloatDigits(const decFloat *df) {
1902 uInt dpd; /* work */
1903 uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
1905 uInt sourmh, sourml;
1909 if (DFISINF(df)) return 1;
1910 /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
1911 /* then the coefficient is full-length */
1912 if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
1915 if (sourhi&0x0003ffff) { /* ends in first */
1916 dpdlenchk(0, sourhi>>8);
1917 sourlo=DFWORD(df, 1);
1918 dpdlendun(1, (sourhi<<2) | (sourlo>>30));
1919 } /* [cannot drop through] */
1920 sourlo=DFWORD(df, 1); /* sourhi not involved now */
1921 if (sourlo&0xfff00000) { /* in one of first two */
1922 dpdlenchk(1, sourlo>>30); /* very rare */
1923 dpdlendun(2, sourlo>>20);
1924 } /* [cannot drop through] */
1925 dpdlenchk(3, sourlo>>10);
1926 dpdlendun(4, sourlo);
1927 /* [cannot drop through] */
1930 if (sourhi&0x00003fff) { /* ends in first */
1931 dpdlenchk(0, sourhi>>4);
1932 sourmh=DFWORD(df, 1);
1933 dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
1934 } /* [cannot drop through] */
1935 sourmh=DFWORD(df, 1);
1937 dpdlenchk(1, sourmh>>26);
1938 dpdlenchk(2, sourmh>>16);
1939 dpdlenchk(3, sourmh>>6);
1940 sourml=DFWORD(df, 2);
1941 dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
1942 } /* [cannot drop through] */
1943 sourml=DFWORD(df, 2);
1945 dpdlenchk(4, sourml>>28);
1946 dpdlenchk(5, sourml>>18);
1947 dpdlenchk(6, sourml>>8);
1948 sourlo=DFWORD(df, 3);
1949 dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
1950 } /* [cannot drop through] */
1951 sourlo=DFWORD(df, 3);
1952 if (sourlo&0xfff00000) { /* in one of first two */
1953 dpdlenchk(7, sourlo>>30); /* very rare */
1954 dpdlendun(8, sourlo>>20);
1955 } /* [cannot drop through] */
1956 dpdlenchk(9, sourlo>>10);
1957 dpdlendun(10, sourlo);
1958 /* [cannot drop through] */
1960 } /* decFloatDigits */
1962 /* ------------------------------------------------------------------ */
1963 /* decFloatDivide -- divide a decFloat by another */
1965 /* result gets the result of dividing dfl by dfr: */
1966 /* dfl is the first decFloat (lhs) */
1967 /* dfr is the second decFloat (rhs) */
1968 /* set is the context */
1969 /* returns result */
1971 /* ------------------------------------------------------------------ */
1972 /* This is just a wrapper. */
1973 decFloat * decFloatDivide(decFloat *result,
1974 const decFloat *dfl, const decFloat *dfr,
1976 return decDivide(result, dfl, dfr, set, DIVIDE);
1977 } /* decFloatDivide */
1979 /* ------------------------------------------------------------------ */
1980 /* decFloatDivideInteger -- integer divide a decFloat by another */
1982 /* result gets the result of dividing dfl by dfr: */
1983 /* dfl is the first decFloat (lhs) */
1984 /* dfr is the second decFloat (rhs) */
1985 /* set is the context */
1986 /* returns result */
1988 /* ------------------------------------------------------------------ */
1989 decFloat * decFloatDivideInteger(decFloat *result,
1990 const decFloat *dfl, const decFloat *dfr,
1992 return decDivide(result, dfl, dfr, set, DIVIDEINT);
1993 } /* decFloatDivideInteger */
1995 /* ------------------------------------------------------------------ */
1996 /* decFloatFMA -- multiply and add three decFloats, fused */
1998 /* result gets the result of (dfl*dfr)+dff with a single rounding */
1999 /* dfl is the first decFloat (lhs) */
2000 /* dfr is the second decFloat (rhs) */
2001 /* dff is the final decFloat (fhs) */
2002 /* set is the context */
2003 /* returns result */
2005 /* ------------------------------------------------------------------ */
2006 decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
2007 const decFloat *dfr, const decFloat *dff,
2010 /* The accumulator has the bytes needed for FiniteMultiply, plus */
2011 /* one byte to the left in case of carry, plus DECPMAX+2 to the */
2012 /* right for the final addition (up to full fhs + round & sticky) */
2013 #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
2014 uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
2015 /* .. and for final result */
2016 bcdnum mul; /* for multiplication result */
2017 bcdnum fin; /* for final operand, expanded */
2018 uByte coe[ROUNDUP4(DECPMAX)]; /* dff coefficient in BCD */
2019 bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
2020 uInt diffsign; /* non-zero if signs differ */
2021 uInt hipad; /* pad digit for hi if needed */
2022 Int padding; /* excess exponent */
2023 uInt carry; /* +1 for ten's complement and during add */
2024 uByte *ub, *uh, *ul; /* work */
2025 uInt uiwork; /* for macros */
2027 /* handle all the special values [any special operand leads to a */
2028 /* special result] */
2029 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
2030 decFloat proxy; /* multiplication result proxy */
2031 /* NaNs are handled as usual, giving priority to sNaNs */
2032 if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2033 if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
2034 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2035 if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
2036 /* One or more of the three is infinite */
2037 /* infinity times zero is bad */
2038 decFloatZero(&proxy);
2040 if (DFISZERO(dfr)) return decInvalid(result, set);
2041 decInfinity(&proxy, &proxy);
2043 else if (DFISINF(dfr)) {
2044 if (DFISZERO(dfl)) return decInvalid(result, set);
2045 decInfinity(&proxy, &proxy);
2047 /* compute sign of multiplication and place in proxy */
2048 DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
2049 if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
2050 /* dff is Infinite */
2051 if (!DFISINF(&proxy)) return decInfinity(result, dff);
2052 /* both sides of addition are infinite; different sign is bad */
2053 if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
2054 return decInvalid(result, set);
2055 return decFloatCopy(result, &proxy);
2058 /* Here when all operands are finite */
2060 /* First multiply dfl*dfr */
2061 decFiniteMultiply(&mul, acc+1, dfl, dfr);
2062 /* The multiply is complete, exact and unbounded, and described in */
2063 /* mul with the coefficient held in acc[1...] */
2065 /* now add in dff; the algorithm is essentially the same as */
2066 /* decFloatAdd, but the code is different because the code there */
2067 /* is highly optimized for adding two numbers of the same size */
2068 fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
2069 fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
2070 diffsign=mul.sign^fin.sign; /* note if signs differ */
2072 fin.lsd=coe+DECPMAX-1;
2073 GETCOEFF(dff, coe); /* extract the coefficient */
2075 /* now set hi and lo so that hi points to whichever of mul and fin */
2076 /* has the higher exponent and lo points to the other [don't care, */
2077 /* if the same]. One coefficient will be in acc, the other in coe. */
2078 if (mul.exponent>=fin.exponent) {
2087 /* remove leading zeros on both operands; this will save time later */
2088 /* and make testing for zero trivial (tests are safe because acc */
2089 /* and coe are rounded up to uInts) */
2090 for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
2091 for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
2092 for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2093 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2095 /* if hi is zero then result will be lo (which has the smaller */
2096 /* exponent), which also may need to be tested for zero for the */
2097 /* weird IEEE 754 sign rules */
2098 if (*hi->msd==0) { /* hi is zero */
2099 /* "When the sum of two operands with opposite signs is */
2100 /* exactly zero, the sign of that sum shall be '+' in all */
2101 /* rounding modes except round toward -Infinity, in which */
2102 /* mode that sign shall be '-'." */
2104 if (*lo->msd==0) { /* lo is zero */
2106 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2107 } /* diffsign && lo=0 */
2109 return decFinalize(result, lo, set); /* may need clamping */
2110 } /* numfl is zero */
2111 /* [here, both are minimal length and hi is non-zero] */
2112 /* (if lo is zero then padding with zeros may be needed, below) */
2114 /* if signs differ, take the ten's complement of hi (zeros to the */
2115 /* right do not matter because the complement of zero is zero); the */
2116 /* +1 is done later, as part of the addition, inserted at the */
2123 /* exactly the correct number of digits must be inverted */
2124 for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
2125 for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
2128 /* ready to add; note that hi has no leading zeros so gap */
2129 /* calculation does not have to be as pessimistic as in decFloatAdd */
2130 /* (this is much more like the arbitrary-precision algorithm in */
2131 /* Rexx and decNumber) */
2133 /* padding is the number of zeros that would need to be added to hi */
2134 /* for its lsd to be aligned with the lsd of lo */
2135 padding=hi->exponent-lo->exponent;
2136 /* printf("FMA pad %ld\n", (LI)padding); */
2138 /* the result of the addition will be built into the accumulator, */
2139 /* starting from the far right; this could be either hi or lo, and */
2140 /* will be aligned */
2141 ub=acc+FMALEN-1; /* where lsd of result will go */
2142 ul=lo->lsd; /* lsd of rhs */
2144 if (padding!=0) { /* unaligned */
2145 /* if the msd of lo is more than DECPMAX+2 digits to the right of */
2146 /* the original msd of hi then it can be reduced to a single */
2147 /* digit at the right place, as it stays clear of hi digits */
2148 /* [it must be DECPMAX+2 because during a subtraction the msd */
2149 /* could become 0 after a borrow from 1.000 to 0.9999...] */
2151 Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
2152 Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
2154 if (hilen+padding-lolen > DECPMAX+2) { /* can reduce lo to single */
2155 /* make sure it is virtually at least DECPMAX from hi->msd, at */
2156 /* least to right of hi->lsd (in case of destructive subtract), */
2157 /* and separated by at least two digits from either of those */
2158 /* (the tricky DOUBLE case is when hi is a 1 that will become a */
2159 /* 0.9999... by subtraction: */
2161 /* lo: .................1000000000000000 E-16 */
2162 /* which for the addition pads to: */
2163 /* hi: 1000000000000000000 E-16 */
2164 /* lo: .................1000000000000000 E-16 */
2165 Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
2167 /* printf("FMA reduce: %ld\n", (LI)reduce); */
2168 lo->lsd=lo->msd; /* to single digit [maybe 0] */
2169 lo->exponent=newexp; /* new lowest exponent */
2170 padding=hi->exponent-lo->exponent; /* recalculate */
2171 ul=lo->lsd; /* .. and repoint */
2174 /* padding is still > 0, but will fit in acc (less leading carry slot) */
2176 if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
2177 if (hilen+padding+1>FMALEN)
2178 printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding);
2179 /* printf("FMA padding: %ld\n", (LI)padding); */
2182 /* padding digits can now be set in the result; one or more of */
2183 /* these will come from lo; others will be zeros in the gap */
2184 for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
2185 UBFROMUI(ub-3, UBTOUI(ul-3)); /* [cannot overlap] */
2187 for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
2188 for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
2191 /* addition now complete to the right of the rightmost digit of hi */
2194 /* dow do the add from hi->lsd to the left */
2195 /* [bytewise, because either operand can run out at any time] */
2196 /* carry was set up depending on ten's complement above */
2197 /* first assume both operands have some digits */
2199 if (uh<hi->msd || ul<lo->msd) break;
2200 *ub=(uByte)(carry+(*uh--)+(*ul--));
2202 if (*ub<10) continue;
2207 if (ul<lo->msd) { /* to left of lo */
2209 if (uh<hi->msd) break;
2210 *ub=(uByte)(carry+(*uh--)); /* [+0] */
2212 if (*ub<10) continue;
2217 else { /* to left of hi */
2219 if (ul<lo->msd) break;
2220 *ub=(uByte)(carry+hipad+(*ul--));
2222 if (*ub<10) continue;
2228 /* addition complete -- now handle carry, borrow, etc. */
2229 /* use lo to set up the num (its exponent is already correct, and */
2230 /* sign usually is) */
2232 lo->lsd=acc+FMALEN-1;
2233 /* decShowNum(lo, "lo"); */
2234 if (!diffsign) { /* same-sign addition */
2235 if (carry) { /* carry out */
2236 *ub=1; /* place the 1 .. */
2237 lo->msd--; /* .. and update */
2240 else { /* signs differed (subtraction) */
2241 if (!carry) { /* no carry out means hi<lo */
2242 /* borrowed -- take ten's complement of the right digits */
2243 lo->sign=hi->sign; /* sign is lhs sign */
2244 for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul));
2245 for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
2246 /* complete the ten's complement by adding 1 [cannot overrun] */
2247 for (ul--; *ul==9; ul--) *ul=0;
2250 else { /* carry out means hi>=lo */
2251 /* sign to use is lo->sign */
2252 /* all done except for the special IEEE 754 exact-zero-result */
2253 /* rule (see above); while testing for zero, strip leading */
2254 /* zeros (which will save decFinalize doing it) */
2255 for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2256 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2257 if (*lo->msd==0) { /* must be true zero (and diffsign) */
2258 lo->sign=0; /* assume + */
2259 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2261 /* [else was not zero, might still have leading zeros] */
2262 } /* subtraction gave positive result */
2266 /* assert no left underrun */
2268 printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
2272 return decFinalize(result, lo, set); /* round, check, and lay out */
2275 /* ------------------------------------------------------------------ */
2276 /* decFloatFromInt -- initialise a decFloat from an Int */
2278 /* result gets the converted Int */
2279 /* n is the Int to convert */
2280 /* returns result */
2282 /* The result is Exact; no errors or exceptions are possible. */
2283 /* ------------------------------------------------------------------ */
2284 decFloat * decFloatFromInt32(decFloat *result, Int n) {
2285 uInt u=(uInt)n; /* copy as bits */
2286 uInt encode; /* work */
2287 DFWORD(result, 0)=ZEROWORD; /* always */
2289 DFWORD(result, 1)=0;
2290 DFWORD(result, 2)=0;
2292 if (n<0) { /* handle -n with care */
2293 /* [This can be done without the test, but is then slightly slower] */
2295 DFWORD(result, 0)|=DECFLOAT_Sign;
2297 /* Since the maximum value of u now is 2**31, only the low word of */
2298 /* result is affected */
2299 encode=BIN2DPD[u%1000];
2301 encode|=BIN2DPD[u%1000]<<10;
2303 encode|=BIN2DPD[u%1000]<<20;
2304 u/=1000; /* now 0, 1, or 2 */
2306 DFWORD(result, DECWORDS-1)=encode;
2308 } /* decFloatFromInt32 */
2310 /* ------------------------------------------------------------------ */
2311 /* decFloatFromUInt -- initialise a decFloat from a uInt */
2313 /* result gets the converted uInt */
2314 /* n is the uInt to convert */
2315 /* returns result */
2317 /* The result is Exact; no errors or exceptions are possible. */
2318 /* ------------------------------------------------------------------ */
2319 decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
2320 uInt encode; /* work */
2321 DFWORD(result, 0)=ZEROWORD; /* always */
2323 DFWORD(result, 1)=0;
2324 DFWORD(result, 2)=0;
2326 encode=BIN2DPD[u%1000];
2328 encode|=BIN2DPD[u%1000]<<10;
2330 encode|=BIN2DPD[u%1000]<<20;
2331 u/=1000; /* now 0 -> 4 */
2333 DFWORD(result, DECWORDS-1)=encode;
2334 DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
2336 } /* decFloatFromUInt32 */
2338 /* ------------------------------------------------------------------ */
2339 /* decFloatInvert -- logical digitwise INVERT of a decFloat */
2341 /* result gets the result of INVERTing df */
2342 /* df is the decFloat to invert */
2343 /* set is the context */
2344 /* returns result, which will be canonical with sign=0 */
2346 /* The operand must be positive, finite with exponent q=0, and */
2347 /* comprise just zeros and ones; if not, Invalid operation results. */
2348 /* ------------------------------------------------------------------ */
2349 decFloat * decFloatInvert(decFloat *result, const decFloat *df,
2351 uInt sourhi=DFWORD(df, 0); /* top word of dfs */
2353 if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
2354 /* the operand is a finite integer (q=0) */
2356 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
2357 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
2359 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
2360 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
2361 DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
2362 DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
2365 } /* decFloatInvert */
2367 /* ------------------------------------------------------------------ */
2368 /* decFloatIs -- decFloat tests (IsSigned, etc.) */
2370 /* df is the decFloat to test */
2371 /* returns 0 or 1 in a uInt */
2373 /* Many of these could be macros, but having them as real functions */
2374 /* is a little cleaner (and they can be referred to here by the */
2375 /* generic names) */
2376 /* ------------------------------------------------------------------ */
2377 uInt decFloatIsCanonical(const decFloat *df) {
2378 if (DFISSPECIAL(df)) {
2380 if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
2381 if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
2385 if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
2386 if (DFISCCZERO(df)) return 1; /* coefficient continuation */
2387 /* drop through to check payload */
2389 { /* declare block */
2391 uInt sourhi=DFWORD(df, 0);
2392 uInt sourlo=DFWORD(df, 1);
2393 if (CANONDPDOFF(sourhi, 8)
2394 && CANONDPDTWO(sourhi, sourlo, 30)
2395 && CANONDPDOFF(sourlo, 20)
2396 && CANONDPDOFF(sourlo, 10)
2397 && CANONDPDOFF(sourlo, 0)) return 1;
2399 uInt sourhi=DFWORD(df, 0);
2400 uInt sourmh=DFWORD(df, 1);
2401 uInt sourml=DFWORD(df, 2);
2402 uInt sourlo=DFWORD(df, 3);
2403 if (CANONDPDOFF(sourhi, 4)
2404 && CANONDPDTWO(sourhi, sourmh, 26)
2405 && CANONDPDOFF(sourmh, 16)
2406 && CANONDPDOFF(sourmh, 6)
2407 && CANONDPDTWO(sourmh, sourml, 28)
2408 && CANONDPDOFF(sourml, 18)
2409 && CANONDPDOFF(sourml, 8)
2410 && CANONDPDTWO(sourml, sourlo, 30)
2411 && CANONDPDOFF(sourlo, 20)
2412 && CANONDPDOFF(sourlo, 10)
2413 && CANONDPDOFF(sourlo, 0)) return 1;
2416 return 0; /* a declet is non-canonical */
2419 uInt decFloatIsFinite(const decFloat *df) {
2420 return !DFISSPECIAL(df);
2422 uInt decFloatIsInfinite(const decFloat *df) {
2425 uInt decFloatIsInteger(const decFloat *df) {
2428 uInt decFloatIsNaN(const decFloat *df) {
2431 uInt decFloatIsNormal(const decFloat *df) {
2432 Int exp; /* exponent */
2433 if (DFISSPECIAL(df)) return 0;
2434 if (DFISZERO(df)) return 0;
2435 /* is finite and non-zero */
2436 exp=GETEXPUN(df) /* get unbiased exponent .. */
2437 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2438 return (exp>=DECEMIN); /* < DECEMIN is subnormal */
2440 uInt decFloatIsSignaling(const decFloat *df) {
2441 return DFISSNAN(df);
2443 uInt decFloatIsSignalling(const decFloat *df) {
2444 return DFISSNAN(df);
2446 uInt decFloatIsSigned(const decFloat *df) {
2447 return DFISSIGNED(df);
2449 uInt decFloatIsSubnormal(const decFloat *df) {
2450 if (DFISSPECIAL(df)) return 0;
2452 if (decFloatIsNormal(df)) return 0;
2453 /* it is <Nmin, but could be zero */
2454 if (DFISZERO(df)) return 0;
2455 return 1; /* is subnormal */
2457 uInt decFloatIsZero(const decFloat *df) {
2458 return DFISZERO(df);
2459 } /* decFloatIs... */
2461 /* ------------------------------------------------------------------ */
2462 /* decFloatLogB -- return adjusted exponent, by 754 rules */
2464 /* result gets the adjusted exponent as an integer, or a NaN etc. */
2465 /* df is the decFloat to be examined */
2466 /* set is the context */
2467 /* returns result */
2469 /* Notable cases: */
2470 /* A<0 -> Use |A| */
2471 /* A=0 -> -Infinity (Division by zero) */
2472 /* A=Infinite -> +Infinity (Exact) */
2473 /* A=1 exactly -> 0 (Exact) */
2474 /* NaNs are propagated as usual */
2475 /* ------------------------------------------------------------------ */
2476 decFloat * decFloatLogB(decFloat *result, const decFloat *df,
2478 Int ae; /* adjusted exponent */
2479 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2481 DFWORD(result, 0)=0; /* need +ve */
2482 return decInfinity(result, result); /* canonical +Infinity */
2485 set->status|=DEC_Division_by_zero; /* as per 754 */
2486 DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
2487 return decInfinity(result, result); /* canonical -Infinity */
2489 ae=GETEXPUN(df) /* get unbiased exponent .. */
2490 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2491 /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
2492 /* it is worth using a special case of decFloatFromInt32 */
2493 DFWORD(result, 0)=ZEROWORD; /* always */
2495 DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
2499 DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
2501 DFWORD(result, 1)=0;
2502 DFWORD(result, 2)=0;
2503 DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
2504 DFWORD(result, 3)|=BIN2DPD[ae%1000];
2507 } /* decFloatLogB */
2509 /* ------------------------------------------------------------------ */
2510 /* decFloatMax -- return maxnum of two operands */
2512 /* result gets the chosen decFloat */
2513 /* dfl is the first decFloat (lhs) */
2514 /* dfr is the second decFloat (rhs) */
2515 /* set is the context */
2516 /* returns result */
2518 /* If just one operand is a quiet NaN it is ignored. */
2519 /* ------------------------------------------------------------------ */
2520 decFloat * decFloatMax(decFloat *result,
2521 const decFloat *dfl, const decFloat *dfr,
2525 /* sNaN or both NaNs leads to normal NaN processing */
2526 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2527 return decCanonical(result, dfr); /* RHS is numeric */
2530 /* sNaN leads to normal NaN processing (both NaN handled above) */
2531 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2532 return decCanonical(result, dfl); /* LHS is numeric */
2534 /* Both operands are numeric; numeric comparison needed -- use */
2535 /* total order for a well-defined choice (and +0 > -0) */
2536 comp=decNumCompare(dfl, dfr, 1);
2537 if (comp>=0) return decCanonical(result, dfl);
2538 return decCanonical(result, dfr);
2541 /* ------------------------------------------------------------------ */
2542 /* decFloatMaxMag -- return maxnummag of two operands */
2544 /* result gets the chosen decFloat */
2545 /* dfl is the first decFloat (lhs) */
2546 /* dfr is the second decFloat (rhs) */
2547 /* set is the context */
2548 /* returns result */
2550 /* Returns according to the magnitude comparisons if both numeric and */
2551 /* unequal, otherwise returns maxnum */
2552 /* ------------------------------------------------------------------ */
2553 decFloat * decFloatMaxMag(decFloat *result,
2554 const decFloat *dfl, const decFloat *dfr,
2557 decFloat absl, absr;
2558 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
2560 decFloatCopyAbs(&absl, dfl);
2561 decFloatCopyAbs(&absr, dfr);
2562 comp=decNumCompare(&absl, &absr, 0);
2563 if (comp>0) return decCanonical(result, dfl);
2564 if (comp<0) return decCanonical(result, dfr);
2565 return decFloatMax(result, dfl, dfr, set);
2566 } /* decFloatMaxMag */
2568 /* ------------------------------------------------------------------ */
2569 /* decFloatMin -- return minnum of two operands */
2571 /* result gets the chosen decFloat */
2572 /* dfl is the first decFloat (lhs) */
2573 /* dfr is the second decFloat (rhs) */
2574 /* set is the context */
2575 /* returns result */
2577 /* If just one operand is a quiet NaN it is ignored. */
2578 /* ------------------------------------------------------------------ */
2579 decFloat * decFloatMin(decFloat *result,
2580 const decFloat *dfl, const decFloat *dfr,
2584 /* sNaN or both NaNs leads to normal NaN processing */
2585 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2586 return decCanonical(result, dfr); /* RHS is numeric */
2589 /* sNaN leads to normal NaN processing (both NaN handled above) */
2590 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2591 return decCanonical(result, dfl); /* LHS is numeric */
2593 /* Both operands are numeric; numeric comparison needed -- use */
2594 /* total order for a well-defined choice (and +0 > -0) */
2595 comp=decNumCompare(dfl, dfr, 1);
2596 if (comp<=0) return decCanonical(result, dfl);
2597 return decCanonical(result, dfr);
2600 /* ------------------------------------------------------------------ */
2601 /* decFloatMinMag -- return minnummag of two operands */
2603 /* result gets the chosen decFloat */
2604 /* dfl is the first decFloat (lhs) */
2605 /* dfr is the second decFloat (rhs) */
2606 /* set is the context */
2607 /* returns result */
2609 /* Returns according to the magnitude comparisons if both numeric and */
2610 /* unequal, otherwise returns minnum */
2611 /* ------------------------------------------------------------------ */
2612 decFloat * decFloatMinMag(decFloat *result,
2613 const decFloat *dfl, const decFloat *dfr,
2616 decFloat absl, absr;
2617 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
2619 decFloatCopyAbs(&absl, dfl);
2620 decFloatCopyAbs(&absr, dfr);
2621 comp=decNumCompare(&absl, &absr, 0);
2622 if (comp<0) return decCanonical(result, dfl);
2623 if (comp>0) return decCanonical(result, dfr);
2624 return decFloatMin(result, dfl, dfr, set);
2625 } /* decFloatMinMag */
2627 /* ------------------------------------------------------------------ */
2628 /* decFloatMinus -- negate value, heeding NaNs, etc. */
2630 /* result gets the canonicalized 0-df */
2631 /* df is the decFloat to minus */
2632 /* set is the context */
2633 /* returns result */
2635 /* This has the same effect as 0-df where the exponent of the zero is */
2636 /* the same as that of df (if df is finite). */
2637 /* The effect is also the same as decFloatCopyNegate except that NaNs */
2638 /* are handled normally (the sign of a NaN is not affected, and an */
2639 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2640 /* ------------------------------------------------------------------ */
2641 decFloat * decFloatMinus(decFloat *result, const decFloat *df,
2643 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2644 decCanonical(result, df); /* copy and check */
2645 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2646 else DFBYTE(result, 0)^=0x80; /* flip sign bit */
2648 } /* decFloatMinus */
2650 /* ------------------------------------------------------------------ */
2651 /* decFloatMultiply -- multiply two decFloats */
2653 /* result gets the result of multiplying dfl and dfr: */
2654 /* dfl is the first decFloat (lhs) */
2655 /* dfr is the second decFloat (rhs) */
2656 /* set is the context */
2657 /* returns result */
2659 /* ------------------------------------------------------------------ */
2660 decFloat * decFloatMultiply(decFloat *result,
2661 const decFloat *dfl, const decFloat *dfr,
2663 bcdnum num; /* for final conversion */
2664 uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
2666 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
2667 /* NaNs are handled as usual */
2668 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2669 /* infinity times zero is bad */
2670 if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
2671 if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
2672 /* both infinite; return canonical infinity with computed sign */
2673 DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
2674 return decInfinity(result, result);
2677 /* Here when both operands are finite */
2678 decFiniteMultiply(&num, bcdacc, dfl, dfr);
2679 return decFinalize(result, &num, set); /* round, check, and lay out */
2680 } /* decFloatMultiply */
2682 /* ------------------------------------------------------------------ */
2683 /* decFloatNextMinus -- next towards -Infinity */
2685 /* result gets the next lesser decFloat */
2686 /* dfl is the decFloat to start with */
2687 /* set is the context */
2688 /* returns result */
2690 /* This is 754 nextdown; Invalid is the only status possible (from */
2692 /* ------------------------------------------------------------------ */
2693 decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
2695 decFloat delta; /* tiny increment */
2696 uInt savestat; /* saves status */
2697 enum rounding saveround; /* .. and mode */
2699 /* +Infinity is the special case */
2700 if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
2702 return result; /* [no status to set] */
2704 /* other cases are effected by sutracting a tiny delta -- this */
2705 /* should be done in a wider format as the delta is unrepresentable */
2706 /* here (but can be done with normal add if the sign of zero is */
2707 /* treated carefully, because no Inexactitude is interesting); */
2708 /* rounding to -Infinity then pushes the result to next below */
2709 decFloatZero(&delta); /* set up tiny delta */
2710 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2711 DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
2712 /* set up for the directional round */
2713 saveround=set->round; /* save mode */
2714 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2715 savestat=set->status; /* save status */
2716 decFloatAdd(result, dfl, &delta, set);
2717 /* Add rules mess up the sign when going from +Ntiny to 0 */
2718 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2719 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2720 set->status|=savestat; /* restore pending flags */
2721 set->round=saveround; /* .. and mode */
2723 } /* decFloatNextMinus */
2725 /* ------------------------------------------------------------------ */
2726 /* decFloatNextPlus -- next towards +Infinity */
2728 /* result gets the next larger decFloat */
2729 /* dfl is the decFloat to start with */
2730 /* set is the context */
2731 /* returns result */
2733 /* This is 754 nextup; Invalid is the only status possible (from */
2735 /* ------------------------------------------------------------------ */
2736 decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
2738 uInt savestat; /* saves status */
2739 enum rounding saveround; /* .. and mode */
2740 decFloat delta; /* tiny increment */
2742 /* -Infinity is the special case */
2743 if (DFISINF(dfl) && DFISSIGNED(dfl)) {
2745 DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
2746 return result; /* [no status to set] */
2748 /* other cases are effected by sutracting a tiny delta -- this */
2749 /* should be done in a wider format as the delta is unrepresentable */
2750 /* here (but can be done with normal add if the sign of zero is */
2751 /* treated carefully, because no Inexactitude is interesting); */
2752 /* rounding to +Infinity then pushes the result to next above */
2753 decFloatZero(&delta); /* set up tiny delta */
2754 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2755 DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
2756 /* set up for the directional round */
2757 saveround=set->round; /* save mode */
2758 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2759 savestat=set->status; /* save status */
2760 decFloatAdd(result, dfl, &delta, set);
2761 /* Add rules mess up the sign when going from -Ntiny to -0 */
2762 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2763 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2764 set->status|=savestat; /* restore pending flags */
2765 set->round=saveround; /* .. and mode */
2767 } /* decFloatNextPlus */
2769 /* ------------------------------------------------------------------ */
2770 /* decFloatNextToward -- next towards a decFloat */
2772 /* result gets the next decFloat */
2773 /* dfl is the decFloat to start with */
2774 /* dfr is the decFloat to move toward */
2775 /* set is the context */
2776 /* returns result */
2778 /* This is 754-1985 nextafter, as modified during revision (dropped */
2779 /* from 754-2008); status may be set unless the result is a normal */
2781 /* ------------------------------------------------------------------ */
2782 decFloat * decFloatNextToward(decFloat *result,
2783 const decFloat *dfl, const decFloat *dfr,
2785 decFloat delta; /* tiny increment or decrement */
2786 decFloat pointone; /* 1e-1 */
2787 uInt savestat; /* saves status */
2788 enum rounding saveround; /* .. and mode */
2789 uInt deltatop; /* top word for delta */
2790 Int comp; /* work */
2792 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2793 /* Both are numeric, so Invalid no longer a possibility */
2794 comp=decNumCompare(dfl, dfr, 0);
2795 if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
2796 /* unequal; do NextPlus or NextMinus but with different status rules */
2798 if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
2799 if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
2801 DFWORD(result, 0)|=DECFLOAT_Sign;
2804 saveround=set->round; /* save mode */
2805 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2806 deltatop=0; /* positive delta */
2808 else { /* lhs>rhs, do NextMinus, see above for commentary */
2809 if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
2813 saveround=set->round; /* save mode */
2814 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2815 deltatop=DECFLOAT_Sign; /* negative delta */
2817 savestat=set->status; /* save status */
2818 /* Here, Inexact is needed where appropriate (and hence Underflow, */
2819 /* etc.). Therefore the tiny delta which is otherwise */
2820 /* unrepresentable (see NextPlus and NextMinus) is constructed */
2821 /* using the multiplication of FMA. */
2822 decFloatZero(&delta); /* set up tiny delta */
2823 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2824 DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
2825 decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
2826 decFloatFMA(result, &delta, &pointone, dfl, set);
2827 /* [Delta is truly tiny, so no need to correct sign of zero] */
2828 /* use new status unless the result is normal */
2829 if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
2830 set->round=saveround; /* restore mode */
2832 } /* decFloatNextToward */
2834 /* ------------------------------------------------------------------ */
2835 /* decFloatOr -- logical digitwise OR of two decFloats */
2837 /* result gets the result of ORing dfl and dfr */
2838 /* dfl is the first decFloat (lhs) */
2839 /* dfr is the second decFloat (rhs) */
2840 /* set is the context */
2841 /* returns result, which will be canonical with sign=0 */
2843 /* The operands must be positive, finite with exponent q=0, and */
2844 /* comprise just zeros and ones; if not, Invalid operation results. */
2845 /* ------------------------------------------------------------------ */
2846 decFloat * decFloatOr(decFloat *result,
2847 const decFloat *dfl, const decFloat *dfr,
2849 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
2850 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
2851 /* the operands are positive finite integers (q=0) with just 0s and 1s */
2853 DFWORD(result, 0)=ZEROWORD
2854 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
2855 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
2857 DFWORD(result, 0)=ZEROWORD
2858 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
2859 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
2860 DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
2861 DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
2866 /* ------------------------------------------------------------------ */
2867 /* decFloatPlus -- add value to 0, heeding NaNs, etc. */
2869 /* result gets the canonicalized 0+df */
2870 /* df is the decFloat to plus */
2871 /* set is the context */
2872 /* returns result */
2874 /* This has the same effect as 0+df where the exponent of the zero is */
2875 /* the same as that of df (if df is finite). */
2876 /* The effect is also the same as decFloatCopy except that NaNs */
2877 /* are handled normally (the sign of a NaN is not affected, and an */
2878 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2879 /* ------------------------------------------------------------------ */
2880 decFloat * decFloatPlus(decFloat *result, const decFloat *df,
2882 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2883 decCanonical(result, df); /* copy and check */
2884 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2886 } /* decFloatPlus */
2888 /* ------------------------------------------------------------------ */
2889 /* decFloatQuantize -- quantize a decFloat */
2891 /* result gets the result of quantizing dfl to match dfr */
2892 /* dfl is the first decFloat (lhs) */
2893 /* dfr is the second decFloat (rhs), which sets the exponent */
2894 /* set is the context */
2895 /* returns result */
2897 /* Unless there is an error or the result is infinite, the exponent */
2898 /* of result is guaranteed to be the same as that of dfr. */
2899 /* ------------------------------------------------------------------ */
2900 decFloat * decFloatQuantize(decFloat *result,
2901 const decFloat *dfl, const decFloat *dfr,
2903 Int explb, exprb; /* left and right biased exponents */
2904 uByte *ulsd; /* local LSD pointer */
2905 uByte *ub, *uc; /* work */
2908 uInt encode; /* encoding accumulator */
2909 uInt sourhil, sourhir; /* top words from source decFloats */
2910 uInt uiwork; /* for macros */
2912 uShort uswork; /* .. */
2914 /* the following buffer holds the coefficient for manipulation */
2915 uByte buf[4+DECPMAX*3+2*QUAD]; /* + space for zeros to left or right */
2917 bcdnum num; /* for trace displays */
2920 /* Start decoding the arguments */
2921 sourhil=DFWORD(dfl, 0); /* LHS top word */
2922 explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
2923 sourhir=DFWORD(dfr, 0); /* RHS top word */
2924 exprb=DECCOMBEXP[sourhir>>26];
2926 if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
2927 /* NaNs are handled as usual */
2928 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2929 /* one infinity but not both is bad */
2930 if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
2931 /* both infinite; return canonical infinity with sign of LHS */
2932 return decInfinity(result, dfl);
2935 /* Here when both arguments are finite */
2936 /* complete extraction of the exponents [no need to unbias] */
2937 explb+=GETECON(dfl); /* + continuation */
2938 exprb+=GETECON(dfr); /* .. */
2940 /* calculate the number of digits to drop from the coefficient */
2941 drop=exprb-explb; /* 0 if nothing to do */
2942 if (drop==0) return decCanonical(result, dfl); /* return canonical */
2944 /* the coefficient is needed; lay it out into buf, offset so zeros */
2945 /* can be added before or after as needed -- an extra heading is */
2946 /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
2948 #define BUFOFF (buf+4+DECPMAX)
2949 GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
2950 /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
2954 num.lsd=BUFOFF+DECPMAX-1;
2955 num.exponent=explb-DECBIAS;
2956 num.sign=sourhil & DECFLOAT_Sign;
2957 decShowNum(&num, "dfl");
2960 if (drop>0) { /* [most common case] */
2961 /* (this code is very similar to that in decFloatFinalize, but */
2962 /* has many differences so is duplicated here -- so any changes */
2963 /* may need to be made there, too) */
2964 uByte *roundat; /* -> re-round digit */
2965 uByte reround; /* reround value */
2966 /* printf("Rounding; drop=%ld\n", (LI)drop); */
2968 /* there is at least one zero needed to the left, in all but one */
2969 /* exceptional (all-nines) case, so place four zeros now; this is */
2970 /* needed almost always and makes rounding all-nines by fours safe */
2971 UBFROMUI(BUFOFF-4, 0);
2973 /* Three cases here: */
2974 /* 1. new LSD is in coefficient (almost always) */
2975 /* 2. new LSD is digit to left of coefficient (so MSD is */
2976 /* round-for-reround digit) */
2977 /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
2978 /* Note that leading zeros can safely be treated as useful digits */
2980 /* [duplicate check-stickies code to save a test] */
2981 /* [by-digit check for stickies as runs of zeros are rare] */
2982 if (drop<DECPMAX) { /* NB lengths not addresses */
2983 roundat=BUFOFF+DECPMAX-drop;
2985 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2986 if (*ub!=0) { /* non-zero to be discarded */
2987 reround=DECSTICKYTAB[reround]; /* apply sticky bit */
2988 break; /* [remainder don't-care] */
2990 } /* check stickies */
2991 ulsd=roundat-1; /* set LSD */
2993 else { /* edge case */
2994 if (drop==DECPMAX) {