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
) {
3002 for (ub
=roundat
+1; ub
<BUFOFF
+DECPMAX
; ub
++) {
3003 if (*ub
!=0) { /* non-zero to be discarded */
3004 reround
=DECSTICKYTAB
[reround
]; /* apply sticky bit */
3005 break; /* [remainder don't-care] */
3007 } /* check stickies */
3008 *BUFOFF
=0; /* make a coefficient of 0 */
3009 ulsd
=BUFOFF
; /* .. at the MSD place */
3012 if (reround
!=0) { /* discarding non-zero */
3014 set
->status
|=DEC_Inexact
;
3016 /* next decide whether to increment the coefficient */
3017 if (set
->round
==DEC_ROUND_HALF_EVEN
) { /* fastpath slowest case */
3018 if (reround
>5) bump
=1; /* >0.5 goes up */
3019 else if (reround
==5) /* exactly 0.5000 .. */
3020 bump
=*ulsd
& 0x01; /* .. up iff [new] lsd is odd */
3022 else switch (set
->round
) {
3023 case DEC_ROUND_DOWN
: {
3026 case DEC_ROUND_HALF_DOWN
: {
3027 if (reround
>5) bump
=1;
3029 case DEC_ROUND_HALF_UP
: {
3030 if (reround
>=5) bump
=1;
3032 case DEC_ROUND_UP
: {
3033 if (reround
>0) bump
=1;
3035 case DEC_ROUND_CEILING
: {
3036 /* same as _UP for positive numbers, and as _DOWN for negatives */
3037 if (!(sourhil
&DECFLOAT_Sign
) && reround
>0) bump
=1;
3039 case DEC_ROUND_FLOOR
: {
3040 /* same as _UP for negative numbers, and as _DOWN for positive */
3041 /* [negative reround cannot occur on 0] */
3042 if (sourhil
&DECFLOAT_Sign
&& reround
>0) bump
=1;
3044 case DEC_ROUND_05UP
: {
3045 if (reround
>0) { /* anything out there is 'sticky' */
3046 /* bump iff lsd=0 or 5; this cannot carry so it could be */
3047 /* effected immediately with no bump -- but the code */
3048 /* is clearer if this is done the same way as the others */
3049 if (*ulsd
==0 || *ulsd
==5) bump
=1;
3052 default: { /* e.g., DEC_ROUND_MAX */
3053 set
->status
|=DEC_Invalid_context
;
3055 printf("Unknown rounding mode: %ld\n", (LI
)set
->round
);
3058 } /* switch (not r-h-e) */
3059 /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
3061 if (bump
!=0) { /* need increment */
3062 /* increment the coefficient; this could give 1000... (after */
3063 /* the all nines case) */
3065 for (; UBTOUI(ub
-3)==0x09090909; ub
-=4) UBFROMUI(ub
-3, 0);
3066 /* now at most 3 digits left to non-9 (usually just the one) */
3067 for (; *ub
==9; ub
--) *ub
=0;
3069 /* [the all-nines case will have carried one digit to the */
3070 /* left of the original MSD -- just where it is needed] */
3072 } /* inexact rounding */
3074 /* now clear zeros to the left so exactly DECPMAX digits will be */
3075 /* available in the coefficent -- the first word to the left was */
3076 /* cleared earlier for safe carry; now add any more needed */
3078 UBFROMUI(BUFOFF
-8, 0); /* must be at least 5 */
3079 for (uc
=BUFOFF
-12; uc
>ulsd
-DECPMAX
-3; uc
-=4) UBFROMUI(uc
, 0);
3081 } /* need round (drop>0) */
3083 else { /* drop<0; padding with -drop digits is needed */
3084 /* This is the case where an error can occur if the padded */
3085 /* coefficient will not fit; checking for this can be done in the */
3086 /* same loop as padding for zeros if the no-hope and zero cases */
3087 /* are checked first */
3088 if (-drop
>DECPMAX
-1) { /* cannot fit unless 0 */
3089 if (!ISCOEFFZERO(BUFOFF
)) return decInvalid(result
, set
);
3090 /* a zero can have any exponent; just drop through and use it */
3091 ulsd
=BUFOFF
+DECPMAX
-1;
3093 else { /* padding will fit (but may still be too long) */
3094 /* final-word mask depends on endianess */
3096 static const uInt dmask
[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
3098 static const uInt dmask
[]={0, 0xff000000, 0xffff0000, 0xffffff00};
3100 /* note that here zeros to the right are added by fours, so in */
3101 /* the Quad case this could write 36 zeros if the coefficient has */
3102 /* fewer than three significant digits (hence the +2*QUAD for buf) */
3103 for (uc
=BUFOFF
+DECPMAX
;; uc
+=4) {
3105 if (UBTOUI(uc
-DECPMAX
)!=0) { /* could be bad */
3106 /* if all four digits should be zero, definitely bad */
3107 if (uc
<=BUFOFF
+DECPMAX
+(-drop
)-4)
3108 return decInvalid(result
, set
);
3109 /* must be a 1- to 3-digit sequence; check more carefully */
3110 if ((UBTOUI(uc
-DECPMAX
)&dmask
[(-drop
)%4])!=0)
3111 return decInvalid(result
, set
);
3112 break; /* no need for loop end test */
3114 if (uc
>=BUFOFF
+DECPMAX
+(-drop
)-4) break; /* done */
3116 ulsd
=BUFOFF
+DECPMAX
+(-drop
)-1;
3117 } /* pad and check leading zeros */
3121 num
.msd
=ulsd
-DECPMAX
+1;
3123 num
.exponent
=explb
-DECBIAS
;
3124 num
.sign
=sourhil
& DECFLOAT_Sign
;
3125 decShowNum(&num
, "res");
3128 /*------------------------------------------------------------------*/
3129 /* At this point the result is DECPMAX digits, ending at ulsd, so */
3130 /* fits the encoding exactly; there is no possibility of error */
3131 /*------------------------------------------------------------------*/
3132 encode
=((exprb
>>DECECONL
)<<4) + *(ulsd
-DECPMAX
+1); /* make index */
3133 encode
=DECCOMBFROM
[encode
]; /* indexed by (0-2)*16+msd */
3134 /* the exponent continuation can be extracted from the original RHS */
3135 encode
|=sourhir
& ECONMASK
;
3136 encode
|=sourhil
&DECFLOAT_Sign
; /* add the sign from LHS */
3138 /* finally encode the coefficient */
3139 /* private macro to encode a declet; this version can be used */
3140 /* because all coefficient digits exist */
3141 #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
3142 dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
3145 getDPD3q(dpd
, 4); encode
|=dpd
<<8;
3146 getDPD3q(dpd
, 3); encode
|=dpd
>>2;
3147 DFWORD(result
, 0)=encode
;
3149 getDPD3q(dpd
, 2); encode
|=dpd
<<20;
3150 getDPD3q(dpd
, 1); encode
|=dpd
<<10;
3151 getDPD3q(dpd
, 0); encode
|=dpd
;
3152 DFWORD(result
, 1)=encode
;
3155 getDPD3q(dpd
,10); encode
|=dpd
<<4;
3156 getDPD3q(dpd
, 9); encode
|=dpd
>>6;
3157 DFWORD(result
, 0)=encode
;
3159 getDPD3q(dpd
, 8); encode
|=dpd
<<16;
3160 getDPD3q(dpd
, 7); encode
|=dpd
<<6;
3161 getDPD3q(dpd
, 6); encode
|=dpd
>>4;
3162 DFWORD(result
, 1)=encode
;
3164 getDPD3q(dpd
, 5); encode
|=dpd
<<18;
3165 getDPD3q(dpd
, 4); encode
|=dpd
<<8;
3166 getDPD3q(dpd
, 3); encode
|=dpd
>>2;
3167 DFWORD(result
, 2)=encode
;
3169 getDPD3q(dpd
, 2); encode
|=dpd
<<20;
3170 getDPD3q(dpd
, 1); encode
|=dpd
<<10;
3171 getDPD3q(dpd
, 0); encode
|=dpd
;
3172 DFWORD(result
, 3)=encode
;
3175 } /* decFloatQuantize */
3177 /* ------------------------------------------------------------------ */
3178 /* decFloatReduce -- reduce finite coefficient to minimum length */
3180 /* result gets the reduced decFloat */
3181 /* df is the source decFloat */
3182 /* set is the context */
3183 /* returns result, which will be canonical */
3185 /* This removes all possible trailing zeros from the coefficient; */
3186 /* some may remain when the number is very close to Nmax. */
3187 /* Special values are unchanged and no status is set unless df=sNaN. */
3188 /* Reduced zero has an exponent q=0. */
3189 /* ------------------------------------------------------------------ */
3190 decFloat
* decFloatReduce(decFloat
*result
, const decFloat
*df
,
3192 bcdnum num
; /* work */
3193 uByte buf
[DECPMAX
], *ub
; /* coefficient and pointer */
3194 if (df
!=result
) *result
=*df
; /* copy, if needed */
3195 if (DFISNAN(df
)) return decNaNs(result
, df
, NULL
, set
); /* sNaN */
3196 /* zeros and infinites propagate too */
3197 if (DFISINF(df
)) return decInfinity(result
, df
); /* canonical */
3199 uInt sign
=DFWORD(df
, 0)&DECFLOAT_Sign
;
3200 decFloatZero(result
);
3201 DFWORD(result
, 0)|=sign
;
3202 return result
; /* exponent dropped, sign OK */
3204 /* non-zero finite */
3206 ub
=buf
+DECPMAX
-1; /* -> lsd */
3207 if (*ub
) return result
; /* no trailing zeros */
3208 for (ub
--; *ub
==0;) ub
--; /* terminates because non-zero */
3209 /* *ub is the first non-zero from the right */
3210 num
.sign
=DFWORD(df
, 0)&DECFLOAT_Sign
; /* set up number... */
3211 num
.exponent
=GETEXPUN(df
)+(Int
)(buf
+DECPMAX
-1-ub
); /* adjusted exponent */
3214 return decFinalize(result
, &num
, set
);
3215 } /* decFloatReduce */
3217 /* ------------------------------------------------------------------ */
3218 /* decFloatRemainder -- integer divide and return remainder */
3220 /* result gets the remainder of dividing dfl by dfr: */
3221 /* dfl is the first decFloat (lhs) */
3222 /* dfr is the second decFloat (rhs) */
3223 /* set is the context */
3224 /* returns result */
3226 /* ------------------------------------------------------------------ */
3227 decFloat
* decFloatRemainder(decFloat
*result
,
3228 const decFloat
*dfl
, const decFloat
*dfr
,
3230 return decDivide(result
, dfl
, dfr
, set
, REMAINDER
);
3231 } /* decFloatRemainder */
3233 /* ------------------------------------------------------------------ */
3234 /* decFloatRemainderNear -- integer divide to nearest and remainder */
3236 /* result gets the remainder of dividing dfl by dfr: */
3237 /* dfl is the first decFloat (lhs) */
3238 /* dfr is the second decFloat (rhs) */
3239 /* set is the context */
3240 /* returns result */
3242 /* This is the IEEE remainder, where the nearest integer is used. */
3243 /* ------------------------------------------------------------------ */
3244 decFloat
* decFloatRemainderNear(decFloat
*result
,
3245 const decFloat
*dfl
, const decFloat
*dfr
,
3247 return decDivide(result
, dfl
, dfr
, set
, REMNEAR
);
3248 } /* decFloatRemainderNear */
3250 /* ------------------------------------------------------------------ */
3251 /* decFloatRotate -- rotate the coefficient of a decFloat left/right */
3253 /* result gets the result of rotating dfl */
3254 /* dfl is the source decFloat to rotate */
3255 /* dfr is the count of digits to rotate, an integer (with q=0) */
3256 /* set is the context */
3257 /* returns result */
3259 /* The digits of the coefficient of dfl are rotated to the left (if */
3260 /* dfr is positive) or to the right (if dfr is negative) without */
3261 /* adjusting the exponent or the sign of dfl. */
3263 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3264 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3265 /* dfr must be valid). No status is set unless dfr is invalid or an */
3266 /* operand is an sNaN. The result is canonical. */
3267 /* ------------------------------------------------------------------ */
3268 #define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
3269 decFloat
* decFloatRotate(decFloat
*result
,
3270 const decFloat
*dfl
, const decFloat
*dfr
,
3272 Int rotate
; /* dfr as an Int */
3273 uByte buf
[DECPMAX
+PHALF
]; /* coefficient + half */
3274 uInt digits
, savestat
; /* work */
3275 bcdnum num
; /* .. */
3278 if (DFISNAN(dfl
)||DFISNAN(dfr
)) return decNaNs(result
, dfl
, dfr
, set
);
3279 if (!DFISINT(dfr
)) return decInvalid(result
, set
);
3280 digits
=decFloatDigits(dfr
); /* calculate digits */
3281 if (digits
>2) return decInvalid(result
, set
); /* definitely out of range */
3282 rotate
=DPD2BIN
[DFWORD(dfr
, DECWORDS
-1)&0x3ff]; /* is in bottom declet */
3283 if (rotate
>DECPMAX
) return decInvalid(result
, set
); /* too big */
3284 /* [from here on no error or status change is possible] */
3285 if (DFISINF(dfl
)) return decInfinity(result
, dfl
); /* canonical */
3286 /* handle no-rotate cases */
3287 if (rotate
==0 || rotate
==DECPMAX
) return decCanonical(result
, dfl
);
3288 /* a real rotate is needed: 0 < rotate < DECPMAX */
3289 /* reduce the rotation to no more than half to reduce copying later */
3290 /* (for QUAD in fact half + 2 digits) */
3291 if (DFISSIGNED(dfr
)) rotate
=-rotate
;
3292 if (abs(rotate
)>PHALF
) {
3293 if (rotate
<0) rotate
=DECPMAX
+rotate
;
3294 else rotate
=rotate
-DECPMAX
;
3296 /* now lay out the coefficient, leaving room to the right or the */
3297 /* left depending on the direction of rotation */
3299 if (rotate
<0) ub
+=PHALF
; /* rotate right, so space to left */
3301 /* copy half the digits to left or right, and set num.msd */
3303 memcpy(buf
, buf
+DECPMAX
, PHALF
);
3304 num
.msd
=buf
+PHALF
+rotate
;
3307 memcpy(buf
+DECPMAX
, buf
, PHALF
);
3310 /* fill in rest of num */
3311 num
.lsd
=num
.msd
+DECPMAX
-1;
3312 num
.sign
=DFWORD(dfl
, 0)&DECFLOAT_Sign
;
3313 num
.exponent
=GETEXPUN(dfl
);
3314 savestat
=set
->status
; /* record */
3315 decFinalize(result
, &num
, set
);
3316 set
->status
=savestat
; /* restore */
3318 } /* decFloatRotate */
3320 /* ------------------------------------------------------------------ */
3321 /* decFloatSameQuantum -- test decFloats for same quantum */
3323 /* dfl is the first decFloat (lhs) */
3324 /* dfr is the second decFloat (rhs) */
3325 /* returns 1 if the operands have the same quantum, 0 otherwise */
3327 /* No error is possible and no status results. */
3328 /* ------------------------------------------------------------------ */
3329 uInt
decFloatSameQuantum(const decFloat
*dfl
, const decFloat
*dfr
) {
3330 if (DFISSPECIAL(dfl
) || DFISSPECIAL(dfr
)) {
3331 if (DFISNAN(dfl
) && DFISNAN(dfr
)) return 1;
3332 if (DFISINF(dfl
) && DFISINF(dfr
)) return 1;
3333 return 0; /* any other special mixture gives false */
3335 if (GETEXP(dfl
)==GETEXP(dfr
)) return 1; /* biased exponents match */
3337 } /* decFloatSameQuantum */
3339 /* ------------------------------------------------------------------ */
3340 /* decFloatScaleB -- multiply by a power of 10, as per 754 */
3342 /* result gets the result of the operation */
3343 /* dfl is the first decFloat (lhs) */
3344 /* dfr is the second decFloat (rhs), am integer (with q=0) */
3345 /* set is the context */
3346 /* returns result */
3348 /* This computes result=dfl x 10**dfr where dfr is an integer in the */
3349 /* range +/-2*(emax+pmax), typically resulting from LogB. */
3350 /* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
3352 /* ------------------------------------------------------------------ */
3353 #define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
3354 decFloat
* decFloatScaleB(decFloat
*result
,
3355 const decFloat
*dfl
, const decFloat
*dfr
,
3357 uInt digits
; /* work */
3358 Int expr
; /* dfr as an Int */
3360 if (DFISNAN(dfl
)||DFISNAN(dfr
)) return decNaNs(result
, dfl
, dfr
, set
);
3361 if (!DFISINT(dfr
)) return decInvalid(result
, set
);
3362 digits
=decFloatDigits(dfr
); /* calculate digits */
3365 if (digits
>3) return decInvalid(result
, set
); /* definitely out of range */
3366 expr
=DPD2BIN
[DFWORD(dfr
, 1)&0x3ff]; /* must be in bottom declet */
3368 if (digits
>5) return decInvalid(result
, set
); /* definitely out of range */
3369 expr
=DPD2BIN
[DFWORD(dfr
, 3)&0x3ff] /* in bottom 2 declets .. */
3370 +DPD2BIN
[(DFWORD(dfr
, 3)>>10)&0x3ff]*1000; /* .. */
3372 if (expr
>SCALEBMAX
) return decInvalid(result
, set
); /* oops */
3373 /* [from now on no error possible] */
3374 if (DFISINF(dfl
)) return decInfinity(result
, dfl
); /* canonical */
3375 if (DFISSIGNED(dfr
)) expr
=-expr
;
3376 /* dfl is finite and expr is valid */
3377 *result
=*dfl
; /* copy to target */
3378 return decFloatSetExponent(result
, set
, GETEXPUN(result
)+expr
);
3379 } /* decFloatScaleB */
3381 /* ------------------------------------------------------------------ */
3382 /* decFloatShift -- shift the coefficient of a decFloat left or right */
3384 /* result gets the result of shifting dfl */
3385 /* dfl is the source decFloat to shift */
3386 /* dfr is the count of digits to shift, an integer (with q=0) */
3387 /* set is the context */
3388 /* returns result */
3390 /* The digits of the coefficient of dfl are shifted to the left (if */
3391 /* dfr is positive) or to the right (if dfr is negative) without */
3392 /* adjusting the exponent or the sign of dfl. */
3394 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3395 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3396 /* dfr must be valid). No status is set unless dfr is invalid or an */
3397 /* operand is an sNaN. The result is canonical. */
3398 /* ------------------------------------------------------------------ */
3399 decFloat
* decFloatShift(decFloat
*result
,
3400 const decFloat
*dfl
, const decFloat
*dfr
,
3402 Int shift
; /* dfr as an Int */
3403 uByte buf
[DECPMAX
*2]; /* coefficient + padding */
3404 uInt digits
, savestat
; /* work */
3405 bcdnum num
; /* .. */
3406 uInt uiwork
; /* for macros */
3408 if (DFISNAN(dfl
)||DFISNAN(dfr
)) return decNaNs(result
, dfl
, dfr
, set
);
3409 if (!DFISINT(dfr
)) return decInvalid(result
, set
);
3410 digits
=decFloatDigits(dfr
); /* calculate digits */
3411 if (digits
>2) return decInvalid(result
, set
); /* definitely out of range */
3412 shift
=DPD2BIN
[DFWORD(dfr
, DECWORDS
-1)&0x3ff]; /* is in bottom declet */
3413 if (shift
>DECPMAX
) return decInvalid(result
, set
); /* too big */
3414 /* [from here on no error or status change is possible] */
3416 if (DFISINF(dfl
)) return decInfinity(result
, dfl
); /* canonical */
3417 /* handle no-shift and all-shift (clear to zero) cases */
3418 if (shift
==0) return decCanonical(result
, dfl
);
3419 if (shift
==DECPMAX
) { /* zero with sign */
3420 uByte sign
=(uByte
)(DFBYTE(dfl
, 0)&0x80); /* save sign bit */
3421 decFloatZero(result
); /* make +0 */
3422 DFBYTE(result
, 0)=(uByte
)(DFBYTE(result
, 0)|sign
); /* and set sign */
3423 /* [cannot safely use CopySign] */
3426 /* a real shift is needed: 0 < shift < DECPMAX */
3427 num
.sign
=DFWORD(dfl
, 0)&DECFLOAT_Sign
;
3428 num
.exponent
=GETEXPUN(dfl
);
3431 if (DFISSIGNED(dfr
)) { /* shift right */
3432 /* edge cases are taken care of, so this is easy */
3433 num
.lsd
=buf
+DECPMAX
-shift
-1;
3435 else { /* shift left -- zero padding needed to right */
3436 UBFROMUI(buf
+DECPMAX
, 0); /* 8 will handle most cases */
3437 UBFROMUI(buf
+DECPMAX
+4, 0); /* .. */
3438 if (shift
>8) memset(buf
+DECPMAX
+8, 0, 8+QUAD
*18); /* all other cases */
3440 num
.lsd
=num
.msd
+DECPMAX
-1;
3442 savestat
=set
->status
; /* record */
3443 decFinalize(result
, &num
, set
);
3444 set
->status
=savestat
; /* restore */
3446 } /* decFloatShift */
3448 /* ------------------------------------------------------------------ */
3449 /* decFloatSubtract -- subtract a decFloat from another */
3451 /* result gets the result of subtracting dfr from dfl: */
3452 /* dfl is the first decFloat (lhs) */
3453 /* dfr is the second decFloat (rhs) */
3454 /* set is the context */
3455 /* returns result */
3457 /* ------------------------------------------------------------------ */
3458 decFloat
* decFloatSubtract(decFloat
*result
,
3459 const decFloat
*dfl
, const decFloat
*dfr
,
3462 /* NaNs must propagate without sign change */
3463 if (DFISNAN(dfr
)) return decFloatAdd(result
, dfl
, dfr
, set
);
3464 temp
=*dfr
; /* make a copy */
3465 DFBYTE(&temp
, 0)^=0x80; /* flip sign */
3466 return decFloatAdd(result
, dfl
, &temp
, set
); /* and add to the lhs */
3467 } /* decFloatSubtract */
3469 /* ------------------------------------------------------------------ */
3470 /* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
3472 /* df is the decFloat to round */
3473 /* set is the context */
3474 /* round is the rounding mode to use */
3475 /* returns a uInt or an Int, rounded according to the name */
3477 /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
3478 /* outside the range of the target; Inexact will not be signaled for */
3479 /* simple rounding unless 'Exact' appears in the name. */
3480 /* ------------------------------------------------------------------ */
3481 uInt
decFloatToUInt32(const decFloat
*df
, decContext
*set
,
3482 enum rounding round
) {
3483 return decToInt32(df
, set
, round
, 0, 1);}
3485 uInt
decFloatToUInt32Exact(const decFloat
*df
, decContext
*set
,
3486 enum rounding round
) {
3487 return decToInt32(df
, set
, round
, 1, 1);}
3489 Int
decFloatToInt32(const decFloat
*df
, decContext
*set
,
3490 enum rounding round
) {
3491 return (Int
)decToInt32(df
, set
, round
, 0, 0);}
3493 Int
decFloatToInt32Exact(const decFloat
*df
, decContext
*set
,
3494 enum rounding round
) {
3495 return (Int
)decToInt32(df
, set
, round
, 1, 0);}
3497 /* ------------------------------------------------------------------ */
3498 /* decFloatToIntegral -- round to integral value (two flavours) */
3500 /* result gets the result */
3501 /* df is the decFloat to round */
3502 /* set is the context */
3503 /* round is the rounding mode to use */
3504 /* returns result */
3506 /* No exceptions, even Inexact, are raised except for sNaN input, or */
3507 /* if 'Exact' appears in the name. */
3508 /* ------------------------------------------------------------------ */
3509 decFloat
* decFloatToIntegralValue(decFloat
*result
, const decFloat
*df
,
3510 decContext
*set
, enum rounding round
) {
3511 return decToIntegral(result
, df
, set
, round
, 0);}
3513 decFloat
* decFloatToIntegralExact(decFloat
*result
, const decFloat
*df
,
3515 return decToIntegral(result
, df
, set
, set
->round
, 1);}
3517 /* ------------------------------------------------------------------ */
3518 /* decFloatXor -- logical digitwise XOR of two decFloats */
3520 /* result gets the result of XORing dfl and dfr */
3521 /* dfl is the first decFloat (lhs) */
3522 /* dfr is the second decFloat (rhs) */
3523 /* set is the context */
3524 /* returns result, which will be canonical with sign=0 */
3526 /* The operands must be positive, finite with exponent q=0, and */
3527 /* comprise just zeros and ones; if not, Invalid operation results. */
3528 /* ------------------------------------------------------------------ */
3529 decFloat
* decFloatXor(decFloat
*result
,
3530 const decFloat
*dfl
, const decFloat
*dfr
,
3532 if (!DFISUINT01(dfl
) || !DFISUINT01(dfr
)
3533 || !DFISCC01(dfl
) || !DFISCC01(dfr
)) return decInvalid(result
, set
);
3534 /* the operands are positive finite integers (q=0) with just 0s and 1s */
3536 DFWORD(result
, 0)=ZEROWORD
3537 |((DFWORD(dfl
, 0) ^ DFWORD(dfr
, 0))&0x04009124);
3538 DFWORD(result
, 1)=(DFWORD(dfl
, 1) ^ DFWORD(dfr
, 1))&0x49124491;
3540 DFWORD(result
, 0)=ZEROWORD
3541 |((DFWORD(dfl
, 0) ^ DFWORD(dfr
, 0))&0x04000912);
3542 DFWORD(result
, 1)=(DFWORD(dfl
, 1) ^ DFWORD(dfr
, 1))&0x44912449;
3543 DFWORD(result
, 2)=(DFWORD(dfl
, 2) ^ DFWORD(dfr
, 2))&0x12449124;
3544 DFWORD(result
, 3)=(DFWORD(dfl
, 3) ^ DFWORD(dfr
, 3))&0x49124491;
3549 /* ------------------------------------------------------------------ */
3550 /* decInvalid -- set Invalid_operation result */
3552 /* result gets a canonical NaN */
3553 /* set is the context */
3554 /* returns result */
3556 /* status has Invalid_operation added */
3557 /* ------------------------------------------------------------------ */
3558 static decFloat
*decInvalid(decFloat
*result
, decContext
*set
) {
3559 decFloatZero(result
);
3560 DFWORD(result
, 0)=DECFLOAT_qNaN
;
3561 set
->status
|=DEC_Invalid_operation
;
3565 /* ------------------------------------------------------------------ */
3566 /* decInfinity -- set canonical Infinity with sign from a decFloat */
3568 /* result gets a canonical Infinity */
3569 /* df is source decFloat (only the sign is used) */
3570 /* returns result */
3572 /* df may be the same as result */
3573 /* ------------------------------------------------------------------ */
3574 static decFloat
*decInfinity(decFloat
*result
, const decFloat
*df
) {
3575 uInt sign
=DFWORD(df
, 0); /* save source signword */
3576 decFloatZero(result
); /* clear everything */
3577 DFWORD(result
, 0)=DECFLOAT_Inf
| (sign
& DECFLOAT_Sign
);
3581 /* ------------------------------------------------------------------ */
3582 /* decNaNs -- handle NaN argument(s) */
3584 /* result gets the result of handling dfl and dfr, one or both of */
3585 /* which is a NaN */
3586 /* dfl is the first decFloat (lhs) */
3587 /* dfr is the second decFloat (rhs) -- may be NULL for a single- */
3588 /* operand operation */
3589 /* set is the context */
3590 /* returns result */
3592 /* Called when one or both operands is a NaN, and propagates the */
3593 /* appropriate result to res. When an sNaN is found, it is changed */
3594 /* to a qNaN and Invalid operation is set. */
3595 /* ------------------------------------------------------------------ */
3596 static decFloat
*decNaNs(decFloat
*result
,
3597 const decFloat
*dfl
, const decFloat
*dfr
,
3599 /* handle sNaNs first */
3600 if (dfr
!=NULL
&& DFISSNAN(dfr
) && !DFISSNAN(dfl
)) dfl
=dfr
; /* use RHS */
3601 if (DFISSNAN(dfl
)) {
3602 decCanonical(result
, dfl
); /* propagate canonical sNaN */
3603 DFWORD(result
, 0)&=~(DECFLOAT_qNaN
^ DECFLOAT_sNaN
); /* quiet */
3604 set
->status
|=DEC_Invalid_operation
;
3607 /* one or both is a quiet NaN */
3608 if (!DFISNAN(dfl
)) dfl
=dfr
; /* RHS must be NaN, use it */
3609 return decCanonical(result
, dfl
); /* propagate canonical qNaN */
3612 /* ------------------------------------------------------------------ */
3613 /* decNumCompare -- numeric comparison of two decFloats */
3615 /* dfl is the left-hand decFloat, which is not a NaN */
3616 /* dfr is the right-hand decFloat, which is not a NaN */
3617 /* tot is 1 for total order compare, 0 for simple numeric */
3618 /* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
3620 /* No error is possible; status and mode are unchanged. */
3621 /* ------------------------------------------------------------------ */
3622 static Int
decNumCompare(const decFloat
*dfl
, const decFloat
*dfr
, Flag tot
) {
3623 Int sigl
, sigr
; /* LHS and RHS non-0 signums */
3624 Int shift
; /* shift needed to align operands */
3625 uByte
*ub
, *uc
; /* work */
3626 uInt uiwork
; /* for macros */
3627 /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
3628 uByte bufl
[DECPMAX
*2+QUAD
*2+4]; /* for LHS coefficient + padding */
3629 uByte bufr
[DECPMAX
*2+QUAD
*2+4]; /* for RHS coefficient + padding */
3632 if (DFISSIGNED(dfl
)) {
3633 if (!DFISSIGNED(dfr
)) { /* -LHS +RHS */
3634 if (DFISZERO(dfl
) && DFISZERO(dfr
) && !tot
) return 0;
3635 return -1; /* RHS wins */
3639 if (DFISSIGNED(dfr
)) {
3640 if (!DFISSIGNED(dfl
)) { /* +LHS -RHS */
3641 if (DFISZERO(dfl
) && DFISZERO(dfr
) && !tot
) return 0;
3642 return +1; /* LHS wins */
3646 /* signs are the same; operand(s) could be zero */
3647 sigr
=-sigl
; /* sign to return if abs(RHS) wins */
3650 if (DFISINF(dfr
)) return 0; /* both infinite & same sign */
3651 return sigl
; /* inf > n */
3653 if (DFISINF(dfr
)) return sigr
; /* n < inf [dfl is finite] */
3655 /* here, both are same sign and finite; calculate their offset */
3656 shift
=GETEXP(dfl
)-GETEXP(dfr
); /* [0 means aligned] */
3657 /* [bias can be ignored -- the absolute exponent is not relevant] */
3659 if (DFISZERO(dfl
)) {
3660 if (!DFISZERO(dfr
)) return sigr
; /* LHS=0, RHS!=0 */
3661 /* both are zero, return 0 if both same exponent or numeric compare */
3662 if (shift
==0 || !tot
) return 0;
3663 if (shift
>0) return sigl
;
3664 return sigr
; /* [shift<0] */
3667 if (DFISZERO(dfr
)) return sigl
; /* LHS!=0, RHS=0 */
3669 /* both are known to be non-zero at this point */
3671 /* if the exponents are so different that the coefficients do not */
3672 /* overlap (by even one digit) then a full comparison is not needed */
3673 if (abs(shift
)>=DECPMAX
) { /* no overlap */
3674 /* coefficients are known to be non-zero */
3675 if (shift
>0) return sigl
;
3676 return sigr
; /* [shift<0] */
3679 /* decode the coefficients */
3680 /* (shift both right two if Quad to make a multiple of four) */
3685 GETCOEFF(dfl
, bufl
+QUAD
*2); /* decode from decFloat */
3686 GETCOEFF(dfr
, bufr
+QUAD
*2); /* .. */
3687 if (shift
==0) { /* aligned; common and easy */
3688 /* all multiples of four, here */
3689 for (ub
=bufl
, uc
=bufr
; ub
<bufl
+DECPMAX
+QUAD
*2; ub
+=4, uc
+=4) {
3691 if (ui
==UBTOUI(uc
)) continue; /* so far so same */
3692 /* about to find a winner; go by bytes in case little-endian */
3693 for (;; ub
++, uc
++) {
3694 if (*ub
>*uc
) return sigl
; /* difference found */
3695 if (*ub
<*uc
) return sigr
; /* .. */
3699 else if (shift
>0) { /* lhs to left */
3700 ub
=bufl
; /* RHS pointer */
3701 /* pad bufl so right-aligned; most shifts will fit in 8 */
3702 UBFROMUI(bufl
+DECPMAX
+QUAD
*2, 0); /* add eight zeros */
3703 UBFROMUI(bufl
+DECPMAX
+QUAD
*2+4, 0); /* .. */
3705 /* more than eight; fill the rest, and also worth doing the */
3706 /* lead-in by fours */
3707 uByte
*up
; /* work */
3708 uByte
*upend
=bufl
+DECPMAX
+QUAD
*2+shift
;
3709 for (up
=bufl
+DECPMAX
+QUAD
*2+8; up
<upend
; up
+=4) UBFROMUI(up
, 0);
3710 /* [pads up to 36 in all for Quad] */
3712 if (UBTOUI(ub
)!=0) return sigl
;
3713 if (ub
+4>bufl
+shift
-4) break;
3716 /* check remaining leading digits */
3717 for (; ub
<bufl
+shift
; ub
++) if (*ub
!=0) return sigl
;
3718 /* now start the overlapped part; bufl has been padded, so the */
3719 /* comparison can go for the full length of bufr, which is a */
3720 /* multiple of 4 bytes */
3721 for (uc
=bufr
; ; uc
+=4, ub
+=4) {
3723 if (ui
!=UBTOUI(uc
)) { /* mismatch found */
3724 for (;; uc
++, ub
++) { /* check from left [little-endian?] */
3725 if (*ub
>*uc
) return sigl
; /* difference found */
3726 if (*ub
<*uc
) return sigr
; /* .. */
3729 if (uc
==bufr
+QUAD
*2+DECPMAX
-4) break; /* all checked */
3733 else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
3734 uc
=bufr
; /* RHS pointer */
3735 /* pad bufr so right-aligned; most shifts will fit in 8 */
3736 UBFROMUI(bufr
+DECPMAX
+QUAD
*2, 0); /* add eight zeros */
3737 UBFROMUI(bufr
+DECPMAX
+QUAD
*2+4, 0); /* .. */
3739 /* more than eight; fill the rest, and also worth doing the */
3740 /* lead-in by fours */
3741 uByte
*up
; /* work */
3742 uByte
*upend
=bufr
+DECPMAX
+QUAD
*2-shift
;
3743 for (up
=bufr
+DECPMAX
+QUAD
*2+8; up
<upend
; up
+=4) UBFROMUI(up
, 0);
3744 /* [pads up to 36 in all for Quad] */
3746 if (UBTOUI(uc
)!=0) return sigr
;
3747 if (uc
+4>bufr
-shift
-4) break;
3750 /* check remaining leading digits */
3751 for (; uc
<bufr
-shift
; uc
++) if (*uc
!=0) return sigr
;
3752 /* now start the overlapped part; bufr has been padded, so the */
3753 /* comparison can go for the full length of bufl, which is a */
3754 /* multiple of 4 bytes */
3755 for (ub
=bufl
; ; ub
+=4, uc
+=4) {
3757 if (ui
!=UBTOUI(uc
)) { /* mismatch found */
3758 for (;; ub
++, uc
++) { /* check from left [little-endian?] */
3759 if (*ub
>*uc
) return sigl
; /* difference found */
3760 if (*ub
<*uc
) return sigr
; /* .. */
3763 if (ub
==bufl
+QUAD
*2+DECPMAX
-4) break; /* all checked */
3767 /* Here when compare equal */
3768 if (!tot
) return 0; /* numerically equal */
3769 /* total ordering .. exponent matters */
3770 if (shift
>0) return sigl
; /* total order by exponent */
3771 if (shift
<0) return sigr
; /* .. */
3773 } /* decNumCompare */
3775 /* ------------------------------------------------------------------ */
3776 /* decToInt32 -- local routine to effect ToInteger conversions */
3778 /* df is the decFloat to convert */
3779 /* set is the context */
3780 /* rmode is the rounding mode to use */
3781 /* exact is 1 if Inexact should be signalled */
3782 /* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
3783 /* returns 32-bit result as a uInt */
3785 /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
3786 /* these cases 0 is returned. */
3787 /* ------------------------------------------------------------------ */
3788 static uInt
decToInt32(const decFloat
*df
, decContext
*set
,
3789 enum rounding rmode
, Flag exact
, Flag unsign
) {
3790 Int exp
; /* exponent */
3791 uInt sourhi
, sourpen
, sourlo
; /* top word from source decFloat .. */
3792 uInt hi
, lo
; /* .. penultimate, least, etc. */
3793 decFloat zero
, result
; /* work */
3796 /* Start decoding the argument */
3797 sourhi
=DFWORD(df
, 0); /* top word */
3798 exp
=DECCOMBEXP
[sourhi
>>26]; /* get exponent high bits (in place) */
3799 if (EXPISSPECIAL(exp
)) { /* is special? */
3800 set
->status
|=DEC_Invalid_operation
; /* signal */
3804 /* Here when the argument is finite */
3805 if (GETEXPUN(df
)==0) result
=*df
; /* already a true integer */
3806 else { /* need to round to integer */
3807 enum rounding saveround
; /* saver */
3808 uInt savestatus
; /* .. */
3809 saveround
=set
->round
; /* save rounding mode .. */
3810 savestatus
=set
->status
; /* .. and status */
3811 set
->round
=rmode
; /* set mode */
3812 decFloatZero(&zero
); /* make 0E+0 */
3813 set
->status
=0; /* clear */
3814 decFloatQuantize(&result
, df
, &zero
, set
); /* [this may fail] */
3815 set
->round
=saveround
; /* restore rounding mode .. */
3816 if (exact
) set
->status
|=savestatus
; /* include Inexact */
3817 else set
->status
=savestatus
; /* .. or just original status */
3820 /* only the last four declets of the coefficient can contain */
3821 /* non-zero; check for others (and also NaN or Infinity from the */
3822 /* Quantize) first (see DFISZERO for explanation): */
3823 /* decFloatShow(&result, "sofar"); */
3825 if ((DFWORD(&result
, 0)&0x1c03ff00)!=0
3826 || (DFWORD(&result
, 0)&0x60000000)==0x60000000) {
3828 if ((DFWORD(&result
, 2)&0xffffff00)!=0
3829 || DFWORD(&result
, 1)!=0
3830 || (DFWORD(&result
, 0)&0x1c003fff)!=0
3831 || (DFWORD(&result
, 0)&0x60000000)==0x60000000) {
3833 set
->status
|=DEC_Invalid_operation
; /* Invalid or out of range */
3836 /* get last twelve digits of the coefficent into hi & ho, base */
3837 /* 10**9 (see GETCOEFFBILL): */
3838 sourlo
=DFWORD(&result
, DECWORDS
-1);
3839 lo
=DPD2BIN0
[sourlo
&0x3ff]
3840 +DPD2BINK
[(sourlo
>>10)&0x3ff]
3841 +DPD2BINM
[(sourlo
>>20)&0x3ff];
3842 sourpen
=DFWORD(&result
, DECWORDS
-2);
3843 hi
=DPD2BIN0
[((sourpen
<<2) | (sourlo
>>30))&0x3ff];
3845 /* according to request, check range carefully */
3847 if (hi
>4 || (hi
==4 && lo
>294967295) || (hi
+lo
!=0 && DFISSIGNED(&result
))) {
3848 set
->status
|=DEC_Invalid_operation
; /* out of range */
3851 return hi
*BILLION
+lo
;
3854 if (hi
>2 || (hi
==2 && lo
>147483647)) {
3855 /* handle the usual edge case */
3856 if (lo
==147483648 && hi
==2 && DFISSIGNED(&result
)) return 0x80000000;
3857 set
->status
|=DEC_Invalid_operation
; /* truly out of range */
3861 if (DFISSIGNED(&result
)) i
=-i
;
3865 /* ------------------------------------------------------------------ */
3866 /* decToIntegral -- local routine to effect ToIntegral value */
3868 /* result gets the result */
3869 /* df is the decFloat to round */
3870 /* set is the context */
3871 /* rmode is the rounding mode to use */
3872 /* exact is 1 if Inexact should be signalled */
3873 /* returns result */
3874 /* ------------------------------------------------------------------ */
3875 static decFloat
* decToIntegral(decFloat
*result
, const decFloat
*df
,
3876 decContext
*set
, enum rounding rmode
,
3878 Int exp
; /* exponent */
3879 uInt sourhi
; /* top word from source decFloat */
3880 enum rounding saveround
; /* saver */
3881 uInt savestatus
; /* .. */
3882 decFloat zero
; /* work */
3884 /* Start decoding the argument */
3885 sourhi
=DFWORD(df
, 0); /* top word */
3886 exp
=DECCOMBEXP
[sourhi
>>26]; /* get exponent high bits (in place) */
3888 if (EXPISSPECIAL(exp
)) { /* is special? */
3889 /* NaNs are handled as usual */
3890 if (DFISNAN(df
)) return decNaNs(result
, df
, NULL
, set
);
3891 /* must be infinite; return canonical infinity with sign of df */
3892 return decInfinity(result
, df
);
3895 /* Here when the argument is finite */
3896 /* complete extraction of the exponent */
3897 exp
+=GETECON(df
)-DECBIAS
; /* .. + continuation and unbias */
3899 if (exp
>=0) return decCanonical(result
, df
); /* already integral */
3901 saveround
=set
->round
; /* save rounding mode .. */
3902 savestatus
=set
->status
; /* .. and status */
3903 set
->round
=rmode
; /* set mode */
3904 decFloatZero(&zero
); /* make 0E+0 */
3905 decFloatQuantize(result
, df
, &zero
, set
); /* 'integrate'; cannot fail */
3906 set
->round
=saveround
; /* restore rounding mode .. */
3907 if (!exact
) set
->status
=savestatus
; /* .. and status, unless exact */
3909 } /* decToIntegral */
This page took 0.148085 seconds and 4 git commands to generate.