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