]> git.ipfire.org Git - thirdparty/gcc.git/blob - libdecnumber/decBasic.c
projects / thirdparty / gcc.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Temporarily revert revision 145222.
[thirdparty/gcc.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 ub = bufl; /* avoid type-pun violation */
1664 USHORTAT(ub)=0;
1665 uc = bufr; /* avoid type-pun violation */
1666 USHORTAT(uc)=0;
1667 #endif
1668 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
1669 GETCOEFF(dfr, bufr+QUAD*2); /* .. */
1670 /* all multiples of four, here */
1671 comp=0; /* assume equal */
1672 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
1673 if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
1674 /* about to find a winner; go by bytes in case little-endian */
1675 for (;; ub++, uc++) {
1676 if (*ub==*uc) continue;
1677 if (*ub>*uc) comp=sigl; /* difference found */
1678 else comp=-sigl; /* .. */
1679 break;
1680 }
1681 }
1682 } /* same NaN type and sign */
1683 }
1684 else {
1685 /* numeric comparison needed */
1686 comp=decNumCompare(dfl, dfr, 1); /* total ordering */
1687 }
1688 decFloatZero(result);
1689 if (comp==0) return result;
1690 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
1691 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
1692 return result;
1693 } /* decFloatCompareTotal */
1694
1695 /* ------------------------------------------------------------------ */
1696 /* decFloatCompareTotalMag -- compare magnitudes with total ordering */
1697 /* */
1698 /* result gets the result of comparing abs(dfl) and abs(dfr) */
1699 /* dfl is the first decFloat (lhs) */
1700 /* dfr is the second decFloat (rhs) */
1701 /* returns result, which may be -1, 0, or 1 */
1702 /* ------------------------------------------------------------------ */
1703 decFloat * decFloatCompareTotalMag(decFloat *result,
1704 const decFloat *dfl, const decFloat *dfr) {
1705 decFloat a, b; /* for copy if needed */
1706 /* copy and redirect signed operand(s) */
1707 if (DFISSIGNED(dfl)) {
1708 decFloatCopyAbs(&a, dfl);
1709 dfl=&a;
1710 }
1711 if (DFISSIGNED(dfr)) {
1712 decFloatCopyAbs(&b, dfr);
1713 dfr=&b;
1714 }
1715 return decFloatCompareTotal(result, dfl, dfr);
1716 } /* decFloatCompareTotalMag */
1717
1718 /* ------------------------------------------------------------------ */
1719 /* decFloatCopy -- copy a decFloat as-is */
1720 /* */
1721 /* result gets the copy of dfl */
1722 /* dfl is the decFloat to copy */
1723 /* returns result */
1724 /* */
1725 /* This is a bitwise operation; no errors or exceptions are possible. */
1726 /* ------------------------------------------------------------------ */
1727 decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
1728 if (dfl!=result) *result=*dfl; /* copy needed */
1729 return result;
1730 } /* decFloatCopy */
1731
1732 /* ------------------------------------------------------------------ */
1733 /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
1734 /* */
1735 /* result gets the copy of dfl with sign bit 0 */
1736 /* dfl is the decFloat to copy */
1737 /* returns result */
1738 /* */
1739 /* This is a bitwise operation; no errors or exceptions are possible. */
1740 /* ------------------------------------------------------------------ */
1741 decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
1742 if (dfl!=result) *result=*dfl; /* copy needed */
1743 DFBYTE(result, 0)&=~0x80; /* zero sign bit */
1744 return result;
1745 } /* decFloatCopyAbs */
1746
1747 /* ------------------------------------------------------------------ */
1748 /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
1749 /* */
1750 /* result gets the copy of dfl with sign bit inverted */
1751 /* dfl is the decFloat to copy */
1752 /* returns result */
1753 /* */
1754 /* This is a bitwise operation; no errors or exceptions are possible. */
1755 /* ------------------------------------------------------------------ */
1756 decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
1757 if (dfl!=result) *result=*dfl; /* copy needed */
1758 DFBYTE(result, 0)^=0x80; /* invert sign bit */
1759 return result;
1760 } /* decFloatCopyNegate */
1761
1762 /* ------------------------------------------------------------------ */
1763 /* decFloatCopySign -- copy a decFloat with the sign of another */
1764 /* */
1765 /* result gets the result of copying dfl with the sign of dfr */
1766 /* dfl is the first decFloat (lhs) */
1767 /* dfr is the second decFloat (rhs) */
1768 /* returns result */
1769 /* */
1770 /* This is a bitwise operation; no errors or exceptions are possible. */
1771 /* ------------------------------------------------------------------ */
1772 decFloat * decFloatCopySign(decFloat *result,
1773 const decFloat *dfl, const decFloat *dfr) {
1774 uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
1775 if (dfl!=result) *result=*dfl; /* copy needed */
1776 DFBYTE(result, 0)&=~0x80; /* clear sign .. */
1777 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
1778 return result;
1779 } /* decFloatCopySign */
1780
1781 /* ------------------------------------------------------------------ */
1782 /* decFloatDigits -- return the number of digits in a decFloat */
1783 /* */
1784 /* df is the decFloat to investigate */
1785 /* returns the number of significant digits in the decFloat; a */
1786 /* zero coefficient returns 1 as does an infinity (a NaN returns */
1787 /* the number of digits in the payload) */
1788 /* ------------------------------------------------------------------ */
1789 /* private macro to extract a declet according to provided formula */
1790 /* (form), and if it is non-zero then return the calculated digits */
1791 /* depending on the declet number (n), where n=0 for the most */
1792 /* significant declet; uses uInt dpd for work */
1793 #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
1794 if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1795 /* next one is used when it is known that the declet must be */
1796 /* non-zero, or is the final zero declet */
1797 #define dpdlendun(n, form) {dpd=(form)&0x3ff; \
1798 if (dpd==0) return 1; \
1799 return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1800
1801 uInt decFloatDigits(const decFloat *df) {
1802 uInt dpd; /* work */
1803 uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
1804 #if QUAD
1805 uInt sourmh, sourml;
1806 #endif
1807 uInt sourlo;
1808
1809 if (DFISINF(df)) return 1;
1810 /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
1811 /* then the coefficient is full-length */
1812 if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
1813
1814 #if DOUBLE
1815 if (sourhi&0x0003ffff) { /* ends in first */
1816 dpdlenchk(0, sourhi>>8);
1817 sourlo=DFWORD(df, 1);
1818 dpdlendun(1, (sourhi<<2) | (sourlo>>30));
1819 } /* [cannot drop through] */
1820 sourlo=DFWORD(df, 1); /* sourhi not involved now */
1821 if (sourlo&0xfff00000) { /* in one of first two */
1822 dpdlenchk(1, sourlo>>30); /* very rare */
1823 dpdlendun(2, sourlo>>20);
1824 } /* [cannot drop through] */
1825 dpdlenchk(3, sourlo>>10);
1826 dpdlendun(4, sourlo);
1827 /* [cannot drop through] */
1828
1829 #elif QUAD
1830 if (sourhi&0x00003fff) { /* ends in first */
1831 dpdlenchk(0, sourhi>>4);
1832 sourmh=DFWORD(df, 1);
1833 dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
1834 } /* [cannot drop through] */
1835 sourmh=DFWORD(df, 1);
1836 if (sourmh) {
1837 dpdlenchk(1, sourmh>>26);
1838 dpdlenchk(2, sourmh>>16);
1839 dpdlenchk(3, sourmh>>6);
1840 sourml=DFWORD(df, 2);
1841 dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
1842 } /* [cannot drop through] */
1843 sourml=DFWORD(df, 2);
1844 if (sourml) {
1845 dpdlenchk(4, sourml>>28);
1846 dpdlenchk(5, sourml>>18);
1847 dpdlenchk(6, sourml>>8);
1848 sourlo=DFWORD(df, 3);
1849 dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
1850 } /* [cannot drop through] */
1851 sourlo=DFWORD(df, 3);
1852 if (sourlo&0xfff00000) { /* in one of first two */
1853 dpdlenchk(7, sourlo>>30); /* very rare */
1854 dpdlendun(8, sourlo>>20);
1855 } /* [cannot drop through] */
1856 dpdlenchk(9, sourlo>>10);
1857 dpdlendun(10, sourlo);
1858 /* [cannot drop through] */
1859 #endif
1860 } /* decFloatDigits */
1861
1862 /* ------------------------------------------------------------------ */
1863 /* decFloatDivide -- divide a decFloat by another */
1864 /* */
1865 /* result gets the result of dividing dfl by dfr: */
1866 /* dfl is the first decFloat (lhs) */
1867 /* dfr is the second decFloat (rhs) */
1868 /* set is the context */
1869 /* returns result */
1870 /* */
1871 /* ------------------------------------------------------------------ */
1872 /* This is just a wrapper. */
1873 decFloat * decFloatDivide(decFloat *result,
1874 const decFloat *dfl, const decFloat *dfr,
1875 decContext *set) {
1876 return decDivide(result, dfl, dfr, set, DIVIDE);
1877 } /* decFloatDivide */
1878
1879 /* ------------------------------------------------------------------ */
1880 /* decFloatDivideInteger -- integer divide a decFloat by another */
1881 /* */
1882 /* result gets the result of dividing dfl by dfr: */
1883 /* dfl is the first decFloat (lhs) */
1884 /* dfr is the second decFloat (rhs) */
1885 /* set is the context */
1886 /* returns result */
1887 /* */
1888 /* ------------------------------------------------------------------ */
1889 decFloat * decFloatDivideInteger(decFloat *result,
1890 const decFloat *dfl, const decFloat *dfr,
1891 decContext *set) {
1892 return decDivide(result, dfl, dfr, set, DIVIDEINT);
1893 } /* decFloatDivideInteger */
1894
1895 /* ------------------------------------------------------------------ */
1896 /* decFloatFMA -- multiply and add three decFloats, fused */
1897 /* */
1898 /* result gets the result of (dfl*dfr)+dff with a single rounding */
1899 /* dfl is the first decFloat (lhs) */
1900 /* dfr is the second decFloat (rhs) */
1901 /* dff is the final decFloat (fhs) */
1902 /* set is the context */
1903 /* returns result */
1904 /* */
1905 /* ------------------------------------------------------------------ */
1906 decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
1907 const decFloat *dfr, const decFloat *dff,
1908 decContext *set) {
1909 /* The accumulator has the bytes needed for FiniteMultiply, plus */
1910 /* one byte to the left in case of carry, plus DECPMAX+2 to the */
1911 /* right for the final addition (up to full fhs + round & sticky) */
1912 #define FMALEN (1+ (DECPMAX9*18) +DECPMAX+2)
1913 uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
1914 /* .. and for final result */
1915 bcdnum mul; /* for multiplication result */
1916 bcdnum fin; /* for final operand, expanded */
1917 uByte coe[DECPMAX]; /* dff coefficient in BCD */
1918 bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
1919 uInt diffsign; /* non-zero if signs differ */
1920 uInt hipad; /* pad digit for hi if needed */
1921 Int padding; /* excess exponent */
1922 uInt carry; /* +1 for ten's complement and during add */
1923 uByte *ub, *uh, *ul; /* work */
1924
1925 /* handle all the special values [any special operand leads to a */
1926 /* special result] */
1927 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
1928 decFloat proxy; /* multiplication result proxy */
1929 /* NaNs are handled as usual, giving priority to sNaNs */
1930 if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1931 if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
1932 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1933 if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
1934 /* One or more of the three is infinite */
1935 /* infinity times zero is bad */
1936 decFloatZero(&proxy);
1937 if (DFISINF(dfl)) {
1938 if (DFISZERO(dfr)) return decInvalid(result, set);
1939 decInfinity(&proxy, &proxy);
1940 }
1941 else if (DFISINF(dfr)) {
1942 if (DFISZERO(dfl)) return decInvalid(result, set);
1943 decInfinity(&proxy, &proxy);
1944 }
1945 /* compute sign of multiplication and place in proxy */
1946 DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
1947 if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
1948 /* dff is Infinite */
1949 if (!DFISINF(&proxy)) return decInfinity(result, dff);
1950 /* both sides of addition are infinite; different sign is bad */
1951 if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
1952 return decInvalid(result, set);
1953 return decFloatCopy(result, &proxy);
1954 }
1955
1956 /* Here when all operands are finite */
1957
1958 /* First multiply dfl*dfr */
1959 decFiniteMultiply(&mul, acc+1, dfl, dfr);
1960 /* The multiply is complete, exact and unbounded, and described in */
1961 /* mul with the coefficient held in acc[1...] */
1962
1963 /* now add in dff; the algorithm is essentially the same as */
1964 /* decFloatAdd, but the code is different because the code there */
1965 /* is highly optimized for adding two numbers of the same size */
1966 fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
1967 fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
1968 diffsign=mul.sign^fin.sign; /* note if signs differ */
1969 fin.msd=coe;
1970 fin.lsd=coe+DECPMAX-1;
1971 GETCOEFF(dff, coe); /* extract the coefficient */
1972
1973 /* now set hi and lo so that hi points to whichever of mul and fin */
1974 /* has the higher exponent and lo point to the other [don't care if */
1975 /* the same] */
1976 if (mul.exponent>=fin.exponent) {
1977 hi=&mul;
1978 lo=&fin;
1979 }
1980 else {
1981 hi=&fin;
1982 lo=&mul;
1983 }
1984
1985 /* remove leading zeros on both operands; this will save time later */
1986 /* and make testing for zero trivial */
1987 for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
1988 for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
1989 for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
1990 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
1991
1992 /* if hi is zero then result will be lo (which has the smaller */
1993 /* exponent), which also may need to be tested for zero for the */
1994 /* weird IEEE 754 sign rules */
1995 if (*hi->msd==0 && hi->msd==hi->lsd) { /* hi is zero */
1996 /* "When the sum of two operands with opposite signs is */
1997 /* exactly zero, the sign of that sum shall be '+' in all */
1998 /* rounding modes except round toward -Infinity, in which */
1999 /* mode that sign shall be '-'." */
2000 if (diffsign) {
2001 if (*lo->msd==0 && lo->msd==lo->lsd) { /* lo is zero */
2002 lo->sign=0;
2003 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2004 } /* diffsign && lo=0 */
2005 } /* diffsign */
2006 return decFinalize(result, lo, set); /* may need clamping */
2007 } /* numfl is zero */
2008 /* [here, both are minimal length and hi is non-zero] */
2009
2010 /* if signs differ, take the ten's complement of hi (zeros to the */
2011 /* right do not matter because the complement of zero is zero); */
2012 /* the +1 is done later, as part of the addition, inserted at the */
2013 /* correct digit */
2014 hipad=0;
2015 carry=0;
2016 if (diffsign) {
2017 hipad=9;
2018 carry=1;
2019 /* exactly the correct number of digits must be inverted */
2020 for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UINTAT(uh)=0x09090909-UINTAT(uh);
2021 for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
2022 }
2023
2024 /* ready to add; note that hi has no leading zeros so gap */
2025 /* calculation does not have to be as pessimistic as in decFloatAdd */
2026 /* (this is much more like the arbitrary-precision algorithm in */
2027 /* Rexx and decNumber) */
2028
2029 /* padding is the number of zeros that would need to be added to hi */
2030 /* for its lsd to be aligned with the lsd of lo */
2031 padding=hi->exponent-lo->exponent;
2032 /* printf("FMA pad %ld\n", (LI)padding); */
2033
2034 /* the result of the addition will be built into the accumulator, */
2035 /* starting from the far right; this could be either hi or lo */
2036 ub=acc+FMALEN-1; /* where lsd of result will go */
2037 ul=lo->lsd; /* lsd of rhs */
2038
2039 if (padding!=0) { /* unaligned */
2040 /* if the msd of lo is more than DECPMAX+2 digits to the right of */
2041 /* the original msd of hi then it can be reduced to a single */
2042 /* digit at the right place, as it stays clear of hi digits */
2043 /* [it must be DECPMAX+2 because during a subtraction the msd */
2044 /* could become 0 after a borrow from 1.000 to 0.9999...] */
2045 Int hilen=(Int)(hi->lsd-hi->msd+1); /* lengths */
2046 Int lolen=(Int)(lo->lsd-lo->msd+1); /* .. */
2047 Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
2048 Int reduce=newexp-lo->exponent;
2049 if (reduce>0) { /* [= case gives reduce=0 nop] */
2050 /* printf("FMA reduce: %ld\n", (LI)reduce); */
2051 if (reduce>=lolen) { /* eating all */
2052 lo->lsd=lo->msd; /* reduce to single digit */
2053 lo->exponent=newexp; /* [known to be non-zero] */
2054 }
2055 else { /* < */
2056 uByte *up=lo->lsd;
2057 lo->lsd=lo->lsd-reduce;
2058 if (*lo->lsd==0) /* could need sticky bit */
2059 for (; up>lo->lsd; up--) { /* search discarded digits */
2060 if (*up!=0) { /* found one... */
2061 *lo->lsd=1; /* set sticky bit */
2062 break;
2063 }
2064 }
2065 lo->exponent+=reduce;
2066 }
2067 padding=hi->exponent-lo->exponent; /* recalculate */
2068 ul=lo->lsd; /* .. */
2069 } /* maybe reduce */
2070 /* padding is now <= DECPMAX+2 but still > 0; tricky DOUBLE case */
2071 /* is when hi is a 1 that will become a 0.9999... by subtraction: */
2072 /* hi: 1 E+16 */
2073 /* lo: .................1000000000000000 E-16 */
2074 /* which for the addition pads and reduces to: */
2075 /* hi: 1000000000000000000 E-2 */
2076 /* lo: .................1 E-2 */
2077 #if DECCHECK
2078 if (padding>DECPMAX+2) printf("FMA excess padding: %ld\n", (LI)padding);
2079 if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
2080 /* printf("FMA padding: %ld\n", (LI)padding); */
2081 #endif
2082 /* padding digits can now be set in the result; one or more of */
2083 /* these will come from lo; others will be zeros in the gap */
2084 for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
2085 for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
2086 }
2087
2088 /* addition now complete to the right of the rightmost digit of hi */
2089 uh=hi->lsd;
2090
2091 /* carry was set up depending on ten's complement above; do the add... */
2092 for (;; ub--) {
2093 uInt hid, lod;
2094 if (uh<hi->msd) {
2095 if (ul<lo->msd) break;
2096 hid=hipad;
2097 }
2098 else hid=*uh--;
2099 if (ul<lo->msd) lod=0;
2100 else lod=*ul--;
2101 *ub=(uByte)(carry+hid+lod);
2102 if (*ub<10) carry=0;
2103 else {
2104 *ub-=10;
2105 carry=1;
2106 }
2107 } /* addition loop */
2108
2109 /* addition complete -- now handle carry, borrow, etc. */
2110 /* use lo to set up the num (its exponent is already correct, and */
2111 /* sign usually is) */
2112 lo->msd=ub+1;
2113 lo->lsd=acc+FMALEN-1;
2114 /* decShowNum(lo, "lo"); */
2115 if (!diffsign) { /* same-sign addition */
2116 if (carry) { /* carry out */
2117 *ub=1; /* place the 1 .. */
2118 lo->msd--; /* .. and update */
2119 }
2120 } /* same sign */
2121 else { /* signs differed (subtraction) */
2122 if (!carry) { /* no carry out means hi<lo */
2123 /* borrowed -- take ten's complement of the right digits */
2124 lo->sign=hi->sign; /* sign is lhs sign */
2125 for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UINTAT(ul)=0x09090909-UINTAT(ul);
2126 for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
2127 /* complete the ten's complement by adding 1 [cannot overrun] */
2128 for (ul--; *ul==9; ul--) *ul=0;
2129 *ul+=1;
2130 } /* borrowed */
2131 else { /* carry out means hi>=lo */
2132 /* sign to use is lo->sign */
2133 /* all done except for the special IEEE 754 exact-zero-result */
2134 /* rule (see above); while testing for zero, strip leading */
2135 /* zeros (which will save decFinalize doing it) */
2136 for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2137 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2138 if (*lo->msd==0) { /* must be true zero (and diffsign) */
2139 lo->sign=0; /* assume + */
2140 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2141 }
2142 /* [else was not zero, might still have leading zeros] */
2143 } /* subtraction gave positive result */
2144 } /* diffsign */
2145
2146 return decFinalize(result, lo, set); /* round, check, and lay out */
2147 } /* decFloatFMA */
2148
2149 /* ------------------------------------------------------------------ */
2150 /* decFloatFromInt -- initialise a decFloat from an Int */
2151 /* */
2152 /* result gets the converted Int */
2153 /* n is the Int to convert */
2154 /* returns result */
2155 /* */
2156 /* The result is Exact; no errors or exceptions are possible. */
2157 /* ------------------------------------------------------------------ */
2158 decFloat * decFloatFromInt32(decFloat *result, Int n) {
2159 uInt u=(uInt)n; /* copy as bits */
2160 uInt encode; /* work */
2161 DFWORD(result, 0)=ZEROWORD; /* always */
2162 #if QUAD
2163 DFWORD(result, 1)=0;
2164 DFWORD(result, 2)=0;
2165 #endif
2166 if (n<0) { /* handle -n with care */
2167 /* [This can be done without the test, but is then slightly slower] */
2168 u=(~u)+1;
2169 DFWORD(result, 0)|=DECFLOAT_Sign;
2170 }
2171 /* Since the maximum value of u now is 2**31, only the low word of */
2172 /* result is affected */
2173 encode=BIN2DPD[u%1000];
2174 u/=1000;
2175 encode|=BIN2DPD[u%1000]<<10;
2176 u/=1000;
2177 encode|=BIN2DPD[u%1000]<<20;
2178 u/=1000; /* now 0, 1, or 2 */
2179 encode|=u<<30;
2180 DFWORD(result, DECWORDS-1)=encode;
2181 return result;
2182 } /* decFloatFromInt32 */
2183
2184 /* ------------------------------------------------------------------ */
2185 /* decFloatFromUInt -- initialise a decFloat from a uInt */
2186 /* */
2187 /* result gets the converted uInt */
2188 /* n is the uInt to convert */
2189 /* returns result */
2190 /* */
2191 /* The result is Exact; no errors or exceptions are possible. */
2192 /* ------------------------------------------------------------------ */
2193 decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
2194 uInt encode; /* work */
2195 DFWORD(result, 0)=ZEROWORD; /* always */
2196 #if QUAD
2197 DFWORD(result, 1)=0;
2198 DFWORD(result, 2)=0;
2199 #endif
2200 encode=BIN2DPD[u%1000];
2201 u/=1000;
2202 encode|=BIN2DPD[u%1000]<<10;
2203 u/=1000;
2204 encode|=BIN2DPD[u%1000]<<20;
2205 u/=1000; /* now 0 -> 4 */
2206 encode|=u<<30;
2207 DFWORD(result, DECWORDS-1)=encode;
2208 DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
2209 return result;
2210 } /* decFloatFromUInt32 */
2211
2212 /* ------------------------------------------------------------------ */
2213 /* decFloatInvert -- logical digitwise INVERT of a decFloat */
2214 /* */
2215 /* result gets the result of INVERTing df */
2216 /* df is the decFloat to invert */
2217 /* set is the context */
2218 /* returns result, which will be canonical with sign=0 */
2219 /* */
2220 /* The operand must be positive, finite with exponent q=0, and */
2221 /* comprise just zeros and ones; if not, Invalid operation results. */
2222 /* ------------------------------------------------------------------ */
2223 decFloat * decFloatInvert(decFloat *result, const decFloat *df,
2224 decContext *set) {
2225 uInt sourhi=DFWORD(df, 0); /* top word of dfs */
2226
2227 if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
2228 /* the operand is a finite integer (q=0) */
2229 #if DOUBLE
2230 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
2231 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
2232 #elif QUAD
2233 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
2234 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
2235 DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
2236 DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
2237 #endif
2238 return result;
2239 } /* decFloatInvert */
2240
2241 /* ------------------------------------------------------------------ */
2242 /* decFloatIs -- decFloat tests (IsSigned, etc.) */
2243 /* */
2244 /* df is the decFloat to test */
2245 /* returns 0 or 1 in an int32_t */
2246 /* */
2247 /* Many of these could be macros, but having them as real functions */
2248 /* is a bit cleaner (and they can be referred to here by the generic */
2249 /* names) */
2250 /* ------------------------------------------------------------------ */
2251 uInt decFloatIsCanonical(const decFloat *df) {
2252 if (DFISSPECIAL(df)) {
2253 if (DFISINF(df)) {
2254 if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
2255 if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
2256 return 1;
2257 }
2258 /* is a NaN */
2259 if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
2260 if (DFISCCZERO(df)) return 1; /* coefficient continuation */
2261 /* drop through to check payload */
2262 }
2263 { /* declare block */
2264 #if DOUBLE
2265 uInt sourhi=DFWORD(df, 0);
2266 uInt sourlo=DFWORD(df, 1);
2267 if (CANONDPDOFF(sourhi, 8)
2268 && CANONDPDTWO(sourhi, sourlo, 30)
2269 && CANONDPDOFF(sourlo, 20)
2270 && CANONDPDOFF(sourlo, 10)
2271 && CANONDPDOFF(sourlo, 0)) return 1;
2272 #elif QUAD
2273 uInt sourhi=DFWORD(df, 0);
2274 uInt sourmh=DFWORD(df, 1);
2275 uInt sourml=DFWORD(df, 2);
2276 uInt sourlo=DFWORD(df, 3);
2277 if (CANONDPDOFF(sourhi, 4)
2278 && CANONDPDTWO(sourhi, sourmh, 26)
2279 && CANONDPDOFF(sourmh, 16)
2280 && CANONDPDOFF(sourmh, 6)
2281 && CANONDPDTWO(sourmh, sourml, 28)
2282 && CANONDPDOFF(sourml, 18)
2283 && CANONDPDOFF(sourml, 8)
2284 && CANONDPDTWO(sourml, sourlo, 30)
2285 && CANONDPDOFF(sourlo, 20)
2286 && CANONDPDOFF(sourlo, 10)
2287 && CANONDPDOFF(sourlo, 0)) return 1;
2288 #endif
2289 } /* block */
2290 return 0; /* a declet is non-canonical */
2291 }
2292
2293 uInt decFloatIsFinite(const decFloat *df) {
2294 return !DFISSPECIAL(df);
2295 }
2296 uInt decFloatIsInfinite(const decFloat *df) {
2297 return DFISINF(df);
2298 }
2299 uInt decFloatIsInteger(const decFloat *df) {
2300 return DFISINT(df);
2301 }
2302 uInt decFloatIsNaN(const decFloat *df) {
2303 return DFISNAN(df);
2304 }
2305 uInt decFloatIsNormal(const decFloat *df) {
2306 Int exp; /* exponent */
2307 if (DFISSPECIAL(df)) return 0;
2308 if (DFISZERO(df)) return 0;
2309 /* is finite and non-zero */
2310 exp=GETEXPUN(df) /* get unbiased exponent .. */
2311 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2312 return (exp>=DECEMIN); /* < DECEMIN is subnormal */
2313 }
2314 uInt decFloatIsSignaling(const decFloat *df) {
2315 return DFISSNAN(df);
2316 }
2317 uInt decFloatIsSignalling(const decFloat *df) {
2318 return DFISSNAN(df);
2319 }
2320 uInt decFloatIsSigned(const decFloat *df) {
2321 return DFISSIGNED(df);
2322 }
2323 uInt decFloatIsSubnormal(const decFloat *df) {
2324 if (DFISSPECIAL(df)) return 0;
2325 /* is finite */
2326 if (decFloatIsNormal(df)) return 0;
2327 /* it is <Nmin, but could be zero */
2328 if (DFISZERO(df)) return 0;
2329 return 1; /* is subnormal */
2330 }
2331 uInt decFloatIsZero(const decFloat *df) {
2332 return DFISZERO(df);
2333 } /* decFloatIs... */
2334
2335 /* ------------------------------------------------------------------ */
2336 /* decFloatLogB -- return adjusted exponent, by 754r rules */
2337 /* */
2338 /* result gets the adjusted exponent as an integer, or a NaN etc. */
2339 /* df is the decFloat to be examined */
2340 /* set is the context */
2341 /* returns result */
2342 /* */
2343 /* Notable cases: */
2344 /* A<0 -> Use |A| */
2345 /* A=0 -> -Infinity (Division by zero) */
2346 /* A=Infinite -> +Infinity (Exact) */
2347 /* A=1 exactly -> 0 (Exact) */
2348 /* NaNs are propagated as usual */
2349 /* ------------------------------------------------------------------ */
2350 decFloat * decFloatLogB(decFloat *result, const decFloat *df,
2351 decContext *set) {
2352 Int ae; /* adjusted exponent */
2353 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2354 if (DFISINF(df)) {
2355 DFWORD(result, 0)=0; /* need +ve */
2356 return decInfinity(result, result); /* canonical +Infinity */
2357 }
2358 if (DFISZERO(df)) {
2359 set->status|=DEC_Division_by_zero; /* as per 754r */
2360 DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
2361 return decInfinity(result, result); /* canonical -Infinity */
2362 }
2363 ae=GETEXPUN(df) /* get unbiased exponent .. */
2364 +decFloatDigits(df)-1; /* .. and make adjusted exponent */
2365 /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
2366 /* it is worth using a special case of decFloatFromInt32 */
2367 DFWORD(result, 0)=ZEROWORD; /* always */
2368 if (ae<0) {
2369 DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
2370 ae=-ae;
2371 }
2372 #if DOUBLE
2373 DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
2374 #elif QUAD
2375 DFWORD(result, 1)=0;
2376 DFWORD(result, 2)=0;
2377 DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
2378 DFWORD(result, 3)|=BIN2DPD[ae%1000];
2379 #endif
2380 return result;
2381 } /* decFloatLogB */
2382
2383 /* ------------------------------------------------------------------ */
2384 /* decFloatMax -- return maxnum of two operands */
2385 /* */
2386 /* result gets the chosen decFloat */
2387 /* dfl is the first decFloat (lhs) */
2388 /* dfr is the second decFloat (rhs) */
2389 /* set is the context */
2390 /* returns result */
2391 /* */
2392 /* If just one operand is a quiet NaN it is ignored. */
2393 /* ------------------------------------------------------------------ */
2394 decFloat * decFloatMax(decFloat *result,
2395 const decFloat *dfl, const decFloat *dfr,
2396 decContext *set) {
2397 Int comp;
2398 if (DFISNAN(dfl)) {
2399 /* sNaN or both NaNs leads to normal NaN processing */
2400 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2401 return decCanonical(result, dfr); /* RHS is numeric */
2402 }
2403 if (DFISNAN(dfr)) {
2404 /* sNaN leads to normal NaN processing (both NaN handled above) */
2405 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2406 return decCanonical(result, dfl); /* LHS is numeric */
2407 }
2408 /* Both operands are numeric; numeric comparison needed -- use */
2409 /* total order for a well-defined choice (and +0 > -0) */
2410 comp=decNumCompare(dfl, dfr, 1);
2411 if (comp>=0) return decCanonical(result, dfl);
2412 return decCanonical(result, dfr);
2413 } /* decFloatMax */
2414
2415 /* ------------------------------------------------------------------ */
2416 /* decFloatMaxMag -- return maxnummag of two operands */
2417 /* */
2418 /* result gets the chosen decFloat */
2419 /* dfl is the first decFloat (lhs) */
2420 /* dfr is the second decFloat (rhs) */
2421 /* set is the context */
2422 /* returns result */
2423 /* */
2424 /* Returns according to the magnitude comparisons if both numeric and */
2425 /* unequal, otherwise returns maxnum */
2426 /* ------------------------------------------------------------------ */
2427 decFloat * decFloatMaxMag(decFloat *result,
2428 const decFloat *dfl, const decFloat *dfr,
2429 decContext *set) {
2430 Int comp;
2431 decFloat absl, absr;
2432 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
2433
2434 decFloatCopyAbs(&absl, dfl);
2435 decFloatCopyAbs(&absr, dfr);
2436 comp=decNumCompare(&absl, &absr, 0);
2437 if (comp>0) return decCanonical(result, dfl);
2438 if (comp<0) return decCanonical(result, dfr);
2439 return decFloatMax(result, dfl, dfr, set);
2440 } /* decFloatMaxMag */
2441
2442 /* ------------------------------------------------------------------ */
2443 /* decFloatMin -- return minnum of two operands */
2444 /* */
2445 /* result gets the chosen decFloat */
2446 /* dfl is the first decFloat (lhs) */
2447 /* dfr is the second decFloat (rhs) */
2448 /* set is the context */
2449 /* returns result */
2450 /* */
2451 /* If just one operand is a quiet NaN it is ignored. */
2452 /* ------------------------------------------------------------------ */
2453 decFloat * decFloatMin(decFloat *result,
2454 const decFloat *dfl, const decFloat *dfr,
2455 decContext *set) {
2456 Int comp;
2457 if (DFISNAN(dfl)) {
2458 /* sNaN or both NaNs leads to normal NaN processing */
2459 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2460 return decCanonical(result, dfr); /* RHS is numeric */
2461 }
2462 if (DFISNAN(dfr)) {
2463 /* sNaN leads to normal NaN processing (both NaN handled above) */
2464 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2465 return decCanonical(result, dfl); /* LHS is numeric */
2466 }
2467 /* Both operands are numeric; numeric comparison needed -- use */
2468 /* total order for a well-defined choice (and +0 > -0) */
2469 comp=decNumCompare(dfl, dfr, 1);
2470 if (comp<=0) return decCanonical(result, dfl);
2471 return decCanonical(result, dfr);
2472 } /* decFloatMin */
2473
2474 /* ------------------------------------------------------------------ */
2475 /* decFloatMinMag -- return minnummag of two operands */
2476 /* */
2477 /* result gets the chosen decFloat */
2478 /* dfl is the first decFloat (lhs) */
2479 /* dfr is the second decFloat (rhs) */
2480 /* set is the context */
2481 /* returns result */
2482 /* */
2483 /* Returns according to the magnitude comparisons if both numeric and */
2484 /* unequal, otherwise returns minnum */
2485 /* ------------------------------------------------------------------ */
2486 decFloat * decFloatMinMag(decFloat *result,
2487 const decFloat *dfl, const decFloat *dfr,
2488 decContext *set) {
2489 Int comp;
2490 decFloat absl, absr;
2491 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
2492
2493 decFloatCopyAbs(&absl, dfl);
2494 decFloatCopyAbs(&absr, dfr);
2495 comp=decNumCompare(&absl, &absr, 0);
2496 if (comp<0) return decCanonical(result, dfl);
2497 if (comp>0) return decCanonical(result, dfr);
2498 return decFloatMin(result, dfl, dfr, set);
2499 } /* decFloatMinMag */
2500
2501 /* ------------------------------------------------------------------ */
2502 /* decFloatMinus -- negate value, heeding NaNs, etc. */
2503 /* */
2504 /* result gets the canonicalized 0-df */
2505 /* df is the decFloat to minus */
2506 /* set is the context */
2507 /* returns result */
2508 /* */
2509 /* This has the same effect as 0-df where the exponent of the zero is */
2510 /* the same as that of df (if df is finite). */
2511 /* The effect is also the same as decFloatCopyNegate except that NaNs */
2512 /* are handled normally (the sign of a NaN is not affected, and an */
2513 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2514 /* ------------------------------------------------------------------ */
2515 decFloat * decFloatMinus(decFloat *result, const decFloat *df,
2516 decContext *set) {
2517 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2518 decCanonical(result, df); /* copy and check */
2519 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2520 else DFBYTE(result, 0)^=0x80; /* flip sign bit */
2521 return result;
2522 } /* decFloatMinus */
2523
2524 /* ------------------------------------------------------------------ */
2525 /* decFloatMultiply -- multiply two decFloats */
2526 /* */
2527 /* result gets the result of multiplying dfl and dfr: */
2528 /* dfl is the first decFloat (lhs) */
2529 /* dfr is the second decFloat (rhs) */
2530 /* set is the context */
2531 /* returns result */
2532 /* */
2533 /* ------------------------------------------------------------------ */
2534 decFloat * decFloatMultiply(decFloat *result,
2535 const decFloat *dfl, const decFloat *dfr,
2536 decContext *set) {
2537 bcdnum num; /* for final conversion */
2538 uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
2539
2540 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
2541 /* NaNs are handled as usual */
2542 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2543 /* infinity times zero is bad */
2544 if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
2545 if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
2546 /* both infinite; return canonical infinity with computed sign */
2547 DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
2548 return decInfinity(result, result);
2549 }
2550
2551 /* Here when both operands are finite */
2552 decFiniteMultiply(&num, bcdacc, dfl, dfr);
2553 return decFinalize(result, &num, set); /* round, check, and lay out */
2554 } /* decFloatMultiply */
2555
2556 /* ------------------------------------------------------------------ */
2557 /* decFloatNextMinus -- next towards -Infinity */
2558 /* */
2559 /* result gets the next lesser decFloat */
2560 /* dfl is the decFloat to start with */
2561 /* set is the context */
2562 /* returns result */
2563 /* */
2564 /* This is 754r nextdown; Invalid is the only status possible (from */
2565 /* an sNaN). */
2566 /* ------------------------------------------------------------------ */
2567 decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
2568 decContext *set) {
2569 decFloat delta; /* tiny increment */
2570 uInt savestat; /* saves status */
2571 enum rounding saveround; /* .. and mode */
2572
2573 /* +Infinity is the special case */
2574 if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
2575 DFSETNMAX(result);
2576 return result; /* [no status to set] */
2577 }
2578 /* other cases are effected by sutracting a tiny delta -- this */
2579 /* should be done in a wider format as the delta is unrepresentable */
2580 /* here (but can be done with normal add if the sign of zero is */
2581 /* treated carefully, because no Inexactitude is interesting); */
2582 /* rounding to -Infinity then pushes the result to next below */
2583 decFloatZero(&delta); /* set up tiny delta */
2584 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2585 DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
2586 /* set up for the directional round */
2587 saveround=set->round; /* save mode */
2588 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2589 savestat=set->status; /* save status */
2590 decFloatAdd(result, dfl, &delta, set);
2591 /* Add rules mess up the sign when going from +Ntiny to 0 */
2592 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2593 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2594 set->status|=savestat; /* restore pending flags */
2595 set->round=saveround; /* .. and mode */
2596 return result;
2597 } /* decFloatNextMinus */
2598
2599 /* ------------------------------------------------------------------ */
2600 /* decFloatNextPlus -- next towards +Infinity */
2601 /* */
2602 /* result gets the next larger decFloat */
2603 /* dfl is the decFloat to start with */
2604 /* set is the context */
2605 /* returns result */
2606 /* */
2607 /* This is 754r nextup; Invalid is the only status possible (from */
2608 /* an sNaN). */
2609 /* ------------------------------------------------------------------ */
2610 decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
2611 decContext *set) {
2612 uInt savestat; /* saves status */
2613 enum rounding saveround; /* .. and mode */
2614 decFloat delta; /* tiny increment */
2615
2616 /* -Infinity is the special case */
2617 if (DFISINF(dfl) && DFISSIGNED(dfl)) {
2618 DFSETNMAX(result);
2619 DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
2620 return result; /* [no status to set] */
2621 }
2622 /* other cases are effected by sutracting a tiny delta -- this */
2623 /* should be done in a wider format as the delta is unrepresentable */
2624 /* here (but can be done with normal add if the sign of zero is */
2625 /* treated carefully, because no Inexactitude is interesting); */
2626 /* rounding to +Infinity then pushes the result to next above */
2627 decFloatZero(&delta); /* set up tiny delta */
2628 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2629 DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
2630 /* set up for the directional round */
2631 saveround=set->round; /* save mode */
2632 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2633 savestat=set->status; /* save status */
2634 decFloatAdd(result, dfl, &delta, set);
2635 /* Add rules mess up the sign when going from -Ntiny to -0 */
2636 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2637 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
2638 set->status|=savestat; /* restore pending flags */
2639 set->round=saveround; /* .. and mode */
2640 return result;
2641 } /* decFloatNextPlus */
2642
2643 /* ------------------------------------------------------------------ */
2644 /* decFloatNextToward -- next towards a decFloat */
2645 /* */
2646 /* result gets the next decFloat */
2647 /* dfl is the decFloat to start with */
2648 /* dfr is the decFloat to move toward */
2649 /* set is the context */
2650 /* returns result */
2651 /* */
2652 /* This is 754r nextafter; status may be set unless the result is a */
2653 /* normal number. */
2654 /* ------------------------------------------------------------------ */
2655 decFloat * decFloatNextToward(decFloat *result,
2656 const decFloat *dfl, const decFloat *dfr,
2657 decContext *set) {
2658 decFloat delta; /* tiny increment or decrement */
2659 decFloat pointone; /* 1e-1 */
2660 uInt savestat; /* saves status */
2661 enum rounding saveround; /* .. and mode */
2662 uInt deltatop; /* top word for delta */
2663 Int comp; /* work */
2664
2665 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2666 /* Both are numeric, so Invalid no longer a possibility */
2667 comp=decNumCompare(dfl, dfr, 0);
2668 if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
2669 /* unequal; do NextPlus or NextMinus but with different status rules */
2670
2671 if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
2672 if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
2673 DFSETNMAX(result);
2674 DFWORD(result, 0)|=DECFLOAT_Sign;
2675 return result;
2676 }
2677 saveround=set->round; /* save mode */
2678 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
2679 deltatop=0; /* positive delta */
2680 }
2681 else { /* lhs>rhs, do NextMinus, see above for commentary */
2682 if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
2683 DFSETNMAX(result);
2684 return result;
2685 }
2686 saveround=set->round; /* save mode */
2687 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
2688 deltatop=DECFLOAT_Sign; /* negative delta */
2689 }
2690 savestat=set->status; /* save status */
2691 /* Here, Inexact is needed where appropriate (and hence Underflow, */
2692 /* etc.). Therefore the tiny delta which is otherwise */
2693 /* unrepresentable (see NextPlus and NextMinus) is constructed */
2694 /* using the multiplication of FMA. */
2695 decFloatZero(&delta); /* set up tiny delta */
2696 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
2697 DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
2698 decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
2699 decFloatFMA(result, &delta, &pointone, dfl, set);
2700 /* [Delta is truly tiny, so no need to correct sign of zero] */
2701 /* use new status unless the result is normal */
2702 if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
2703 set->round=saveround; /* restore mode */
2704 return result;
2705 } /* decFloatNextToward */
2706
2707 /* ------------------------------------------------------------------ */
2708 /* decFloatOr -- logical digitwise OR of two decFloats */
2709 /* */
2710 /* result gets the result of ORing dfl and dfr */
2711 /* dfl is the first decFloat (lhs) */
2712 /* dfr is the second decFloat (rhs) */
2713 /* set is the context */
2714 /* returns result, which will be canonical with sign=0 */
2715 /* */
2716 /* The operands must be positive, finite with exponent q=0, and */
2717 /* comprise just zeros and ones; if not, Invalid operation results. */
2718 /* ------------------------------------------------------------------ */
2719 decFloat * decFloatOr(decFloat *result,
2720 const decFloat *dfl, const decFloat *dfr,
2721 decContext *set) {
2722 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
2723 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
2724 /* the operands are positive finite integers (q=0) with just 0s and 1s */
2725 #if DOUBLE
2726 DFWORD(result, 0)=ZEROWORD
2727 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
2728 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
2729 #elif QUAD
2730 DFWORD(result, 0)=ZEROWORD
2731 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
2732 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
2733 DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
2734 DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
2735 #endif
2736 return result;
2737 } /* decFloatOr */
2738
2739 /* ------------------------------------------------------------------ */
2740 /* decFloatPlus -- add value to 0, heeding NaNs, etc. */
2741 /* */
2742 /* result gets the canonicalized 0+df */
2743 /* df is the decFloat to plus */
2744 /* set is the context */
2745 /* returns result */
2746 /* */
2747 /* This has the same effect as 0+df where the exponent of the zero is */
2748 /* the same as that of df (if df is finite). */
2749 /* The effect is also the same as decFloatCopy except that NaNs */
2750 /* are handled normally (the sign of a NaN is not affected, and an */
2751 /* sNaN will signal), the result is canonical, and zero gets sign 0. */
2752 /* ------------------------------------------------------------------ */
2753 decFloat * decFloatPlus(decFloat *result, const decFloat *df,
2754 decContext *set) {
2755 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2756 decCanonical(result, df); /* copy and check */
2757 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
2758 return result;
2759 } /* decFloatPlus */
2760
2761 /* ------------------------------------------------------------------ */
2762 /* decFloatQuantize -- quantize a decFloat */
2763 /* */
2764 /* result gets the result of quantizing dfl to match dfr */
2765 /* dfl is the first decFloat (lhs) */
2766 /* dfr is the second decFloat (rhs), which sets the exponent */
2767 /* set is the context */
2768 /* returns result */
2769 /* */
2770 /* Unless there is an error or the result is infinite, the exponent */
2771 /* of result is guaranteed to be the same as that of dfr. */
2772 /* ------------------------------------------------------------------ */
2773 decFloat * decFloatQuantize(decFloat *result,
2774 const decFloat *dfl, const decFloat *dfr,
2775 decContext *set) {
2776 Int explb, exprb; /* left and right biased exponents */
2777 uByte *ulsd; /* local LSD pointer */
2778 uInt *ui; /* work */
2779 uByte *ub; /* .. */
2780 Int drop; /* .. */
2781 uInt dpd; /* .. */
2782 uInt encode; /* encoding accumulator */
2783 uInt sourhil, sourhir; /* top words from source decFloats */
2784 /* the following buffer holds the coefficient for manipulation */
2785 uByte buf[4+DECPMAX*3]; /* + space for zeros to left or right */
2786 #if DECTRACE
2787 bcdnum num; /* for trace displays */
2788 #endif
2789
2790 /* Start decoding the arguments */
2791 sourhil=DFWORD(dfl, 0); /* LHS top word */
2792 explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
2793 sourhir=DFWORD(dfr, 0); /* RHS top word */
2794 exprb=DECCOMBEXP[sourhir>>26];
2795
2796 if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
2797 /* NaNs are handled as usual */
2798 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2799 /* one infinity but not both is bad */
2800 if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
2801 /* both infinite; return canonical infinity with sign of LHS */
2802 return decInfinity(result, dfl);
2803 }
2804
2805 /* Here when both arguments are finite */
2806 /* complete extraction of the exponents [no need to unbias] */
2807 explb+=GETECON(dfl); /* + continuation */
2808 exprb+=GETECON(dfr); /* .. */
2809
2810 /* calculate the number of digits to drop from the coefficient */
2811 drop=exprb-explb; /* 0 if nothing to do */
2812 if (drop==0) return decCanonical(result, dfl); /* return canonical */
2813
2814 /* the coefficient is needed; lay it out into buf, offset so zeros */
2815 /* can be added before or after as needed -- an extra heading is */
2816 /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
2817 /* fours */
2818 #define BUFOFF (buf+4+DECPMAX)
2819 GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
2820 /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
2821
2822 #if DECTRACE
2823 num.msd=BUFOFF;
2824 num.lsd=BUFOFF+DECPMAX-1;
2825 num.exponent=explb-DECBIAS;
2826 num.sign=sourhil & DECFLOAT_Sign;
2827 decShowNum(&num, "dfl");
2828 #endif
2829
2830 if (drop>0) { /* [most common case] */
2831 /* (this code is very similar to that in decFloatFinalize, but */
2832 /* has many differences so is duplicated here -- so any changes */
2833 /* may need to be made there, too) */
2834 uByte *roundat; /* -> re-round digit */
2835 uByte reround; /* reround value */
2836 /* printf("Rounding; drop=%ld\n", (LI)drop); */
2837
2838 /* there is at least one zero needed to the left, in all but one */
2839 /* exceptional (all-nines) case, so place four zeros now; this is */
2840 /* needed almost always and makes rounding all-nines by fours safe */
2841 UINTAT(BUFOFF-4)=0;
2842
2843 /* Three cases here: */
2844 /* 1. new LSD is in coefficient (almost always) */
2845 /* 2. new LSD is digit to left of coefficient (so MSD is */
2846 /* round-for-reround digit) */
2847 /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
2848 /* Note that leading zeros can safely be treated as useful digits */
2849
2850 /* [duplicate check-stickies code to save a test] */
2851 /* [by-digit check for stickies as runs of zeros are rare] */
2852 if (drop<DECPMAX) { /* NB lengths not addresses */
2853 roundat=BUFOFF+DECPMAX-drop;
2854 reround=*roundat;
2855 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2856 if (*ub!=0) { /* non-zero to be discarded */
2857 reround=DECSTICKYTAB[reround]; /* apply sticky bit */
2858 break; /* [remainder don't-care] */
2859 }
2860 } /* check stickies */
2861 ulsd=roundat-1; /* set LSD */
2862 }
2863 else { /* edge case */
2864 if (drop==DECPMAX) {
2865 roundat=BUFOFF;
2866 reround=*roundat;
2867 }
2868 else {
2869 roundat=BUFOFF-1;
2870 reround=0;
2871 }
2872 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2873 if (*ub!=0) { /* non-zero to be discarded */
2874 reround=DECSTICKYTAB[reround]; /* apply sticky bit */
2875 break; /* [remainder don't-care] */
2876 }
2877 } /* check stickies */
2878 *BUFOFF=0; /* make a coefficient of 0 */
2879 ulsd=BUFOFF; /* .. at the MSD place */
2880 }
2881
2882 if (reround!=0) { /* discarding non-zero */
2883 uInt bump=0;
2884 set->status|=DEC_Inexact;
2885
2886 /* next decide whether to increment the coefficient */
2887 if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
2888 if (reround>5) bump=1; /* >0.5 goes up */
2889 else if (reround==5) /* exactly 0.5000 .. */
2890 bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */
2891 } /* r-h-e */
2892 else switch (set->round) {
2893 case DEC_ROUND_DOWN: {
2894 /* no change */
2895 break;} /* r-d */
2896 case DEC_ROUND_HALF_DOWN: {
2897 if (reround>5) bump=1;
2898 break;} /* r-h-d */
2899 case DEC_ROUND_HALF_UP: {
2900 if (reround>=5) bump=1;
2901 break;} /* r-h-u */
2902 case DEC_ROUND_UP: {
2903 if (reround>0) bump=1;
2904 break;} /* r-u */
2905 case DEC_ROUND_CEILING: {
2906 /* same as _UP for positive numbers, and as _DOWN for negatives */
2907 if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
2908 break;} /* r-c */
2909 case DEC_ROUND_FLOOR: {
2910 /* same as _UP for negative numbers, and as _DOWN for positive */
2911 /* [negative reround cannot occur on 0] */
2912 if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
2913 break;} /* r-f */
2914 case DEC_ROUND_05UP: {
2915 if (reround>0) { /* anything out there is 'sticky' */
2916 /* bump iff lsd=0 or 5; this cannot carry so it could be */
2917 /* effected immediately with no bump -- but the code */
2918 /* is clearer if this is done the same way as the others */
2919 if (*ulsd==0 || *ulsd==5) bump=1;
2920 }
2921 break;} /* r-r */
2922 default: { /* e.g., DEC_ROUND_MAX */
2923 set->status|=DEC_Invalid_context;
2924 #if DECCHECK
2925 printf("Unknown rounding mode: %ld\n", (LI)set->round);
2926 #endif
2927 break;}
2928 } /* switch (not r-h-e) */
2929 /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
2930
2931 if (bump!=0) { /* need increment */
2932 /* increment the coefficient; this could give 1000... (after */
2933 /* the all nines case) */
2934 ub=ulsd;
2935 for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
2936 /* now at most 3 digits left to non-9 (usually just the one) */
2937 for (; *ub==9; ub--) *ub=0;
2938 *ub+=1;
2939 /* [the all-nines case will have carried one digit to the */
2940 /* left of the original MSD -- just where it is needed] */
2941 } /* bump needed */
2942 } /* inexact rounding */
2943
2944 /* now clear zeros to the left so exactly DECPMAX digits will be */
2945 /* available in the coefficent -- the first word to the left was */
2946 /* cleared earlier for safe carry; now add any more needed */
2947 if (drop>4) {
2948 UINTAT(BUFOFF-8)=0; /* must be at least 5 */
2949 for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
2950 }
2951 } /* need round (drop>0) */
2952
2953 else { /* drop<0; padding with -drop digits is needed */
2954 /* This is the case where an error can occur if the padded */
2955 /* coefficient will not fit; checking for this can be done in the */
2956 /* same loop as padding for zeros if the no-hope and zero cases */
2957 /* are checked first */
2958 if (-drop>DECPMAX-1) { /* cannot fit unless 0 */
2959 if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
2960 /* a zero can have any exponent; just drop through and use it */
2961 ulsd=BUFOFF+DECPMAX-1;
2962 }
2963 else { /* padding will fit (but may still be too long) */
2964 /* final-word mask depends on endianess */
2965 #if DECLITEND
2966 static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
2967 #else
2968 static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
2969 #endif
2970 for (ui=&UINTAT(BUFOFF+DECPMAX);; ui++) {
2971 *ui=0;
2972 if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
2973 /* if all four digits should be zero, definitely bad */
2974 if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
2975 return decInvalid(result, set);
2976 /* must be a 1- to 3-digit sequence; check more carefully */
2977 if ((UINTAT(&UBYTEAT(ui)-DECPMAX)&dmask[(-drop)%4])!=0)
2978 return decInvalid(result, set);
2979 break; /* no need for loop end test */
2980 }
2981 if (ui>=&UINTAT(BUFOFF+DECPMAX+(-drop)-4)) break; /* done */
2982 }
2983 ulsd=BUFOFF+DECPMAX+(-drop)-1;
2984 } /* pad and check leading zeros */
2985 } /* drop<0 */
2986
2987 #if DECTRACE
2988 num.msd=ulsd-DECPMAX+1;
2989 num.lsd=ulsd;
2990 num.exponent=explb-DECBIAS;
2991 num.sign=sourhil & DECFLOAT_Sign;
2992 decShowNum(&num, "res");
2993 #endif
2994
2995 /*------------------------------------------------------------------*/
2996 /* At this point the result is DECPMAX digits, ending at ulsd, so */
2997 /* fits the encoding exactly; there is no possibility of error */
2998 /*------------------------------------------------------------------*/
2999 encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
3000 encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */
3001 /* the exponent continuation can be extracted from the original RHS */
3002 encode|=sourhir & ECONMASK;
3003 encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */
3004
3005 /* finally encode the coefficient */
3006 /* private macro to encode a declet; this version can be used */
3007 /* because all coefficient digits exist */
3008 #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
3009 dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
3010
3011 #if DOUBLE
3012 getDPD3q(dpd, 4); encode|=dpd<<8;
3013 getDPD3q(dpd, 3); encode|=dpd>>2;
3014 DFWORD(result, 0)=encode;
3015 encode=dpd<<30;
3016 getDPD3q(dpd, 2); encode|=dpd<<20;
3017 getDPD3q(dpd, 1); encode|=dpd<<10;
3018 getDPD3q(dpd, 0); encode|=dpd;
3019 DFWORD(result, 1)=encode;
3020
3021 #elif QUAD
3022 getDPD3q(dpd,10); encode|=dpd<<4;
3023 getDPD3q(dpd, 9); encode|=dpd>>6;
3024 DFWORD(result, 0)=encode;
3025 encode=dpd<<26;
3026 getDPD3q(dpd, 8); encode|=dpd<<16;
3027 getDPD3q(dpd, 7); encode|=dpd<<6;
3028 getDPD3q(dpd, 6); encode|=dpd>>4;
3029 DFWORD(result, 1)=encode;
3030 encode=dpd<<28;
3031 getDPD3q(dpd, 5); encode|=dpd<<18;
3032 getDPD3q(dpd, 4); encode|=dpd<<8;
3033 getDPD3q(dpd, 3); encode|=dpd>>2;
3034 DFWORD(result, 2)=encode;
3035 encode=dpd<<30;
3036 getDPD3q(dpd, 2); encode|=dpd<<20;
3037 getDPD3q(dpd, 1); encode|=dpd<<10;
3038 getDPD3q(dpd, 0); encode|=dpd;
3039 DFWORD(result, 3)=encode;
3040 #endif
3041 return result;
3042 } /* decFloatQuantize */
3043
3044 /* ------------------------------------------------------------------ */
3045 /* decFloatReduce -- reduce finite coefficient to minimum length */
3046 /* */
3047 /* result gets the reduced decFloat */
3048 /* df is the source decFloat */
3049 /* set is the context */
3050 /* returns result, which will be canonical */
3051 /* */
3052 /* This removes all possible trailing zeros from the coefficient; */
3053 /* some may remain when the number is very close to Nmax. */
3054 /* Special values are unchanged and no status is set unless df=sNaN. */
3055 /* Reduced zero has an exponent q=0. */
3056 /* ------------------------------------------------------------------ */
3057 decFloat * decFloatReduce(decFloat *result, const decFloat *df,
3058 decContext *set) {
3059 bcdnum num; /* work */
3060 uByte buf[DECPMAX], *ub; /* coefficient and pointer */
3061 if (df!=result) *result=*df; /* copy, if needed */
3062 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */
3063 /* zeros and infinites propagate too */
3064 if (DFISINF(df)) return decInfinity(result, df); /* canonical */
3065 if (DFISZERO(df)) {
3066 uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
3067 decFloatZero(result);
3068 DFWORD(result, 0)|=sign;
3069 return result; /* exponent dropped, sign OK */
3070 }
3071 /* non-zero finite */
3072 GETCOEFF(df, buf);
3073 ub=buf+DECPMAX-1; /* -> lsd */
3074 if (*ub) return result; /* no trailing zeros */
3075 for (ub--; *ub==0;) ub--; /* terminates because non-zero */
3076 /* *ub is the first non-zero from the right */
3077 num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
3078 num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
3079 num.msd=buf;
3080 num.lsd=ub;
3081 return decFinalize(result, &num, set);
3082 } /* decFloatReduce */
3083
3084 /* ------------------------------------------------------------------ */
3085 /* decFloatRemainder -- integer divide and return remainder */
3086 /* */
3087 /* result gets the remainder of dividing dfl by dfr: */
3088 /* dfl is the first decFloat (lhs) */
3089 /* dfr is the second decFloat (rhs) */
3090 /* set is the context */
3091 /* returns result */
3092 /* */
3093 /* ------------------------------------------------------------------ */
3094 decFloat * decFloatRemainder(decFloat *result,
3095 const decFloat *dfl, const decFloat *dfr,
3096 decContext *set) {
3097 return decDivide(result, dfl, dfr, set, REMAINDER);
3098 } /* decFloatRemainder */
3099
3100 /* ------------------------------------------------------------------ */
3101 /* decFloatRemainderNear -- integer divide to nearest and remainder */
3102 /* */
3103 /* result gets the remainder of dividing dfl by dfr: */
3104 /* dfl is the first decFloat (lhs) */
3105 /* dfr is the second decFloat (rhs) */
3106 /* set is the context */
3107 /* returns result */
3108 /* */
3109 /* This is the IEEE remainder, where the nearest integer is used. */
3110 /* ------------------------------------------------------------------ */
3111 decFloat * decFloatRemainderNear(decFloat *result,
3112 const decFloat *dfl, const decFloat *dfr,
3113 decContext *set) {
3114 return decDivide(result, dfl, dfr, set, REMNEAR);
3115 } /* decFloatRemainderNear */
3116
3117 /* ------------------------------------------------------------------ */
3118 /* decFloatRotate -- rotate the coefficient of a decFloat left/right */
3119 /* */
3120 /* result gets the result of rotating dfl */
3121 /* dfl is the source decFloat to rotate */
3122 /* dfr is the count of digits to rotate, an integer (with q=0) */
3123 /* set is the context */
3124 /* returns result */
3125 /* */
3126 /* The digits of the coefficient of dfl are rotated to the left (if */
3127 /* dfr is positive) or to the right (if dfr is negative) without */
3128 /* adjusting the exponent or the sign of dfl. */
3129 /* */
3130 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3131 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3132 /* dfr must be valid). No status is set unless dfr is invalid or an */
3133 /* operand is an sNaN. The result is canonical. */
3134 /* ------------------------------------------------------------------ */
3135 #define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
3136 decFloat * decFloatRotate(decFloat *result,
3137 const decFloat *dfl, const decFloat *dfr,
3138 decContext *set) {
3139 Int rotate; /* dfr as an Int */
3140 uByte buf[DECPMAX+PHALF]; /* coefficient + half */
3141 uInt digits, savestat; /* work */
3142 bcdnum num; /* .. */
3143 uByte *ub; /* .. */
3144
3145 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3146 if (!DFISINT(dfr)) return decInvalid(result, set);
3147 digits=decFloatDigits(dfr); /* calculate digits */
3148 if (digits>2) return decInvalid(result, set); /* definitely out of range */
3149 rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3150 if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
3151 /* [from here on no error or status change is possible] */
3152 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3153 /* handle no-rotate cases */
3154 if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
3155 /* a real rotate is needed: 0 < rotate < DECPMAX */
3156 /* reduce the rotation to no more than half to reduce copying later */
3157 /* (for QUAD in fact half + 2 digits) */
3158 if (DFISSIGNED(dfr)) rotate=-rotate;
3159 if (abs(rotate)>PHALF) {
3160 if (rotate<0) rotate=DECPMAX+rotate;
3161 else rotate=rotate-DECPMAX;
3162 }
3163 /* now lay out the coefficient, leaving room to the right or the */
3164 /* left depending on the direction of rotation */
3165 ub=buf;
3166 if (rotate<0) ub+=PHALF; /* rotate right, so space to left */
3167 GETCOEFF(dfl, ub);
3168 /* copy half the digits to left or right, and set num.msd */
3169 if (rotate<0) {
3170 memcpy(buf, buf+DECPMAX, PHALF);
3171 num.msd=buf+PHALF+rotate;
3172 }
3173 else {
3174 memcpy(buf+DECPMAX, buf, PHALF);
3175 num.msd=buf+rotate;
3176 }
3177 /* fill in rest of num */
3178 num.lsd=num.msd+DECPMAX-1;
3179 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3180 num.exponent=GETEXPUN(dfl);
3181 savestat=set->status; /* record */
3182 decFinalize(result, &num, set);
3183 set->status=savestat; /* restore */
3184 return result;
3185 } /* decFloatRotate */
3186
3187 /* ------------------------------------------------------------------ */
3188 /* decFloatSameQuantum -- test decFloats for same quantum */
3189 /* */
3190 /* dfl is the first decFloat (lhs) */
3191 /* dfr is the second decFloat (rhs) */
3192 /* returns 1 if the operands have the same quantum, 0 otherwise */
3193 /* */
3194 /* No error is possible and no status results. */
3195 /* ------------------------------------------------------------------ */
3196 uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
3197 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
3198 if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
3199 if (DFISINF(dfl) && DFISINF(dfr)) return 1;
3200 return 0; /* any other special mixture gives false */
3201 }
3202 if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
3203 return 0;
3204 } /* decFloatSameQuantum */
3205
3206 /* ------------------------------------------------------------------ */
3207 /* decFloatScaleB -- multiply by a power of 10, as per 754r */
3208 /* */
3209 /* result gets the result of the operation */
3210 /* dfl is the first decFloat (lhs) */
3211 /* dfr is the second decFloat (rhs), am integer (with q=0) */
3212 /* set is the context */
3213 /* returns result */
3214 /* */
3215 /* This computes result=dfl x 10**dfr where dfr is an integer in the */
3216 /* range +/-2*(emax+pmax), typically resulting from LogB. */
3217 /* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
3218 /* as usual. */
3219 /* ------------------------------------------------------------------ */
3220 #define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
3221 decFloat * decFloatScaleB(decFloat *result,
3222 const decFloat *dfl, const decFloat *dfr,
3223 decContext *set) {
3224 uInt digits; /* work */
3225 Int expr; /* dfr as an Int */
3226
3227 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3228 if (!DFISINT(dfr)) return decInvalid(result, set);
3229 digits=decFloatDigits(dfr); /* calculate digits */
3230
3231 #if DOUBLE
3232 if (digits>3) return decInvalid(result, set); /* definitely out of range */
3233 expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */
3234 #elif QUAD
3235 if (digits>5) return decInvalid(result, set); /* definitely out of range */
3236 expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */
3237 +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */
3238 #endif
3239 if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */
3240 /* [from now on no error possible] */
3241 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3242 if (DFISSIGNED(dfr)) expr=-expr;
3243 /* dfl is finite and expr is valid */
3244 *result=*dfl; /* copy to target */
3245 return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
3246 } /* decFloatScaleB */
3247
3248 /* ------------------------------------------------------------------ */
3249 /* decFloatShift -- shift the coefficient of a decFloat left or right */
3250 /* */
3251 /* result gets the result of shifting dfl */
3252 /* dfl is the source decFloat to shift */
3253 /* dfr is the count of digits to shift, an integer (with q=0) */
3254 /* set is the context */
3255 /* returns result */
3256 /* */
3257 /* The digits of the coefficient of dfl are shifted to the left (if */
3258 /* dfr is positive) or to the right (if dfr is negative) without */
3259 /* adjusting the exponent or the sign of dfl. */
3260 /* */
3261 /* dfr must be in the range -DECPMAX through +DECPMAX. */
3262 /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
3263 /* dfr must be valid). No status is set unless dfr is invalid or an */
3264 /* operand is an sNaN. The result is canonical. */
3265 /* ------------------------------------------------------------------ */
3266 decFloat * decFloatShift(decFloat *result,
3267 const decFloat *dfl, const decFloat *dfr,
3268 decContext *set) {
3269 Int shift; /* dfr as an Int */
3270 uByte buf[DECPMAX*2]; /* coefficient + padding */
3271 uInt digits, savestat; /* work */
3272 bcdnum num; /* .. */
3273
3274 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3275 if (!DFISINT(dfr)) return decInvalid(result, set);
3276 digits=decFloatDigits(dfr); /* calculate digits */
3277 if (digits>2) return decInvalid(result, set); /* definitely out of range */
3278 shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3279 if (shift>DECPMAX) return decInvalid(result, set); /* too big */
3280 /* [from here on no error or status change is possible] */
3281
3282 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3283 /* handle no-shift and all-shift (clear to zero) cases */
3284 if (shift==0) return decCanonical(result, dfl);
3285 if (shift==DECPMAX) { /* zero with sign */
3286 uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
3287 decFloatZero(result); /* make +0 */
3288 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
3289 /* [cannot safely use CopySign] */
3290 return result;
3291 }
3292 /* a real shift is needed: 0 < shift < DECPMAX */
3293 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3294 num.exponent=GETEXPUN(dfl);
3295 num.msd=buf;
3296 GETCOEFF(dfl, buf);
3297 if (DFISSIGNED(dfr)) { /* shift right */
3298 /* edge cases are taken care of, so this is easy */
3299 num.lsd=buf+DECPMAX-shift-1;
3300 }
3301 else { /* shift left -- zero padding needed to right */
3302 UINTAT(buf+DECPMAX)=0; /* 8 will handle most cases */
3303 UINTAT(buf+DECPMAX+4)=0; /* .. */
3304 if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
3305 num.msd+=shift;
3306 num.lsd=num.msd+DECPMAX-1;
3307 }
3308 savestat=set->status; /* record */
3309 decFinalize(result, &num, set);
3310 set->status=savestat; /* restore */
3311 return result;
3312 } /* decFloatShift */
3313
3314 /* ------------------------------------------------------------------ */
3315 /* decFloatSubtract -- subtract a decFloat from another */
3316 /* */
3317 /* result gets the result of subtracting dfr from dfl: */
3318 /* dfl is the first decFloat (lhs) */
3319 /* dfr is the second decFloat (rhs) */
3320 /* set is the context */
3321 /* returns result */
3322 /* */
3323 /* ------------------------------------------------------------------ */
3324 decFloat * decFloatSubtract(decFloat *result,
3325 const decFloat *dfl, const decFloat *dfr,
3326 decContext *set) {
3327 decFloat temp;
3328 /* NaNs must propagate without sign change */
3329 if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
3330 temp=*dfr; /* make a copy */
3331 DFBYTE(&temp, 0)^=0x80; /* flip sign */
3332 return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
3333 } /* decFloatSubtract */
3334
3335 /* ------------------------------------------------------------------ */
3336 /* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
3337 /* */
3338 /* df is the decFloat to round */
3339 /* set is the context */
3340 /* round is the rounding mode to use */
3341 /* returns a uInt or an Int, rounded according to the name */
3342 /* */
3343 /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
3344 /* outside the range of the target; Inexact will not be signaled for */
3345 /* simple rounding unless 'Exact' appears in the name. */
3346 /* ------------------------------------------------------------------ */
3347 uInt decFloatToUInt32(const decFloat *df, decContext *set,
3348 enum rounding round) {
3349 return decToInt32(df, set, round, 0, 1);}
3350
3351 uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
3352 enum rounding round) {
3353 return decToInt32(df, set, round, 1, 1);}
3354
3355 Int decFloatToInt32(const decFloat *df, decContext *set,
3356 enum rounding round) {
3357 return (Int)decToInt32(df, set, round, 0, 0);}
3358
3359 Int decFloatToInt32Exact(const decFloat *df, decContext *set,
3360 enum rounding round) {
3361 return (Int)decToInt32(df, set, round, 1, 0);}
3362
3363 /* ------------------------------------------------------------------ */
3364 /* decFloatToIntegral -- round to integral value (two flavours) */
3365 /* */
3366 /* result gets the result */
3367 /* df is the decFloat to round */
3368 /* set is the context */
3369 /* round is the rounding mode to use */
3370 /* returns result */
3371 /* */
3372 /* No exceptions, even Inexact, are raised except for sNaN input, or */
3373 /* if 'Exact' appears in the name. */
3374 /* ------------------------------------------------------------------ */
3375 decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
3376 decContext *set, enum rounding round) {
3377 return decToIntegral(result, df, set, round, 0);}
3378
3379 decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
3380 decContext *set) {
3381 return decToIntegral(result, df, set, set->round, 1);}
3382
3383 /* ------------------------------------------------------------------ */
3384 /* decFloatXor -- logical digitwise XOR of two decFloats */
3385 /* */
3386 /* result gets the result of XORing dfl and dfr */
3387 /* dfl is the first decFloat (lhs) */
3388 /* dfr is the second decFloat (rhs) */
3389 /* set is the context */
3390 /* returns result, which will be canonical with sign=0 */
3391 /* */
3392 /* The operands must be positive, finite with exponent q=0, and */
3393 /* comprise just zeros and ones; if not, Invalid operation results. */
3394 /* ------------------------------------------------------------------ */
3395 decFloat * decFloatXor(decFloat *result,
3396 const decFloat *dfl, const decFloat *dfr,
3397 decContext *set) {
3398 if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
3399 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
3400 /* the operands are positive finite integers (q=0) with just 0s and 1s */
3401 #if DOUBLE
3402 DFWORD(result, 0)=ZEROWORD
3403 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
3404 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
3405 #elif QUAD
3406 DFWORD(result, 0)=ZEROWORD
3407 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
3408 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
3409 DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
3410 DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
3411 #endif
3412 return result;
3413 } /* decFloatXor */
3414
3415 /* ------------------------------------------------------------------ */
3416 /* decInvalid -- set Invalid_operation result */
3417 /* */
3418 /* result gets a canonical NaN */
3419 /* set is the context */
3420 /* returns result */
3421 /* */
3422 /* status has Invalid_operation added */
3423 /* ------------------------------------------------------------------ */
3424 static decFloat *decInvalid(decFloat *result, decContext *set) {
3425 decFloatZero(result);
3426 DFWORD(result, 0)=DECFLOAT_qNaN;
3427 set->status|=DEC_Invalid_operation;
3428 return result;
3429 } /* decInvalid */
3430
3431 /* ------------------------------------------------------------------ */
3432 /* decInfinity -- set canonical Infinity with sign from a decFloat */
3433 /* */
3434 /* result gets a canonical Infinity */
3435 /* df is source decFloat (only the sign is used) */
3436 /* returns result */
3437 /* */
3438 /* df may be the same as result */
3439 /* ------------------------------------------------------------------ */
3440 static decFloat *decInfinity(decFloat *result, const decFloat *df) {
3441 uInt sign=DFWORD(df, 0); /* save source signword */
3442 decFloatZero(result); /* clear everything */
3443 DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
3444 return result;
3445 } /* decInfinity */
3446
3447 /* ------------------------------------------------------------------ */
3448 /* decNaNs -- handle NaN argument(s) */
3449 /* */
3450 /* result gets the result of handling dfl and dfr, one or both of */
3451 /* which is a NaN */
3452 /* dfl is the first decFloat (lhs) */
3453 /* dfr is the second decFloat (rhs) -- may be NULL for a single- */
3454 /* operand operation */
3455 /* set is the context */
3456 /* returns result */
3457 /* */
3458 /* Called when one or both operands is a NaN, and propagates the */
3459 /* appropriate result to res. When an sNaN is found, it is changed */
3460 /* to a qNaN and Invalid operation is set. */
3461 /* ------------------------------------------------------------------ */
3462 static decFloat *decNaNs(decFloat *result,
3463 const decFloat *dfl, const decFloat *dfr,
3464 decContext *set) {
3465 /* handle sNaNs first */
3466 if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
3467 if (DFISSNAN(dfl)) {
3468 decCanonical(result, dfl); /* propagate canonical sNaN */
3469 DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
3470 set->status|=DEC_Invalid_operation;
3471 return result;
3472 }
3473 /* one or both is a quiet NaN */
3474 if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */
3475 return decCanonical(result, dfl); /* propagate canonical qNaN */
3476 } /* decNaNs */
3477
3478 /* ------------------------------------------------------------------ */
3479 /* decNumCompare -- numeric comparison of two decFloats */
3480 /* */
3481 /* dfl is the left-hand decFloat, which is not a NaN */
3482 /* dfr is the right-hand decFloat, which is not a NaN */
3483 /* tot is 1 for total order compare, 0 for simple numeric */
3484 /* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
3485 /* */
3486 /* No error is possible; status and mode are unchanged. */
3487 /* ------------------------------------------------------------------ */
3488 static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
3489 Int sigl, sigr; /* LHS and RHS non-0 signums */
3490 Int shift; /* shift needed to align operands */
3491 uByte *ub, *uc; /* work */
3492 /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
3493 uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */
3494 uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */
3495
3496 sigl=1;
3497 if (DFISSIGNED(dfl)) {
3498 if (!DFISSIGNED(dfr)) { /* -LHS +RHS */
3499 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3500 return -1; /* RHS wins */
3501 }
3502 sigl=-1;
3503 }
3504 if (DFISSIGNED(dfr)) {
3505 if (!DFISSIGNED(dfl)) { /* +LHS -RHS */
3506 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3507 return +1; /* LHS wins */
3508 }
3509 }
3510
3511 /* signs are the same; operand(s) could be zero */
3512 sigr=-sigl; /* sign to return if abs(RHS) wins */
3513
3514 if (DFISINF(dfl)) {
3515 if (DFISINF(dfr)) return 0; /* both infinite & same sign */
3516 return sigl; /* inf > n */
3517 }
3518 if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
3519
3520 /* here, both are same sign and finite; calculate their offset */
3521 shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */
3522 /* [bias can be ignored -- the absolute exponent is not relevant] */
3523
3524 if (DFISZERO(dfl)) {
3525 if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */
3526 /* both are zero, return 0 if both same exponent or numeric compare */
3527 if (shift==0 || !tot) return 0;
3528 if (shift>0) return sigl;
3529 return sigr; /* [shift<0] */
3530 }
3531 else { /* LHS!=0 */
3532 if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */
3533 }
3534 /* both are known to be non-zero at this point */
3535
3536 /* if the exponents are so different that the coefficients do not */
3537 /* overlap (by even one digit) then a full comparison is not needed */
3538 if (abs(shift)>=DECPMAX) { /* no overlap */
3539 /* coefficients are known to be non-zero */
3540 if (shift>0) return sigl;
3541 return sigr; /* [shift<0] */
3542 }
3543
3544 /* decode the coefficients */
3545 /* (shift both right two if Quad to make a multiple of four) */
3546 #if QUAD
3547 ub=bufl; /* avoid type-pun violation */
3548 UINTAT(ub)=0;
3549 uc=bufr; /* avoid type-pun violation */
3550 UINTAT(uc)=0;
3551 #endif
3552 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
3553 GETCOEFF(dfr, bufr+QUAD*2); /* .. */
3554 if (shift==0) { /* aligned; common and easy */
3555 /* all multiples of four, here */
3556 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
3557 if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
3558 /* about to find a winner; go by bytes in case little-endian */
3559 for (;; ub++, uc++) {
3560 if (*ub>*uc) return sigl; /* difference found */
3561 if (*ub<*uc) return sigr; /* .. */
3562 }
3563 }
3564 } /* aligned */
3565 else if (shift>0) { /* lhs to left */
3566 ub=bufl; /* RHS pointer */
3567 /* pad bufl so right-aligned; most shifts will fit in 8 */
3568 UINTAT(bufl+DECPMAX+QUAD*2)=0; /* add eight zeros */
3569 UINTAT(bufl+DECPMAX+QUAD*2+4)=0; /* .. */
3570 if (shift>8) {
3571 /* more than eight; fill the rest, and also worth doing the */
3572 /* lead-in by fours */
3573 uByte *up; /* work */
3574 uByte *upend=bufl+DECPMAX+QUAD*2+shift;
3575 for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3576 /* [pads up to 36 in all for Quad] */
3577 for (;; ub+=4) {
3578 if (UINTAT(ub)!=0) return sigl;
3579 if (ub+4>bufl+shift-4) break;
3580 }
3581 }
3582 /* check remaining leading digits */
3583 for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
3584 /* now start the overlapped part; bufl has been padded, so the */
3585 /* comparison can go for the full length of bufr, which is a */
3586 /* multiple of 4 bytes */
3587 for (uc=bufr; ; uc+=4, ub+=4) {
3588 if (UINTAT(uc)!=UINTAT(ub)) { /* mismatch found */
3589 for (;; uc++, ub++) { /* check from left [little-endian?] */
3590 if (*ub>*uc) return sigl; /* difference found */
3591 if (*ub<*uc) return sigr; /* .. */
3592 }
3593 } /* mismatch */
3594 if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
3595 }
3596 } /* shift>0 */
3597
3598 else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
3599 uc=bufr; /* RHS pointer */
3600 /* pad bufr so right-aligned; most shifts will fit in 8 */
3601 UINTAT(bufr+DECPMAX+QUAD*2)=0; /* add eight zeros */
3602 UINTAT(bufr+DECPMAX+QUAD*2+4)=0; /* .. */
3603 if (shift<-8) {
3604 /* more than eight; fill the rest, and also worth doing the */
3605 /* lead-in by fours */
3606 uByte *up; /* work */
3607 uByte *upend=bufr+DECPMAX+QUAD*2-shift;
3608 for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3609 /* [pads up to 36 in all for Quad] */
3610 for (;; uc+=4) {
3611 if (UINTAT(uc)!=0) return sigr;
3612 if (uc+4>bufr-shift-4) break;
3613 }
3614 }
3615 /* check remaining leading digits */
3616 for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
3617 /* now start the overlapped part; bufr has been padded, so the */
3618 /* comparison can go for the full length of bufl, which is a */
3619 /* multiple of 4 bytes */
3620 for (ub=bufl; ; ub+=4, uc+=4) {
3621 if (UINTAT(ub)!=UINTAT(uc)) { /* mismatch found */
3622 for (;; ub++, uc++) { /* check from left [little-endian?] */
3623 if (*ub>*uc) return sigl; /* difference found */
3624 if (*ub<*uc) return sigr; /* .. */
3625 }
3626 } /* mismatch */
3627 if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
3628 }
3629 } /* shift<0 */
3630
3631 /* Here when compare equal */
3632 if (!tot) return 0; /* numerically equal */
3633 /* total ordering .. exponent matters */
3634 if (shift>0) return sigl; /* total order by exponent */
3635 if (shift<0) return sigr; /* .. */
3636 return 0;
3637 } /* decNumCompare */
3638
3639 /* ------------------------------------------------------------------ */
3640 /* decToInt32 -- local routine to effect ToInteger conversions */
3641 /* */
3642 /* df is the decFloat to convert */
3643 /* set is the context */
3644 /* rmode is the rounding mode to use */
3645 /* exact is 1 if Inexact should be signalled */
3646 /* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
3647 /* returns 32-bit result as a uInt */
3648 /* */
3649 /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
3650 /* these cases 0 is returned. */
3651 /* ------------------------------------------------------------------ */
3652 static uInt decToInt32(const decFloat *df, decContext *set,
3653 enum rounding rmode, Flag exact, Flag unsign) {
3654 Int exp; /* exponent */
3655 uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */
3656 uInt hi, lo; /* .. penultimate, least, etc. */
3657 decFloat zero, result; /* work */
3658 Int i; /* .. */
3659
3660 /* Start decoding the argument */
3661 sourhi=DFWORD(df, 0); /* top word */
3662 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
3663 if (EXPISSPECIAL(exp)) { /* is special? */
3664 set->status|=DEC_Invalid_operation; /* signal */
3665 return 0;
3666 }
3667
3668 /* Here when the argument is finite */
3669 if (GETEXPUN(df)==0) result=*df; /* already a true integer */
3670 else { /* need to round to integer */
3671 enum rounding saveround; /* saver */
3672 uInt savestatus; /* .. */
3673 saveround=set->round; /* save rounding mode .. */
3674 savestatus=set->status; /* .. and status */
3675 set->round=rmode; /* set mode */
3676 decFloatZero(&zero); /* make 0E+0 */
3677 set->status=0; /* clear */
3678 decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
3679 set->round=saveround; /* restore rounding mode .. */
3680 if (exact) set->status|=savestatus; /* include Inexact */
3681 else set->status=savestatus; /* .. or just original status */
3682 }
3683
3684 /* only the last four declets of the coefficient can contain */
3685 /* non-zero; check for others (and also NaN or Infinity from the */
3686 /* Quantize) first (see DFISZERO for explanation): */
3687 /* decFloatShow(&result, "sofar"); */
3688 #if DOUBLE
3689 if ((DFWORD(&result, 0)&0x1c03ff00)!=0
3690 || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3691 #elif QUAD
3692 if ((DFWORD(&result, 2)&0xffffff00)!=0
3693 || DFWORD(&result, 1)!=0
3694 || (DFWORD(&result, 0)&0x1c003fff)!=0
3695 || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3696 #endif
3697 set->status|=DEC_Invalid_operation; /* Invalid or out of range */
3698 return 0;
3699 }
3700 /* get last twelve digits of the coefficent into hi & ho, base */
3701 /* 10**9 (see GETCOEFFBILL): */
3702 sourlo=DFWORD(&result, DECWORDS-1);
3703 lo=DPD2BIN0[sourlo&0x3ff]
3704 +DPD2BINK[(sourlo>>10)&0x3ff]
3705 +DPD2BINM[(sourlo>>20)&0x3ff];
3706 sourpen=DFWORD(&result, DECWORDS-2);
3707 hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
3708
3709 /* according to request, check range carefully */
3710 if (unsign) {
3711 if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
3712 set->status|=DEC_Invalid_operation; /* out of range */
3713 return 0;
3714 }
3715 return hi*BILLION+lo;
3716 }
3717 /* signed */
3718 if (hi>2 || (hi==2 && lo>147483647)) {
3719 /* handle the usual edge case */
3720 if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
3721 set->status|=DEC_Invalid_operation; /* truly out of range */
3722 return 0;
3723 }
3724 i=hi*BILLION+lo;
3725 if (DFISSIGNED(&result)) i=-i;
3726 return (uInt)i;
3727 } /* decToInt32 */
3728
3729 /* ------------------------------------------------------------------ */
3730 /* decToIntegral -- local routine to effect ToIntegral value */
3731 /* */
3732 /* result gets the result */
3733 /* df is the decFloat to round */
3734 /* set is the context */
3735 /* rmode is the rounding mode to use */
3736 /* exact is 1 if Inexact should be signalled */
3737 /* returns result */
3738 /* ------------------------------------------------------------------ */
3739 static decFloat * decToIntegral(decFloat *result, const decFloat *df,
3740 decContext *set, enum rounding rmode,
3741 Flag exact) {
3742 Int exp; /* exponent */
3743 uInt sourhi; /* top word from source decFloat */
3744 enum rounding saveround; /* saver */
3745 uInt savestatus; /* .. */
3746 decFloat zero; /* work */
3747
3748 /* Start decoding the argument */
3749 sourhi=DFWORD(df, 0); /* top word */
3750 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
3751
3752 if (EXPISSPECIAL(exp)) { /* is special? */
3753 /* NaNs are handled as usual */
3754 if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
3755 /* must be infinite; return canonical infinity with sign of df */
3756 return decInfinity(result, df);
3757 }
3758
3759 /* Here when the argument is finite */
3760 /* complete extraction of the exponent */
3761 exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */
3762
3763 if (exp>=0) return decCanonical(result, df); /* already integral */
3764
3765 saveround=set->round; /* save rounding mode .. */
3766 savestatus=set->status; /* .. and status */
3767 set->round=rmode; /* set mode */
3768 decFloatZero(&zero); /* make 0E+0 */
3769 decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
3770 set->round=saveround; /* restore rounding mode .. */
3771 if (!exact) set->status=savestatus; /* .. and status, unless exact */
3772 return result;
3773 } /* decToIntegral */
Mirror of https://gcc.gnu.org/git/gcc.git