projects / deliverable / binutils-gdb.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
merge from gcc
[deliverable/binutils-gdb.git] / libdecnumber / decBasic.c
1 /* Common base code for the decNumber C Library.
2 Copyright (C) 2007 Free Software Foundation, Inc.
3 Contributed by IBM Corporation. Author Mike Cowlishaw.
4
5 This file is part of GCC.
6
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 2, or (at your option) any later
10 version.
11
12 In addition to the permissions in the GNU General Public License,
13 the Free Software Foundation gives you unlimited permission to link
14 the compiled version of this file into combinations with other
15 programs, and to distribute those combinations without any
16 restriction coming from the use of this file. (The General Public
17 License restrictions do apply in other respects; for example, they
18 cover modification of the file, and distribution when not linked
19 into a combine executable.)
20
21 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
22 WARRANTY; without even the implied warranty of MERCHANTABILITY or
23 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24 for more details.
25
26 You should have received a copy of the GNU General Public License
27 along with GCC; see the file COPYING. If not, write to the Free
28 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
29 02110-1301, USA. */
30
31 /* ------------------------------------------------------------------ */
32 /* decBasic.c -- common base code for Basic decimal types */
33 /* ------------------------------------------------------------------ */
34 /* This module comprises code that is shared between decDouble and */
35 /* decQuad (but not decSingle). The main arithmetic operations are */
36 /* here (Add, Subtract, Multiply, FMA, and Division operators). */
37 /* */
38 /* Unlike decNumber, parameterization takes place at compile time */
39 /* rather than at runtime. The parameters are set in the decDouble.c */
40 /* (etc.) files, which then include this one to produce the compiled */
41 /* code. The functions here, therefore, are code shared between */
42 /* multiple formats. */
43 /* */
44 /* This must be included after decCommon.c. */
45 /* ------------------------------------------------------------------ */
46 /* Names here refer to decFloat rather than to decDouble, etc., and */
47 /* the functions are in strict alphabetical order. */
48
49 /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
50 /* decCommon.c */
51 #if !defined(QUAD)
52 #error decBasic.c must be included after decCommon.c
53 #endif
54 #if SINGLE
55 #error Routines in decBasic.c are for decDouble and decQuad only
56 #endif
57
58 /* Private constants */
59 #define DIVIDE 0x80000000 /* Divide operations [as flags] */
60 #define REMAINDER 0x40000000 /* .. */
61 #define DIVIDEINT 0x20000000 /* .. */
62 #define REMNEAR 0x10000000 /* .. */
63
64 /* Private functions (local, used only by routines in this module) */
65 static decFloat *decDivide(decFloat *, const decFloat *,
66 const decFloat *, decContext *, uInt);
67 static decFloat *decCanonical(decFloat *, const decFloat *);
68 static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
69 const decFloat *);
70 static decFloat *decInfinity(decFloat *, const decFloat *);
71 static decFloat *decInvalid(decFloat *, decContext *);
72 static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
73 decContext *);
74 static Int decNumCompare(const decFloat *, const decFloat *, Flag);
75 static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
76 enum rounding, Flag);
77 static uInt decToInt32(const decFloat *, decContext *, enum rounding,
78 Flag, Flag);
79
80 /* ------------------------------------------------------------------ */
81 /* decCanonical -- copy a decFloat, making canonical */
82 /* */
83 /* result gets the canonicalized df */
84 /* df is the decFloat to copy and make canonical */
85 /* returns result */
86 /* */
87 /* This is exposed via decFloatCanonical for Double and Quad only. */
88 /* This works on specials, too; no error or exception is possible. */
89 /* ------------------------------------------------------------------ */
90 static decFloat * decCanonical(decFloat *result, const decFloat *df) {
91 uInt encode, precode, dpd; /* work */
92 uInt inword, uoff, canon; /* .. */
93 Int n; /* counter (down) */
94 if (df!=result) *result=*df; /* effect copy if needed */
95 if (DFISSPECIAL(result)) {
96 if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
97 /* is a NaN */
98 DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */
99 if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */
100 /* drop through to check payload */
101 }
102 /* return quickly if the coefficient continuation is canonical */
103 { /* declare block */
104 #if DOUBLE
105 uInt sourhi=DFWORD(df, 0);
106 uInt sourlo=DFWORD(df, 1);
107 if (CANONDPDOFF(sourhi, 8)
108 && CANONDPDTWO(sourhi, sourlo, 30)
109 && CANONDPDOFF(sourlo, 20)
110 && CANONDPDOFF(sourlo, 10)
111 && CANONDPDOFF(sourlo, 0)) return result;
112 #elif QUAD
113 uInt sourhi=DFWORD(df, 0);
114 uInt sourmh=DFWORD(df, 1);
115 uInt sourml=DFWORD(df, 2);
116 uInt sourlo=DFWORD(df, 3);
117 if (CANONDPDOFF(sourhi, 4)
118 && CANONDPDTWO(sourhi, sourmh, 26)
119 && CANONDPDOFF(sourmh, 16)
120 && CANONDPDOFF(sourmh, 6)
121 && CANONDPDTWO(sourmh, sourml, 28)
122 && CANONDPDOFF(sourml, 18)
123 && CANONDPDOFF(sourml, 8)
124 && CANONDPDTWO(sourml, sourlo, 30)
125 && CANONDPDOFF(sourlo, 20)
126 && CANONDPDOFF(sourlo, 10)
127 && CANONDPDOFF(sourlo, 0)) return result;
128 #endif
129 } /* block */
130
131 /* Loop to repair a non-canonical coefficent, as needed */
132 inword=DECWORDS-1; /* current input word */
133 uoff=0; /* bit offset of declet */
134 encode=DFWORD(result, inword);
135 for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */
136 dpd=encode>>uoff;
137 uoff+=10;
138 if (uoff>32) { /* crossed uInt boundary */
139 inword--;
140 encode=DFWORD(result, inword);
141 uoff-=32;
142 dpd|=encode<<(10-uoff); /* get pending bits */
143 }
144 dpd&=0x3ff; /* clear uninteresting bits */
145 if (dpd<0x16e) continue; /* must be canonical */
146 canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */
147 if (canon==dpd) continue; /* have canonical declet */
148 /* need to replace declet */
149 if (uoff>=10) { /* all within current word */
150 encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
151 encode|=canon<<(uoff-10); /* insert the canonical form */
152 DFWORD(result, inword)=encode; /* .. and save */
153 continue;
154 }
155 /* straddled words */
156 precode=DFWORD(result, inword+1); /* get previous */
157 precode&=0xffffffff>>(10-uoff); /* clear top bits */
158 DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
159 encode&=0xffffffff<<uoff; /* clear bottom bits */
160 encode|=canon>>(10-uoff); /* insert canonical */
161 DFWORD(result, inword)=encode; /* .. and save */
162 } /* n */
163 return result;
164 } /* decCanonical */
165
166 /* ------------------------------------------------------------------ */
167 /* decDivide -- divide operations */
168 /* */
169 /* result gets the result of dividing dfl by dfr: */
170 /* dfl is the first decFloat (lhs) */
171 /* dfr is the second decFloat (rhs) */
172 /* set is the context */
173 /* op is the operation selector */
174 /* returns result */
175 /* */
176 /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
177 /* ------------------------------------------------------------------ */
178 #define DIVCOUNT 0 /* 1 to instrument subtractions counter */
179 #define DIVBASE BILLION /* the base used for divide */
180 #define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
181 #define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */
182 static decFloat * decDivide(decFloat *result, const decFloat *dfl,
183 const decFloat *dfr, decContext *set, uInt op) {
184 decFloat quotient; /* for remainders */
185 bcdnum num; /* for final conversion */
186 uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
187 uInt div[DIVOPLEN]; /* divisor in base-billion .. */
188 uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */
189 uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
190 uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */
191 Int divunits, accunits; /* lengths */
192 Int quodigits; /* digits in quotient */
193 uInt *lsua, *lsuq; /* -> current acc and quo lsus */
194 Int length, multiplier; /* work */
195 uInt carry, sign; /* .. */
196 uInt *ua, *ud, *uq; /* .. */
197 uByte *ub; /* .. */
198 uInt divtop; /* top unit of div adjusted for estimating */
199 #if DIVCOUNT
200 static uInt maxcount=0; /* worst-seen subtractions count */
201 uInt divcount=0; /* subtractions count [this divide] */
202 #endif
203
204 /* calculate sign */
205 num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
206
207 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
208 /* NaNs are handled as usual */
209 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
210 /* one or two infinities */
211 if (DFISINF(dfl)) {
212 if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
213 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
214 /* Infinity/x is infinite and quiet, even if x=0 */
215 DFWORD(result, 0)=num.sign;
216 return decInfinity(result, result);
217 }
218 /* must be x/Infinity -- remainders are lhs */
219 if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
220 /* divides: return zero with correct sign and exponent depending */
221 /* on op (Etiny for divide, 0 for divideInt) */
222 decFloatZero(result);
223 if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
224 else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */
225 return result;
226 }
227 /* next, handle zero operands (x/0 and 0/x) */
228 if (DFISZERO(dfr)) { /* x/0 */
229 if (DFISZERO(dfl)) { /* 0/0 is undefined */
230 decFloatZero(result);
231 DFWORD(result, 0)=DECFLOAT_qNaN;
232 set->status|=DEC_Division_undefined;
233 return result;
234 }
235 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
236 set->status|=DEC_Division_by_zero;
237 DFWORD(result, 0)=num.sign;
238 return decInfinity(result, result); /* x/0 -> signed Infinity */
239 }
240 num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
241 if (DFISZERO(dfl)) { /* 0/x (x!=0) */
242 /* if divide, result is 0 with ideal exponent; divideInt has */
243 /* exponent=0, remainders give zero with lower exponent */
244 if (op&DIVIDEINT) {
245 decFloatZero(result);
246 DFWORD(result, 0)|=num.sign; /* add sign */
247 return result;
248 }
249 if (!(op&DIVIDE)) { /* a remainder */
250 /* exponent is the minimum of the operands */
251 num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
252 /* if the result is zero the sign shall be sign of dfl */
253 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
254 }
255 bcdacc[0]=0;
256 num.msd=bcdacc; /* -> 0 */
257 num.lsd=bcdacc; /* .. */
258 return decFinalize(result, &num, set); /* [divide may clamp exponent] */
259 } /* 0/x */
260 /* [here, both operands are known to be finite and non-zero] */
261
262 /* extract the operand coefficents into 'units' which are */
263 /* base-billion; the lhs is high-aligned in acc and the msu of both */
264 /* acc and div is at the right-hand end of array (offset length-1); */
265 /* the quotient can need one more unit than the operands as digits */
266 /* in it are not necessarily aligned neatly; further, the quotient */
267 /* may not start accumulating until after the end of the initial */
268 /* operand in acc if that is small (e.g., 1) so the accumulator */
269 /* must have at least that number of units extra (at the ls end) */
270 GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
271 GETCOEFFBILL(dfr, div);
272 /* zero the low uInts of acc */
273 acc[0]=0;
274 acc[1]=0;
275 acc[2]=0;
276 acc[3]=0;
277 #if DOUBLE
278 #if DIVOPLEN!=2
279 #error Unexpected Double DIVOPLEN
280 #endif
281 #elif QUAD
282 acc[4]=0;
283 acc[5]=0;
284 acc[6]=0;
285 acc[7]=0;
286 #if DIVOPLEN!=4
287 #error Unexpected Quad DIVOPLEN
288 #endif
289 #endif
290
291 /* set msu and lsu pointers */
292 msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
293 msuq=quo+DIVOPLEN;
294 /*[loop for div will terminate because operands are non-zero] */
295 for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
296 /* the initial least-significant unit of acc is set so acc appears */
297 /* to have the same length as div. */
298 /* This moves one position towards the least possible for each */
299 /* iteration */
300 divunits=(Int)(msud-div+1); /* precalculate */
301 lsua=msua-divunits+1; /* initial working lsu of acc */
302 lsuq=msuq; /* and of quo */
303
304 /* set up the estimator for the multiplier; this is the msu of div, */
305 /* plus two bits from the unit below (if any) rounded up by one if */
306 /* there are any non-zero bits or units below that [the extra two */
307 /* bits makes for a much better estimate when the top unit is small] */
308 divtop=*msud<<2;
309 if (divunits>1) {
310 uInt *um=msud-1;
311 uInt d=*um;
312 if (d>=750000000) {divtop+=3; d-=750000000;}
313 else if (d>=500000000) {divtop+=2; d-=500000000;}
314 else if (d>=250000000) {divtop++; d-=250000000;}
315 if (d) divtop++;
316 else for (um--; um>=div; um--) if (*um) {
317 divtop++;
318 break;
319 }
320 } /* >1 unit */
321
322 #if DECTRACE
323 {Int i;
324 printf("----- div=");
325 for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
326 printf("\n");}
327 #endif
328
329 /* now collect up to DECPMAX+1 digits in the quotient (this may */
330 /* need OPLEN+1 uInts if unaligned) */
331 quodigits=0; /* no digits yet */
332 for (;; lsua--) { /* outer loop -- each input position */
333 #if DECCHECK
334 if (lsua<acc) {
335 printf("Acc underrun...\n");
336 break;
337 }
338 #endif
339 #if DECTRACE
340 printf("Outer: quodigits=%ld acc=", (LI)quodigits);
341 for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
342 printf("\n");
343 #endif
344 *lsuq=0; /* default unit result is 0 */
345 for (;;) { /* inner loop -- calculate quotient unit */
346 /* strip leading zero units from acc (either there initially or */
347 /* from subtraction below); this may strip all if exactly 0 */
348 for (; *msua==0 && msua>=lsua;) msua--;
349 accunits=(Int)(msua-lsua+1); /* [maybe 0] */
350 /* subtraction is only necessary and possible if there are as */
351 /* least as many units remaining in acc for this iteration as */
352 /* there are in div */
353 if (accunits<divunits) {
354 if (accunits==0) msua++; /* restore */
355 break;
356 }
357
358 /* If acc is longer than div then subtraction is definitely */
359 /* possible (as msu of both is non-zero), but if they are the */
360 /* same length a comparison is needed. */
361 /* If a subtraction is needed then a good estimate of the */
362 /* multiplier for the subtraction is also needed in order to */
363 /* minimise the iterations of this inner loop because the */
364 /* subtractions needed dominate division performance. */
365 if (accunits==divunits) {
366 /* compare the high divunits of acc and div: */
367 /* acc<div: this quotient unit is unchanged; subtraction */
368 /* will be possible on the next iteration */
369 /* acc==div: quotient gains 1, set acc=0 */
370 /* acc>div: subtraction necessary at this position */
371 for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
372 /* [now at first mismatch or lsu] */
373 if (*ud>*ua) break; /* next time... */
374 if (*ud==*ua) { /* all compared equal */
375 *lsuq+=1; /* increment result */
376 msua=lsua; /* collapse acc units */
377 *msua=0; /* .. to a zero */
378 break;
379 }
380
381 /* subtraction necessary; estimate multiplier [see above] */
382 /* if both *msud and *msua are small it is cost-effective to */
383 /* bring in part of the following units (if any) to get a */
384 /* better estimate (assume some other non-zero in div) */
385 #define DIVLO 1000000U
386 #define DIVHI (DIVBASE/DIVLO)
387 #if DECUSE64
388 if (divunits>1) {
389 /* there cannot be a *(msud-2) for DECDOUBLE so next is */
390 /* an exact calculation unless DECQUAD (which needs to */
391 /* assume bits out there if divunits>2) */
392 uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
393 uLong div=(uLong)*msud * DIVBASE + *(msud-1);
394 #if QUAD
395 if (divunits>2) div++;
396 #endif
397 mul/=div;
398 multiplier=(Int)mul;
399 }
400 else multiplier=*msua/(*msud);
401 #else
402 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
403 multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
404 /(*msud*DIVHI + *(msud-1)/DIVLO +1);
405 }
406 else multiplier=(*msua<<2)/divtop;
407 #endif
408 }
409 else { /* accunits>divunits */
410 /* msud is one unit 'lower' than msua, so estimate differently */
411 #if DECUSE64
412 uLong mul;
413 /* as before, bring in extra digits if possible */
414 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
415 mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
416 + *(msua-2)/DIVLO;
417 mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
418 }
419 else if (divunits==1) {
420 mul=(uLong)*msua * DIVBASE + *(msua-1);
421 mul/=*msud; /* no more to the right */
422 }
423 else {
424 mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) + (*(msua-1)<<2);
425 mul/=divtop; /* [divtop already allows for sticky bits] */
426 }
427 multiplier=(Int)mul;
428 #else
429 multiplier=*msua * ((DIVBASE<<2)/divtop);
430 #endif
431 }
432 if (multiplier==0) multiplier=1; /* marginal case */
433 *lsuq+=multiplier;
434
435 #if DIVCOUNT
436 /* printf("Multiplier: %ld\n", (LI)multiplier); */
437 divcount++;
438 #endif
439
440 /* Carry out the subtraction acc-(div*multiplier); for each */
441 /* unit in div, do the multiply, split to units (see */
442 /* decFloatMultiply for the algorithm), and subtract from acc */
443 #define DIVMAGIC 2305843009U /* 2**61/10**9 */
444 #define DIVSHIFTA 29
445 #define DIVSHIFTB 32
446 carry=0;
447 for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
448 uInt lo, hop;
449 #if DECUSE64
450 uLong sub=(uLong)multiplier*(*ud)+carry;
451 if (sub<DIVBASE) {
452 carry=0;
453 lo=(uInt)sub;
454 }
455 else {
456 hop=(uInt)(sub>>DIVSHIFTA);
457 carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
458 /* the estimate is now in hi; now calculate sub-hi*10**9 */
459 /* to get the remainder (which will be <DIVBASE)) */
460 lo=(uInt)sub;
461 lo-=carry*DIVBASE; /* low word of result */
462 if (lo>=DIVBASE) {
463 lo-=DIVBASE; /* correct by +1 */
464 carry++;
465 }
466 }
467 #else /* 32-bit */
468 uInt hi;
469 /* calculate multiplier*(*ud) into hi and lo */
470 LONGMUL32HI(hi, *ud, multiplier); /* get the high word */
471 lo=multiplier*(*ud); /* .. and the low */
472 lo+=carry; /* add the old hi */
473 carry=hi+(lo<carry); /* .. with any carry */
474 if (carry || lo>=DIVBASE) { /* split is needed */
475 hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */
476 LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */
477 /* [DIVSHIFTB is 32, so carry can be used directly] */
478 /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
479 /* happily the top word of the result is irrelevant because it */
480 /* will always be zero so this needs only one multiplication */
481 lo-=(carry*DIVBASE);
482 /* the correction here will be at most +1; do it */
483 if (lo>=DIVBASE) {
484 lo-=DIVBASE;
485 carry++;
486 }
487 }
488 #endif
489 if (lo>*ua) { /* borrow needed */
490 *ua+=DIVBASE;
491 carry++;
492 }
493 *ua-=lo;
494 } /* ud loop */
495 if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */
496 } /* inner loop */
497
498 /* the outer loop terminates when there is either an exact result */
499 /* or enough digits; first update the quotient digit count and */
500 /* pointer (if any significant digits) */
501 #if DECTRACE
502 if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
503 #endif
504 if (quodigits) {
505 quodigits+=9; /* had leading unit earlier */
506 lsuq--;
507 if (quodigits>DECPMAX+1) break; /* have enough */
508 }
509 else if (*lsuq) { /* first quotient digits */
510 const uInt *pow;
511 for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
512 lsuq--;
513 /* [cannot have >DECPMAX+1 on first unit] */
514 }
515
516 if (*msua!=0) continue; /* not an exact result */
517 /* acc is zero iff used all of original units and zero down to lsua */
518 /* (must also continue to original lsu for correct quotient length) */
519 if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
520 for (; msua>lsua && *msua==0;) msua--;
521 if (*msua==0 && msua==lsua) break;
522 } /* outer loop */
523
524 /* all of the original operand in acc has been covered at this point */
525 /* quotient now has at least DECPMAX+2 digits */
526 /* *msua is now non-0 if inexact and sticky bits */
527 /* lsuq is one below the last uint of the quotient */
528 lsuq++; /* set -> true lsu of quo */
529 if (*msua) *lsuq|=1; /* apply sticky bit */
530
531 /* quo now holds the (unrounded) quotient in base-billion; one */
532 /* base-billion 'digit' per uInt. */
533 #if DECTRACE
534 printf("DivQuo:");
535 for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
536 printf("\n");
537 #endif
538
539 /* Now convert to BCD for rounding and cleanup, starting from the */
540 /* most significant end [offset by one into bcdacc to leave room */
541 /* for a possible carry digit if rounding for REMNEAR is needed] */
542 for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
543 uInt top, mid, rem; /* work */
544 if (*uq==0) { /* no split needed */
545 UINTAT(ub)=0; /* clear 9 BCD8s */
546 UINTAT(ub+4)=0; /* .. */
547 *(ub+8)=0; /* .. */
548 continue;
549 }
550 /* *uq is non-zero -- split the base-billion digit into */
551 /* hi, mid, and low three-digits */
552 #define divsplit9 1000000 /* divisor */
553 #define divsplit6 1000 /* divisor */
554 /* The splitting is done by simple divides and remainders, */
555 /* assuming the compiler will optimize these [GCC does] */
556 top=*uq/divsplit9;
557 rem=*uq%divsplit9;
558 mid=rem/divsplit6;
559 rem=rem%divsplit6;
560 /* lay out the nine BCD digits (plus one unwanted byte) */
561 UINTAT(ub) =UINTAT(&BIN2BCD8[top*4]);
562 UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
563 UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
564 } /* BCD conversion loop */
565 ub--; /* -> lsu */
566
567 /* complete the bcdnum; quodigits is correct, so the position of */
568 /* the first non-zero is known */
569 num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
570 num.lsd=ub;
571
572 /* make exponent adjustments, etc */
573 if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */
574 num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
575 /* if the result was exact then there may be up to 8 extra */
576 /* trailing zeros in the overflowed quotient final unit */
577 if (*msua==0) {
578 for (; *ub==0;) ub--; /* drop zeros */
579 num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */
580 num.lsd=ub;
581 }
582 } /* adjustment needed */
583
584 #if DIVCOUNT
585 if (divcount>maxcount) { /* new high-water nark */
586 maxcount=divcount;
587 printf("DivNewMaxCount: %ld\n", (LI)maxcount);
588 }
589 #endif
590
591 if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
592
593 /* Is DIVIDEINT or a remainder; there is more to do -- first form */
594 /* the integer (this is done 'after the fact', unlike as in */
595 /* decNumber, so as not to tax DIVIDE) */
596
597 /* The first non-zero digit will be in the first 9 digits, known */
598 /* from quodigits and num.msd, so there is always space for DECPMAX */
599 /* digits */
600
601 length=(Int)(num.lsd-num.msd+1);
602 /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
603
604 if (length+num.exponent>DECPMAX) { /* cannot fit */
605 decFloatZero(result);
606 DFWORD(result, 0)=DECFLOAT_qNaN;
607 set->status|=DEC_Division_impossible;
608 return result;
609 }
610
611 if (num.exponent>=0) { /* already an int, or need pad zeros */
612 for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
613 num.lsd+=num.exponent;
614 }
615 else { /* too long: round or truncate needed */
616 Int drop=-num.exponent;
617 if (!(op&REMNEAR)) { /* simple truncate */
618 num.lsd-=drop;
619 if (num.lsd<num.msd) { /* truncated all */
620 num.lsd=num.msd; /* make 0 */
621 *num.lsd=0; /* .. [sign still relevant] */
622 }
623 }
624 else { /* round to nearest even [sigh] */
625 /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
626 /* (this is a special case of Quantize -- q.v. for commentary) */
627 uByte *roundat; /* -> re-round digit */
628 uByte reround; /* reround value */
629 *(num.msd-1)=0; /* in case of left carry, or make 0 */
630 if (drop<length) roundat=num.lsd-drop+1;
631 else if (drop==length) roundat=num.msd;
632 else roundat=num.msd-1; /* [-> 0] */
633 reround=*roundat;
634 for (ub=roundat+1; ub<=num.lsd; ub++) {
635 if (*ub!=0) {
636 reround=DECSTICKYTAB[reround];
637 break;
638 }
639 } /* check stickies */
640 if (roundat>num.msd) num.lsd=roundat-1;
641 else {
642 num.msd--; /* use the 0 .. */
643 num.lsd=num.msd; /* .. at the new MSD place */
644 }
645 if (reround!=0) { /* discarding non-zero */
646 uInt bump=0;
647 /* rounding is DEC_ROUND_HALF_EVEN always */
648 if (reround>5) bump=1; /* >0.5 goes up */
649 else if (reround==5) /* exactly 0.5000 .. */
650 bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */
651 if (bump!=0) { /* need increment */
652 /* increment the coefficient; this might end up with 1000... */
653 ub=num.lsd;
654 for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
655 for (; *ub==9; ub--) *ub=0; /* at most 3 more */
656 *ub+=1;
657 if (ub<num.msd) num.msd--; /* carried */
658 } /* bump needed */
659 } /* reround!=0 */
660 } /* remnear */
661 } /* round or truncate needed */
662 num.exponent=0; /* all paths */
663 /*decShowNum(&num, "int"); */
664
665 if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
666
667 /* Have a remainder to calculate */
668 decFinalize(&quotient, &num, set); /* lay out the integer so far */
669 DFWORD(&quotient, 0)^=DECFLOAT_Sign; /* negate it */
670 sign=DFWORD(dfl, 0); /* save sign of dfl */
671 decFloatFMA(result, &quotient, dfr, dfl, set);
672 if (!DFISZERO(result)) return result;
673 /* if the result is zero the sign shall be sign of dfl */
674 DFWORD(&quotient, 0)=sign; /* construct decFloat of sign */
675 return decFloatCopySign(result, result, &quotient);
676 } /* decDivide */
677
678 /* ------------------------------------------------------------------ */
679 /* decFiniteMultiply -- multiply two finite decFloats */
680 /* */
681 /* num gets the result of multiplying dfl and dfr */
682 /* bcdacc .. with the coefficient in this array */
683 /* dfl is the first decFloat (lhs) */
684 /* dfr is the second decFloat (rhs) */
685 /* */
686 /* This effects the multiplication of two decFloats, both known to be */
687 /* finite, leaving the result in a bcdnum ready for decFinalize (for */
688 /* use in Multiply) or in a following addition (FMA). */
689 /* */
690 /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
691 /* No error is possible and no status is set. */
692 /* ------------------------------------------------------------------ */
693 /* This routine has two separate implementations of the core */
694 /* multiplication; both using base-billion. One uses only 32-bit */
695 /* variables (Ints and uInts) or smaller; the other uses uLongs (for */
696 /* multiplication and addition only). Both implementations cover */
697 /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
698 /* comparisons. In any one compilation only one implementation for */
699 /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
700 /* version is forced. */
701 /* */
702 /* Historical note: an earlier version of this code also supported the */
703 /* 256-bit format and has been preserved. That is somewhat trickier */
704 /* during lazy carry splitting because the initial quotient estimate */
705 /* (est) can exceed 32 bits. */
706
707 #define MULTBASE BILLION /* the base used for multiply */
708 #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
709 #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */
710 #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
711
712 /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
713 #if DECEMAXD>9
714 #error Exponent may overflow when doubled for Multiply
715 #endif
716 #if MULACCLEN!=(MULACCLEN/4)*4
717 /* This assumption is used below only for initialization */
718 #error MULACCLEN is not a multiple of 4
719 #endif
720
721 static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
722 const decFloat *dfl, const decFloat *dfr) {
723 uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
724 uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
725 uInt *ui, *uj; /* work */
726 uByte *ub; /* .. */
727
728 #if DECUSE64
729 uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
730 uLong *pl; /* work -> lazy accumulator */
731 uInt acc[MULACCLEN]; /* coefficent in base-billion .. */
732 #else
733 uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */
734 #endif
735 uInt *pa; /* work -> accumulator */
736 /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
737
738 /* Calculate sign and exponent */
739 num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
740 num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
741
742 /* Extract the coefficients and prepare the accumulator */
743 /* the coefficients of the operands are decoded into base-billion */
744 /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
745 /* appropriate size. */
746 GETCOEFFBILL(dfl, bufl);
747 GETCOEFFBILL(dfr, bufr);
748 #if DECTRACE && 0
749 printf("CoeffbL:");
750 for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
751 printf("\n");
752 printf("CoeffbR:");
753 for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
754 printf("\n");
755 #endif
756
757 /* start the 64-bit/32-bit differing paths... */
758 #if DECUSE64
759
760 /* zero the accumulator */
761 #if MULACCLEN==4
762 accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
763 #else /* use a loop */
764 /* MULACCLEN is a multiple of four, asserted above */
765 for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
766 *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
767 } /* pl */
768 #endif
769
770 /* Effect the multiplication */
771 /* The multiplcation proceeds using MFC's lazy-carry resolution */
772 /* algorithm from decNumber. First, the multiplication is */
773 /* effected, allowing accumulation of the partial products (which */
774 /* are in base-billion at each column position) into 64 bits */
775 /* without resolving back to base=billion after each addition. */
776 /* These 64-bit numbers (which may contain up to 19 decimal digits) */
777 /* are then split using the Clark & Cowlishaw algorithm (see below). */
778 /* [Testing for 0 in the inner loop is not really a 'win'] */
779 for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
780 if (*ui==0) continue; /* product cannot affect result */
781 pl=accl+(ui-bufr); /* where to add the lhs */
782 for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
783 /* if (*uj==0) continue; // product cannot affect result */
784 *pl+=((uLong)*ui)*(*uj);
785 } /* uj */
786 } /* ui */
787
788 /* The 64-bit carries must now be resolved; this means that a */
789 /* quotient/remainder has to be calculated for base-billion (1E+9). */
790 /* For this, Clark & Cowlishaw's quotient estimation approach (also */
791 /* used in decNumber) is needed, because 64-bit divide is generally */
792 /* extremely slow on 32-bit machines, and may be slower than this */
793 /* approach even on 64-bit machines. This algorithm splits X */
794 /* using: */
795 /* */
796 /* magic=2**(A+B)/1E+9; // 'magic number' */
797 /* hop=X/2**A; // high order part of X (by shift) */
798 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
799 /* */
800 /* A and B are quite constrained; hop and magic must fit in 32 bits, */
801 /* and 2**(A+B) must be as large as possible (which is 2**61 if */
802 /* magic is to fit). Further, maxX increases with the length of */
803 /* the operands (and hence the number of partial products */
804 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
805 /* */
806 /* It can be shown that when OPLEN is 2 then the maximum error in */
807 /* the estimated quotient is <1, but for larger maximum x the */
808 /* maximum error is above 1 so a correction that is >1 may be */
809 /* needed. Values of A and B are chosen to satisfy the constraints */
810 /* just mentioned while minimizing the maximum error (and hence the */
811 /* maximum correction), as shown in the following table: */
812 /* */
813 /* Type OPLEN A B maxX maxError maxCorrection */
814 /* --------------------------------------------------------- */
815 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
816 /* QUAD 4 30 31 <4*10**18 1.17 2 */
817 /* */
818 /* In the OPLEN==2 case there is most choice, but the value for B */
819 /* of 32 has a big advantage as then the calculation of the */
820 /* estimate requires no shifting; the compiler can extract the high */
821 /* word directly after multiplying magic*hop. */
822 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
823 #if DOUBLE
824 #define MULSHIFTA 29
825 #define MULSHIFTB 32
826 #elif QUAD
827 #define MULSHIFTA 30
828 #define MULSHIFTB 31
829 #else
830 #error Unexpected type
831 #endif
832
833 #if DECTRACE
834 printf("MulAccl:");
835 for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
836 printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
837 printf("\n");
838 #endif
839
840 for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
841 uInt lo, hop; /* work */
842 uInt est; /* cannot exceed 4E+9 */
843 if (*pl>MULTBASE) {
844 /* *pl holds a binary number which needs to be split */
845 hop=(uInt)(*pl>>MULSHIFTA);
846 est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
847 /* the estimate is now in est; now calculate hi:lo-est*10**9; */
848 /* happily the top word of the result is irrelevant because it */
849 /* will always be zero so this needs only one multiplication */
850 lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */
851 /* If QUAD, the correction here could be +2 */
852 if (lo>=MULTBASE) {
853 lo-=MULTBASE; /* correct by +1 */
854 est++;
855 #if QUAD
856 /* may need to correct by +2 */
857 if (lo>=MULTBASE) {
858 lo-=MULTBASE;
859 est++;
860 }
861 #endif
862 }
863 /* finally place lo as the new coefficient 'digit' and add est to */
864 /* the next place up [this is safe because this path is never */
865 /* taken on the final iteration as *pl will fit] */
866 *pa=lo;
867 *(pl+1)+=est;
868 } /* *pl needed split */
869 else { /* *pl<MULTBASE */
870 *pa=(uInt)*pl; /* just copy across */
871 }
872 } /* pl loop */
873
874 #else /* 32-bit */
875 for (pa=acc;; pa+=4) { /* zero the accumulator */
876 *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */
877 if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */
878 } /* pa */
879
880 /* Effect the multiplication */
881 /* uLongs are not available (and in particular, there is no uLong */
882 /* divide) but it is still possible to use MFC's lazy-carry */
883 /* resolution algorithm from decNumber. First, the multiplication */
884 /* is effected, allowing accumulation of the partial products */
885 /* (which are in base-billion at each column position) into 64 bits */
886 /* [with the high-order 32 bits in each position being held at */
887 /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
888 /* These 64-bit numbers (which may contain up to 19 decimal digits) */
889 /* are then split using the Clark & Cowlishaw algorithm (see */
890 /* below). */
891 for (ui=bufr;; ui++) { /* over each item in rhs */
892 uInt hi, lo; /* words of exact multiply result */
893 pa=acc+(ui-bufr); /* where to add the lhs */
894 for (uj=bufl;; uj++, pa++) { /* over each item in lhs */
895 LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */
896 lo=(*ui)*(*uj); /* .. */
897 *pa+=lo; /* accumulate low bits and .. */
898 *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */
899 if (uj==bufl+MULOPLEN-1) break;
900 }
901 if (ui==bufr+MULOPLEN-1) break;
902 }
903
904 /* The 64-bit carries must now be resolved; this means that a */
905 /* quotient/remainder has to be calculated for base-billion (1E+9). */
906 /* For this, Clark & Cowlishaw's quotient estimation approach (also */
907 /* used in decNumber) is needed, because 64-bit divide is generally */
908 /* extremely slow on 32-bit machines. This algorithm splits X */
909 /* using: */
910 /* */
911 /* magic=2**(A+B)/1E+9; // 'magic number' */
912 /* hop=X/2**A; // high order part of X (by shift) */
913 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
914 /* */
915 /* A and B are quite constrained; hop and magic must fit in 32 bits, */
916 /* and 2**(A+B) must be as large as possible (which is 2**61 if */
917 /* magic is to fit). Further, maxX increases with the length of */
918 /* the operands (and hence the number of partial products */
919 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
920 /* */
921 /* It can be shown that when OPLEN is 2 then the maximum error in */
922 /* the estimated quotient is <1, but for larger maximum x the */
923 /* maximum error is above 1 so a correction that is >1 may be */
924 /* needed. Values of A and B are chosen to satisfy the constraints */
925 /* just mentioned while minimizing the maximum error (and hence the */
926 /* maximum correction), as shown in the following table: */
927 /* */
928 /* Type OPLEN A B maxX maxError maxCorrection */
929 /* --------------------------------------------------------- */
930 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
931 /* QUAD 4 30 31 <4*10**18 1.17 2 */
932 /* */
933 /* In the OPLEN==2 case there is most choice, but the value for B */
934 /* of 32 has a big advantage as then the calculation of the */
935 /* estimate requires no shifting; the high word is simply */
936 /* calculated from multiplying magic*hop. */
937 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
938 #if DOUBLE
939 #define MULSHIFTA 29
940 #define MULSHIFTB 32
941 #elif QUAD
942 #define MULSHIFTA 30
943 #define MULSHIFTB 31
944 #else
945 #error Unexpected type
946 #endif
947
948 #if DECTRACE
949 printf("MulHiLo:");
950 for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
951 printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
952 printf("\n");
953 #endif
954
955 for (pa=acc;; pa++) { /* each low uInt */
956 uInt hi, lo; /* words of exact multiply result */
957 uInt hop, estlo; /* work */
958 #if QUAD
959 uInt esthi; /* .. */
960 #endif
961
962 lo=*pa;
963 hi=*(pa+MULACCLEN); /* top 32 bits */
964 /* hi and lo now hold a binary number which needs to be split */
965
966 #if DOUBLE
967 hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */
968 LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
969 /* [MULSHIFTB is 32, so estlo can be used directly] */
970 /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
971 /* happily the top word of the result is irrelevant because it */
972 /* will always be zero so this needs only one multiplication */
973 lo-=(estlo*MULTBASE);
974 /* esthi=0; // high word is ignored below */
975 /* the correction here will be at most +1; do it */
976 if (lo>=MULTBASE) {
977 lo-=MULTBASE;
978 estlo++;
979 }
980 #elif QUAD
981 hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */
982 LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
983 estlo=hop*MULMAGIC; /* .. so low word needed */
984 estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
985 /* esthi=0; // high word is ignored below */
986 lo-=(estlo*MULTBASE); /* as above */
987 /* the correction here could be +1 or +2 */
988 if (lo>=MULTBASE) {
989 lo-=MULTBASE;
990 estlo++;
991 }
992 if (lo>=MULTBASE) {
993 lo-=MULTBASE;
994 estlo++;
995 }
996 #else
997 #error Unexpected type
998 #endif
999
1000 /* finally place lo as the new accumulator digit and add est to */
1001 /* the next place up; this latter add could cause a carry of 1 */
1002 /* to the high word of the next place */
1003 *pa=lo;
1004 *(pa+1)+=estlo;
1005 /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
1006 /* *(pa+1+MULACCLEN)+=esthi; */
1007 if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
1008 if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */
1009 } /* pa loop */
1010 #endif
1011
1012 /* At this point, whether using the 64-bit or the 32-bit paths, the */
1013 /* accumulator now holds the (unrounded) result in base-billion; */
1014 /* one base-billion 'digit' per uInt. */
1015 #if DECTRACE
1016 printf("MultAcc:");
1017 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
1018 printf("\n");
1019 #endif
1020
1021 /* Now convert to BCD for rounding and cleanup, starting from the */
1022 /* most significant end */
1023 pa=acc+MULACCLEN-1;
1024 if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
1025 else { /* >=1 word of leading zeros */
1026 num->msd=bcdacc; /* known leading zeros are gone */
1027 pa--; /* skip first word .. */
1028 for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
1029 }
1030 for (ub=bcdacc;; pa--, ub+=9) {
1031 if (*pa!=0) { /* split(s) needed */
1032 uInt top, mid, rem; /* work */
1033 /* *pa is non-zero -- split the base-billion acc digit into */
1034 /* hi, mid, and low three-digits */
1035 #define mulsplit9 1000000 /* divisor */
1036 #define mulsplit6 1000 /* divisor */
1037 /* The splitting is done by simple divides and remainders, */
1038 /* assuming the compiler will optimize these where useful */
1039 /* [GCC does] */
1040 top=*pa/mulsplit9;
1041 rem=*pa%mulsplit9;
1042 mid=rem/mulsplit6;
1043 rem=rem%mulsplit6;
1044 /* lay out the nine BCD digits (plus one unwanted byte) */
1045 UINTAT(ub) =UINTAT(&BIN2BCD8[top*4]);
1046 UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
1047 UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
1048 }
1049 else { /* *pa==0 */
1050 UINTAT(ub)=0; /* clear 9 BCD8s */
1051 UINTAT(ub+4)=0; /* .. */
1052 *(ub+8)=0; /* .. */
1053 }
1054 if (pa==acc) break;
1055 } /* BCD conversion loop */
1056
1057 num->lsd=ub+8; /* complete the bcdnum .. */
1058
1059 #if DECTRACE
1060 decShowNum(num, "postmult");
1061 decFloatShow(dfl, "dfl");
1062 decFloatShow(dfr, "dfr");
1063 #endif
1064 return;
1065 } /* decFiniteMultiply */
1066
1067 /* ------------------------------------------------------------------ */
1068 /* decFloatAbs -- absolute value, heeding NaNs, etc. */
1069 /* */
1070 /* result gets the canonicalized df with sign 0 */
1071 /* df is the decFloat to abs */
1072 /* set is the context */
1073 /* returns result */
1074 /* */
1075 /* This has the same effect as decFloatPlus unless df is negative, */
1076 /* in which case it has the same effect as decFloatMinus. The */
1077 /* effect is also the same as decFloatCopyAbs except that NaNs are */
1078 /* handled normally (the sign of a NaN is not affected, and an sNaN */
1079 /* will signal) and the result will be canonical. */
1080 /* ------------------------------------------------------------------ */
1081 decFloat * decFloatAbs(decFloat *result, const decFloat *df,
1082 decContext *set) {
1083 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
1084 decCanonical(result, df); /* copy and check */
1085 DFBYTE(result, 0)&=~0x80; /* zero sign bit */
1086 return result;
1087 } /* decFloatAbs */
1088
1089 /* ------------------------------------------------------------------ */
1090 /* decFloatAdd -- add two decFloats */
1091 /* */
1092 /* result gets the result of adding dfl and dfr: */
1093 /* dfl is the first decFloat (lhs) */
1094 /* dfr is the second decFloat (rhs) */
1095 /* set is the context */
1096 /* returns result */
1097 /* */
1098 /* ------------------------------------------------------------------ */
1099 decFloat * decFloatAdd(decFloat *result,
1100 const decFloat *dfl, const decFloat *dfr,
1101 decContext *set) {
1102 bcdnum num; /* for final conversion */
1103 Int expl, expr; /* left and right exponents */
1104 uInt *ui, *uj; /* work */
1105 uByte *ub; /* .. */
1106
1107 uInt sourhil, sourhir; /* top words from source decFloats */
1108 /* [valid only until specials */
1109 /* handled or exponents decoded] */
1110 uInt diffsign; /* non-zero if signs differ */
1111 uInt carry; /* carry: 0 or 1 before add loop */
1112 Int overlap; /* coefficient overlap (if full) */
1113 /* the following buffers hold coefficients with various alignments */
1114 /* (see commentary and diagrams below) */
1115 uByte acc[4+2+DECPMAX*3+8];
1116 uByte buf[4+2+DECPMAX*2];
1117 uByte *umsd, *ulsd; /* local MSD and LSD pointers */
1118
1119 #if DECLITEND
1120 #define CARRYPAT 0x01000000 /* carry=1 pattern */
1121 #else
1122 #define CARRYPAT 0x00000001 /* carry=1 pattern */
1123 #endif
1124
1125 /* Start decoding the arguments */
1126 /* the initial exponents are placed into the opposite Ints to */
1127 /* that which might be expected; there are two sets of data to */
1128 /* keep track of (each decFloat and the corresponding exponent), */
1129 /* and this scheme means that at the swap point (after comparing */
1130 /* exponents) only one pair of words needs to be swapped */
1131 /* whichever path is taken (thereby minimising worst-case path) */
1132 sourhil=DFWORD(dfl, 0); /* LHS top word */
1133 expr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
1134 sourhir=DFWORD(dfr, 0); /* RHS top word */
1135 expl=DECCOMBEXP[sourhir>>26];
1136
1137 diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
1138
1139 if (EXPISSPECIAL(expl | expr)) { /* either is special? */
1140 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1141 /* one or two infinities */
1142 /* two infinities with different signs is invalid */
1143 if (diffsign && DFISINF(dfl) && DFISINF(dfr))
1144 return decInvalid(result, set);
1145 if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
1146 return decInfinity(result, dfr); /* RHS must be Infinite */
1147 }
1148
1149 /* Here when both arguments are finite */
1150
1151 /* complete exponent gathering (keeping swapped) */
1152 expr+=GETECON(dfl)-DECBIAS; /* .. + continuation and unbias */
1153 expl+=GETECON(dfr)-DECBIAS;
1154 /* here expr has exponent from lhs, and vice versa */
1155
1156 /* now swap either exponents or argument pointers */
1157 if (expl<=expr) {
1158 /* original left is bigger */
1159 Int expswap=expl;
1160 expl=expr;
1161 expr=expswap;
1162 /* printf("left bigger\n"); */
1163 }
1164 else {
1165 const decFloat *dfswap=dfl;
1166 dfl=dfr;
1167 dfr=dfswap;
1168 /* printf("right bigger\n"); */
1169 }
1170 /* [here dfl and expl refer to the datum with the larger exponent, */
1171 /* of if the exponents are equal then the original LHS argument] */
1172
1173 /* if lhs is zero then result will be the rhs (now known to have */
1174 /* the smaller exponent), which also may need to be tested for zero */
1175 /* for the weird IEEE 754 sign rules */
1176 if (DFISZERO(dfl)) {
1177 decCanonical(result, dfr); /* clean copy */
1178 /* "When the sum of two operands with opposite signs is */
1179 /* exactly zero, the sign of that sum shall be '+' in all */
1180 /* rounding modes except round toward -Infinity, in which */
1181 /* mode that sign shall be '-'." */
1182 if (diffsign && DFISZERO(result)) {
1183 DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */
1184 if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
1185 }
1186 return result;
1187 } /* numfl is zero */
1188 /* [here, LHS is non-zero; code below assumes that] */
1189
1190 /* Coefficients layout during the calculations to follow: */
1191 /* */
1192 /* Overlap case: */
1193 /* +------------------------------------------------+ */
1194 /* acc: |0000| coeffa | tail B | | */
1195 /* +------------------------------------------------+ */
1196 /* buf: |0000| pad0s | coeffb | | */
1197 /* +------------------------------------------------+ */
1198 /* */
1199 /* Touching coefficients or gap: */
1200 /* +------------------------------------------------+ */
1201 /* acc: |0000| coeffa | gap | coeffb | */
1202 /* +------------------------------------------------+ */
1203 /* [buf not used or needed; gap clamped to Pmax] */
1204
1205 /* lay out lhs coefficient into accumulator; this starts at acc+4 */
1206 /* for decDouble or acc+6 for decQuad so the LSD is word- */
1207 /* aligned; the top word gap is there only in case a carry digit */
1208 /* is prefixed after the add -- it does not need to be zeroed */
1209 #if DOUBLE
1210 #define COFF 4 /* offset into acc */
1211 #elif QUAD
1212 USHORTAT(acc+4)=0; /* prefix 00 */
1213 #define COFF 6 /* offset into acc */
1214 #endif
1215
1216 GETCOEFF(dfl, acc+COFF); /* decode from decFloat */
1217 ulsd=acc+COFF+DECPMAX-1;
1218 umsd=acc+4; /* [having this here avoids */
1219 /* weird GCC optimizer failure] */
1220 #if DECTRACE
1221 {bcdnum tum;
1222 tum.msd=umsd;
1223 tum.lsd=ulsd;
1224 tum.exponent=expl;
1225 tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1226 decShowNum(&tum, "dflx");}
1227 #endif
1228
1229 /* if signs differ, take ten's complement of lhs (here the */
1230 /* coefficient is subtracted from all-nines; the 1 is added during */
1231 /* the later add cycle -- zeros to the right do not matter because */
1232 /* the complement of zero is zero); these are fixed-length inverts */
1233 /* where the lsd is known to be at a 4-byte boundary (so no borrow */
1234 /* possible) */
1235 carry=0; /* assume no carry */
1236 if (diffsign) {
1237 carry=CARRYPAT; /* for +1 during add */
1238 UINTAT(acc+ 4)=0x09090909-UINTAT(acc+ 4);
1239 UINTAT(acc+ 8)=0x09090909-UINTAT(acc+ 8);
1240 UINTAT(acc+12)=0x09090909-UINTAT(acc+12);
1241 UINTAT(acc+16)=0x09090909-UINTAT(acc+16);
1242 #if QUAD
1243 UINTAT(acc+20)=0x09090909-UINTAT(acc+20);
1244 UINTAT(acc+24)=0x09090909-UINTAT(acc+24);
1245 UINTAT(acc+28)=0x09090909-UINTAT(acc+28);
1246 UINTAT(acc+32)=0x09090909-UINTAT(acc+32);
1247 UINTAT(acc+36)=0x09090909-UINTAT(acc+36);
1248 #endif
1249 } /* diffsign */
1250
1251 /* now process the rhs coefficient; if it cannot overlap lhs then */
1252 /* it can be put straight into acc (with an appropriate gap, if */
1253 /* needed) because no actual addition will be needed (except */
1254 /* possibly to complete ten's complement) */
1255 overlap=DECPMAX-(expl-expr);
1256 #if DECTRACE
1257 printf("exps: %ld %ld\n", (LI)expl, (LI)expr);
1258 printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
1259 #endif
1260
1261 if (overlap<=0) { /* no overlap possible */
1262 uInt gap; /* local work */
1263 /* since a full addition is not needed, a ten's complement */
1264 /* calculation started above may need to be completed */
1265 if (carry) {
1266 for (ub=ulsd; *ub==9; ub--) *ub=0;
1267 *ub+=1;
1268 carry=0; /* taken care of */
1269 }
1270 /* up to DECPMAX-1 digits of the final result can extend down */
1271 /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
1272 /* rhs will be simply sticky bits. In this case the gap is */
1273 /* clamped to DECPMAX and the exponent adjusted to suit [this is */
1274 /* safe because the lhs is non-zero]. */
1275 gap=-overlap;
1276 if (gap>DECPMAX) {
1277 expr+=gap-1;
1278 gap=DECPMAX;
1279 }
1280 ub=ulsd+gap+1; /* where MSD will go */
1281 /* Fill the gap with 0s; note that there is no addition to do */
1282 ui=&UINTAT(acc+COFF+DECPMAX); /* start of gap */
1283 for (; ui<&UINTAT(ub); ui++) *ui=0; /* mind the gap */
1284 if (overlap<-DECPMAX) { /* gap was > DECPMAX */
1285 *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
1286 }
1287 else { /* need full coefficient */
1288 GETCOEFF(dfr, ub); /* decode from decFloat */
1289 ub+=DECPMAX-1; /* new LSD... */
1290 }
1291 ulsd=ub; /* save new LSD */
1292 } /* no overlap possible */
1293
1294 else { /* overlap>0 */
1295 /* coefficients overlap (perhaps completely, although also */
1296 /* perhaps only where zeros) */
1297 ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */
1298 /* Fill the prefix gap with 0s; 8 will cover most common */
1299 /* unalignments, so start with direct assignments (a loop is */
1300 /* then used for any remaining -- the loop (and the one in a */
1301 /* moment) is not then on the critical path because the number */
1302 /* of additions is reduced by (at least) two in this case) */
1303 UINTAT(buf+4)=0; /* [clears decQuad 00 too] */
1304 UINTAT(buf+8)=0;
1305 if (ub>buf+12) {
1306 ui=&UINTAT(buf+12); /* start of any remaining */
1307 for (; ui<&UINTAT(ub); ui++) *ui=0; /* fill them */
1308 }
1309 GETCOEFF(dfr, ub); /* decode from decFloat */
1310
1311 /* now move tail of rhs across to main acc; again use direct */
1312 /* assignment for 8 digits-worth */
1313 UINTAT(acc+COFF+DECPMAX)=UINTAT(buf+COFF+DECPMAX);
1314 UINTAT(acc+COFF+DECPMAX+4)=UINTAT(buf+COFF+DECPMAX+4);
1315 if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
1316 uj=&UINTAT(buf+COFF+DECPMAX+8); /* source */
1317 ui=&UINTAT(acc+COFF+DECPMAX+8); /* target */
1318 for (; uj<&UINTAT(ub+DECPMAX); ui++, uj++) *ui=*uj;
1319 }
1320
1321 ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
1322
1323 /* now do the add of the non-tail; this is all nicely aligned, */
1324 /* and is over a multiple of four digits (because for Quad two */
1325 /* two 0 digits were added on the left); words in both acc and */
1326 /* buf (buf especially) will often be zero */
1327 /* [byte-by-byte add, here, is about 15% slower than the by-fours] */
1328
1329 /* Now effect the add; this is harder on a little-endian */
1330 /* machine as the inter-digit carry cannot use the usual BCD */
1331 /* addition trick because the bytes are loaded in the wrong order */
1332 /* [this loop could be unrolled, but probably scarcely worth it] */
1333
1334 ui=&UINTAT(acc+COFF+DECPMAX-4); /* target LSW (acc) */
1335 uj=&UINTAT(buf+COFF+DECPMAX-4); /* source LSW (buf, to add to acc) */
1336
1337 #if !DECLITEND
1338 for (; ui>=&UINTAT(acc+4); ui--, uj--) {
1339 /* bcd8 add */
1340 carry+=*uj; /* rhs + carry */
1341 if (carry==0) continue; /* no-op */
1342 carry+=*ui; /* lhs */
1343 /* Big-endian BCD adjust (uses internal carry) */
1344 carry+=0x76f6f6f6; /* note top nibble not all bits */
1345 *ui=(carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4); /* BCD adjust */
1346 carry>>=31; /* true carry was at far left */
1347 } /* add loop */
1348 #else
1349 for (; ui>=&UINTAT(acc+4); ui--, uj--) {
1350 /* bcd8 add */
1351 carry+=*uj; /* rhs + carry */
1352 if (carry==0) continue; /* no-op [common if unaligned] */
1353 carry+=*ui; /* lhs */
1354 /* Little-endian BCD adjust; inter-digit carry must be manual */
1355 /* because the lsb from the array will be in the most-significant */
1356 /* byte of carry */
1357 carry+=0x76767676; /* note no inter-byte carries */
1358 carry+=(carry & 0x80000000)>>15;
1359 carry+=(carry & 0x00800000)>>15;
1360 carry+=(carry & 0x00008000)>>15;
1361 carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
1362 *ui=carry & 0x0f0f0f0f; /* clear debris and save */
1363 /* here, final carry-out bit is at 0x00000080; move it ready */
1364 /* for next word-add (i.e., to 0x01000000) */
1365 carry=(carry & 0x00000080)<<17;
1366 } /* add loop */
1367 #endif
1368 #if DECTRACE
1369 {bcdnum tum;
1370 printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
1371 tum.msd=umsd; /* acc+4; */
1372 tum.lsd=ulsd;
1373 tum.exponent=0;
1374 tum.sign=0;
1375 decShowNum(&tum, "dfadd");}
1376 #endif
1377 } /* overlap possible */
1378
1379 /* ordering here is a little strange in order to have slowest path */
1380 /* first in GCC asm listing */
1381 if (diffsign) { /* subtraction */
1382 if (!carry) { /* no carry out means RHS<LHS */
1383 /* borrowed -- take ten's complement */
1384 /* sign is lhs sign */
1385 num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1386
1387 /* invert the coefficient first by fours, then add one; space */
1388 /* at the end of the buffer ensures the by-fours is always */
1389 /* safe, but lsd+1 must be cleared to prevent a borrow */
1390 /* if big-endian */
1391 #if !DECLITEND
1392 *(ulsd+1)=0;
1393 #endif
1394 /* there are always at least four coefficient words */
1395 UINTAT(umsd) =0x09090909-UINTAT(umsd);
1396 UINTAT(umsd+4) =0x09090909-UINTAT(umsd+4);
1397 UINTAT(umsd+8) =0x09090909-UINTAT(umsd+8);
1398 UINTAT(umsd+12)=0x09090909-UINTAT(umsd+12);
1399 #if DOUBLE
1400 #define BNEXT 16
1401 #elif QUAD
1402 UINTAT(umsd+16)=0x09090909-UINTAT(umsd+16);
1403 UINTAT(umsd+20)=0x09090909-UINTAT(umsd+20);
1404 UINTAT(umsd+24)=0x09090909-UINTAT(umsd+24);
1405 UINTAT(umsd+28)=0x09090909-UINTAT(umsd+28);
1406 UINTAT(umsd+32)=0x09090909-UINTAT(umsd+32);
1407 #define BNEXT 36
1408 #endif
1409 if (ulsd>=umsd+BNEXT) { /* unaligned */
1410 /* eight will handle most unaligments for Double; 16 for Quad */
1411 UINTAT(umsd+BNEXT)=0x09090909-UINTAT(umsd+BNEXT);
1412 UINTAT(umsd+BNEXT+4)=0x09090909-UINTAT(umsd+BNEXT+4);
1413 #if DOUBLE
1414 #define BNEXTY (BNEXT+8)
1415 #elif QUAD
1416 UINTAT(umsd+BNEXT+8)=0x09090909-UINTAT(umsd+BNEXT+8);
1417 UINTAT(umsd+BNEXT+12)=0x09090909-UINTAT(umsd+BNEXT+12);
1418 #define BNEXTY (BNEXT+16)
1419 #endif
1420 if (ulsd>=umsd+BNEXTY) { /* very unaligned */
1421 ui=&UINTAT(umsd+BNEXTY); /* -> continue */
1422 for (;;ui++) {
1423 *ui=0x09090909-*ui; /* invert four digits */
1424 if (ui>=&UINTAT(ulsd-3)) break; /* all done */
1425 }
1426 }
1427 }
1428 /* complete the ten's complement by adding 1 */
1429 for (ub=ulsd; *ub==9; ub--) *ub=0;
1430 *ub+=1;
1431 } /* borrowed */
1432
1433 else { /* carry out means RHS>=LHS */
1434 num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
1435 /* all done except for the special IEEE 754 exact-zero-result */
1436 /* rule (see above); while testing for zero, strip leading */
1437 /* zeros (which will save decFinalize doing it) (this is in */
1438 /* diffsign path, so carry impossible and true umsd is */
1439 /* acc+COFF) */
1440
1441 /* Check the initial coefficient area using the fast macro; */
1442 /* this will often be all that needs to be done (as on the */
1443 /* worst-case path when the subtraction was aligned and */
1444 /* full-length) */
1445 if (ISCOEFFZERO(acc+COFF)) {
1446 umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
1447 if (ulsd>umsd) { /* more to check */
1448 umsd++; /* to align after checked area */
1449 for (; UINTAT(umsd)==0 && umsd+3<ulsd;) umsd+=4;
1450 for (; *umsd==0 && umsd<ulsd;) umsd++;
1451 }
1452 if (*umsd==0) { /* must be true zero (and diffsign) */
1453 num.sign=0; /* assume + */
1454 if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
1455 }
1456 }
1457 /* [else was not zero, might still have leading zeros] */
1458 } /* subtraction gave positive result */
1459 } /* diffsign */
1460
1461 else { /* same-sign addition */
1462 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
1463 #if DOUBLE
1464 if (carry) { /* only possible with decDouble */
1465 *(acc+3)=1; /* [Quad has leading 00] */
1466 umsd=acc+3;
1467 }
1468 #endif
1469 } /* same sign */
1470
1471 num.msd=umsd; /* set MSD .. */
1472 num.lsd=ulsd; /* .. and LSD */
1473 num.exponent=expr; /* set exponent to smaller */
1474
1475 #if DECTRACE
1476 decFloatShow(dfl, "dfl");
1477 decFloatShow(dfr, "dfr");
1478 decShowNum(&num, "postadd");
1479 #endif
1480 return decFinalize(result, &num, set); /* round, check, and lay out */
1481 } /* decFloatAdd */
1482
1483 /* ------------------------------------------------------------------ */
1484 /* decFloatAnd -- logical digitwise AND of two decFloats */
1485 /* */
1486 /* result gets the result of ANDing dfl and dfr */
1487 /* dfl is the first decFloat (lhs) */
1488 /* dfr is the second decFloat (rhs) */
1489 /* set is the context */
1490 /* returns result, which will be canonical with sign=0 */
1491 /* */
1492 /* The operands must be positive, finite with exponent q=0, and */
1493 /* comprise just zeros and ones; if not, Invalid operation results. */
1494 /* ------------------------------------------------------------------ */
1495 decFloat * decFloatAnd(decFloat *result,
1496 const decFloat *dfl, const decFloat *dfr,
1497 decContext *set) {
1498 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
1499 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
1500 /* the operands are positive finite integers (q=0) with just 0s and 1s */
1501 #if DOUBLE
1502 DFWORD(result, 0)=ZEROWORD
1503 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
1504 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
1505 #elif QUAD
1506 DFWORD(result, 0)=ZEROWORD
1507 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
1508 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
1509 DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
1510 DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
1511 #endif
1512 return result;
1513 } /* decFloatAnd */
1514
1515 /* ------------------------------------------------------------------ */
1516 /* decFloatCanonical -- copy a decFloat, making canonical */
1517 /* */
1518 /* result gets the canonicalized df */
1519 /* df is the decFloat to copy and make canonical */
1520 /* returns result */
1521 /* */
1522 /* This works on specials, too; no error or exception is possible. */
1523 /* ------------------------------------------------------------------ */
1524 decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
1525 return decCanonical(result, df);
1526 } /* decFloatCanonical */
1527
1528 /* ------------------------------------------------------------------ */
1529 /* decFloatClass -- return the class of a decFloat */
1530 /* */
1531 /* df is the decFloat to test */
1532 /* returns the decClass that df falls into */
1533 /* ------------------------------------------------------------------ */
1534 enum decClass decFloatClass(const decFloat *df) {
1535 Int exp; /* exponent */
1536 if (DFISSPECIAL(df)) {
1537 if (DFISQNAN(df)) return DEC_CLASS_QNAN;
1538 if (DFISSNAN(df)) return DEC_CLASS_SNAN;
1539 /* must be an infinity */
1540 if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
1541 return DEC_CLASS_POS_INF;
1542 }
1543 if (DFISZERO(df)) { /* quite common */
1544 if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
1545 return DEC_CLASS_POS_ZERO;
1546 }
1547 /* is finite and non-zero; similar code to decFloatIsNormal, here */
1548 /* [this could be speeded up slightly by in-lining decFloatDigits] */
1549 exp=GETEXPUN(df) /* get unbiased exponent .. */
1550 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
1551 if (exp>=DECEMIN) { /* is normal */
1552 if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
1553 return DEC_CLASS_POS_NORMAL;
1554 }
1555 /* is subnormal */
1556 if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
1557 return DEC_CLASS_POS_SUBNORMAL;
1558 } /* decFloatClass */
1559
1560 /* ------------------------------------------------------------------ */
1561 /* decFloatClassString -- return the class of a decFloat as a string */
1562 /* */
1563 /* df is the decFloat to test */
1564 /* returns a constant string describing the class df falls into */
1565 /* ------------------------------------------------------------------ */
1566 const char *decFloatClassString(const decFloat *df) {
1567 enum decClass eclass=decFloatClass(df);
1568 if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
1569 if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
1570 if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
1571 if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
1572 if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
1573 if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
1574 if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
1575 if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
1576 if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
1577 if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
1578 return DEC_ClassString_UN; /* Unknown */
1579 } /* decFloatClassString */
1580
1581 /* ------------------------------------------------------------------ */
1582 /* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
1583 /* */
1584 /* result gets the result of comparing dfl and dfr */
1585 /* dfl is the first decFloat (lhs) */
1586 /* dfr is the second decFloat (rhs) */
1587 /* set is the context */
1588 /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
1589 /* ------------------------------------------------------------------ */
1590 decFloat * decFloatCompare(decFloat *result,
1591 const decFloat *dfl, const decFloat *dfr,
1592 decContext *set) {
1593 Int comp; /* work */
1594 /* NaNs are handled as usual */
1595 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1596 /* numeric comparison needed */
1597 comp=decNumCompare(dfl, dfr, 0);
1598 decFloatZero(result);
1599 if (comp==0) return result;
1600 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1601 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1602 return result;
1603 } /* decFloatCompare */
1604
1605 /* ------------------------------------------------------------------ */
1606 /* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
1607 /* */
1608 /* result gets the result of comparing dfl and dfr */
1609 /* dfl is the first decFloat (lhs) */
1610 /* dfr is the second decFloat (rhs) */
1611 /* set is the context */
1612 /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
1613 /* ------------------------------------------------------------------ */
1614 decFloat * decFloatCompareSignal(decFloat *result,
1615 const decFloat *dfl, const decFloat *dfr,
1616 decContext *set) {
1617 Int comp; /* work */
1618 /* NaNs are handled as usual, except that all NaNs signal */
1619 if (DFISNAN(dfl) || DFISNAN(dfr)) {
1620 set->status|=DEC_Invalid_operation;
1621 return decNaNs(result, dfl, dfr, set);
1622 }
1623 /* numeric comparison needed */
1624 comp=decNumCompare(dfl, dfr, 0);
1625 decFloatZero(result);
1626 if (comp==0) return result;
1627 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1628 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1629 return result;
1630 } /* decFloatCompareSignal */
1631
1632 /* ------------------------------------------------------------------ */
1633 /* decFloatCompareTotal -- compare two decFloats with total ordering */
1634 /* */
1635 /* result gets the result of comparing dfl and dfr */
1636 /* dfl is the first decFloat (lhs) */
1637 /* dfr is the second decFloat (rhs) */
1638 /* returns result, which may be -1, 0, or 1 */
1639 /* ------------------------------------------------------------------ */
1640 decFloat * decFloatCompareTotal(decFloat *result,
1641 const decFloat *dfl, const decFloat *dfr) {
1642 Int comp; /* work */
1643 if (DFISNAN(dfl) || DFISNAN(dfr)) {
1644 Int nanl, nanr; /* work */
1645 /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
1646 nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
1647 if (DFISSIGNED(dfl)) nanl=-nanl;
1648 nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
1649 if (DFISSIGNED(dfr)) nanr=-nanr;
1650 if (nanl>nanr) comp=+1;
1651 else if (nanl<nanr) comp=-1;
1652 else { /* NaNs are the same type and sign .. must compare payload */
1653 /* buffers need +2 for QUAD */
1654 uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */
1655 uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */
1656 uByte *ub, *uc; /* work */
1657 Int sigl; /* signum of LHS */
1658 sigl=(DFISSIGNED(dfl) ? -1 : +1);
1659
1660 /* decode the coefficients */
1661 /* (shift both right two if Quad to make a multiple of four) */
1662 #if QUAD
1663 USHORTAT(bufl)=0;
1664 USHORTAT(bufr)=0;
1665 #endif
1666 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
1667 GETCOEFF(dfr, bufr+QUAD*2); /* .. */
1668 /* all multiples of four, here */
1669 comp=0; /* assume equal */
1670 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
1671 if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
1672 /* about to find a winner; go by bytes in case little-endian */
1673 for (;; ub++, uc++) {
1674 if (*ub==*uc) continue;
1675 if (*ub>*uc) comp=sigl; /* difference found */
1676 else comp=-sigl; /* .. */
1677 break;
1678 }
1679 }
1680 } /* same NaN type and sign */
1681 }
1682 else {
1683 /* numeric comparison needed */
1684 comp=decNumCompare(dfl, dfr, 1); /* total ordering */
1685 }
1686 decFloatZero(result);
1687 if (comp==0) return result;
1688 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1689 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1690 return result;
1691 } /* decFloatCompareTotal */
1692
1693 /* ------------------------------------------------------------------ */
1694 /* decFloatCompareTotalMag -- compare magnitudes with total ordering */
1695 /* */
1696 /* result gets the result of comparing abs(dfl) and abs(dfr) */
1697 /* dfl is the first decFloat (lhs) */
1698 /* dfr is the second decFloat (rhs) */
1699 /* returns result, which may be -1, 0, or 1 */
1700 /* ------------------------------------------------------------------ */
1701 decFloat * decFloatCompareTotalMag(decFloat *result,
1702 const decFloat *dfl, const decFloat *dfr) {
1703 decFloat a, b; /* for copy if needed */
1704 /* copy and redirect signed operand(s) */
1705 if (DFISSIGNED(dfl)) {
1706 decFloatCopyAbs(&a, dfl);
1707 dfl=&a;
1708 }
1709 if (DFISSIGNED(dfr)) {
1710 decFloatCopyAbs(&b, dfr);
1711 dfr=&b;
1712 }
1713 return decFloatCompareTotal(result, dfl, dfr);
1714 } /* decFloatCompareTotalMag */
1715
1716 /* ------------------------------------------------------------------ */
1717 /* decFloatCopy -- copy a decFloat as-is */
1718 /* */
1719 /* result gets the copy of dfl */
1720 /* dfl is the decFloat to copy */
1721 /* returns result */
1722 /* */
1723 /* This is a bitwise operation; no errors or exceptions are possible. */
1724 /* ------------------------------------------------------------------ */
1725 decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
1726 if (dfl!=result) *result=*dfl; /* copy needed */
1727 return result;
1728 } /* decFloatCopy */
1729
1730 /* ------------------------------------------------------------------ */
1731 /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
1732 /* */
1733 /* result gets the copy of dfl with sign bit 0 */
1734 /* dfl is the decFloat to copy */
1735 /* returns result */
1736 /* */
1737 /* This is a bitwise operation; no errors or exceptions are possible. */
1738 /* ------------------------------------------------------------------ */
1739 decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
1740 if (dfl!=result) *result=*dfl; /* copy needed */
1741 DFBYTE(result, 0)&=~0x80; /* zero sign bit */
1742 return result;
1743 } /* decFloatCopyAbs */
1744
1745 /* ------------------------------------------------------------------ */
1746 /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
1747 /* */
1748 /* result gets the copy of dfl with sign bit inverted */
1749 /* dfl is the decFloat to copy */
1750 /* returns result */
1751 /* */
1752 /* This is a bitwise operation; no errors or exceptions are possible. */
1753 /* ------------------------------------------------------------------ */
1754 decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
1755 if (dfl!=result) *result=*dfl; /* copy needed */
1756 DFBYTE(result, 0)^=0x80; /* invert sign bit */
1757 return result;
1758 } /* decFloatCopyNegate */
1759
1760 /* ------------------------------------------------------------------ */
1761 /* decFloatCopySign -- copy a decFloat with the sign of another */
1762 /* */
1763 /* result gets the result of copying dfl with the sign of dfr */
1764 /* dfl is the first decFloat (lhs) */
1765 /* dfr is the second decFloat (rhs) */
1766 /* returns result */
1767 /* */
1768 /* This is a bitwise operation; no errors or exceptions are possible. */
1769 /* ------------------------------------------------------------------ */
1770 decFloat * decFloatCopySign(decFloat *result,
1771 const decFloat *dfl, const decFloat *dfr) {
1772 uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
1773 if (dfl!=result) *result=*dfl; /* copy needed */
1774 DFBYTE(result, 0)&=~0x80; /* clear sign .. */
1775 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
1776 return result;
1777 } /* decFloatCopySign */
1778
1779 /* ------------------------------------------------------------------ */
1780 /* decFloatDigits -- return the number of digits in a decFloat */
1781 /* */
1782 /* df is the decFloat to investigate */
1783 /* returns the number of significant digits in the decFloat; a */
1784 /* zero coefficient returns 1 as does an infinity (a NaN returns */
1785 /* the number of digits in the payload) */
1786 /* ------------------------------------------------------------------ */
1787 /* private macro to extract a declet according to provided formula */
1788 /* (form), and if it is non-zero then return the calculated digits */
1789 /* depending on the declet number (n), where n=0 for the most */
1790 /* significant declet; uses uInt dpd for work */
1791 #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
1792 if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1793 /* next one is used when it is known that the declet must be */
1794 /* non-zero, or is the final zero declet */
1795 #define dpdlendun(n, form) {dpd=(form)&0x3ff; \
1796 if (dpd==0) return 1; \
1797 return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1798
1799 uInt decFloatDigits(const decFloat *df) {
1800 uInt dpd; /* work */
1801 uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
1802 #if QUAD
1803 uInt sourmh, sourml;
1804 #endif
1805 uInt sourlo;
1806
1807 if (DFISINF(df)) return 1;
1808 /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
1809 /* then the coefficient is full-length */
1810 if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
1811
1812 #if DOUBLE
1813 if (sourhi&0x0003ffff) { /* ends in first */
1814 dpdlenchk(0, sourhi>>8);
1815 sourlo=DFWORD(df, 1);
1816 dpdlendun(1, (sourhi<<2) | (sourlo>>30));
1817 } /* [cannot drop through] */
1818 sourlo=DFWORD(df, 1); /* sourhi not involved now */
1819 if (sourlo&0xfff00000) { /* in one of first two */
1820 dpdlenchk(1, sourlo>>30); /* very rare */
1821 dpdlendun(2, sourlo>>20);
1822 } /* [cannot drop through] */
1823 dpdlenchk(3, sourlo>>10);
1824 dpdlendun(4, sourlo);
1825 /* [cannot drop through] */
1826
1827 #elif QUAD
1828 if (sourhi&0x00003fff) { /* ends in first */
1829 dpdlenchk(0, sourhi>>4);
1830 sourmh=DFWORD(df, 1);
1831 dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
1832 } /* [cannot drop through] */
1833 sourmh=DFWORD(df, 1);
1834 if (sourmh) {
1835 dpdlenchk(1, sourmh>>26);
1836 dpdlenchk(2, sourmh>>16);
1837 dpdlenchk(3, sourmh>>6);
1838 sourml=DFWORD(df, 2);
1839 dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
1840 } /* [cannot drop through] */
1841 sourml=DFWORD(df, 2);
1842 if (sourml) {
1843 dpdlenchk(4, sourml>>28);
1844 dpdlenchk(5, sourml>>18);
1845 dpdlenchk(6, sourml>>8);
1846 sourlo=DFWORD(df, 3);
1847 dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
1848 } /* [cannot drop through] */
1849 sourlo=DFWORD(df, 3);
1850 if (sourlo&0xfff00000) { /* in one of first two */
1851 dpdlenchk(7, sourlo>>30); /* very rare */
1852 dpdlendun(8, sourlo>>20);
1853 } /* [cannot drop through] */
1854 dpdlenchk(9, sourlo>>10);
1855 dpdlendun(10, sourlo);
1856 /* [cannot drop through] */
1857 #endif
1858 } /* decFloatDigits */
1859
1860 /* ------------------------------------------------------------------ */
1861 /* decFloatDivide -- divide a decFloat by another */
1862 /* */
1863 /* result gets the result of dividing dfl by dfr: */
1864 /* dfl is the first decFloat (lhs) */
1865 /* dfr is the second decFloat (rhs) */
1866 /* set is the context */
1867 /* returns result */
1868 /* */
1869 /* ------------------------------------------------------------------ */
1870 /* This is just a wrapper. */
1871 decFloat * decFloatDivide(decFloat *result,
1872 const decFloat *dfl, const decFloat *dfr,
1873 decContext *set) {
1874 return decDivide(result, dfl, dfr, set, DIVIDE);
1875 } /* decFloatDivide */
1876
1877 /* ------------------------------------------------------------------ */
1878 /* decFloatDivideInteger -- integer divide a decFloat by another */
1879 /* */
1880 /* result gets the result of dividing dfl by dfr: */
1881 /* dfl is the first decFloat (lhs) */
1882 /* dfr is the second decFloat (rhs) */
1883 /* set is the context */
1884 /* returns result */
1885 /* */
1886 /* ------------------------------------------------------------------ */
1887 decFloat * decFloatDivideInteger(decFloat *result,
1888 const decFloat *dfl, const decFloat *dfr,
1889 decContext *set) {
1890 return decDivide(result, dfl, dfr, set, DIVIDEINT);
1891 } /* decFloatDivideInteger */
1892
1893 /* ------------------------------------------------------------------ */
1894 /* decFloatFMA -- multiply and add three decFloats, fused */
1895 /* */
1896 /* result gets the result of (dfl*dfr)+dff with a single rounding */
1897 /* dfl is the first decFloat (lhs) */
1898 /* dfr is the second decFloat (rhs) */
1899 /* dff is the final decFloat (fhs) */
1900 /* set is the context */
1901 /* returns result */
1902 /* */
1903 /* ------------------------------------------------------------------ */
1904 decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
1905 const decFloat *dfr, const decFloat *dff,
1906 decContext *set) {
1907 /* The accumulator has the bytes needed for FiniteMultiply, plus */
1908 /* one byte to the left in case of carry, plus DECPMAX+2 to the */
1909 /* right for the final addition (up to full fhs + round & sticky) */
1910 #define FMALEN (1+ (DECPMAX9*18) +DECPMAX+2)
1911 uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
1912 /* .. and for final result */
1913 bcdnum mul; /* for multiplication result */
1914 bcdnum fin; /* for final operand, expanded */
1915 uByte coe[DECPMAX]; /* dff coefficient in BCD */
1916 bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
1917 uInt diffsign; /* non-zero if signs differ */
1918 uInt hipad; /* pad digit for hi if needed */
1919 Int padding; /* excess exponent */
1920 uInt carry; /* +1 for ten's complement and during add */
1921 uByte *ub, *uh, *ul; /* work */
1922
1923 /* handle all the special values [any special operand leads to a */
1924 /* special result] */
1925 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
1926 decFloat proxy; /* multiplication result proxy */
1927 /* NaNs are handled as usual, giving priority to sNaNs */
1928 if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1929 if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
1930 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1931 if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
1932 /* One or more of the three is infinite */
1933 /* infinity times zero is bad */
1934 decFloatZero(&proxy);
1935 if (DFISINF(dfl)) {
1936 if (DFISZERO(dfr)) return decInvalid(result, set);
1937 decInfinity(&proxy, &proxy);
1938 }
1939 else if (DFISINF(dfr)) {
1940 if (DFISZERO(dfl)) return decInvalid(result, set);
1941 decInfinity(&proxy, &proxy);
1942 }
1943 /* compute sign of multiplication and place in proxy */
1944 DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
1945 if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
1946 /* dff is Infinite */
1947 if (!DFISINF(&proxy)) return decInfinity(result, dff);
1948 /* both sides of addition are infinite; different sign is bad */
1949 if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
1950 return decInvalid(result, set);
1951 return decFloatCopy(result, &proxy);
1952 }
1953
1954 /* Here when all operands are finite */
1955
1956 /* First multiply dfl*dfr */
1957 decFiniteMultiply(&mul, acc+1, dfl, dfr);
1958 /* The multiply is complete, exact and unbounded, and described in */
1959 /* mul with the coefficient held in acc[1...] */
1960
1961 /* now add in dff; the algorithm is essentially the same as */
1962 /* decFloatAdd, but the code is different because the code there */
1963 /* is highly optimized for adding two numbers of the same size */
1964 fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
1965 fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
1966 diffsign=mul.sign^fin.sign; /* note if signs differ */
1967 fin.msd=coe;
1968 fin.lsd=coe+DECPMAX-1;
1969 GETCOEFF(dff, coe); /* extract the coefficient */
1970
1971 /* now set hi and lo so that hi points to whichever of mul and fin */
1972 /* has the higher exponent and lo point to the other [don't care if */
1973 /* the same] */
1974 if (mul.exponent>=fin.exponent) {
1975 hi=&mul;
1976 lo=&fin;
1977 }
1978 else {
1979 hi=&fin;
1980 lo=&mul;
1981 }
1982
1983 /* remove leading zeros on both operands; this will save time later */
1984 /* and make testing for zero trivial */
1985 for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
1986 for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
1987 for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
1988 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
1989
1990 /* if hi is zero then result will be lo (which has the smaller */
1991 /* exponent), which also may need to be tested for zero for the */
1992 /* weird IEEE 754 sign rules */
1993 if (*hi->msd==0 && hi->msd==hi->lsd) { /* hi is zero */
1994 /* "When the sum of two operands with opposite signs is */
1995 /* exactly zero, the sign of that sum shall be '+' in all */
1996 /* rounding modes except round toward -Infinity, in which */
1997 /* mode that sign shall be '-'." */
1998 if (diffsign) {
1999 if (*lo->msd==0 && lo->msd==lo->lsd) { /* lo is zero */
2000 lo->sign=0;
2001 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2002 } /* diffsign && lo=0 */
2003 } /* diffsign */
2004 return decFinalize(result, lo, set); /* may need clamping */
2005 } /* numfl is zero */
2006 /* [here, both are minimal length and hi is non-zero] */
2007
2008 /* if signs differ, take the ten's complement of hi (zeros to the */
2009 /* right do not matter because the complement of zero is zero); */
2010 /* the +1 is done later, as part of the addition, inserted at the */
2011 /* correct digit */
2012 hipad=0;
2013 carry=0;
2014 if (diffsign) {
2015 hipad=9;
2016 carry=1;
2017 /* exactly the correct number of digits must be inverted */
2018 for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UINTAT(uh)=0x09090909-UINTAT(uh);
2019 for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
2020 }
2021
2022 /* ready to add; note that hi has no leading zeros so gap */
2023 /* calculation does not have to be as pessimistic as in decFloatAdd */
2024 /* (this is much more like the arbitrary-precision algorithm in */
2025 /* Rexx and decNumber) */
2026
2027 /* padding is the number of zeros that would need to be added to hi */
2028 /* for its lsd to be aligned with the lsd of lo */
2029 padding=hi->exponent-lo->exponent;
2030 /* printf("FMA pad %ld\n", (LI)padding); */
2031
2032 /* the result of the addition will be built into the accumulator, */
2033 /* starting from the far right; this could be either hi or lo */
2034 ub=acc+FMALEN-1; /* where lsd of result will go */
2035 ul=lo->lsd; /* lsd of rhs */
2036
2037 if (padding!=0) { /* unaligned */
2038 /* if the msd of lo is more than DECPMAX+2 digits to the right of */
2039 /* the original msd of hi then it can be reduced to a single */
2040 /* digit at the right place, as it stays clear of hi digits */
2041 /* [it must be DECPMAX+2 because during a subtraction the msd */
2042 /* could become 0 after a borrow from 1.000 to 0.9999...] */
2043 Int hilen=(Int)(hi->lsd-hi->msd+1); /* lengths */
2044 Int lolen=(Int)(lo->lsd-lo->msd+1); /* .. */
2045 Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
2046 Int reduce=newexp-lo->exponent;
2047 if (reduce>0) { /* [= case gives reduce=0 nop] */
2048 /* printf("FMA reduce: %ld\n", (LI)reduce); */
2049 if (reduce>=lolen) { /* eating all */
2050 lo->lsd=lo->msd; /* reduce to single digit */
2051 lo->exponent=newexp; /* [known to be non-zero] */
2052 }
2053 else { /* < */
2054 uByte *up=lo->lsd;
2055 lo->lsd=lo->lsd-reduce;
2056 if (*lo->lsd==0) /* could need sticky bit */
2057 for (; up>lo->lsd; up--) { /* search discarded digits */
2058 if (*up!=0) { /* found one... */
2059 *lo->lsd=1; /* set sticky bit */
2060 break;
2061 }
2062 }
2063 lo->exponent+=reduce;
2064 }
2065 padding=hi->exponent-lo->exponent; /* recalculate */
2066 ul=lo->lsd; /* .. */
2067 } /* maybe reduce */
2068 /* padding is now <= DECPMAX+2 but still > 0; tricky DOUBLE case */
2069 /* is when hi is a 1 that will become a 0.9999... by subtraction: */
2070 /* hi: 1 E+16 */
2071 /* lo: .................1000000000000000 E-16 */
2072 /* which for the addition pads and reduces to: */
2073 /* hi: 1000000000000000000 E-2 */
2074 /* lo: .................1 E-2 */
2075 #if DECCHECK
2076 if (padding>DECPMAX+2) printf("FMA excess padding: %ld\n", (LI)padding);
2077 if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
2078 /* printf("FMA padding: %ld\n", (LI)padding); */
2079 #endif
2080 /* padding digits can now be set in the result; one or more of */
2081 /* these will come from lo; others will be zeros in the gap */
2082 for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
2083 for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
2084 }
2085
2086 /* addition now complete to the right of the rightmost digit of hi */
2087 uh=hi->lsd;
2088
2089 /* carry was set up depending on ten's complement above; do the add... */
2090 for (;; ub--) {
2091 uInt hid, lod;
2092 if (uh<hi->msd) {
2093 if (ul<lo->msd) break;
2094 hid=hipad;
2095 }
2096 else hid=*uh--;
2097 if (ul<lo->msd) lod=0;
2098 else lod=*ul--;
2099 *ub=(uByte)(carry+hid+lod);
2100 if (*ub<10) carry=0;
2101 else {
2102 *ub-=10;
2103 carry=1;
2104 }
2105 } /* addition loop */
2106
2107 /* addition complete -- now handle carry, borrow, etc. */
2108 /* use lo to set up the num (its exponent is already correct, and */
2109 /* sign usually is) */
2110 lo->msd=ub+1;
2111 lo->lsd=acc+FMALEN-1;
2112 /* decShowNum(lo, "lo"); */
2113 if (!diffsign) { /* same-sign addition */
2114 if (carry) { /* carry out */
2115 *ub=1; /* place the 1 .. */
2116 lo->msd--; /* .. and update */
2117 }
2118 } /* same sign */
2119 else { /* signs differed (subtraction) */
2120 if (!carry) { /* no carry out means hi<lo */
2121 /* borrowed -- take ten's complement of the right digits */
2122 lo->sign=hi->sign; /* sign is lhs sign */
2123 for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UINTAT(ul)=0x09090909-UINTAT(ul);
2124 for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
2125 /* complete the ten's complement by adding 1 [cannot overrun] */
2126 for (ul--; *ul==9; ul--) *ul=0;
2127 *ul+=1;
2128 } /* borrowed */
2129 else { /* carry out means hi>=lo */
2130 /* sign to use is lo->sign */
2131 /* all done except for the special IEEE 754 exact-zero-result */
2132 /* rule (see above); while testing for zero, strip leading */
2133 /* zeros (which will save decFinalize doing it) */
2134 for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2135 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2136 if (*lo->msd==0) { /* must be true zero (and diffsign) */
2137 lo->sign=0; /* assume + */
2138 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2139 }
2140 /* [else was not zero, might still have leading zeros] */
2141 } /* subtraction gave positive result */
2142 } /* diffsign */
2143
2144 return decFinalize(result, lo, set); /* round, check, and lay out */
2145 } /* decFloatFMA */
2146
2147 /* ------------------------------------------------------------------ */
2148 /* decFloatFromInt -- initialise a decFloat from an Int */
2149 /* */
2150 /* result gets the converted Int */
2151 /* n is the Int to convert */
2152 /* returns result */
2153 /* */
2154 /* The result is Exact; no errors or exceptions are possible. */
2155 /* ------------------------------------------------------------------ */
2156 decFloat * decFloatFromInt32(decFloat *result, Int n) {
2157 uInt u=(uInt)n; /* copy as bits */
2158 uInt encode; /* work */
2159 DFWORD(result, 0)=ZEROWORD; /* always */
2160 #if QUAD
2161 DFWORD(result, 1)=0;
2162 DFWORD(result, 2)=0;
2163 #endif
2164 if (n<0) { /* handle -n with care */
2165 /* [This can be done without the test, but is then slightly slower] */
2166 u=(~u)+1;
2167 DFWORD(result, 0)|=DECFLOAT_Sign;
2168 }
2169 /* Since the maximum value of u now is 2**31, only the low word of */
2170 /* result is affected */
2171 encode=BIN2DPD[u%1000];
2172 u/=1000;
2173 encode|=BIN2DPD[u%1000]<<10;
2174 u/=1000;
2175 encode|=BIN2DPD[u%1000]<<20;
2176 u/=1000; /* now 0, 1, or 2 */
2177 encode|=u<<30;
2178 DFWORD(result, DECWORDS-1)=encode;
2179 return result;
2180 } /* decFloatFromInt32 */
2181
2182 /* ------------------------------------------------------------------ */
2183 /* decFloatFromUInt -- initialise a decFloat from a uInt */
2184 /* */
2185 /* result gets the converted uInt */
2186 /* n is the uInt to convert */
2187 /* returns result */
2188 /* */
2189 /* The result is Exact; no errors or exceptions are possible. */
2190 /* ------------------------------------------------------------------ */
2191 decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
2192 uInt encode; /* work */
2193 DFWORD(result, 0)=ZEROWORD; /* always */
2194 #if QUAD
2195 DFWORD(result, 1)=0;
2196 DFWORD(result, 2)=0;
2197 #endif
2198 encode=BIN2DPD[u%1000];
2199 u/=1000;
2200 encode|=BIN2DPD[u%1000]<<10;
2201 u/=1000;
2202 encode|=BIN2DPD[u%1000]<<20;
2203 u/=1000; /* now 0 -> 4 */
2204 encode|=u<<30;
2205 DFWORD(result, DECWORDS-1)=encode;
2206 DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
2207 return result;
2208 } /* decFloatFromUInt32 */
2209
2210 /* ------------------------------------------------------------------ */
2211 /* decFloatInvert -- logical digitwise INVERT of a decFloat */
2212 /* */
2213 /* result gets the result of INVERTing df */
2214 /* df is the decFloat to invert */
2215 /* set is the context */
2216 /* returns result, which will be canonical with sign=0 */
2217 /* */
2218 /* The operand must be positive, finite with exponent q=0, and */
2219 /* comprise just zeros and ones; if not, Invalid operation results. */
2220 /* ------------------------------------------------------------------ */
2221 decFloat * decFloatInvert(decFloat *result, const decFloat *df,
2222 decContext *set) {
2223 uInt sourhi=DFWORD(df, 0); /* top word of dfs */
2224
2225 if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
2226 /* the operand is a finite integer (q=0) */
2227 #if DOUBLE
2228 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
2229 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
2230 #elif QUAD
2231 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
2232 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
2233 DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
2234 DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
2235 #endif
2236 return result;
2237 } /* decFloatInvert */
2238
2239 /* ------------------------------------------------------------------ */
2240 /* decFloatIs -- decFloat tests (IsSigned, etc.) */
2241 /* */
2242 /* df is the decFloat to test */
2243 /* returns 0 or 1 in an int32_t */
2244 /* */
2245 /* Many of these could be macros, but having them as real functions */
2246 /* is a bit cleaner (and they can be referred to here by the generic */
2247 /* names) */
2248 /* ------------------------------------------------------------------ */
2249 uInt decFloatIsCanonical(const decFloat *df) {
2250 if (DFISSPECIAL(df)) {
2251 if (DFISINF(df)) {
2252 if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
2253 if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
2254 return 1;
2255 }
2256 /* is a NaN */
2257 if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
2258 if (DFISCCZERO(df)) return 1; /* coefficient continuation */
2259 /* drop through to check payload */
2260 }
2261 { /* declare block */
2262 #if DOUBLE
2263 uInt sourhi=DFWORD(df, 0);
2264 uInt sourlo=DFWORD(df, 1);
2265 if (CANONDPDOFF(sourhi, 8)
2266 && CANONDPDTWO(sourhi, sourlo, 30)
2267 && CANONDPDOFF(sourlo, 20)
2268 && CANONDPDOFF(sourlo, 10)
2269 && CANONDPDOFF(sourlo, 0)) return 1;
2270 #elif QUAD
2271 uInt sourhi=DFWORD(df, 0);
2272 uInt sourmh=DFWORD(df, 1);
2273 uInt sourml=DFWORD(df, 2);
2274 uInt sourlo=DFWORD(df, 3);
2275 if (CANONDPDOFF(sourhi, 4)
2276 && CANONDPDTWO(sourhi, sourmh, 26)
2277 && CANONDPDOFF(sourmh, 16)
2278 && CANONDPDOFF(sourmh, 6)
2279 && CANONDPDTWO(sourmh, sourml, 28)
2280 && CANONDPDOFF(sourml, 18)
2281 && CANONDPDOFF(sourml, 8)
2282 && CANONDPDTWO(sourml, sourlo, 30)
2283 && CANONDPDOFF(sourlo, 20)
2284 && CANONDPDOFF(sourlo, 10)
2285 && CANONDPDOFF(sourlo, 0)) return 1;
2286 #endif
2287 } /* block */
2288 return 0; /* a declet is non-canonical */
2289 }
2290
2291 uInt decFloatIsFinite(const decFloat *df) {
2292 return !DFISSPECIAL(df);
2293 }
2294 uInt decFloatIsInfinite(const decFloat *df) {
2295 return DFISINF(df);
2296 }
2297 uInt decFloatIsInteger(const decFloat *df) {
2298 return DFISINT(df);
2299 }
2300 uInt decFloatIsNaN(const decFloat *df) {
2301 return DFISNAN(df);
2302 }
2303 uInt decFloatIsNormal(const decFloat *df) {
2304 Int exp; /* exponent */
2305 if (DFISSPECIAL(df)) return 0;
2306 if (DFISZERO(df)) return 0;
2307 /* is finite and non-zero */
2308 exp=GETEXPUN(df) /* get unbiased exponent .. */
2309 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2310 return (exp>=DECEMIN); /* < DECEMIN is subnormal */
2311 }
2312 uInt decFloatIsSignaling(const decFloat *df) {
2313 return DFISSNAN(df);
2314 }
2315 uInt decFloatIsSignalling(const decFloat *df) {
2316 return DFISSNAN(df);
2317 }
2318 uInt decFloatIsSigned(const decFloat *df) {
2319 return DFISSIGNED(df);
2320 }
2321 uInt decFloatIsSubnormal(const decFloat *df) {
2322 if (DFISSPECIAL(df)) return 0;
2323 /* is finite */
2324 if (decFloatIsNormal(df)) return 0;
2325 /* it is <Nmin, but could be zero */
2326 if (DFISZERO(df)) return 0;
2327 return 1; /* is subnormal */
2328 }
2329 uInt decFloatIsZero(const decFloat *df) {
2330 return DFISZERO(df);
2331 } /* decFloatIs... */
2332
2333 /* ------------------------------------------------------------------ */
2334 /* decFloatLogB -- return adjusted exponent, by 754r rules */
2335 /* */
2336 /* result gets the adjusted exponent as an integer, or a NaN etc. */
2337 /* df is the decFloat to be examined */
2338 /* set is the context */
2339 /* returns result */
2340 /* */
2341 /* Notable cases: */
2342 /* A<0 -> Use |A| */
2343 /* A=0 -> -Infinity (Division by zero) */
2344 /* A=Infinite -> +Infinity (Exact) */
2345 /* A=1 exactly -> 0 (Exact) */
2346 /* NaNs are propagated as usual */
2347 /* ------------------------------------------------------------------ */
2348 decFloat * decFloatLogB(decFloat *result, const decFloat *df,
2349 decContext *set) {
2350 Int ae; /* adjusted exponent */
2351 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2352 if (DFISINF(df)) {
2353 DFWORD(result, 0)=0; /* need +ve */
2354 return decInfinity(result, result); /* canonical +Infinity */
2355 }
2356 if (DFISZERO(df)) {
2357 set->status|=DEC_Division_by_zero; /* as per 754r */
2358 DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
2359 return decInfinity(result, result); /* canonical -Infinity */
2360 }
2361 ae=GETEXPUN(df) /* get unbiased exponent .. */
2362 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2363 /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
2364 /* it is worth using a special case of decFloatFromInt32 */
2365 DFWORD(result, 0)=ZEROWORD; /* always */
2366 if (ae<0) {
2367 DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
2368 ae=-ae;
2369 }
2370 #if DOUBLE
2371 DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
2372 #elif QUAD
2373 DFWORD(result, 1)=0;
2374 DFWORD(result, 2)=0;
2375 DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
2376 DFWORD(result, 3)|=BIN2DPD[ae%1000];
2377 #endif
2378 return result;
2379 } /* decFloatLogB */
2380
2381 /* ------------------------------------------------------------------ */
2382 /* decFloatMax -- return maxnum of two operands */
2383 /* */
2384 /* result gets the chosen decFloat */
2385 /* dfl is the first decFloat (lhs) */
2386 /* dfr is the second decFloat (rhs) */
2387 /* set is the context */
2388 /* returns result */
2389 /* */
2390 /* If just one operand is a quiet NaN it is ignored. */
2391 /* ------------------------------------------------------------------ */
2392 decFloat * decFloatMax(decFloat *result,
2393 const decFloat *dfl, const decFloat *dfr,
2394 decContext *set) {
2395 Int comp;
2396 if (DFISNAN(dfl)) {
2397 /* sNaN or both NaNs leads to normal NaN processing */
2398 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2399 return decCanonical(result, dfr); /* RHS is numeric */
2400 }
2401 if (DFISNAN(dfr)) {
2402 /* sNaN leads to normal NaN processing (both NaN handled above) */
2403 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2404 return decCanonical(result, dfl); /* LHS is numeric */
2405 }
2406 /* Both operands are numeric; numeric comparison needed -- use */
2407 /* total order for a well-defined choice (and +0 > -0) */
2408 comp=decNumCompare(dfl, dfr, 1);
2409 if (comp>=0) return decCanonical(result, dfl);
2410 return decCanonical(result, dfr);
2411 } /* decFloatMax */
2412
2413 /* ------------------------------------------------------------------ */
2414 /* decFloatMaxMag -- return maxnummag of two operands */
2415 /* */
2416 /* result gets the chosen decFloat */
2417 /* dfl is the first decFloat (lhs) */
2418 /* dfr is the second decFloat (rhs) */
2419 /* set is the context */
2420 /* returns result */
2421 /* */
2422 /* Returns according to the magnitude comparisons if both numeric and */
2423 /* unequal, otherwise returns maxnum */
2424 /* ------------------------------------------------------------------ */
2425 decFloat * decFloatMaxMag(decFloat *result,
2426 const decFloat *dfl, const decFloat *dfr,
2427 decContext *set) {
2428 Int comp;
2429 decFloat absl, absr;
2430 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
2431
2432 decFloatCopyAbs(&absl, dfl);
2433 decFloatCopyAbs(&absr, dfr);
2434 comp=decNumCompare(&absl, &absr, 0);
2435 if (comp>0) return decCanonical(result, dfl);
2436 if (comp<0) return decCanonical(result, dfr);
2437 return decFloatMax(result, dfl, dfr, set);
2438 } /* decFloatMaxMag */
2439
2440 /* ------------------------------------------------------------------ */
2441 /* decFloatMin -- return minnum of two operands */
2442 /* */
2443 /* result gets the chosen decFloat */
2444 /* dfl is the first decFloat (lhs) */
2445 /* dfr is the second decFloat (rhs) */
2446 /* set is the context */
2447 /* returns result */
2448 /* */
2449 /* If just one operand is a quiet NaN it is ignored. */
2450 /* ------------------------------------------------------------------ */
2451 decFloat * decFloatMin(decFloat *result,
2452 const decFloat *dfl, const decFloat *dfr,
2453 decContext *set) {
2454 Int comp;
2455 if (DFISNAN(dfl)) {
2456 /* sNaN or both NaNs leads to normal NaN processing */
2457 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2458 return decCanonical(result, dfr); /* RHS is numeric */
2459 }
2460 if (DFISNAN(dfr)) {
2461 /* sNaN leads to normal NaN processing (both NaN handled above) */
2462 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2463 return decCanonical(result, dfl); /* LHS is numeric */
2464 }
2465 /* Both operands are numeric; numeric comparison needed -- use */
2466 /* total order for a well-defined choice (and +0 > -0) */
2467 comp=decNumCompare(dfl, dfr, 1);
2468 if (comp<=0) return decCanonical(result, dfl);
2469 return decCanonical(result, dfr);
2470 } /* decFloatMin */
2471
2472 /* ------------------------------------------------------------------ */
2473 /* decFloatMinMag -- return minnummag of two operands */
2474 /* */
2475 /* result gets the chosen decFloat */
2476 /* dfl is the first decFloat (lhs) */
2477 /* dfr is the second decFloat (rhs) */
2478 /* set is the context */
2479 /* returns result */
2480 /* */
2481 /* Returns according to the magnitude comparisons if both numeric and */
2482 /* unequal, otherwise returns minnum */
2483 /* ------------------------------------------------------------------ */
2484 decFloat * decFloatMinMag(decFloat *result,
2485 const decFloat *dfl, const decFloat *dfr,
2486 decContext *set) {
2487 Int comp;
2488 decFloat absl, absr;
2489 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
2490
2491 decFloatCopyAbs(&absl, dfl);
2492 decFloatCopyAbs(&absr, dfr);
2493 comp=decNumCompare(&absl, &absr, 0);
2494 if (comp<0) return decCanonical(result, dfl);
2495 if (comp>0) return decCanonical(result, dfr);
2496 return decFloatMin(result, dfl, dfr, set);
2497 } /* decFloatMinMag */
2498
2499 /* ------------------------------------------------------------------ */
2500 /* decFloatMinus -- negate value, heeding NaNs, etc. */
2501 /* */
2502 /* result gets the canonicalized 0-df */
2503 /* df is the decFloat to minus */
2504 /* set is the context */
2505 /* returns result */
2506 /* */
2507 /* This has the same effect as 0-df where the exponent of the zero is */
2508 /* the same as that of df (if df is finite). */
2509 /* The effect is also the same as decFloatCopyNegate except that NaNs */
2510 /* are handled normally (the sign of a NaN is not affected, and an */
2511 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2512 /* ------------------------------------------------------------------ */
2513 decFloat * decFloatMinus(decFloat *result, const decFloat *df,
2514 decContext *set) {
2515 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2516 decCanonical(result, df); /* copy and check */
2517 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2518 else DFBYTE(result, 0)^=0x80; /* flip sign bit */
2519 return result;
2520 } /* decFloatMinus */
2521
2522 /* ------------------------------------------------------------------ */
2523 /* decFloatMultiply -- multiply two decFloats */
2524 /* */
2525 /* result gets the result of multiplying dfl and dfr: */
2526 /* dfl is the first decFloat (lhs) */
2527 /* dfr is the second decFloat (rhs) */
2528 /* set is the context */
2529 /* returns result */
2530 /* */
2531 /* ------------------------------------------------------------------ */
2532 decFloat * decFloatMultiply(decFloat *result,
2533 const decFloat *dfl, const decFloat *dfr,
2534 decContext *set) {
2535 bcdnum num; /* for final conversion */
2536 uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
2537
2538 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
2539 /* NaNs are handled as usual */
2540 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2541 /* infinity times zero is bad */
2542 if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
2543 if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
2544 /* both infinite; return canonical infinity with computed sign */
2545 DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
2546 return decInfinity(result, result);
2547 }
2548
2549 /* Here when both operands are finite */
2550 decFiniteMultiply(&num, bcdacc, dfl, dfr);
2551 return decFinalize(result, &num, set); /* round, check, and lay out */
2552 } /* decFloatMultiply */
2553
2554 /* ------------------------------------------------------------------ */
2555 /* decFloatNextMinus -- next towards -Infinity */
2556 /* */
2557 /* result gets the next lesser decFloat */
2558 /* dfl is the decFloat to start with */
2559 /* set is the context */
2560 /* returns result */
2561 /* */
2562 /* This is 754r nextdown; Invalid is the only status possible (from */
2563 /* an sNaN). */
2564 /* ------------------------------------------------------------------ */
2565 decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
2566 decContext *set) {
2567 decFloat delta; /* tiny increment */
2568 uInt savestat; /* saves status */
2569 enum rounding saveround; /* .. and mode */
2570
2571 /* +Infinity is the special case */
2572 if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
2573 DFSETNMAX(result);
2574 return result; /* [no status to set] */
2575 }
2576 /* other cases are effected by sutracting a tiny delta -- this */
2577 /* should be done in a wider format as the delta is unrepresentable */
2578 /* here (but can be done with normal add if the sign of zero is */
2579 /* treated carefully, because no Inexactitude is interesting); */
2580 /* rounding to -Infinity then pushes the result to next below */
2581 decFloatZero(&delta); /* set up tiny delta */
2582 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2583 DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
2584 /* set up for the directional round */
2585 saveround=set->round; /* save mode */
2586 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2587 savestat=set->status; /* save status */
2588 decFloatAdd(result, dfl, &delta, set);
2589 /* Add rules mess up the sign when going from +Ntiny to 0 */
2590 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2591 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2592 set->status|=savestat; /* restore pending flags */
2593 set->round=saveround; /* .. and mode */
2594 return result;
2595 } /* decFloatNextMinus */
2596
2597 /* ------------------------------------------------------------------ */
2598 /* decFloatNextPlus -- next towards +Infinity */
2599 /* */
2600 /* result gets the next larger decFloat */
2601 /* dfl is the decFloat to start with */
2602 /* set is the context */
2603 /* returns result */
2604 /* */
2605 /* This is 754r nextup; Invalid is the only status possible (from */
2606 /* an sNaN). */
2607 /* ------------------------------------------------------------------ */
2608 decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
2609 decContext *set) {
2610 uInt savestat; /* saves status */
2611 enum rounding saveround; /* .. and mode */
2612 decFloat delta; /* tiny increment */
2613
2614 /* -Infinity is the special case */
2615 if (DFISINF(dfl) && DFISSIGNED(dfl)) {
2616 DFSETNMAX(result);
2617 DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
2618 return result; /* [no status to set] */
2619 }
2620 /* other cases are effected by sutracting a tiny delta -- this */
2621 /* should be done in a wider format as the delta is unrepresentable */
2622 /* here (but can be done with normal add if the sign of zero is */
2623 /* treated carefully, because no Inexactitude is interesting); */
2624 /* rounding to +Infinity then pushes the result to next above */
2625 decFloatZero(&delta); /* set up tiny delta */
2626 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2627 DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
2628 /* set up for the directional round */
2629 saveround=set->round; /* save mode */
2630 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2631 savestat=set->status; /* save status */
2632 decFloatAdd(result, dfl, &delta, set);
2633 /* Add rules mess up the sign when going from -Ntiny to -0 */
2634 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2635 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2636 set->status|=savestat; /* restore pending flags */
2637 set->round=saveround; /* .. and mode */
2638 return result;
2639 } /* decFloatNextPlus */
2640
2641 /* ------------------------------------------------------------------ */
2642 /* decFloatNextToward -- next towards a decFloat */
2643 /* */
2644 /* result gets the next decFloat */
2645 /* dfl is the decFloat to start with */
2646 /* dfr is the decFloat to move toward */
2647 /* set is the context */
2648 /* returns result */
2649 /* */
2650 /* This is 754r nextafter; status may be set unless the result is a */
2651 /* normal number. */
2652 /* ------------------------------------------------------------------ */
2653 decFloat * decFloatNextToward(decFloat *result,
2654 const decFloat *dfl, const decFloat *dfr,
2655 decContext *set) {
2656 decFloat delta; /* tiny increment or decrement */
2657 decFloat pointone; /* 1e-1 */
2658 uInt savestat; /* saves status */
2659 enum rounding saveround; /* .. and mode */
2660 uInt deltatop; /* top word for delta */
2661 Int comp; /* work */
2662
2663 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2664 /* Both are numeric, so Invalid no longer a possibility */
2665 comp=decNumCompare(dfl, dfr, 0);
2666 if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
2667 /* unequal; do NextPlus or NextMinus but with different status rules */
2668
2669 if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
2670 if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
2671 DFSETNMAX(result);
2672 DFWORD(result, 0)|=DECFLOAT_Sign;
2673 return result;
2674 }
2675 saveround=set->round; /* save mode */
2676 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2677 deltatop=0; /* positive delta */
2678 }
2679 else { /* lhs>rhs, do NextMinus, see above for commentary */
2680 if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
2681 DFSETNMAX(result);
2682 return result;
2683 }
2684 saveround=set->round; /* save mode */
2685 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2686 deltatop=DECFLOAT_Sign; /* negative delta */
2687 }
2688 savestat=set->status; /* save status */
2689 /* Here, Inexact is needed where appropriate (and hence Underflow, */
2690 /* etc.). Therefore the tiny delta which is otherwise */
2691 /* unrepresentable (see NextPlus and NextMinus) is constructed */
2692 /* using the multiplication of FMA. */
2693 decFloatZero(&delta); /* set up tiny delta */
2694 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2695 DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
2696 decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
2697 decFloatFMA(result, &delta, &pointone, dfl, set);
2698 /* [Delta is truly tiny, so no need to correct sign of zero] */
2699 /* use new status unless the result is normal */
2700 if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
2701 set->round=saveround; /* restore mode */
2702 return result;
2703 } /* decFloatNextToward */
2704
2705 /* ------------------------------------------------------------------ */
2706 /* decFloatOr -- logical digitwise OR of two decFloats */
2707 /* */
2708 /* result gets the result of ORing dfl and dfr */
2709 /* dfl is the first decFloat (lhs) */
2710 /* dfr is the second decFloat (rhs) */
2711 /* set is the context */
2712 /* returns result, which will be canonical with sign=0 */
2713 /* */
2714 /* The operands must be positive, finite with exponent q=0, and */
2715 /* comprise just zeros and ones; if not, Invalid operation results. */
2716 /* ------------------------------------------------------------------ */
2717 decFloat * decFloatOr(decFloat *result,
2718 const decFloat *dfl, const decFloat *dfr,
2719 decContext *set) {
2720 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
2721 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
2722 /* the operands are positive finite integers (q=0) with just 0s and 1s */
2723 #if DOUBLE
2724 DFWORD(result, 0)=ZEROWORD
2725 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
2726 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
2727 #elif QUAD
2728 DFWORD(result, 0)=ZEROWORD
2729 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
2730 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
2731 DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
2732 DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
2733 #endif
2734 return result;
2735 } /* decFloatOr */
2736
2737 /* ------------------------------------------------------------------ */
2738 /* decFloatPlus -- add value to 0, heeding NaNs, etc. */
2739 /* */
2740 /* result gets the canonicalized 0+df */
2741 /* df is the decFloat to plus */
2742 /* set is the context */
2743 /* returns result */
2744 /* */
2745 /* This has the same effect as 0+df where the exponent of the zero is */
2746 /* the same as that of df (if df is finite). */
2747 /* The effect is also the same as decFloatCopy except that NaNs */
2748 /* are handled normally (the sign of a NaN is not affected, and an */
2749 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2750 /* ------------------------------------------------------------------ */
2751 decFloat * decFloatPlus(decFloat *result, const decFloat *df,
2752 decContext *set) {
2753 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2754 decCanonical(result, df); /* copy and check */
2755 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2756 return result;
2757 } /* decFloatPlus */
2758
2759 /* ------------------------------------------------------------------ */
2760 /* decFloatQuantize -- quantize a decFloat */
2761 /* */
2762 /* result gets the result of quantizing dfl to match dfr */
2763 /* dfl is the first decFloat (lhs) */
2764 /* dfr is the second decFloat (rhs), which sets the exponent */
2765 /* set is the context */
2766 /* returns result */
2767 /* */
2768 /* Unless there is an error or the result is infinite, the exponent */
2769 /* of result is guaranteed to be the same as that of dfr. */
2770 /* ------------------------------------------------------------------ */
2771 decFloat * decFloatQuantize(decFloat *result,
2772 const decFloat *dfl, const decFloat *dfr,
2773 decContext *set) {
2774 Int explb, exprb; /* left and right biased exponents */
2775 uByte *ulsd; /* local LSD pointer */
2776 uInt *ui; /* work */
2777 uByte *ub; /* .. */
2778 Int drop; /* .. */
2779 uInt dpd; /* .. */
2780 uInt encode; /* encoding accumulator */
2781 uInt sourhil, sourhir; /* top words from source decFloats */
2782 /* the following buffer holds the coefficient for manipulation */
2783 uByte buf[4+DECPMAX*3]; /* + space for zeros to left or right */
2784 #if DECTRACE
2785 bcdnum num; /* for trace displays */
2786 #endif
2787
2788 /* Start decoding the arguments */
2789 sourhil=DFWORD(dfl, 0); /* LHS top word */
2790 explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
2791 sourhir=DFWORD(dfr, 0); /* RHS top word */
2792 exprb=DECCOMBEXP[sourhir>>26];
2793
2794 if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
2795 /* NaNs are handled as usual */
2796 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2797 /* one infinity but not both is bad */
2798 if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
2799 /* both infinite; return canonical infinity with sign of LHS */
2800 return decInfinity(result, dfl);
2801 }
2802
2803 /* Here when both arguments are finite */
2804 /* complete extraction of the exponents [no need to unbias] */
2805 explb+=GETECON(dfl); /* + continuation */
2806 exprb+=GETECON(dfr); /* .. */
2807
2808 /* calculate the number of digits to drop from the coefficient */
2809 drop=exprb-explb; /* 0 if nothing to do */
2810 if (drop==0) return decCanonical(result, dfl); /* return canonical */
2811
2812 /* the coefficient is needed; lay it out into buf, offset so zeros */
2813 /* can be added before or after as needed -- an extra heading is */
2814 /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
2815 /* fours */
2816 #define BUFOFF (buf+4+DECPMAX)
2817 GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
2818 /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
2819
2820 #if DECTRACE
2821 num.msd=BUFOFF;
2822 num.lsd=BUFOFF+DECPMAX-1;
2823 num.exponent=explb-DECBIAS;
2824 num.sign=sourhil & DECFLOAT_Sign;
2825 decShowNum(&num, "dfl");
2826 #endif
2827
2828 if (drop>0) { /* [most common case] */
2829 /* (this code is very similar to that in decFloatFinalize, but */
2830 /* has many differences so is duplicated here -- so any changes */
2831 /* may need to be made there, too) */
2832 uByte *roundat; /* -> re-round digit */
2833 uByte reround; /* reround value */
2834 /* printf("Rounding; drop=%ld\n", (LI)drop); */
2835
2836 /* there is at least one zero needed to the left, in all but one */
2837 /* exceptional (all-nines) case, so place four zeros now; this is */
2838 /* needed almost always and makes rounding all-nines by fours safe */
2839 UINTAT(BUFOFF-4)=0;
2840
2841 /* Three cases here: */
2842 /* 1. new LSD is in coefficient (almost always) */
2843 /* 2. new LSD is digit to left of coefficient (so MSD is */
2844 /* round-for-reround digit) */
2845 /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
2846 /* Note that leading zeros can safely be treated as useful digits */
2847
2848 /* [duplicate check-stickies code to save a test] */
2849 /* [by-digit check for stickies as runs of zeros are rare] */
2850 if (drop<DECPMAX) { /* NB lengths not addresses */
2851 roundat=BUFOFF+DECPMAX-drop;
2852 reround=*roundat;
2853 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2854 if (*ub!=0) { /* non-zero to be discarded */
2855 reround=DECSTICKYTAB[reround]; /* apply sticky bit */
2856 break; /* [remainder don't-care] */
2857 }
2858 } /* check stickies */
2859 ulsd=roundat-1; /* set LSD */
2860 }
2861 else { /* edge case */
2862 if (drop==DECPMAX) {
2863 roundat=BUFOFF;
2864 reround=*roundat;
2865 }
2866 else {
2867 roundat=BUFOFF-1;
2868 reround=0;
2869 }
2870 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2871 if (*ub!=0) { /* non-zero to be discarded */
2872 reround=DECSTICKYTAB[reround]; /* apply sticky bit */
2873 break; /* [remainder don't-care] */
2874 }
2875 } /* check stickies */
2876 *BUFOFF=0; /* make a coefficient of 0 */
2877 ulsd=BUFOFF; /* .. at the MSD place */
2878 }
2879
2880 if (reround!=0) { /* discarding non-zero */
2881 uInt bump=0;
2882 set->status|=DEC_Inexact;
2883
2884 /* next decide whether to increment the coefficient */
2885 if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
2886 if (reround>5) bump=1; /* >0.5 goes up */
2887 else if (reround==5) /* exactly 0.5000 .. */
2888 bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */
2889 } /* r-h-e */
2890 else switch (set->round) {
2891 case DEC_ROUND_DOWN: {
2892 /* no change */
2893 break;} /* r-d */
2894 case DEC_ROUND_HALF_DOWN: {
2895 if (reround>5) bump=1;
2896 break;} /* r-h-d */
2897 case DEC_ROUND_HALF_UP: {
2898 if (reround>=5) bump=1;
2899 break;} /* r-h-u */
2900 case DEC_ROUND_UP: {
2901 if (reround>0) bump=1;
2902 break;} /* r-u */
2903 case DEC_ROUND_CEILING: {
2904 /* same as _UP for positive numbers, and as _DOWN for negatives */
2905 if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
2906 break;} /* r-c */
2907 case DEC_ROUND_FLOOR: {
2908 /* same as _UP for negative numbers, and as _DOWN for positive */
2909 /* [negative reround cannot occur on 0] */
2910 if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
2911 break;} /* r-f */
2912 case DEC_ROUND_05UP: {
2913 if (reround>0) { /* anything out there is 'sticky' */
2914 /* bump iff lsd=0 or 5; this cannot carry so it could be */
2915 /* effected immediately with no bump -- but the code */
2916 /* is clearer if this is done the same way as the others */
2917 if (*ulsd==0 || *ulsd==5) bump=1;
2918 }
2919 break;} /* r-r */
2920 default: { /* e.g., DEC_ROUND_MAX */
2921 set->status|=DEC_Invalid_context;
2922 #if DECCHECK
2923 printf("Unknown rounding mode: %ld\n", (LI)set->round);
2924 #endif
2925 break;}
2926 } /* switch (not r-h-e) */
2927 /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
2928
2929 if (bump!=0) { /* need increment */
2930 /* increment the coefficient; this could give 1000... (after */
2931 /* the all nines case) */
2932 ub=ulsd;
2933 for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
2934 /* now at most 3 digits left to non-9 (usually just the one) */
2935 for (; *ub==9; ub--) *ub=0;
2936 *ub+=1;
2937 /* [the all-nines case will have carried one digit to the */
2938 /* left of the original MSD -- just where it is needed] */
2939 } /* bump needed */
2940 } /* inexact rounding */
2941
2942 /* now clear zeros to the left so exactly DECPMAX digits will be */
2943 /* available in the coefficent -- the first word to the left was */
2944 /* cleared earlier for safe carry; now add any more needed */
2945 if (drop>4) {
2946 UINTAT(BUFOFF-8)=0; /* must be at least 5 */
2947 for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
2948 }
2949 } /* need round (drop>0) */
2950
2951 else { /* drop<0; padding with -drop digits is needed */
2952 /* This is the case where an error can occur if the padded */
2953 /* coefficient will not fit; checking for this can be done in the */
2954 /* same loop as padding for zeros if the no-hope and zero cases */
2955 /* are checked first */
2956 if (-drop>DECPMAX-1) { /* cannot fit unless 0 */
2957 if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
2958 /* a zero can have any exponent; just drop through and use it */
2959 ulsd=BUFOFF+DECPMAX-1;
2960 }
2961 else { /* padding will fit (but may still be too long) */
2962 /* final-word mask depends on endianess */
2963 #if DECLITEND
2964 static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
2965 #else
2966 static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
2967 #endif
2968 for (ui=&UINTAT(BUFOFF+DECPMAX);; ui++) {
2969 *ui=0;
2970 if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
2971 /* if all four digits should be zero, definitely bad */
2972 if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
2973 return decInvalid(result, set);
2974 /* must be a 1- to 3-digit sequence; check more carefully */
2975 if ((UINTAT(&UBYTEAT(ui)-DECPMAX)&dmask[(-drop)%4])!=0)
2976 return decInvalid(result, set);
2977 break; /* no need for loop end test */
2978 }
2979 if (ui>=&UINTAT(BUFOFF+DECPMAX+(-drop)-4)) break; /* done */
2980 }
2981 ulsd=BUFOFF+DECPMAX+(-drop)-1;
2982 } /* pad and check leading zeros */
2983 } /* drop<0 */
2984
2985 #if DECTRACE
2986 num.msd=ulsd-DECPMAX+1;
2987 num.lsd=ulsd;
2988 num.exponent=explb-DECBIAS;
2989 num.sign=sourhil & DECFLOAT_Sign;
2990 decShowNum(&num, "res");
2991 #endif
2992
2993 /*------------------------------------------------------------------*/
2994 /* At this point the result is DECPMAX digits, ending at ulsd, so */
2995 /* fits the encoding exactly; there is no possibility of error */
2996 /*------------------------------------------------------------------*/
2997 encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
2998 encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */
2999 /* the exponent continuation can be extracted from the original RHS */
3000 encode|=sourhir & ECONMASK;
3001 encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */
3002
3003 /* finally encode the coefficient */
3004 /* private macro to encode a declet; this version can be used */
3005 /* because all coefficient digits exist */
3006 #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
3007 dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
3008
3009 #if DOUBLE
3010 getDPD3q(dpd, 4); encode|=dpd<<8;
3011 getDPD3q(dpd, 3); encode|=dpd>>2;
3012 DFWORD(result, 0)=encode;
3013 encode=dpd<<30;
3014 getDPD3q(dpd, 2); encode|=dpd<<20;
3015 getDPD3q(dpd, 1); encode|=dpd<<10;
3016 getDPD3q(dpd, 0); encode|=dpd;
3017 DFWORD(result, 1)=encode;
3018
3019 #elif QUAD
3020 getDPD3q(dpd,10); encode|=dpd<<4;
3021 getDPD3q(dpd, 9); encode|=dpd>>6;
3022 DFWORD(result, 0)=encode;
3023 encode=dpd<<26;
3024 getDPD3q(dpd, 8); encode|=dpd<<16;
3025 getDPD3q(dpd, 7); encode|=dpd<<6;
3026 getDPD3q(dpd, 6); encode|=dpd>>4;
3027 DFWORD(result, 1)=encode;
3028 encode=dpd<<28;
3029 getDPD3q(dpd, 5); encode|=dpd<<18;
3030 getDPD3q(dpd, 4); encode|=dpd<<8;
3031 getDPD3q(dpd, 3); encode|=dpd>>2;
3032 DFWORD(result, 2)=encode;
3033 encode=dpd<<30;
3034 getDPD3q(dpd, 2); encode|=dpd<<20;
3035 getDPD3q(dpd, 1); encode|=dpd<<10;
3036 getDPD3q(dpd, 0); encode|=dpd;
3037 DFWORD(result, 3)=encode;
3038 #endif
3039 return result;
3040 } /* decFloatQuantize */
3041
3042 /* ------------------------------------------------------------------ */
3043 /* decFloatReduce -- reduce finite coefficient to minimum length */
3044 /* */
3045 /* result gets the reduced decFloat */
3046 /* df is the source decFloat */
3047 /* set is the context */
3048 /* returns result, which will be canonical */
3049 /* */
3050 /* This removes all possible trailing zeros from the coefficient; */
3051 /* some may remain when the number is very close to Nmax. */
3052 /* Special values are unchanged and no status is set unless df=sNaN. */
3053 /* Reduced zero has an exponent q=0. */
3054 /* ------------------------------------------------------------------ */
3055 decFloat * decFloatReduce(decFloat *result, const decFloat *df,
3056 decContext *set) {
3057 bcdnum num; /* work */
3058 uByte buf[DECPMAX], *ub; /* coefficient and pointer */
3059 if (df!=result) *result=*df; /* copy, if needed */
3060 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */
3061 /* zeros and infinites propagate too */
3062 if (DFISINF(df)) return decInfinity(result, df); /* canonical */
3063 if (DFISZERO(df)) {
3064 uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
3065 decFloatZero(result);
3066 DFWORD(result, 0)|=sign;
3067 return result; /* exponent dropped, sign OK */
3068 }
3069 /* non-zero finite */
3070 GETCOEFF(df, buf);
3071 ub=buf+DECPMAX-1; /* -> lsd */
3072 if (*ub) return result; /* no trailing zeros */
3073 for (ub--; *ub==0;) ub--; /* terminates because non-zero */
3074 /* *ub is the first non-zero from the right */
3075 num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
3076 num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
3077 num.msd=buf;
3078 num.lsd=ub;
3079 return decFinalize(result, &num, set);
3080 } /* decFloatReduce */
3081
3082 /* ------------------------------------------------------------------ */
3083 /* decFloatRemainder -- integer divide and return remainder */
3084 /* */
3085 /* result gets the remainder of dividing dfl by dfr: */
3086 /* dfl is the first decFloat (lhs) */
3087 /* dfr is the second decFloat (rhs) */
3088 /* set is the context */
3089 /* returns result */
3090 /* */
3091 /* ------------------------------------------------------------------ */
3092 decFloat * decFloatRemainder(decFloat *result,
3093 const decFloat *dfl, const decFloat *dfr,
3094 decContext *set) {
3095 return decDivide(result, dfl, dfr, set, REMAINDER);
3096 } /* decFloatRemainder */
3097
3098 /* ------------------------------------------------------------------ */
3099 /* decFloatRemainderNear -- integer divide to nearest and remainder */
3100 /* */
3101 /* result gets the remainder of dividing dfl by dfr: */
3102 /* dfl is the first decFloat (lhs) */
3103 /* dfr is the second decFloat (rhs) */
3104 /* set is the context */
3105 /* returns result */
3106 /* */
3107 /* This is the IEEE remainder, where the nearest integer is used. */
3108 /* ------------------------------------------------------------------ */
3109 decFloat * decFloatRemainderNear(decFloat *result,
3110 const decFloat *dfl, const decFloat *dfr,
3111 decContext *set) {
3112 return decDivide(result, dfl, dfr, set, REMNEAR);
3113 } /* decFloatRemainderNear */
3114
3115 /* ------------------------------------------------------------------ */
3116 /* decFloatRotate -- rotate the coefficient of a decFloat left/right */
3117 /* */
3118 /* result gets the result of rotating dfl */
3119 /* dfl is the source decFloat to rotate */
3120 /* dfr is the count of digits to rotate, an integer (with q=0) */
3121 /* set is the context */
3122 /* returns result */
3123 /* */
3124 /* The digits of the coefficient of dfl are rotated to the left (if */
3125 /* dfr is positive) or to the right (if dfr is negative) without */
3126 /* adjusting the exponent or the sign of dfl. */
3127 /* */
3128 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3129 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3130 /* dfr must be valid). No status is set unless dfr is invalid or an */
3131 /* operand is an sNaN. The result is canonical. */
3132 /* ------------------------------------------------------------------ */
3133 #define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
3134 decFloat * decFloatRotate(decFloat *result,
3135 const decFloat *dfl, const decFloat *dfr,
3136 decContext *set) {
3137 Int rotate; /* dfr as an Int */
3138 uByte buf[DECPMAX+PHALF]; /* coefficient + half */
3139 uInt digits, savestat; /* work */
3140 bcdnum num; /* .. */
3141 uByte *ub; /* .. */
3142
3143 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3144 if (!DFISINT(dfr)) return decInvalid(result, set);
3145 digits=decFloatDigits(dfr); /* calculate digits */
3146 if (digits>2) return decInvalid(result, set); /* definitely out of range */
3147 rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3148 if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
3149 /* [from here on no error or status change is possible] */
3150 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3151 /* handle no-rotate cases */
3152 if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
3153 /* a real rotate is needed: 0 < rotate < DECPMAX */
3154 /* reduce the rotation to no more than half to reduce copying later */
3155 /* (for QUAD in fact half + 2 digits) */
3156 if (DFISSIGNED(dfr)) rotate=-rotate;
3157 if (abs(rotate)>PHALF) {
3158 if (rotate<0) rotate=DECPMAX+rotate;
3159 else rotate=rotate-DECPMAX;
3160 }
3161 /* now lay out the coefficient, leaving room to the right or the */
3162 /* left depending on the direction of rotation */
3163 ub=buf;
3164 if (rotate<0) ub+=PHALF; /* rotate right, so space to left */
3165 GETCOEFF(dfl, ub);
3166 /* copy half the digits to left or right, and set num.msd */
3167 if (rotate<0) {
3168 memcpy(buf, buf+DECPMAX, PHALF);
3169 num.msd=buf+PHALF+rotate;
3170 }
3171 else {
3172 memcpy(buf+DECPMAX, buf, PHALF);
3173 num.msd=buf+rotate;
3174 }
3175 /* fill in rest of num */
3176 num.lsd=num.msd+DECPMAX-1;
3177 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3178 num.exponent=GETEXPUN(dfl);
3179 savestat=set->status; /* record */
3180 decFinalize(result, &num, set);
3181 set->status=savestat; /* restore */
3182 return result;
3183 } /* decFloatRotate */
3184
3185 /* ------------------------------------------------------------------ */
3186 /* decFloatSameQuantum -- test decFloats for same quantum */
3187 /* */
3188 /* dfl is the first decFloat (lhs) */
3189 /* dfr is the second decFloat (rhs) */
3190 /* returns 1 if the operands have the same quantum, 0 otherwise */
3191 /* */
3192 /* No error is possible and no status results. */
3193 /* ------------------------------------------------------------------ */
3194 uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
3195 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
3196 if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
3197 if (DFISINF(dfl) && DFISINF(dfr)) return 1;
3198 return 0; /* any other special mixture gives false */
3199 }
3200 if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
3201 return 0;
3202 } /* decFloatSameQuantum */
3203
3204 /* ------------------------------------------------------------------ */
3205 /* decFloatScaleB -- multiply by a power of 10, as per 754r */
3206 /* */
3207 /* result gets the result of the operation */
3208 /* dfl is the first decFloat (lhs) */
3209 /* dfr is the second decFloat (rhs), am integer (with q=0) */
3210 /* set is the context */
3211 /* returns result */
3212 /* */
3213 /* This computes result=dfl x 10**dfr where dfr is an integer in the */
3214 /* range +/-2*(emax+pmax), typically resulting from LogB. */
3215 /* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
3216 /* as usual. */
3217 /* ------------------------------------------------------------------ */
3218 #define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
3219 decFloat * decFloatScaleB(decFloat *result,
3220 const decFloat *dfl, const decFloat *dfr,
3221 decContext *set) {
3222 uInt digits; /* work */
3223 Int expr; /* dfr as an Int */
3224
3225 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3226 if (!DFISINT(dfr)) return decInvalid(result, set);
3227 digits=decFloatDigits(dfr); /* calculate digits */
3228
3229 #if DOUBLE
3230 if (digits>3) return decInvalid(result, set); /* definitely out of range */
3231 expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */
3232 #elif QUAD
3233 if (digits>5) return decInvalid(result, set); /* definitely out of range */
3234 expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */
3235 +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */
3236 #endif
3237 if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */
3238 /* [from now on no error possible] */
3239 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3240 if (DFISSIGNED(dfr)) expr=-expr;
3241 /* dfl is finite and expr is valid */
3242 *result=*dfl; /* copy to target */
3243 return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
3244 } /* decFloatScaleB */
3245
3246 /* ------------------------------------------------------------------ */
3247 /* decFloatShift -- shift the coefficient of a decFloat left or right */
3248 /* */
3249 /* result gets the result of shifting dfl */
3250 /* dfl is the source decFloat to shift */
3251 /* dfr is the count of digits to shift, an integer (with q=0) */
3252 /* set is the context */
3253 /* returns result */
3254 /* */
3255 /* The digits of the coefficient of dfl are shifted to the left (if */
3256 /* dfr is positive) or to the right (if dfr is negative) without */
3257 /* adjusting the exponent or the sign of dfl. */
3258 /* */
3259 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3260 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3261 /* dfr must be valid). No status is set unless dfr is invalid or an */
3262 /* operand is an sNaN. The result is canonical. */
3263 /* ------------------------------------------------------------------ */
3264 decFloat * decFloatShift(decFloat *result,
3265 const decFloat *dfl, const decFloat *dfr,
3266 decContext *set) {
3267 Int shift; /* dfr as an Int */
3268 uByte buf[DECPMAX*2]; /* coefficient + padding */
3269 uInt digits, savestat; /* work */
3270 bcdnum num; /* .. */
3271
3272 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3273 if (!DFISINT(dfr)) return decInvalid(result, set);
3274 digits=decFloatDigits(dfr); /* calculate digits */
3275 if (digits>2) return decInvalid(result, set); /* definitely out of range */
3276 shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3277 if (shift>DECPMAX) return decInvalid(result, set); /* too big */
3278 /* [from here on no error or status change is possible] */
3279
3280 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3281 /* handle no-shift and all-shift (clear to zero) cases */
3282 if (shift==0) return decCanonical(result, dfl);
3283 if (shift==DECPMAX) { /* zero with sign */
3284 uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
3285 decFloatZero(result); /* make +0 */
3286 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
3287 /* [cannot safely use CopySign] */
3288 return result;
3289 }
3290 /* a real shift is needed: 0 < shift < DECPMAX */
3291 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3292 num.exponent=GETEXPUN(dfl);
3293 num.msd=buf;
3294 GETCOEFF(dfl, buf);
3295 if (DFISSIGNED(dfr)) { /* shift right */
3296 /* edge cases are taken care of, so this is easy */
3297 num.lsd=buf+DECPMAX-shift-1;
3298 }
3299 else { /* shift left -- zero padding needed to right */
3300 UINTAT(buf+DECPMAX)=0; /* 8 will handle most cases */
3301 UINTAT(buf+DECPMAX+4)=0; /* .. */
3302 if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
3303 num.msd+=shift;
3304 num.lsd=num.msd+DECPMAX-1;
3305 }
3306 savestat=set->status; /* record */
3307 decFinalize(result, &num, set);
3308 set->status=savestat; /* restore */
3309 return result;
3310 } /* decFloatShift */
3311
3312 /* ------------------------------------------------------------------ */
3313 /* decFloatSubtract -- subtract a decFloat from another */
3314 /* */
3315 /* result gets the result of subtracting dfr from dfl: */
3316 /* dfl is the first decFloat (lhs) */
3317 /* dfr is the second decFloat (rhs) */
3318 /* set is the context */
3319 /* returns result */
3320 /* */
3321 /* ------------------------------------------------------------------ */
3322 decFloat * decFloatSubtract(decFloat *result,
3323 const decFloat *dfl, const decFloat *dfr,
3324 decContext *set) {
3325 decFloat temp;
3326 /* NaNs must propagate without sign change */
3327 if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
3328 temp=*dfr; /* make a copy */
3329 DFBYTE(&temp, 0)^=0x80; /* flip sign */
3330 return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
3331 } /* decFloatSubtract */
3332
3333 /* ------------------------------------------------------------------ */
3334 /* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
3335 /* */
3336 /* df is the decFloat to round */
3337 /* set is the context */
3338 /* round is the rounding mode to use */
3339 /* returns a uInt or an Int, rounded according to the name */
3340 /* */
3341 /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
3342 /* outside the range of the target; Inexact will not be signaled for */
3343 /* simple rounding unless 'Exact' appears in the name. */
3344 /* ------------------------------------------------------------------ */
3345 uInt decFloatToUInt32(const decFloat *df, decContext *set,
3346 enum rounding round) {
3347 return decToInt32(df, set, round, 0, 1);}
3348
3349 uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
3350 enum rounding round) {
3351 return decToInt32(df, set, round, 1, 1);}
3352
3353 Int decFloatToInt32(const decFloat *df, decContext *set,
3354 enum rounding round) {
3355 return (Int)decToInt32(df, set, round, 0, 0);}
3356
3357 Int decFloatToInt32Exact(const decFloat *df, decContext *set,
3358 enum rounding round) {
3359 return (Int)decToInt32(df, set, round, 1, 0);}
3360
3361 /* ------------------------------------------------------------------ */
3362 /* decFloatToIntegral -- round to integral value (two flavours) */
3363 /* */
3364 /* result gets the result */
3365 /* df is the decFloat to round */
3366 /* set is the context */
3367 /* round is the rounding mode to use */
3368 /* returns result */
3369 /* */
3370 /* No exceptions, even Inexact, are raised except for sNaN input, or */
3371 /* if 'Exact' appears in the name. */
3372 /* ------------------------------------------------------------------ */
3373 decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
3374 decContext *set, enum rounding round) {
3375 return decToIntegral(result, df, set, round, 0);}
3376
3377 decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
3378 decContext *set) {
3379 return decToIntegral(result, df, set, set->round, 1);}
3380
3381 /* ------------------------------------------------------------------ */
3382 /* decFloatXor -- logical digitwise XOR of two decFloats */
3383 /* */
3384 /* result gets the result of XORing dfl and dfr */
3385 /* dfl is the first decFloat (lhs) */
3386 /* dfr is the second decFloat (rhs) */
3387 /* set is the context */
3388 /* returns result, which will be canonical with sign=0 */
3389 /* */
3390 /* The operands must be positive, finite with exponent q=0, and */
3391 /* comprise just zeros and ones; if not, Invalid operation results. */
3392 /* ------------------------------------------------------------------ */
3393 decFloat * decFloatXor(decFloat *result,
3394 const decFloat *dfl, const decFloat *dfr,
3395 decContext *set) {
3396 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
3397 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
3398 /* the operands are positive finite integers (q=0) with just 0s and 1s */
3399 #if DOUBLE
3400 DFWORD(result, 0)=ZEROWORD
3401 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
3402 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
3403 #elif QUAD
3404 DFWORD(result, 0)=ZEROWORD
3405 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
3406 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
3407 DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
3408 DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
3409 #endif
3410 return result;
3411 } /* decFloatXor */
3412
3413 /* ------------------------------------------------------------------ */
3414 /* decInvalid -- set Invalid_operation result */
3415 /* */
3416 /* result gets a canonical NaN */
3417 /* set is the context */
3418 /* returns result */
3419 /* */
3420 /* status has Invalid_operation added */
3421 /* ------------------------------------------------------------------ */
3422 static decFloat *decInvalid(decFloat *result, decContext *set) {
3423 decFloatZero(result);
3424 DFWORD(result, 0)=DECFLOAT_qNaN;
3425 set->status|=DEC_Invalid_operation;
3426 return result;
3427 } /* decInvalid */
3428
3429 /* ------------------------------------------------------------------ */
3430 /* decInfinity -- set canonical Infinity with sign from a decFloat */
3431 /* */
3432 /* result gets a canonical Infinity */
3433 /* df is source decFloat (only the sign is used) */
3434 /* returns result */
3435 /* */
3436 /* df may be the same as result */
3437 /* ------------------------------------------------------------------ */
3438 static decFloat *decInfinity(decFloat *result, const decFloat *df) {
3439 uInt sign=DFWORD(df, 0); /* save source signword */
3440 decFloatZero(result); /* clear everything */
3441 DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
3442 return result;
3443 } /* decInfinity */
3444
3445 /* ------------------------------------------------------------------ */
3446 /* decNaNs -- handle NaN argument(s) */
3447 /* */
3448 /* result gets the result of handling dfl and dfr, one or both of */
3449 /* which is a NaN */
3450 /* dfl is the first decFloat (lhs) */
3451 /* dfr is the second decFloat (rhs) -- may be NULL for a single- */
3452 /* operand operation */
3453 /* set is the context */
3454 /* returns result */
3455 /* */
3456 /* Called when one or both operands is a NaN, and propagates the */
3457 /* appropriate result to res. When an sNaN is found, it is changed */
3458 /* to a qNaN and Invalid operation is set. */
3459 /* ------------------------------------------------------------------ */
3460 static decFloat *decNaNs(decFloat *result,
3461 const decFloat *dfl, const decFloat *dfr,
3462 decContext *set) {
3463 /* handle sNaNs first */
3464 if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
3465 if (DFISSNAN(dfl)) {
3466 decCanonical(result, dfl); /* propagate canonical sNaN */
3467 DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
3468 set->status|=DEC_Invalid_operation;
3469 return result;
3470 }
3471 /* one or both is a quiet NaN */
3472 if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */
3473 return decCanonical(result, dfl); /* propagate canonical qNaN */
3474 } /* decNaNs */
3475
3476 /* ------------------------------------------------------------------ */
3477 /* decNumCompare -- numeric comparison of two decFloats */
3478 /* */
3479 /* dfl is the left-hand decFloat, which is not a NaN */
3480 /* dfr is the right-hand decFloat, which is not a NaN */
3481 /* tot is 1 for total order compare, 0 for simple numeric */
3482 /* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
3483 /* */
3484 /* No error is possible; status and mode are unchanged. */
3485 /* ------------------------------------------------------------------ */
3486 static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
3487 Int sigl, sigr; /* LHS and RHS non-0 signums */
3488 Int shift; /* shift needed to align operands */
3489 uByte *ub, *uc; /* work */
3490 /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
3491 uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */
3492 uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */
3493
3494 sigl=1;
3495 if (DFISSIGNED(dfl)) {
3496 if (!DFISSIGNED(dfr)) { /* -LHS +RHS */
3497 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3498 return -1; /* RHS wins */
3499 }
3500 sigl=-1;
3501 }
3502 if (DFISSIGNED(dfr)) {
3503 if (!DFISSIGNED(dfl)) { /* +LHS -RHS */
3504 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3505 return +1; /* LHS wins */
3506 }
3507 }
3508
3509 /* signs are the same; operand(s) could be zero */
3510 sigr=-sigl; /* sign to return if abs(RHS) wins */
3511
3512 if (DFISINF(dfl)) {
3513 if (DFISINF(dfr)) return 0; /* both infinite & same sign */
3514 return sigl; /* inf > n */
3515 }
3516 if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
3517
3518 /* here, both are same sign and finite; calculate their offset */
3519 shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */
3520 /* [bias can be ignored -- the absolute exponent is not relevant] */
3521
3522 if (DFISZERO(dfl)) {
3523 if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */
3524 /* both are zero, return 0 if both same exponent or numeric compare */
3525 if (shift==0 || !tot) return 0;
3526 if (shift>0) return sigl;
3527 return sigr; /* [shift<0] */
3528 }
3529 else { /* LHS!=0 */
3530 if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */
3531 }
3532 /* both are known to be non-zero at this point */
3533
3534 /* if the exponents are so different that the coefficients do not */
3535 /* overlap (by even one digit) then a full comparison is not needed */
3536 if (abs(shift)>=DECPMAX) { /* no overlap */
3537 /* coefficients are known to be non-zero */
3538 if (shift>0) return sigl;
3539 return sigr; /* [shift<0] */
3540 }
3541
3542 /* decode the coefficients */
3543 /* (shift both right two if Quad to make a multiple of four) */
3544 #if QUAD
3545 UINTAT(bufl)=0;
3546 UINTAT(bufr)=0;
3547 #endif
3548 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
3549 GETCOEFF(dfr, bufr+QUAD*2); /* .. */
3550 if (shift==0) { /* aligned; common and easy */
3551 /* all multiples of four, here */
3552 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
3553 if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
3554 /* about to find a winner; go by bytes in case little-endian */
3555 for (;; ub++, uc++) {
3556 if (*ub>*uc) return sigl; /* difference found */
3557 if (*ub<*uc) return sigr; /* .. */
3558 }
3559 }
3560 } /* aligned */
3561 else if (shift>0) { /* lhs to left */
3562 ub=bufl; /* RHS pointer */
3563 /* pad bufl so right-aligned; most shifts will fit in 8 */
3564 UINTAT(bufl+DECPMAX+QUAD*2)=0; /* add eight zeros */
3565 UINTAT(bufl+DECPMAX+QUAD*2+4)=0; /* .. */
3566 if (shift>8) {
3567 /* more than eight; fill the rest, and also worth doing the */
3568 /* lead-in by fours */
3569 uByte *up; /* work */
3570 uByte *upend=bufl+DECPMAX+QUAD*2+shift;
3571 for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3572 /* [pads up to 36 in all for Quad] */
3573 for (;; ub+=4) {
3574 if (UINTAT(ub)!=0) return sigl;
3575 if (ub+4>bufl+shift-4) break;
3576 }
3577 }
3578 /* check remaining leading digits */
3579 for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
3580 /* now start the overlapped part; bufl has been padded, so the */
3581 /* comparison can go for the full length of bufr, which is a */
3582 /* multiple of 4 bytes */
3583 for (uc=bufr; ; uc+=4, ub+=4) {
3584 if (UINTAT(uc)!=UINTAT(ub)) { /* mismatch found */
3585 for (;; uc++, ub++) { /* check from left [little-endian?] */
3586 if (*ub>*uc) return sigl; /* difference found */
3587 if (*ub<*uc) return sigr; /* .. */
3588 }
3589 } /* mismatch */
3590 if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
3591 }
3592 } /* shift>0 */
3593
3594 else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
3595 uc=bufr; /* RHS pointer */
3596 /* pad bufr so right-aligned; most shifts will fit in 8 */
3597 UINTAT(bufr+DECPMAX+QUAD*2)=0; /* add eight zeros */
3598 UINTAT(bufr+DECPMAX+QUAD*2+4)=0; /* .. */
3599 if (shift<-8) {
3600 /* more than eight; fill the rest, and also worth doing the */
3601 /* lead-in by fours */
3602 uByte *up; /* work */
3603 uByte *upend=bufr+DECPMAX+QUAD*2-shift;
3604 for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3605 /* [pads up to 36 in all for Quad] */
3606 for (;; uc+=4) {
3607 if (UINTAT(uc)!=0) return sigr;
3608 if (uc+4>bufr-shift-4) break;
3609 }
3610 }
3611 /* check remaining leading digits */
3612 for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
3613 /* now start the overlapped part; bufr has been padded, so the */
3614 /* comparison can go for the full length of bufl, which is a */
3615 /* multiple of 4 bytes */
3616 for (ub=bufl; ; ub+=4, uc+=4) {
3617 if (UINTAT(ub)!=UINTAT(uc)) { /* mismatch found */
3618 for (;; ub++, uc++) { /* check from left [little-endian?] */
3619 if (*ub>*uc) return sigl; /* difference found */
3620 if (*ub<*uc) return sigr; /* .. */
3621 }
3622 } /* mismatch */
3623 if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
3624 }
3625 } /* shift<0 */
3626
3627 /* Here when compare equal */
3628 if (!tot) return 0; /* numerically equal */
3629 /* total ordering .. exponent matters */
3630 if (shift>0) return sigl; /* total order by exponent */
3631 if (shift<0) return sigr; /* .. */
3632 return 0;
3633 } /* decNumCompare */
3634
3635 /* ------------------------------------------------------------------ */
3636 /* decToInt32 -- local routine to effect ToInteger conversions */
3637 /* */
3638 /* df is the decFloat to convert */
3639 /* set is the context */
3640 /* rmode is the rounding mode to use */
3641 /* exact is 1 if Inexact should be signalled */
3642 /* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
3643 /* returns 32-bit result as a uInt */
3644 /* */
3645 /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
3646 /* these cases 0 is returned. */
3647 /* ------------------------------------------------------------------ */
3648 static uInt decToInt32(const decFloat *df, decContext *set,
3649 enum rounding rmode, Flag exact, Flag unsign) {
3650 Int exp; /* exponent */
3651 uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */
3652 uInt hi, lo; /* .. penultimate, least, etc. */
3653 decFloat zero, result; /* work */
3654 Int i; /* .. */
3655
3656 /* Start decoding the argument */
3657 sourhi=DFWORD(df, 0); /* top word */
3658 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
3659 if (EXPISSPECIAL(exp)) { /* is special? */
3660 set->status|=DEC_Invalid_operation; /* signal */
3661 return 0;
3662 }
3663
3664 /* Here when the argument is finite */
3665 if (GETEXPUN(df)==0) result=*df; /* already a true integer */
3666 else { /* need to round to integer */
3667 enum rounding saveround; /* saver */
3668 uInt savestatus; /* .. */
3669 saveround=set->round; /* save rounding mode .. */
3670 savestatus=set->status; /* .. and status */
3671 set->round=rmode; /* set mode */
3672 decFloatZero(&zero); /* make 0E+0 */
3673 set->status=0; /* clear */
3674 decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
3675 set->round=saveround; /* restore rounding mode .. */
3676 if (exact) set->status|=savestatus; /* include Inexact */
3677 else set->status=savestatus; /* .. or just original status */
3678 }
3679
3680 /* only the last four declets of the coefficient can contain */
3681 /* non-zero; check for others (and also NaN or Infinity from the */
3682 /* Quantize) first (see DFISZERO for explanation): */
3683 /* decFloatShow(&result, "sofar"); */
3684 #if DOUBLE
3685 if ((DFWORD(&result, 0)&0x1c03ff00)!=0
3686 || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3687 #elif QUAD
3688 if ((DFWORD(&result, 2)&0xffffff00)!=0
3689 || DFWORD(&result, 1)!=0
3690 || (DFWORD(&result, 0)&0x1c003fff)!=0
3691 || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3692 #endif
3693 set->status|=DEC_Invalid_operation; /* Invalid or out of range */
3694 return 0;
3695 }
3696 /* get last twelve digits of the coefficent into hi & ho, base */
3697 /* 10**9 (see GETCOEFFBILL): */
3698 sourlo=DFWORD(&result, DECWORDS-1);
3699 lo=DPD2BIN0[sourlo&0x3ff]
3700 +DPD2BINK[(sourlo>>10)&0x3ff]
3701 +DPD2BINM[(sourlo>>20)&0x3ff];
3702 sourpen=DFWORD(&result, DECWORDS-2);
3703 hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
3704
3705 /* according to request, check range carefully */
3706 if (unsign) {
3707 if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
3708 set->status|=DEC_Invalid_operation; /* out of range */
3709 return 0;
3710 }
3711 return hi*BILLION+lo;
3712 }
3713 /* signed */
3714 if (hi>2 || (hi==2 && lo>147483647)) {
3715 /* handle the usual edge case */
3716 if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
3717 set->status|=DEC_Invalid_operation; /* truly out of range */
3718 return 0;
3719 }
3720 i=hi*BILLION+lo;
3721 if (DFISSIGNED(&result)) i=-i;
3722 return (uInt)i;
3723 } /* decToInt32 */
3724
3725 /* ------------------------------------------------------------------ */
3726 /* decToIntegral -- local routine to effect ToIntegral value */
3727 /* */
3728 /* result gets the result */
3729 /* df is the decFloat to round */
3730 /* set is the context */
3731 /* rmode is the rounding mode to use */
3732 /* exact is 1 if Inexact should be signalled */
3733 /* returns result */
3734 /* ------------------------------------------------------------------ */
3735 static decFloat * decToIntegral(decFloat *result, const decFloat *df,
3736 decContext *set, enum rounding rmode,
3737 Flag exact) {
3738 Int exp; /* exponent */
3739 uInt sourhi; /* top word from source decFloat */
3740 enum rounding saveround; /* saver */
3741 uInt savestatus; /* .. */
3742 decFloat zero; /* work */
3743
3744 /* Start decoding the argument */
3745 sourhi=DFWORD(df, 0); /* top word */
3746 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
3747
3748 if (EXPISSPECIAL(exp)) { /* is special? */
3749 /* NaNs are handled as usual */
3750 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
3751 /* must be infinite; return canonical infinity with sign of df */
3752 return decInfinity(result, df);
3753 }
3754
3755 /* Here when the argument is finite */
3756 /* complete extraction of the exponent */
3757 exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */
3758
3759 if (exp>=0) return decCanonical(result, df); /* already integral */
3760
3761 saveround=set->round; /* save rounding mode .. */
3762 savestatus=set->status; /* .. and status */
3763 set->round=rmode; /* set mode */
3764 decFloatZero(&zero); /* make 0E+0 */
3765 decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
3766 set->round=saveround; /* restore rounding mode .. */
3767 if (!exact) set->status=savestatus; /* .. and status, unless exact */
3768 return result;
3769 } /* decToIntegral */
GDB Binutils - Deliverables
This page took 0.135232 seconds and 4 git commands to generate.