Mercurial > hg > CbC > CbC_gcc
comparison libgcc/config/libbid/bid_internal.h @ 0:a06113de4d67
first commit
author | kent <kent@cr.ie.u-ryukyu.ac.jp> |
---|---|
date | Fri, 17 Jul 2009 14:47:48 +0900 |
parents | |
children | 04ced10e8804 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:a06113de4d67 |
---|---|
1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc. | |
2 | |
3 This file is part of GCC. | |
4 | |
5 GCC is free software; you can redistribute it and/or modify it under | |
6 the terms of the GNU General Public License as published by the Free | |
7 Software Foundation; either version 3, or (at your option) any later | |
8 version. | |
9 | |
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY | |
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
13 for more details. | |
14 | |
15 Under Section 7 of GPL version 3, you are granted additional | |
16 permissions described in the GCC Runtime Library Exception, version | |
17 3.1, as published by the Free Software Foundation. | |
18 | |
19 You should have received a copy of the GNU General Public License and | |
20 a copy of the GCC Runtime Library Exception along with this program; | |
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see | |
22 <http://www.gnu.org/licenses/>. */ | |
23 | |
24 #ifndef __BIDECIMAL_H | |
25 #define __BIDECIMAL_H | |
26 | |
27 #include "bid_conf.h" | |
28 #include "bid_functions.h" | |
29 | |
30 #define __BID_INLINE__ static __inline | |
31 | |
32 /********************************************************************* | |
33 * | |
34 * Logical Shift Macros | |
35 * | |
36 *********************************************************************/ | |
37 | |
38 #define __shr_128(Q, A, k) \ | |
39 { \ | |
40 (Q).w[0] = (A).w[0] >> k; \ | |
41 (Q).w[0] |= (A).w[1] << (64-k); \ | |
42 (Q).w[1] = (A).w[1] >> k; \ | |
43 } | |
44 | |
45 #define __shr_128_long(Q, A, k) \ | |
46 { \ | |
47 if((k)<64) { \ | |
48 (Q).w[0] = (A).w[0] >> k; \ | |
49 (Q).w[0] |= (A).w[1] << (64-k); \ | |
50 (Q).w[1] = (A).w[1] >> k; \ | |
51 } \ | |
52 else { \ | |
53 (Q).w[0] = (A).w[1]>>((k)-64); \ | |
54 (Q).w[1] = 0; \ | |
55 } \ | |
56 } | |
57 | |
58 #define __shl_128_long(Q, A, k) \ | |
59 { \ | |
60 if((k)<64) { \ | |
61 (Q).w[1] = (A).w[1] << k; \ | |
62 (Q).w[1] |= (A).w[0] >> (64-k); \ | |
63 (Q).w[0] = (A).w[0] << k; \ | |
64 } \ | |
65 else { \ | |
66 (Q).w[1] = (A).w[0]<<((k)-64); \ | |
67 (Q).w[0] = 0; \ | |
68 } \ | |
69 } | |
70 | |
71 #define __low_64(Q) (Q).w[0] | |
72 /********************************************************************* | |
73 * | |
74 * String Macros | |
75 * | |
76 *********************************************************************/ | |
77 #define tolower_macro(x) (((unsigned char)((x)-'A')<=('Z'-'A'))?((x)-'A'+'a'):(x)) | |
78 /********************************************************************* | |
79 * | |
80 * Compare Macros | |
81 * | |
82 *********************************************************************/ | |
83 // greater than | |
84 // return 0 if A<=B | |
85 // non-zero if A>B | |
86 #define __unsigned_compare_gt_128(A, B) \ | |
87 ((A.w[1]>B.w[1]) || ((A.w[1]==B.w[1]) && (A.w[0]>B.w[0]))) | |
88 // greater-or-equal | |
89 #define __unsigned_compare_ge_128(A, B) \ | |
90 ((A.w[1]>B.w[1]) || ((A.w[1]==B.w[1]) && (A.w[0]>=B.w[0]))) | |
91 #define __test_equal_128(A, B) (((A).w[1]==(B).w[1]) && ((A).w[0]==(B).w[0])) | |
92 // tighten exponent range | |
93 #define __tight_bin_range_128(bp, P, bin_expon) \ | |
94 { \ | |
95 UINT64 M; \ | |
96 M = 1; \ | |
97 (bp) = (bin_expon); \ | |
98 if((bp)<63) { \ | |
99 M <<= ((bp)+1); \ | |
100 if((P).w[0] >= M) (bp)++; } \ | |
101 else if((bp)>64) { \ | |
102 M <<= ((bp)+1-64); \ | |
103 if(((P).w[1]>M) ||((P).w[1]==M && (P).w[0]))\ | |
104 (bp)++; } \ | |
105 else if((P).w[1]) (bp)++; \ | |
106 } | |
107 /********************************************************************* | |
108 * | |
109 * Add/Subtract Macros | |
110 * | |
111 *********************************************************************/ | |
112 // add 64-bit value to 128-bit | |
113 #define __add_128_64(R128, A128, B64) \ | |
114 { \ | |
115 UINT64 R64H; \ | |
116 R64H = (A128).w[1]; \ | |
117 (R128).w[0] = (B64) + (A128).w[0]; \ | |
118 if((R128).w[0] < (B64)) \ | |
119 R64H ++; \ | |
120 (R128).w[1] = R64H; \ | |
121 } | |
122 // subtract 64-bit value from 128-bit | |
123 #define __sub_128_64(R128, A128, B64) \ | |
124 { \ | |
125 UINT64 R64H; \ | |
126 R64H = (A128).w[1]; \ | |
127 if((A128).w[0] < (B64)) \ | |
128 R64H --; \ | |
129 (R128).w[1] = R64H; \ | |
130 (R128).w[0] = (A128).w[0] - (B64); \ | |
131 } | |
132 // add 128-bit value to 128-bit | |
133 // assume no carry-out | |
134 #define __add_128_128(R128, A128, B128) \ | |
135 { \ | |
136 UINT128 Q128; \ | |
137 Q128.w[1] = (A128).w[1]+(B128).w[1]; \ | |
138 Q128.w[0] = (B128).w[0] + (A128).w[0]; \ | |
139 if(Q128.w[0] < (B128).w[0]) \ | |
140 Q128.w[1] ++; \ | |
141 (R128).w[1] = Q128.w[1]; \ | |
142 (R128).w[0] = Q128.w[0]; \ | |
143 } | |
144 #define __sub_128_128(R128, A128, B128) \ | |
145 { \ | |
146 UINT128 Q128; \ | |
147 Q128.w[1] = (A128).w[1]-(B128).w[1]; \ | |
148 Q128.w[0] = (A128).w[0] - (B128).w[0]; \ | |
149 if((A128).w[0] < (B128).w[0]) \ | |
150 Q128.w[1] --; \ | |
151 (R128).w[1] = Q128.w[1]; \ | |
152 (R128).w[0] = Q128.w[0]; \ | |
153 } | |
154 #define __add_carry_out(S, CY, X, Y) \ | |
155 { \ | |
156 UINT64 X1=X; \ | |
157 S = X + Y; \ | |
158 CY = (S<X1) ? 1 : 0; \ | |
159 } | |
160 #define __add_carry_in_out(S, CY, X, Y, CI) \ | |
161 { \ | |
162 UINT64 X1; \ | |
163 X1 = X + CI; \ | |
164 S = X1 + Y; \ | |
165 CY = ((S<X1) || (X1<CI)) ? 1 : 0; \ | |
166 } | |
167 #define __sub_borrow_out(S, CY, X, Y) \ | |
168 { \ | |
169 UINT64 X1=X; \ | |
170 S = X - Y; \ | |
171 CY = (S>X1) ? 1 : 0; \ | |
172 } | |
173 #define __sub_borrow_in_out(S, CY, X, Y, CI) \ | |
174 { \ | |
175 UINT64 X1, X0=X; \ | |
176 X1 = X - CI; \ | |
177 S = X1 - Y; \ | |
178 CY = ((S>X1) || (X1>X0)) ? 1 : 0; \ | |
179 } | |
180 // increment C128 and check for rounding overflow: | |
181 // if (C_128) = 10^34 then (C_128) = 10^33 and increment the exponent | |
182 #define INCREMENT(C_128, exp) \ | |
183 { \ | |
184 C_128.w[0]++; \ | |
185 if (C_128.w[0] == 0) C_128.w[1]++; \ | |
186 if (C_128.w[1] == 0x0001ed09bead87c0ull && \ | |
187 C_128.w[0] == 0x378d8e6400000000ull) { \ | |
188 exp++; \ | |
189 C_128.w[1] = 0x0000314dc6448d93ull; \ | |
190 C_128.w[0] = 0x38c15b0a00000000ull; \ | |
191 } \ | |
192 } | |
193 // decrement C128 and check for rounding underflow, but only at the | |
194 // boundary: if C_128 = 10^33 - 1 and exp > 0 then C_128 = 10^34 - 1 | |
195 // and decrement the exponent | |
196 #define DECREMENT(C_128, exp) \ | |
197 { \ | |
198 C_128.w[0]--; \ | |
199 if (C_128.w[0] == 0xffffffffffffffffull) C_128.w[1]--; \ | |
200 if (C_128.w[1] == 0x0000314dc6448d93ull && \ | |
201 C_128.w[0] == 0x38c15b09ffffffffull && exp > 0) { \ | |
202 exp--; \ | |
203 C_128.w[1] = 0x0001ed09bead87c0ull; \ | |
204 C_128.w[0] = 0x378d8e63ffffffffull; \ | |
205 } \ | |
206 } | |
207 | |
208 /********************************************************************* | |
209 * | |
210 * Multiply Macros | |
211 * | |
212 *********************************************************************/ | |
213 #define __mul_64x64_to_64(P64, CX, CY) (P64) = (CX) * (CY) | |
214 /*************************************** | |
215 * Signed, Full 64x64-bit Multiply | |
216 ***************************************/ | |
217 #define __imul_64x64_to_128(P, CX, CY) \ | |
218 { \ | |
219 UINT64 SX, SY; \ | |
220 __mul_64x64_to_128(P, CX, CY); \ | |
221 \ | |
222 SX = ((SINT64)(CX))>>63; \ | |
223 SY = ((SINT64)(CY))>>63; \ | |
224 SX &= CY; SY &= CX; \ | |
225 \ | |
226 (P).w[1] = (P).w[1] - SX - SY; \ | |
227 } | |
228 /*************************************** | |
229 * Signed, Full 64x128-bit Multiply | |
230 ***************************************/ | |
231 #define __imul_64x128_full(Ph, Ql, A, B) \ | |
232 { \ | |
233 UINT128 ALBL, ALBH, QM2, QM; \ | |
234 \ | |
235 __imul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
236 __imul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
237 \ | |
238 (Ql).w[0] = ALBL.w[0]; \ | |
239 QM.w[0] = ALBL.w[1]; \ | |
240 QM.w[1] = ((SINT64)ALBL.w[1])>>63; \ | |
241 __add_128_128(QM2, ALBH, QM); \ | |
242 (Ql).w[1] = QM2.w[0]; \ | |
243 Ph = QM2.w[1]; \ | |
244 } | |
245 /***************************************************** | |
246 * Unsigned Multiply Macros | |
247 *****************************************************/ | |
248 // get full 64x64bit product | |
249 // | |
250 #define __mul_64x64_to_128(P, CX, CY) \ | |
251 { \ | |
252 UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2;\ | |
253 CXH = (CX) >> 32; \ | |
254 CXL = (UINT32)(CX); \ | |
255 CYH = (CY) >> 32; \ | |
256 CYL = (UINT32)(CY); \ | |
257 \ | |
258 PM = CXH*CYL; \ | |
259 PH = CXH*CYH; \ | |
260 PL = CXL*CYL; \ | |
261 PM2 = CXL*CYH; \ | |
262 PH += (PM>>32); \ | |
263 PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \ | |
264 \ | |
265 (P).w[1] = PH + (PM>>32); \ | |
266 (P).w[0] = (PM<<32)+(UINT32)PL; \ | |
267 } | |
268 // get full 64x64bit product | |
269 // Note: | |
270 // This macro is used for CX < 2^61, CY < 2^61 | |
271 // | |
272 #define __mul_64x64_to_128_fast(P, CX, CY) \ | |
273 { \ | |
274 UINT64 CXH, CXL, CYH, CYL, PL, PH, PM; \ | |
275 CXH = (CX) >> 32; \ | |
276 CXL = (UINT32)(CX); \ | |
277 CYH = (CY) >> 32; \ | |
278 CYL = (UINT32)(CY); \ | |
279 \ | |
280 PM = CXH*CYL; \ | |
281 PL = CXL*CYL; \ | |
282 PH = CXH*CYH; \ | |
283 PM += CXL*CYH; \ | |
284 PM += (PL>>32); \ | |
285 \ | |
286 (P).w[1] = PH + (PM>>32); \ | |
287 (P).w[0] = (PM<<32)+(UINT32)PL; \ | |
288 } | |
289 // used for CX< 2^60 | |
290 #define __sqr64_fast(P, CX) \ | |
291 { \ | |
292 UINT64 CXH, CXL, PL, PH, PM; \ | |
293 CXH = (CX) >> 32; \ | |
294 CXL = (UINT32)(CX); \ | |
295 \ | |
296 PM = CXH*CXL; \ | |
297 PL = CXL*CXL; \ | |
298 PH = CXH*CXH; \ | |
299 PM += PM; \ | |
300 PM += (PL>>32); \ | |
301 \ | |
302 (P).w[1] = PH + (PM>>32); \ | |
303 (P).w[0] = (PM<<32)+(UINT32)PL; \ | |
304 } | |
305 // get full 64x64bit product | |
306 // Note: | |
307 // This implementation is used for CX < 2^61, CY < 2^61 | |
308 // | |
309 #define __mul_64x64_to_64_high_fast(P, CX, CY) \ | |
310 { \ | |
311 UINT64 CXH, CXL, CYH, CYL, PL, PH, PM; \ | |
312 CXH = (CX) >> 32; \ | |
313 CXL = (UINT32)(CX); \ | |
314 CYH = (CY) >> 32; \ | |
315 CYL = (UINT32)(CY); \ | |
316 \ | |
317 PM = CXH*CYL; \ | |
318 PL = CXL*CYL; \ | |
319 PH = CXH*CYH; \ | |
320 PM += CXL*CYH; \ | |
321 PM += (PL>>32); \ | |
322 \ | |
323 (P) = PH + (PM>>32); \ | |
324 } | |
325 // get full 64x64bit product | |
326 // | |
327 #define __mul_64x64_to_128_full(P, CX, CY) \ | |
328 { \ | |
329 UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2;\ | |
330 CXH = (CX) >> 32; \ | |
331 CXL = (UINT32)(CX); \ | |
332 CYH = (CY) >> 32; \ | |
333 CYL = (UINT32)(CY); \ | |
334 \ | |
335 PM = CXH*CYL; \ | |
336 PH = CXH*CYH; \ | |
337 PL = CXL*CYL; \ | |
338 PM2 = CXL*CYH; \ | |
339 PH += (PM>>32); \ | |
340 PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \ | |
341 \ | |
342 (P).w[1] = PH + (PM>>32); \ | |
343 (P).w[0] = (PM<<32)+(UINT32)PL; \ | |
344 } | |
345 #define __mul_128x128_high(Q, A, B) \ | |
346 { \ | |
347 UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2; \ | |
348 \ | |
349 __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]); \ | |
350 __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]); \ | |
351 __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]); \ | |
352 __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]); \ | |
353 \ | |
354 __add_128_128(QM, ALBH, AHBL); \ | |
355 __add_128_64(QM2, QM, ALBL.w[1]); \ | |
356 __add_128_64((Q), AHBH, QM2.w[1]); \ | |
357 } | |
358 #define __mul_128x128_full(Qh, Ql, A, B) \ | |
359 { \ | |
360 UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2; \ | |
361 \ | |
362 __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]); \ | |
363 __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]); \ | |
364 __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]); \ | |
365 __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]); \ | |
366 \ | |
367 __add_128_128(QM, ALBH, AHBL); \ | |
368 (Ql).w[0] = ALBL.w[0]; \ | |
369 __add_128_64(QM2, QM, ALBL.w[1]); \ | |
370 __add_128_64((Qh), AHBH, QM2.w[1]); \ | |
371 (Ql).w[1] = QM2.w[0]; \ | |
372 } | |
373 #define __mul_128x128_low(Ql, A, B) \ | |
374 { \ | |
375 UINT128 ALBL; \ | |
376 UINT64 QM64; \ | |
377 \ | |
378 __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]); \ | |
379 QM64 = (B).w[0]*(A).w[1] + (A).w[0]*(B).w[1]; \ | |
380 \ | |
381 (Ql).w[0] = ALBL.w[0]; \ | |
382 (Ql).w[1] = QM64 + ALBL.w[1]; \ | |
383 } | |
384 #define __mul_64x128_low(Ql, A, B) \ | |
385 { \ | |
386 UINT128 ALBL, ALBH, QM2; \ | |
387 __mul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
388 __mul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
389 (Ql).w[0] = ALBL.w[0]; \ | |
390 __add_128_64(QM2, ALBH, ALBL.w[1]); \ | |
391 (Ql).w[1] = QM2.w[0]; \ | |
392 } | |
393 #define __mul_64x128_full(Ph, Ql, A, B) \ | |
394 { \ | |
395 UINT128 ALBL, ALBH, QM2; \ | |
396 \ | |
397 __mul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
398 __mul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
399 \ | |
400 (Ql).w[0] = ALBL.w[0]; \ | |
401 __add_128_64(QM2, ALBH, ALBL.w[1]); \ | |
402 (Ql).w[1] = QM2.w[0]; \ | |
403 Ph = QM2.w[1]; \ | |
404 } | |
405 #define __mul_64x128_to_192(Q, A, B) \ | |
406 { \ | |
407 UINT128 ALBL, ALBH, QM2; \ | |
408 \ | |
409 __mul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
410 __mul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
411 \ | |
412 (Q).w[0] = ALBL.w[0]; \ | |
413 __add_128_64(QM2, ALBH, ALBL.w[1]); \ | |
414 (Q).w[1] = QM2.w[0]; \ | |
415 (Q).w[2] = QM2.w[1]; \ | |
416 } | |
417 #define __mul_64x128_to192(Q, A, B) \ | |
418 { \ | |
419 UINT128 ALBL, ALBH, QM2; \ | |
420 \ | |
421 __mul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
422 __mul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
423 \ | |
424 (Q).w[0] = ALBL.w[0]; \ | |
425 __add_128_64(QM2, ALBH, ALBL.w[1]); \ | |
426 (Q).w[1] = QM2.w[0]; \ | |
427 (Q).w[2] = QM2.w[1]; \ | |
428 } | |
429 #define __mul_128x128_to_256(P256, A, B) \ | |
430 { \ | |
431 UINT128 Qll, Qlh; \ | |
432 UINT64 Phl, Phh, CY1, CY2; \ | |
433 \ | |
434 __mul_64x128_full(Phl, Qll, A.w[0], B); \ | |
435 __mul_64x128_full(Phh, Qlh, A.w[1], B); \ | |
436 (P256).w[0] = Qll.w[0]; \ | |
437 __add_carry_out((P256).w[1],CY1, Qlh.w[0], Qll.w[1]); \ | |
438 __add_carry_in_out((P256).w[2],CY2, Qlh.w[1], Phl, CY1); \ | |
439 (P256).w[3] = Phh + CY2; \ | |
440 } | |
441 // | |
442 // For better performance, will check A.w[1] against 0, | |
443 // but not B.w[1] | |
444 // Use this macro accordingly | |
445 #define __mul_128x128_to_256_check_A(P256, A, B) \ | |
446 { \ | |
447 UINT128 Qll, Qlh; \ | |
448 UINT64 Phl, Phh, CY1, CY2; \ | |
449 \ | |
450 __mul_64x128_full(Phl, Qll, A.w[0], B); \ | |
451 (P256).w[0] = Qll.w[0]; \ | |
452 if(A.w[1]) { \ | |
453 __mul_64x128_full(Phh, Qlh, A.w[1], B); \ | |
454 __add_carry_out((P256).w[1],CY1, Qlh.w[0], Qll.w[1]); \ | |
455 __add_carry_in_out((P256).w[2],CY2, Qlh.w[1], Phl, CY1); \ | |
456 (P256).w[3] = Phh + CY2; } \ | |
457 else { \ | |
458 (P256).w[1] = Qll.w[1]; \ | |
459 (P256).w[2] = Phl; \ | |
460 (P256).w[3] = 0; } \ | |
461 } | |
462 #define __mul_64x192_to_256(lP, lA, lB) \ | |
463 { \ | |
464 UINT128 lP0,lP1,lP2; \ | |
465 UINT64 lC; \ | |
466 __mul_64x64_to_128(lP0, lA, (lB).w[0]); \ | |
467 __mul_64x64_to_128(lP1, lA, (lB).w[1]); \ | |
468 __mul_64x64_to_128(lP2, lA, (lB).w[2]); \ | |
469 (lP).w[0] = lP0.w[0]; \ | |
470 __add_carry_out((lP).w[1],lC,lP1.w[0],lP0.w[1]); \ | |
471 __add_carry_in_out((lP).w[2],lC,lP2.w[0],lP1.w[1],lC); \ | |
472 (lP).w[3] = lP2.w[1] + lC; \ | |
473 } | |
474 #define __mul_64x256_to_320(P, A, B) \ | |
475 { \ | |
476 UINT128 lP0,lP1,lP2,lP3; \ | |
477 UINT64 lC; \ | |
478 __mul_64x64_to_128(lP0, A, (B).w[0]); \ | |
479 __mul_64x64_to_128(lP1, A, (B).w[1]); \ | |
480 __mul_64x64_to_128(lP2, A, (B).w[2]); \ | |
481 __mul_64x64_to_128(lP3, A, (B).w[3]); \ | |
482 (P).w[0] = lP0.w[0]; \ | |
483 __add_carry_out((P).w[1],lC,lP1.w[0],lP0.w[1]); \ | |
484 __add_carry_in_out((P).w[2],lC,lP2.w[0],lP1.w[1],lC); \ | |
485 __add_carry_in_out((P).w[3],lC,lP3.w[0],lP2.w[1],lC); \ | |
486 (P).w[4] = lP3.w[1] + lC; \ | |
487 } | |
488 #define __mul_192x192_to_384(P, A, B) \ | |
489 { \ | |
490 UINT256 P0,P1,P2; \ | |
491 UINT64 CY; \ | |
492 __mul_64x192_to_256(P0, (A).w[0], B); \ | |
493 __mul_64x192_to_256(P1, (A).w[1], B); \ | |
494 __mul_64x192_to_256(P2, (A).w[2], B); \ | |
495 (P).w[0] = P0.w[0]; \ | |
496 __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]); \ | |
497 __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \ | |
498 __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \ | |
499 (P).w[4] = P1.w[3] + CY; \ | |
500 __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]); \ | |
501 __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY); \ | |
502 __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY); \ | |
503 (P).w[5] = P2.w[3] + CY; \ | |
504 } | |
505 #define __mul_64x320_to_384(P, A, B) \ | |
506 { \ | |
507 UINT128 lP0,lP1,lP2,lP3,lP4; \ | |
508 UINT64 lC; \ | |
509 __mul_64x64_to_128(lP0, A, (B).w[0]); \ | |
510 __mul_64x64_to_128(lP1, A, (B).w[1]); \ | |
511 __mul_64x64_to_128(lP2, A, (B).w[2]); \ | |
512 __mul_64x64_to_128(lP3, A, (B).w[3]); \ | |
513 __mul_64x64_to_128(lP4, A, (B).w[4]); \ | |
514 (P).w[0] = lP0.w[0]; \ | |
515 __add_carry_out((P).w[1],lC,lP1.w[0],lP0.w[1]); \ | |
516 __add_carry_in_out((P).w[2],lC,lP2.w[0],lP1.w[1],lC); \ | |
517 __add_carry_in_out((P).w[3],lC,lP3.w[0],lP2.w[1],lC); \ | |
518 __add_carry_in_out((P).w[4],lC,lP4.w[0],lP3.w[1],lC); \ | |
519 (P).w[5] = lP4.w[1] + lC; \ | |
520 } | |
521 // A*A | |
522 // Full 128x128-bit product | |
523 #define __sqr128_to_256(P256, A) \ | |
524 { \ | |
525 UINT128 Qll, Qlh, Qhh; \ | |
526 UINT64 TMP_C1, TMP_C2; \ | |
527 \ | |
528 __mul_64x64_to_128(Qhh, A.w[1], A.w[1]); \ | |
529 __mul_64x64_to_128(Qlh, A.w[0], A.w[1]); \ | |
530 Qhh.w[1] += (Qlh.w[1]>>63); \ | |
531 Qlh.w[1] = (Qlh.w[1]+Qlh.w[1])|(Qlh.w[0]>>63); \ | |
532 Qlh.w[0] += Qlh.w[0]; \ | |
533 __mul_64x64_to_128(Qll, A.w[0], A.w[0]); \ | |
534 \ | |
535 __add_carry_out((P256).w[1],TMP_C1, Qlh.w[0], Qll.w[1]); \ | |
536 (P256).w[0] = Qll.w[0]; \ | |
537 __add_carry_in_out((P256).w[2],TMP_C2, Qlh.w[1], Qhh.w[0], TMP_C1); \ | |
538 (P256).w[3] = Qhh.w[1]+TMP_C2; \ | |
539 } | |
540 #define __mul_128x128_to_256_low_high(PQh, PQl, A, B) \ | |
541 { \ | |
542 UINT128 Qll, Qlh; \ | |
543 UINT64 Phl, Phh, C1, C2; \ | |
544 \ | |
545 __mul_64x128_full(Phl, Qll, A.w[0], B); \ | |
546 __mul_64x128_full(Phh, Qlh, A.w[1], B); \ | |
547 (PQl).w[0] = Qll.w[0]; \ | |
548 __add_carry_out((PQl).w[1],C1, Qlh.w[0], Qll.w[1]); \ | |
549 __add_carry_in_out((PQh).w[0],C2, Qlh.w[1], Phl, C1); \ | |
550 (PQh).w[1] = Phh + C2; \ | |
551 } | |
552 #define __mul_256x256_to_512(P, A, B) \ | |
553 { \ | |
554 UINT512 P0,P1,P2,P3; \ | |
555 UINT64 CY; \ | |
556 __mul_64x256_to_320(P0, (A).w[0], B); \ | |
557 __mul_64x256_to_320(P1, (A).w[1], B); \ | |
558 __mul_64x256_to_320(P2, (A).w[2], B); \ | |
559 __mul_64x256_to_320(P3, (A).w[3], B); \ | |
560 (P).w[0] = P0.w[0]; \ | |
561 __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]); \ | |
562 __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \ | |
563 __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \ | |
564 __add_carry_in_out((P).w[4],CY,P1.w[3],P0.w[4],CY); \ | |
565 (P).w[5] = P1.w[4] + CY; \ | |
566 __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]); \ | |
567 __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY); \ | |
568 __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY); \ | |
569 __add_carry_in_out((P).w[5],CY,P2.w[3],(P).w[5],CY); \ | |
570 (P).w[6] = P2.w[4] + CY; \ | |
571 __add_carry_out((P).w[3],CY,P3.w[0],(P).w[3]); \ | |
572 __add_carry_in_out((P).w[4],CY,P3.w[1],(P).w[4],CY); \ | |
573 __add_carry_in_out((P).w[5],CY,P3.w[2],(P).w[5],CY); \ | |
574 __add_carry_in_out((P).w[6],CY,P3.w[3],(P).w[6],CY); \ | |
575 (P).w[7] = P3.w[4] + CY; \ | |
576 } | |
577 #define __mul_192x256_to_448(P, A, B) \ | |
578 { \ | |
579 UINT512 P0,P1,P2; \ | |
580 UINT64 CY; \ | |
581 __mul_64x256_to_320(P0, (A).w[0], B); \ | |
582 __mul_64x256_to_320(P1, (A).w[1], B); \ | |
583 __mul_64x256_to_320(P2, (A).w[2], B); \ | |
584 (P).w[0] = P0.w[0]; \ | |
585 __add_carry_out((P).w[1],CY,P1.w[0],P0.w[1]); \ | |
586 __add_carry_in_out((P).w[2],CY,P1.w[1],P0.w[2],CY); \ | |
587 __add_carry_in_out((P).w[3],CY,P1.w[2],P0.w[3],CY); \ | |
588 __add_carry_in_out((P).w[4],CY,P1.w[3],P0.w[4],CY); \ | |
589 (P).w[5] = P1.w[4] + CY; \ | |
590 __add_carry_out((P).w[2],CY,P2.w[0],(P).w[2]); \ | |
591 __add_carry_in_out((P).w[3],CY,P2.w[1],(P).w[3],CY); \ | |
592 __add_carry_in_out((P).w[4],CY,P2.w[2],(P).w[4],CY); \ | |
593 __add_carry_in_out((P).w[5],CY,P2.w[3],(P).w[5],CY); \ | |
594 (P).w[6] = P2.w[4] + CY; \ | |
595 } | |
596 #define __mul_320x320_to_640(P, A, B) \ | |
597 { \ | |
598 UINT512 P0,P1,P2,P3; \ | |
599 UINT64 CY; \ | |
600 __mul_256x256_to_512((P), (A), B); \ | |
601 __mul_64x256_to_320(P1, (A).w[4], B); \ | |
602 __mul_64x256_to_320(P2, (B).w[4], A); \ | |
603 __mul_64x64_to_128(P3, (A).w[4], (B).w[4]); \ | |
604 __add_carry_out((P0).w[0],CY,P1.w[0],P2.w[0]); \ | |
605 __add_carry_in_out((P0).w[1],CY,P1.w[1],P2.w[1],CY); \ | |
606 __add_carry_in_out((P0).w[2],CY,P1.w[2],P2.w[2],CY); \ | |
607 __add_carry_in_out((P0).w[3],CY,P1.w[3],P2.w[3],CY); \ | |
608 __add_carry_in_out((P0).w[4],CY,P1.w[4],P2.w[4],CY); \ | |
609 P3.w[1] += CY; \ | |
610 __add_carry_out((P).w[4],CY,(P).w[4],P0.w[0]); \ | |
611 __add_carry_in_out((P).w[5],CY,(P).w[5],P0.w[1],CY); \ | |
612 __add_carry_in_out((P).w[6],CY,(P).w[6],P0.w[2],CY); \ | |
613 __add_carry_in_out((P).w[7],CY,(P).w[7],P0.w[3],CY); \ | |
614 __add_carry_in_out((P).w[8],CY,P3.w[0],P0.w[4],CY); \ | |
615 (P).w[9] = P3.w[1] + CY; \ | |
616 } | |
617 #define __mul_384x384_to_768(P, A, B) \ | |
618 { \ | |
619 UINT512 P0,P1,P2,P3; \ | |
620 UINT64 CY; \ | |
621 __mul_320x320_to_640((P), (A), B); \ | |
622 __mul_64x320_to_384(P1, (A).w[5], B); \ | |
623 __mul_64x320_to_384(P2, (B).w[5], A); \ | |
624 __mul_64x64_to_128(P3, (A).w[5], (B).w[5]); \ | |
625 __add_carry_out((P0).w[0],CY,P1.w[0],P2.w[0]); \ | |
626 __add_carry_in_out((P0).w[1],CY,P1.w[1],P2.w[1],CY); \ | |
627 __add_carry_in_out((P0).w[2],CY,P1.w[2],P2.w[2],CY); \ | |
628 __add_carry_in_out((P0).w[3],CY,P1.w[3],P2.w[3],CY); \ | |
629 __add_carry_in_out((P0).w[4],CY,P1.w[4],P2.w[4],CY); \ | |
630 __add_carry_in_out((P0).w[5],CY,P1.w[5],P2.w[5],CY); \ | |
631 P3.w[1] += CY; \ | |
632 __add_carry_out((P).w[5],CY,(P).w[5],P0.w[0]); \ | |
633 __add_carry_in_out((P).w[6],CY,(P).w[6],P0.w[1],CY); \ | |
634 __add_carry_in_out((P).w[7],CY,(P).w[7],P0.w[2],CY); \ | |
635 __add_carry_in_out((P).w[8],CY,(P).w[8],P0.w[3],CY); \ | |
636 __add_carry_in_out((P).w[9],CY,(P).w[9],P0.w[4],CY); \ | |
637 __add_carry_in_out((P).w[10],CY,P3.w[0],P0.w[5],CY); \ | |
638 (P).w[11] = P3.w[1] + CY; \ | |
639 } | |
640 #define __mul_64x128_short(Ql, A, B) \ | |
641 { \ | |
642 UINT64 ALBH_L; \ | |
643 \ | |
644 __mul_64x64_to_64(ALBH_L, (A),(B).w[1]); \ | |
645 __mul_64x64_to_128((Ql), (A), (B).w[0]); \ | |
646 \ | |
647 (Ql).w[1] += ALBH_L; \ | |
648 } | |
649 #define __scale128_10(D,_TMP) \ | |
650 { \ | |
651 UINT128 _TMP2,_TMP8; \ | |
652 _TMP2.w[1] = (_TMP.w[1]<<1)|(_TMP.w[0]>>63); \ | |
653 _TMP2.w[0] = _TMP.w[0]<<1; \ | |
654 _TMP8.w[1] = (_TMP.w[1]<<3)|(_TMP.w[0]>>61); \ | |
655 _TMP8.w[0] = _TMP.w[0]<<3; \ | |
656 __add_128_128(D, _TMP2, _TMP8); \ | |
657 } | |
658 // 64x64-bit product | |
659 #define __mul_64x64_to_128MACH(P128, CX64, CY64) \ | |
660 { \ | |
661 UINT64 CXH,CXL,CYH,CYL,PL,PH,PM,PM2; \ | |
662 CXH = (CX64) >> 32; \ | |
663 CXL = (UINT32)(CX64); \ | |
664 CYH = (CY64) >> 32; \ | |
665 CYL = (UINT32)(CY64); \ | |
666 PM = CXH*CYL; \ | |
667 PH = CXH*CYH; \ | |
668 PL = CXL*CYL; \ | |
669 PM2 = CXL*CYH; \ | |
670 PH += (PM>>32); \ | |
671 PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \ | |
672 (P128).w[1] = PH + (PM>>32); \ | |
673 (P128).w[0] = (PM<<32)+(UINT32)PL; \ | |
674 } | |
675 // 64x64-bit product | |
676 #define __mul_64x64_to_128HIGH(P64, CX64, CY64) \ | |
677 { \ | |
678 UINT64 CXH,CXL,CYH,CYL,PL,PH,PM,PM2; \ | |
679 CXH = (CX64) >> 32; \ | |
680 CXL = (UINT32)(CX64); \ | |
681 CYH = (CY64) >> 32; \ | |
682 CYL = (UINT32)(CY64); \ | |
683 PM = CXH*CYL; \ | |
684 PH = CXH*CYH; \ | |
685 PL = CXL*CYL; \ | |
686 PM2 = CXL*CYH; \ | |
687 PH += (PM>>32); \ | |
688 PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \ | |
689 P64 = PH + (PM>>32); \ | |
690 } | |
691 #define __mul_128x64_to_128(Q128, A64, B128) \ | |
692 { \ | |
693 UINT64 ALBH_L; \ | |
694 ALBH_L = (A64) * (B128).w[1]; \ | |
695 __mul_64x64_to_128MACH((Q128), (A64), (B128).w[0]); \ | |
696 (Q128).w[1] += ALBH_L; \ | |
697 } | |
698 // might simplify by calculating just QM2.w[0] | |
699 #define __mul_64x128_to_128(Ql, A, B) \ | |
700 { \ | |
701 UINT128 ALBL, ALBH, QM2; \ | |
702 __mul_64x64_to_128(ALBH, (A), (B).w[1]); \ | |
703 __mul_64x64_to_128(ALBL, (A), (B).w[0]); \ | |
704 (Ql).w[0] = ALBL.w[0]; \ | |
705 __add_128_64(QM2, ALBH, ALBL.w[1]); \ | |
706 (Ql).w[1] = QM2.w[0]; \ | |
707 } | |
708 /********************************************************************* | |
709 * | |
710 * BID Pack/Unpack Macros | |
711 * | |
712 *********************************************************************/ | |
713 ///////////////////////////////////////// | |
714 // BID64 definitions | |
715 //////////////////////////////////////// | |
716 #define DECIMAL_MAX_EXPON_64 767 | |
717 #define DECIMAL_EXPONENT_BIAS 398 | |
718 #define MAX_FORMAT_DIGITS 16 | |
719 ///////////////////////////////////////// | |
720 // BID128 definitions | |
721 //////////////////////////////////////// | |
722 #define DECIMAL_MAX_EXPON_128 12287 | |
723 #define DECIMAL_EXPONENT_BIAS_128 6176 | |
724 #define MAX_FORMAT_DIGITS_128 34 | |
725 ///////////////////////////////////////// | |
726 // BID32 definitions | |
727 //////////////////////////////////////// | |
728 #define DECIMAL_MAX_EXPON_32 191 | |
729 #define DECIMAL_EXPONENT_BIAS_32 101 | |
730 #define MAX_FORMAT_DIGITS_32 7 | |
731 //////////////////////////////////////// | |
732 // Constant Definitions | |
733 /////////////////////////////////////// | |
734 #define SPECIAL_ENCODING_MASK64 0x6000000000000000ull | |
735 #define INFINITY_MASK64 0x7800000000000000ull | |
736 #define SINFINITY_MASK64 0xf800000000000000ull | |
737 #define SSNAN_MASK64 0xfc00000000000000ull | |
738 #define NAN_MASK64 0x7c00000000000000ull | |
739 #define SNAN_MASK64 0x7e00000000000000ull | |
740 #define QUIET_MASK64 0xfdffffffffffffffull | |
741 #define LARGE_COEFF_MASK64 0x0007ffffffffffffull | |
742 #define LARGE_COEFF_HIGH_BIT64 0x0020000000000000ull | |
743 #define SMALL_COEFF_MASK64 0x001fffffffffffffull | |
744 #define EXPONENT_MASK64 0x3ff | |
745 #define EXPONENT_SHIFT_LARGE64 51 | |
746 #define EXPONENT_SHIFT_SMALL64 53 | |
747 #define LARGEST_BID64 0x77fb86f26fc0ffffull | |
748 #define SMALLEST_BID64 0xf7fb86f26fc0ffffull | |
749 #define SMALL_COEFF_MASK128 0x0001ffffffffffffull | |
750 #define LARGE_COEFF_MASK128 0x00007fffffffffffull | |
751 #define EXPONENT_MASK128 0x3fff | |
752 #define LARGEST_BID128_HIGH 0x5fffed09bead87c0ull | |
753 #define LARGEST_BID128_LOW 0x378d8e63ffffffffull | |
754 #define SPECIAL_ENCODING_MASK32 0x60000000ul | |
755 #define INFINITY_MASK32 0x78000000ul | |
756 #define LARGE_COEFF_MASK32 0x007ffffful | |
757 #define LARGE_COEFF_HIGH_BIT32 0x00800000ul | |
758 #define SMALL_COEFF_MASK32 0x001ffffful | |
759 #define EXPONENT_MASK32 0xff | |
760 #define LARGEST_BID32 0x77f8967f | |
761 #define NAN_MASK32 0x7c000000 | |
762 #define SNAN_MASK32 0x7e000000 | |
763 #define MASK_BINARY_EXPONENT 0x7ff0000000000000ull | |
764 #define BINARY_EXPONENT_BIAS 0x3ff | |
765 #define UPPER_EXPON_LIMIT 51 | |
766 // data needed for BID pack/unpack macros | |
767 extern UINT64 round_const_table[][19]; | |
768 extern UINT128 reciprocals10_128[]; | |
769 extern int recip_scale[]; | |
770 extern UINT128 power10_table_128[]; | |
771 extern int estimate_decimal_digits[]; | |
772 extern int estimate_bin_expon[]; | |
773 extern UINT64 power10_index_binexp[]; | |
774 extern int short_recip_scale[]; | |
775 extern UINT64 reciprocals10_64[]; | |
776 extern UINT128 power10_index_binexp_128[]; | |
777 extern UINT128 round_const_table_128[][36]; | |
778 | |
779 | |
780 ////////////////////////////////////////////// | |
781 // Status Flag Handling | |
782 ///////////////////////////////////////////// | |
783 #define __set_status_flags(fpsc, status) *(fpsc) |= status | |
784 #define is_inexact(fpsc) ((*(fpsc))&INEXACT_EXCEPTION) | |
785 | |
786 __BID_INLINE__ UINT64 | |
787 unpack_BID64 (UINT64 * psign_x, int *pexponent_x, | |
788 UINT64 * pcoefficient_x, UINT64 x) { | |
789 UINT64 tmp, coeff; | |
790 | |
791 *psign_x = x & 0x8000000000000000ull; | |
792 | |
793 if ((x & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) { | |
794 // special encodings | |
795 // coefficient | |
796 coeff = (x & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64; | |
797 | |
798 if ((x & INFINITY_MASK64) == INFINITY_MASK64) { | |
799 *pexponent_x = 0; | |
800 *pcoefficient_x = x & 0xfe03ffffffffffffull; | |
801 if ((x & 0x0003ffffffffffffull) >= 1000000000000000ull) | |
802 *pcoefficient_x = x & 0xfe00000000000000ull; | |
803 if ((x & NAN_MASK64) == INFINITY_MASK64) | |
804 *pcoefficient_x = x & SINFINITY_MASK64; | |
805 return 0; // NaN or Infinity | |
806 } | |
807 // check for non-canonical values | |
808 if (coeff >= 10000000000000000ull) | |
809 coeff = 0; | |
810 *pcoefficient_x = coeff; | |
811 // get exponent | |
812 tmp = x >> EXPONENT_SHIFT_LARGE64; | |
813 *pexponent_x = (int) (tmp & EXPONENT_MASK64); | |
814 return coeff; | |
815 } | |
816 // exponent | |
817 tmp = x >> EXPONENT_SHIFT_SMALL64; | |
818 *pexponent_x = (int) (tmp & EXPONENT_MASK64); | |
819 // coefficient | |
820 *pcoefficient_x = (x & SMALL_COEFF_MASK64); | |
821 | |
822 return *pcoefficient_x; | |
823 } | |
824 | |
825 // | |
826 // BID64 pack macro (general form) | |
827 // | |
828 __BID_INLINE__ UINT64 | |
829 get_BID64 (UINT64 sgn, int expon, UINT64 coeff, int rmode, | |
830 unsigned *fpsc) { | |
831 UINT128 Stemp, Q_low; | |
832 UINT64 QH, r, mask, C64, remainder_h, CY, carry; | |
833 int extra_digits, amount, amount2; | |
834 unsigned status; | |
835 | |
836 if (coeff > 9999999999999999ull) { | |
837 expon++; | |
838 coeff = 1000000000000000ull; | |
839 } | |
840 // check for possible underflow/overflow | |
841 if (((unsigned) expon) >= 3 * 256) { | |
842 if (expon < 0) { | |
843 // underflow | |
844 if (expon + MAX_FORMAT_DIGITS < 0) { | |
845 #ifdef SET_STATUS_FLAGS | |
846 __set_status_flags (fpsc, | |
847 UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
848 #endif | |
849 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
850 #ifndef IEEE_ROUND_NEAREST | |
851 if (rmode == ROUNDING_DOWN && sgn) | |
852 return 0x8000000000000001ull; | |
853 if (rmode == ROUNDING_UP && !sgn) | |
854 return 1ull; | |
855 #endif | |
856 #endif | |
857 // result is 0 | |
858 return sgn; | |
859 } | |
860 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
861 #ifndef IEEE_ROUND_NEAREST | |
862 if (sgn && (unsigned) (rmode - 1) < 2) | |
863 rmode = 3 - rmode; | |
864 #endif | |
865 #endif | |
866 // get digits to be shifted out | |
867 extra_digits = -expon; | |
868 coeff += round_const_table[rmode][extra_digits]; | |
869 | |
870 // get coeff*(2^M[extra_digits])/10^extra_digits | |
871 __mul_64x128_full (QH, Q_low, coeff, | |
872 reciprocals10_128[extra_digits]); | |
873 | |
874 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 | |
875 amount = recip_scale[extra_digits]; | |
876 | |
877 C64 = QH >> amount; | |
878 | |
879 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
880 #ifndef IEEE_ROUND_NEAREST | |
881 if (rmode == 0) //ROUNDING_TO_NEAREST | |
882 #endif | |
883 if (C64 & 1) { | |
884 // check whether fractional part of initial_P/10^extra_digits is exactly .5 | |
885 | |
886 // get remainder | |
887 amount2 = 64 - amount; | |
888 remainder_h = 0; | |
889 remainder_h--; | |
890 remainder_h >>= amount2; | |
891 remainder_h = remainder_h & QH; | |
892 | |
893 if (!remainder_h | |
894 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
895 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
896 && Q_low.w[0] < | |
897 reciprocals10_128[extra_digits].w[0]))) { | |
898 C64--; | |
899 } | |
900 } | |
901 #endif | |
902 | |
903 #ifdef SET_STATUS_FLAGS | |
904 | |
905 if (is_inexact (fpsc)) | |
906 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
907 else { | |
908 status = INEXACT_EXCEPTION; | |
909 // get remainder | |
910 remainder_h = QH << (64 - amount); | |
911 | |
912 switch (rmode) { | |
913 case ROUNDING_TO_NEAREST: | |
914 case ROUNDING_TIES_AWAY: | |
915 // test whether fractional part is 0 | |
916 if (remainder_h == 0x8000000000000000ull | |
917 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
918 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
919 && Q_low.w[0] < | |
920 reciprocals10_128[extra_digits].w[0]))) | |
921 status = EXACT_STATUS; | |
922 break; | |
923 case ROUNDING_DOWN: | |
924 case ROUNDING_TO_ZERO: | |
925 if (!remainder_h | |
926 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
927 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
928 && Q_low.w[0] < | |
929 reciprocals10_128[extra_digits].w[0]))) | |
930 status = EXACT_STATUS; | |
931 break; | |
932 default: | |
933 // round up | |
934 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], | |
935 reciprocals10_128[extra_digits].w[0]); | |
936 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], | |
937 reciprocals10_128[extra_digits].w[1], CY); | |
938 if ((remainder_h >> (64 - amount)) + carry >= | |
939 (((UINT64) 1) << amount)) | |
940 status = EXACT_STATUS; | |
941 } | |
942 | |
943 if (status != EXACT_STATUS) | |
944 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
945 } | |
946 | |
947 #endif | |
948 | |
949 return sgn | C64; | |
950 } | |
951 while (coeff < 1000000000000000ull && expon >= 3 * 256) { | |
952 expon--; | |
953 coeff = (coeff << 3) + (coeff << 1); | |
954 } | |
955 if (expon > DECIMAL_MAX_EXPON_64) { | |
956 #ifdef SET_STATUS_FLAGS | |
957 __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
958 #endif | |
959 // overflow | |
960 r = sgn | INFINITY_MASK64; | |
961 switch (rmode) { | |
962 case ROUNDING_DOWN: | |
963 if (!sgn) | |
964 r = LARGEST_BID64; | |
965 break; | |
966 case ROUNDING_TO_ZERO: | |
967 r = sgn | LARGEST_BID64; | |
968 break; | |
969 case ROUNDING_UP: | |
970 // round up | |
971 if (sgn) | |
972 r = SMALLEST_BID64; | |
973 } | |
974 return r; | |
975 } | |
976 } | |
977 | |
978 mask = 1; | |
979 mask <<= EXPONENT_SHIFT_SMALL64; | |
980 | |
981 // check whether coefficient fits in 10*5+3 bits | |
982 if (coeff < mask) { | |
983 r = expon; | |
984 r <<= EXPONENT_SHIFT_SMALL64; | |
985 r |= (coeff | sgn); | |
986 return r; | |
987 } | |
988 // special format | |
989 | |
990 // eliminate the case coeff==10^16 after rounding | |
991 if (coeff == 10000000000000000ull) { | |
992 r = expon + 1; | |
993 r <<= EXPONENT_SHIFT_SMALL64; | |
994 r |= (1000000000000000ull | sgn); | |
995 return r; | |
996 } | |
997 | |
998 r = expon; | |
999 r <<= EXPONENT_SHIFT_LARGE64; | |
1000 r |= (sgn | SPECIAL_ENCODING_MASK64); | |
1001 // add coeff, without leading bits | |
1002 mask = (mask >> 2) - 1; | |
1003 coeff &= mask; | |
1004 r |= coeff; | |
1005 | |
1006 return r; | |
1007 } | |
1008 | |
1009 | |
1010 | |
1011 | |
1012 // | |
1013 // No overflow/underflow checking | |
1014 // | |
1015 __BID_INLINE__ UINT64 | |
1016 fast_get_BID64 (UINT64 sgn, int expon, UINT64 coeff) { | |
1017 UINT64 r, mask; | |
1018 | |
1019 mask = 1; | |
1020 mask <<= EXPONENT_SHIFT_SMALL64; | |
1021 | |
1022 // check whether coefficient fits in 10*5+3 bits | |
1023 if (coeff < mask) { | |
1024 r = expon; | |
1025 r <<= EXPONENT_SHIFT_SMALL64; | |
1026 r |= (coeff | sgn); | |
1027 return r; | |
1028 } | |
1029 // special format | |
1030 | |
1031 // eliminate the case coeff==10^16 after rounding | |
1032 if (coeff == 10000000000000000ull) { | |
1033 r = expon + 1; | |
1034 r <<= EXPONENT_SHIFT_SMALL64; | |
1035 r |= (1000000000000000ull | sgn); | |
1036 return r; | |
1037 } | |
1038 | |
1039 r = expon; | |
1040 r <<= EXPONENT_SHIFT_LARGE64; | |
1041 r |= (sgn | SPECIAL_ENCODING_MASK64); | |
1042 // add coeff, without leading bits | |
1043 mask = (mask >> 2) - 1; | |
1044 coeff &= mask; | |
1045 r |= coeff; | |
1046 | |
1047 return r; | |
1048 } | |
1049 | |
1050 | |
1051 // | |
1052 // no underflow checking | |
1053 // | |
1054 __BID_INLINE__ UINT64 | |
1055 fast_get_BID64_check_OF (UINT64 sgn, int expon, UINT64 coeff, int rmode, | |
1056 unsigned *fpsc) { | |
1057 UINT64 r, mask; | |
1058 | |
1059 if (((unsigned) expon) >= 3 * 256 - 1) { | |
1060 if ((expon == 3 * 256 - 1) && coeff == 10000000000000000ull) { | |
1061 expon = 3 * 256; | |
1062 coeff = 1000000000000000ull; | |
1063 } | |
1064 | |
1065 if (((unsigned) expon) >= 3 * 256) { | |
1066 while (coeff < 1000000000000000ull && expon >= 3 * 256) { | |
1067 expon--; | |
1068 coeff = (coeff << 3) + (coeff << 1); | |
1069 } | |
1070 if (expon > DECIMAL_MAX_EXPON_64) { | |
1071 #ifdef SET_STATUS_FLAGS | |
1072 __set_status_flags (fpsc, | |
1073 OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1074 #endif | |
1075 // overflow | |
1076 r = sgn | INFINITY_MASK64; | |
1077 switch (rmode) { | |
1078 case ROUNDING_DOWN: | |
1079 if (!sgn) | |
1080 r = LARGEST_BID64; | |
1081 break; | |
1082 case ROUNDING_TO_ZERO: | |
1083 r = sgn | LARGEST_BID64; | |
1084 break; | |
1085 case ROUNDING_UP: | |
1086 // round up | |
1087 if (sgn) | |
1088 r = SMALLEST_BID64; | |
1089 } | |
1090 return r; | |
1091 } | |
1092 } | |
1093 } | |
1094 | |
1095 mask = 1; | |
1096 mask <<= EXPONENT_SHIFT_SMALL64; | |
1097 | |
1098 // check whether coefficient fits in 10*5+3 bits | |
1099 if (coeff < mask) { | |
1100 r = expon; | |
1101 r <<= EXPONENT_SHIFT_SMALL64; | |
1102 r |= (coeff | sgn); | |
1103 return r; | |
1104 } | |
1105 // special format | |
1106 | |
1107 // eliminate the case coeff==10^16 after rounding | |
1108 if (coeff == 10000000000000000ull) { | |
1109 r = expon + 1; | |
1110 r <<= EXPONENT_SHIFT_SMALL64; | |
1111 r |= (1000000000000000ull | sgn); | |
1112 return r; | |
1113 } | |
1114 | |
1115 r = expon; | |
1116 r <<= EXPONENT_SHIFT_LARGE64; | |
1117 r |= (sgn | SPECIAL_ENCODING_MASK64); | |
1118 // add coeff, without leading bits | |
1119 mask = (mask >> 2) - 1; | |
1120 coeff &= mask; | |
1121 r |= coeff; | |
1122 | |
1123 return r; | |
1124 } | |
1125 | |
1126 | |
1127 // | |
1128 // No overflow/underflow checking | |
1129 // or checking for coefficients equal to 10^16 (after rounding) | |
1130 // | |
1131 __BID_INLINE__ UINT64 | |
1132 very_fast_get_BID64 (UINT64 sgn, int expon, UINT64 coeff) { | |
1133 UINT64 r, mask; | |
1134 | |
1135 mask = 1; | |
1136 mask <<= EXPONENT_SHIFT_SMALL64; | |
1137 | |
1138 // check whether coefficient fits in 10*5+3 bits | |
1139 if (coeff < mask) { | |
1140 r = expon; | |
1141 r <<= EXPONENT_SHIFT_SMALL64; | |
1142 r |= (coeff | sgn); | |
1143 return r; | |
1144 } | |
1145 // special format | |
1146 r = expon; | |
1147 r <<= EXPONENT_SHIFT_LARGE64; | |
1148 r |= (sgn | SPECIAL_ENCODING_MASK64); | |
1149 // add coeff, without leading bits | |
1150 mask = (mask >> 2) - 1; | |
1151 coeff &= mask; | |
1152 r |= coeff; | |
1153 | |
1154 return r; | |
1155 } | |
1156 | |
1157 // | |
1158 // No overflow/underflow checking or checking for coefficients above 2^53 | |
1159 // | |
1160 __BID_INLINE__ UINT64 | |
1161 very_fast_get_BID64_small_mantissa (UINT64 sgn, int expon, UINT64 coeff) { | |
1162 // no UF/OF | |
1163 UINT64 r; | |
1164 | |
1165 r = expon; | |
1166 r <<= EXPONENT_SHIFT_SMALL64; | |
1167 r |= (coeff | sgn); | |
1168 return r; | |
1169 } | |
1170 | |
1171 | |
1172 // | |
1173 // This pack macro is used when underflow is known to occur | |
1174 // | |
1175 __BID_INLINE__ UINT64 | |
1176 get_BID64_UF (UINT64 sgn, int expon, UINT64 coeff, UINT64 R, int rmode, | |
1177 unsigned *fpsc) { | |
1178 UINT128 C128, Q_low, Stemp; | |
1179 UINT64 C64, remainder_h, QH, carry, CY; | |
1180 int extra_digits, amount, amount2; | |
1181 unsigned status; | |
1182 | |
1183 // underflow | |
1184 if (expon + MAX_FORMAT_DIGITS < 0) { | |
1185 #ifdef SET_STATUS_FLAGS | |
1186 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1187 #endif | |
1188 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1189 #ifndef IEEE_ROUND_NEAREST | |
1190 if (rmode == ROUNDING_DOWN && sgn) | |
1191 return 0x8000000000000001ull; | |
1192 if (rmode == ROUNDING_UP && !sgn) | |
1193 return 1ull; | |
1194 #endif | |
1195 #endif | |
1196 // result is 0 | |
1197 return sgn; | |
1198 } | |
1199 // 10*coeff | |
1200 coeff = (coeff << 3) + (coeff << 1); | |
1201 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1202 #ifndef IEEE_ROUND_NEAREST | |
1203 if (sgn && (unsigned) (rmode - 1) < 2) | |
1204 rmode = 3 - rmode; | |
1205 #endif | |
1206 #endif | |
1207 if (R) | |
1208 coeff |= 1; | |
1209 // get digits to be shifted out | |
1210 extra_digits = 1 - expon; | |
1211 C128.w[0] = coeff + round_const_table[rmode][extra_digits]; | |
1212 | |
1213 // get coeff*(2^M[extra_digits])/10^extra_digits | |
1214 __mul_64x128_full (QH, Q_low, C128.w[0], | |
1215 reciprocals10_128[extra_digits]); | |
1216 | |
1217 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 | |
1218 amount = recip_scale[extra_digits]; | |
1219 | |
1220 C64 = QH >> amount; | |
1221 //__shr_128(C128, Q_high, amount); | |
1222 | |
1223 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1224 #ifndef IEEE_ROUND_NEAREST | |
1225 if (rmode == 0) //ROUNDING_TO_NEAREST | |
1226 #endif | |
1227 if (C64 & 1) { | |
1228 // check whether fractional part of initial_P/10^extra_digits is exactly .5 | |
1229 | |
1230 // get remainder | |
1231 amount2 = 64 - amount; | |
1232 remainder_h = 0; | |
1233 remainder_h--; | |
1234 remainder_h >>= amount2; | |
1235 remainder_h = remainder_h & QH; | |
1236 | |
1237 if (!remainder_h | |
1238 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1239 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1240 && Q_low.w[0] < | |
1241 reciprocals10_128[extra_digits].w[0]))) { | |
1242 C64--; | |
1243 } | |
1244 } | |
1245 #endif | |
1246 | |
1247 #ifdef SET_STATUS_FLAGS | |
1248 | |
1249 if (is_inexact (fpsc)) | |
1250 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
1251 else { | |
1252 status = INEXACT_EXCEPTION; | |
1253 // get remainder | |
1254 remainder_h = QH << (64 - amount); | |
1255 | |
1256 switch (rmode) { | |
1257 case ROUNDING_TO_NEAREST: | |
1258 case ROUNDING_TIES_AWAY: | |
1259 // test whether fractional part is 0 | |
1260 if (remainder_h == 0x8000000000000000ull | |
1261 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1262 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1263 && Q_low.w[0] < | |
1264 reciprocals10_128[extra_digits].w[0]))) | |
1265 status = EXACT_STATUS; | |
1266 break; | |
1267 case ROUNDING_DOWN: | |
1268 case ROUNDING_TO_ZERO: | |
1269 if (!remainder_h | |
1270 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1271 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1272 && Q_low.w[0] < | |
1273 reciprocals10_128[extra_digits].w[0]))) | |
1274 status = EXACT_STATUS; | |
1275 break; | |
1276 default: | |
1277 // round up | |
1278 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], | |
1279 reciprocals10_128[extra_digits].w[0]); | |
1280 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], | |
1281 reciprocals10_128[extra_digits].w[1], CY); | |
1282 if ((remainder_h >> (64 - amount)) + carry >= | |
1283 (((UINT64) 1) << amount)) | |
1284 status = EXACT_STATUS; | |
1285 } | |
1286 | |
1287 if (status != EXACT_STATUS) | |
1288 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
1289 } | |
1290 | |
1291 #endif | |
1292 | |
1293 return sgn | C64; | |
1294 | |
1295 } | |
1296 | |
1297 | |
1298 | |
1299 // | |
1300 // This pack macro doesnot check for coefficients above 2^53 | |
1301 // | |
1302 __BID_INLINE__ UINT64 | |
1303 get_BID64_small_mantissa (UINT64 sgn, int expon, UINT64 coeff, | |
1304 int rmode, unsigned *fpsc) { | |
1305 UINT128 C128, Q_low, Stemp; | |
1306 UINT64 r, mask, C64, remainder_h, QH, carry, CY; | |
1307 int extra_digits, amount, amount2; | |
1308 unsigned status; | |
1309 | |
1310 // check for possible underflow/overflow | |
1311 if (((unsigned) expon) >= 3 * 256) { | |
1312 if (expon < 0) { | |
1313 // underflow | |
1314 if (expon + MAX_FORMAT_DIGITS < 0) { | |
1315 #ifdef SET_STATUS_FLAGS | |
1316 __set_status_flags (fpsc, | |
1317 UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1318 #endif | |
1319 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1320 #ifndef IEEE_ROUND_NEAREST | |
1321 if (rmode == ROUNDING_DOWN && sgn) | |
1322 return 0x8000000000000001ull; | |
1323 if (rmode == ROUNDING_UP && !sgn) | |
1324 return 1ull; | |
1325 #endif | |
1326 #endif | |
1327 // result is 0 | |
1328 return sgn; | |
1329 } | |
1330 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1331 #ifndef IEEE_ROUND_NEAREST | |
1332 if (sgn && (unsigned) (rmode - 1) < 2) | |
1333 rmode = 3 - rmode; | |
1334 #endif | |
1335 #endif | |
1336 // get digits to be shifted out | |
1337 extra_digits = -expon; | |
1338 C128.w[0] = coeff + round_const_table[rmode][extra_digits]; | |
1339 | |
1340 // get coeff*(2^M[extra_digits])/10^extra_digits | |
1341 __mul_64x128_full (QH, Q_low, C128.w[0], | |
1342 reciprocals10_128[extra_digits]); | |
1343 | |
1344 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 | |
1345 amount = recip_scale[extra_digits]; | |
1346 | |
1347 C64 = QH >> amount; | |
1348 | |
1349 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1350 #ifndef IEEE_ROUND_NEAREST | |
1351 if (rmode == 0) //ROUNDING_TO_NEAREST | |
1352 #endif | |
1353 if (C64 & 1) { | |
1354 // check whether fractional part of initial_P/10^extra_digits is exactly .5 | |
1355 | |
1356 // get remainder | |
1357 amount2 = 64 - amount; | |
1358 remainder_h = 0; | |
1359 remainder_h--; | |
1360 remainder_h >>= amount2; | |
1361 remainder_h = remainder_h & QH; | |
1362 | |
1363 if (!remainder_h | |
1364 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1365 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1366 && Q_low.w[0] < | |
1367 reciprocals10_128[extra_digits].w[0]))) { | |
1368 C64--; | |
1369 } | |
1370 } | |
1371 #endif | |
1372 | |
1373 #ifdef SET_STATUS_FLAGS | |
1374 | |
1375 if (is_inexact (fpsc)) | |
1376 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
1377 else { | |
1378 status = INEXACT_EXCEPTION; | |
1379 // get remainder | |
1380 remainder_h = QH << (64 - amount); | |
1381 | |
1382 switch (rmode) { | |
1383 case ROUNDING_TO_NEAREST: | |
1384 case ROUNDING_TIES_AWAY: | |
1385 // test whether fractional part is 0 | |
1386 if (remainder_h == 0x8000000000000000ull | |
1387 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1388 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1389 && Q_low.w[0] < | |
1390 reciprocals10_128[extra_digits].w[0]))) | |
1391 status = EXACT_STATUS; | |
1392 break; | |
1393 case ROUNDING_DOWN: | |
1394 case ROUNDING_TO_ZERO: | |
1395 if (!remainder_h | |
1396 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] | |
1397 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
1398 && Q_low.w[0] < | |
1399 reciprocals10_128[extra_digits].w[0]))) | |
1400 status = EXACT_STATUS; | |
1401 break; | |
1402 default: | |
1403 // round up | |
1404 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], | |
1405 reciprocals10_128[extra_digits].w[0]); | |
1406 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], | |
1407 reciprocals10_128[extra_digits].w[1], CY); | |
1408 if ((remainder_h >> (64 - amount)) + carry >= | |
1409 (((UINT64) 1) << amount)) | |
1410 status = EXACT_STATUS; | |
1411 } | |
1412 | |
1413 if (status != EXACT_STATUS) | |
1414 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
1415 } | |
1416 | |
1417 #endif | |
1418 | |
1419 return sgn | C64; | |
1420 } | |
1421 | |
1422 while (coeff < 1000000000000000ull && expon >= 3 * 256) { | |
1423 expon--; | |
1424 coeff = (coeff << 3) + (coeff << 1); | |
1425 } | |
1426 if (expon > DECIMAL_MAX_EXPON_64) { | |
1427 #ifdef SET_STATUS_FLAGS | |
1428 __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1429 #endif | |
1430 // overflow | |
1431 r = sgn | INFINITY_MASK64; | |
1432 switch (rmode) { | |
1433 case ROUNDING_DOWN: | |
1434 if (!sgn) | |
1435 r = LARGEST_BID64; | |
1436 break; | |
1437 case ROUNDING_TO_ZERO: | |
1438 r = sgn | LARGEST_BID64; | |
1439 break; | |
1440 case ROUNDING_UP: | |
1441 // round up | |
1442 if (sgn) | |
1443 r = SMALLEST_BID64; | |
1444 } | |
1445 return r; | |
1446 } else { | |
1447 mask = 1; | |
1448 mask <<= EXPONENT_SHIFT_SMALL64; | |
1449 if (coeff >= mask) { | |
1450 r = expon; | |
1451 r <<= EXPONENT_SHIFT_LARGE64; | |
1452 r |= (sgn | SPECIAL_ENCODING_MASK64); | |
1453 // add coeff, without leading bits | |
1454 mask = (mask >> 2) - 1; | |
1455 coeff &= mask; | |
1456 r |= coeff; | |
1457 return r; | |
1458 } | |
1459 } | |
1460 } | |
1461 | |
1462 r = expon; | |
1463 r <<= EXPONENT_SHIFT_SMALL64; | |
1464 r |= (coeff | sgn); | |
1465 | |
1466 return r; | |
1467 } | |
1468 | |
1469 | |
1470 /***************************************************************************** | |
1471 * | |
1472 * BID128 pack/unpack macros | |
1473 * | |
1474 *****************************************************************************/ | |
1475 | |
1476 // | |
1477 // Macro for handling BID128 underflow | |
1478 // sticky bit given as additional argument | |
1479 // | |
1480 __BID_INLINE__ UINT128 * | |
1481 handle_UF_128_rem (UINT128 * pres, UINT64 sgn, int expon, UINT128 CQ, | |
1482 UINT64 R, unsigned *prounding_mode, unsigned *fpsc) { | |
1483 UINT128 T128, TP128, Qh, Ql, Qh1, Stemp, Tmp, Tmp1, CQ2, CQ8; | |
1484 UINT64 carry, CY; | |
1485 int ed2, amount; | |
1486 unsigned rmode, status; | |
1487 | |
1488 // UF occurs | |
1489 if (expon + MAX_FORMAT_DIGITS_128 < 0) { | |
1490 #ifdef SET_STATUS_FLAGS | |
1491 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1492 #endif | |
1493 pres->w[1] = sgn; | |
1494 pres->w[0] = 0; | |
1495 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1496 #ifndef IEEE_ROUND_NEAREST | |
1497 if ((sgn && *prounding_mode == ROUNDING_DOWN) | |
1498 || (!sgn && *prounding_mode == ROUNDING_UP)) | |
1499 pres->w[0] = 1ull; | |
1500 #endif | |
1501 #endif | |
1502 return pres; | |
1503 } | |
1504 // CQ *= 10 | |
1505 CQ2.w[1] = (CQ.w[1] << 1) | (CQ.w[0] >> 63); | |
1506 CQ2.w[0] = CQ.w[0] << 1; | |
1507 CQ8.w[1] = (CQ.w[1] << 3) | (CQ.w[0] >> 61); | |
1508 CQ8.w[0] = CQ.w[0] << 3; | |
1509 __add_128_128 (CQ, CQ2, CQ8); | |
1510 | |
1511 // add remainder | |
1512 if (R) | |
1513 CQ.w[0] |= 1; | |
1514 | |
1515 ed2 = 1 - expon; | |
1516 // add rounding constant to CQ | |
1517 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1518 #ifndef IEEE_ROUND_NEAREST | |
1519 rmode = *prounding_mode; | |
1520 if (sgn && (unsigned) (rmode - 1) < 2) | |
1521 rmode = 3 - rmode; | |
1522 #else | |
1523 rmode = 0; | |
1524 #endif | |
1525 #else | |
1526 rmode = 0; | |
1527 #endif | |
1528 T128 = round_const_table_128[rmode][ed2]; | |
1529 __add_carry_out (CQ.w[0], carry, T128.w[0], CQ.w[0]); | |
1530 CQ.w[1] = CQ.w[1] + T128.w[1] + carry; | |
1531 | |
1532 TP128 = reciprocals10_128[ed2]; | |
1533 __mul_128x128_full (Qh, Ql, CQ, TP128); | |
1534 amount = recip_scale[ed2]; | |
1535 | |
1536 if (amount >= 64) { | |
1537 CQ.w[0] = Qh.w[1] >> (amount - 64); | |
1538 CQ.w[1] = 0; | |
1539 } else { | |
1540 __shr_128 (CQ, Qh, amount); | |
1541 } | |
1542 | |
1543 expon = 0; | |
1544 | |
1545 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1546 #ifndef IEEE_ROUND_NEAREST | |
1547 if (!(*prounding_mode)) | |
1548 #endif | |
1549 if (CQ.w[0] & 1) { | |
1550 // check whether fractional part of initial_P/10^ed1 is exactly .5 | |
1551 | |
1552 // get remainder | |
1553 __shl_128_long (Qh1, Qh, (128 - amount)); | |
1554 | |
1555 if (!Qh1.w[1] && !Qh1.w[0] | |
1556 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1557 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1558 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) { | |
1559 CQ.w[0]--; | |
1560 } | |
1561 } | |
1562 #endif | |
1563 | |
1564 #ifdef SET_STATUS_FLAGS | |
1565 | |
1566 if (is_inexact (fpsc)) | |
1567 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
1568 else { | |
1569 status = INEXACT_EXCEPTION; | |
1570 // get remainder | |
1571 __shl_128_long (Qh1, Qh, (128 - amount)); | |
1572 | |
1573 switch (rmode) { | |
1574 case ROUNDING_TO_NEAREST: | |
1575 case ROUNDING_TIES_AWAY: | |
1576 // test whether fractional part is 0 | |
1577 if (Qh1.w[1] == 0x8000000000000000ull && (!Qh1.w[0]) | |
1578 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1579 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1580 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) | |
1581 status = EXACT_STATUS; | |
1582 break; | |
1583 case ROUNDING_DOWN: | |
1584 case ROUNDING_TO_ZERO: | |
1585 if ((!Qh1.w[1]) && (!Qh1.w[0]) | |
1586 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1587 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1588 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) | |
1589 status = EXACT_STATUS; | |
1590 break; | |
1591 default: | |
1592 // round up | |
1593 __add_carry_out (Stemp.w[0], CY, Ql.w[0], | |
1594 reciprocals10_128[ed2].w[0]); | |
1595 __add_carry_in_out (Stemp.w[1], carry, Ql.w[1], | |
1596 reciprocals10_128[ed2].w[1], CY); | |
1597 __shr_128_long (Qh, Qh1, (128 - amount)); | |
1598 Tmp.w[0] = 1; | |
1599 Tmp.w[1] = 0; | |
1600 __shl_128_long (Tmp1, Tmp, amount); | |
1601 Qh.w[0] += carry; | |
1602 if (Qh.w[0] < carry) | |
1603 Qh.w[1]++; | |
1604 if (__unsigned_compare_ge_128 (Qh, Tmp1)) | |
1605 status = EXACT_STATUS; | |
1606 } | |
1607 | |
1608 if (status != EXACT_STATUS) | |
1609 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
1610 } | |
1611 | |
1612 #endif | |
1613 | |
1614 pres->w[1] = sgn | CQ.w[1]; | |
1615 pres->w[0] = CQ.w[0]; | |
1616 | |
1617 return pres; | |
1618 | |
1619 } | |
1620 | |
1621 | |
1622 // | |
1623 // Macro for handling BID128 underflow | |
1624 // | |
1625 __BID_INLINE__ UINT128 * | |
1626 handle_UF_128 (UINT128 * pres, UINT64 sgn, int expon, UINT128 CQ, | |
1627 unsigned *prounding_mode, unsigned *fpsc) { | |
1628 UINT128 T128, TP128, Qh, Ql, Qh1, Stemp, Tmp, Tmp1; | |
1629 UINT64 carry, CY; | |
1630 int ed2, amount; | |
1631 unsigned rmode, status; | |
1632 | |
1633 // UF occurs | |
1634 if (expon + MAX_FORMAT_DIGITS_128 < 0) { | |
1635 #ifdef SET_STATUS_FLAGS | |
1636 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1637 #endif | |
1638 pres->w[1] = sgn; | |
1639 pres->w[0] = 0; | |
1640 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1641 #ifndef IEEE_ROUND_NEAREST | |
1642 if ((sgn && *prounding_mode == ROUNDING_DOWN) | |
1643 || (!sgn && *prounding_mode == ROUNDING_UP)) | |
1644 pres->w[0] = 1ull; | |
1645 #endif | |
1646 #endif | |
1647 return pres; | |
1648 } | |
1649 | |
1650 ed2 = 0 - expon; | |
1651 // add rounding constant to CQ | |
1652 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1653 #ifndef IEEE_ROUND_NEAREST | |
1654 rmode = *prounding_mode; | |
1655 if (sgn && (unsigned) (rmode - 1) < 2) | |
1656 rmode = 3 - rmode; | |
1657 #else | |
1658 rmode = 0; | |
1659 #endif | |
1660 #else | |
1661 rmode = 0; | |
1662 #endif | |
1663 | |
1664 T128 = round_const_table_128[rmode][ed2]; | |
1665 __add_carry_out (CQ.w[0], carry, T128.w[0], CQ.w[0]); | |
1666 CQ.w[1] = CQ.w[1] + T128.w[1] + carry; | |
1667 | |
1668 TP128 = reciprocals10_128[ed2]; | |
1669 __mul_128x128_full (Qh, Ql, CQ, TP128); | |
1670 amount = recip_scale[ed2]; | |
1671 | |
1672 if (amount >= 64) { | |
1673 CQ.w[0] = Qh.w[1] >> (amount - 64); | |
1674 CQ.w[1] = 0; | |
1675 } else { | |
1676 __shr_128 (CQ, Qh, amount); | |
1677 } | |
1678 | |
1679 expon = 0; | |
1680 | |
1681 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1682 #ifndef IEEE_ROUND_NEAREST | |
1683 if (!(*prounding_mode)) | |
1684 #endif | |
1685 if (CQ.w[0] & 1) { | |
1686 // check whether fractional part of initial_P/10^ed1 is exactly .5 | |
1687 | |
1688 // get remainder | |
1689 __shl_128_long (Qh1, Qh, (128 - amount)); | |
1690 | |
1691 if (!Qh1.w[1] && !Qh1.w[0] | |
1692 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1693 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1694 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) { | |
1695 CQ.w[0]--; | |
1696 } | |
1697 } | |
1698 #endif | |
1699 | |
1700 #ifdef SET_STATUS_FLAGS | |
1701 | |
1702 if (is_inexact (fpsc)) | |
1703 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
1704 else { | |
1705 status = INEXACT_EXCEPTION; | |
1706 // get remainder | |
1707 __shl_128_long (Qh1, Qh, (128 - amount)); | |
1708 | |
1709 switch (rmode) { | |
1710 case ROUNDING_TO_NEAREST: | |
1711 case ROUNDING_TIES_AWAY: | |
1712 // test whether fractional part is 0 | |
1713 if (Qh1.w[1] == 0x8000000000000000ull && (!Qh1.w[0]) | |
1714 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1715 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1716 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) | |
1717 status = EXACT_STATUS; | |
1718 break; | |
1719 case ROUNDING_DOWN: | |
1720 case ROUNDING_TO_ZERO: | |
1721 if ((!Qh1.w[1]) && (!Qh1.w[0]) | |
1722 && (Ql.w[1] < reciprocals10_128[ed2].w[1] | |
1723 || (Ql.w[1] == reciprocals10_128[ed2].w[1] | |
1724 && Ql.w[0] < reciprocals10_128[ed2].w[0]))) | |
1725 status = EXACT_STATUS; | |
1726 break; | |
1727 default: | |
1728 // round up | |
1729 __add_carry_out (Stemp.w[0], CY, Ql.w[0], | |
1730 reciprocals10_128[ed2].w[0]); | |
1731 __add_carry_in_out (Stemp.w[1], carry, Ql.w[1], | |
1732 reciprocals10_128[ed2].w[1], CY); | |
1733 __shr_128_long (Qh, Qh1, (128 - amount)); | |
1734 Tmp.w[0] = 1; | |
1735 Tmp.w[1] = 0; | |
1736 __shl_128_long (Tmp1, Tmp, amount); | |
1737 Qh.w[0] += carry; | |
1738 if (Qh.w[0] < carry) | |
1739 Qh.w[1]++; | |
1740 if (__unsigned_compare_ge_128 (Qh, Tmp1)) | |
1741 status = EXACT_STATUS; | |
1742 } | |
1743 | |
1744 if (status != EXACT_STATUS) | |
1745 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
1746 } | |
1747 | |
1748 #endif | |
1749 | |
1750 pres->w[1] = sgn | CQ.w[1]; | |
1751 pres->w[0] = CQ.w[0]; | |
1752 | |
1753 return pres; | |
1754 | |
1755 } | |
1756 | |
1757 | |
1758 | |
1759 // | |
1760 // BID128 unpack, input passed by value | |
1761 // | |
1762 __BID_INLINE__ UINT64 | |
1763 unpack_BID128_value (UINT64 * psign_x, int *pexponent_x, | |
1764 UINT128 * pcoefficient_x, UINT128 x) { | |
1765 UINT128 coeff, T33, T34; | |
1766 UINT64 ex; | |
1767 | |
1768 *psign_x = (x.w[1]) & 0x8000000000000000ull; | |
1769 | |
1770 // special encodings | |
1771 if ((x.w[1] & INFINITY_MASK64) >= SPECIAL_ENCODING_MASK64) { | |
1772 if ((x.w[1] & INFINITY_MASK64) < INFINITY_MASK64) { | |
1773 // non-canonical input | |
1774 pcoefficient_x->w[0] = 0; | |
1775 pcoefficient_x->w[1] = 0; | |
1776 ex = (x.w[1]) >> 47; | |
1777 *pexponent_x = ((int) ex) & EXPONENT_MASK128; | |
1778 return 0; | |
1779 } | |
1780 // 10^33 | |
1781 T33 = power10_table_128[33]; | |
1782 /*coeff.w[0] = x.w[0]; | |
1783 coeff.w[1] = (x.w[1]) & LARGE_COEFF_MASK128; | |
1784 pcoefficient_x->w[0] = x.w[0]; | |
1785 pcoefficient_x->w[1] = x.w[1]; | |
1786 if (__unsigned_compare_ge_128 (coeff, T33)) // non-canonical | |
1787 pcoefficient_x->w[1] &= (~LARGE_COEFF_MASK128); */ | |
1788 | |
1789 pcoefficient_x->w[0] = x.w[0]; | |
1790 pcoefficient_x->w[1] = (x.w[1]) & 0x00003fffffffffffull; | |
1791 if (__unsigned_compare_ge_128 ((*pcoefficient_x), T33)) // non-canonical | |
1792 { | |
1793 pcoefficient_x->w[1] = (x.w[1]) & 0xfe00000000000000ull; | |
1794 pcoefficient_x->w[0] = 0; | |
1795 } else | |
1796 pcoefficient_x->w[1] = (x.w[1]) & 0xfe003fffffffffffull; | |
1797 if ((x.w[1] & NAN_MASK64) == INFINITY_MASK64) { | |
1798 pcoefficient_x->w[0] = 0; | |
1799 pcoefficient_x->w[1] = x.w[1] & SINFINITY_MASK64; | |
1800 } | |
1801 *pexponent_x = 0; | |
1802 return 0; // NaN or Infinity | |
1803 } | |
1804 | |
1805 coeff.w[0] = x.w[0]; | |
1806 coeff.w[1] = (x.w[1]) & SMALL_COEFF_MASK128; | |
1807 | |
1808 // 10^34 | |
1809 T34 = power10_table_128[34]; | |
1810 // check for non-canonical values | |
1811 if (__unsigned_compare_ge_128 (coeff, T34)) | |
1812 coeff.w[0] = coeff.w[1] = 0; | |
1813 | |
1814 pcoefficient_x->w[0] = coeff.w[0]; | |
1815 pcoefficient_x->w[1] = coeff.w[1]; | |
1816 | |
1817 ex = (x.w[1]) >> 49; | |
1818 *pexponent_x = ((int) ex) & EXPONENT_MASK128; | |
1819 | |
1820 return coeff.w[0] | coeff.w[1]; | |
1821 } | |
1822 | |
1823 | |
1824 // | |
1825 // BID128 unpack, input pased by reference | |
1826 // | |
1827 __BID_INLINE__ UINT64 | |
1828 unpack_BID128 (UINT64 * psign_x, int *pexponent_x, | |
1829 UINT128 * pcoefficient_x, UINT128 * px) { | |
1830 UINT128 coeff, T33, T34; | |
1831 UINT64 ex; | |
1832 | |
1833 *psign_x = (px->w[1]) & 0x8000000000000000ull; | |
1834 | |
1835 // special encodings | |
1836 if ((px->w[1] & INFINITY_MASK64) >= SPECIAL_ENCODING_MASK64) { | |
1837 if ((px->w[1] & INFINITY_MASK64) < INFINITY_MASK64) { | |
1838 // non-canonical input | |
1839 pcoefficient_x->w[0] = 0; | |
1840 pcoefficient_x->w[1] = 0; | |
1841 ex = (px->w[1]) >> 47; | |
1842 *pexponent_x = ((int) ex) & EXPONENT_MASK128; | |
1843 return 0; | |
1844 } | |
1845 // 10^33 | |
1846 T33 = power10_table_128[33]; | |
1847 coeff.w[0] = px->w[0]; | |
1848 coeff.w[1] = (px->w[1]) & LARGE_COEFF_MASK128; | |
1849 pcoefficient_x->w[0] = px->w[0]; | |
1850 pcoefficient_x->w[1] = px->w[1]; | |
1851 if (__unsigned_compare_ge_128 (coeff, T33)) { // non-canonical | |
1852 pcoefficient_x->w[1] &= (~LARGE_COEFF_MASK128); | |
1853 pcoefficient_x->w[0] = 0; | |
1854 } | |
1855 *pexponent_x = 0; | |
1856 return 0; // NaN or Infinity | |
1857 } | |
1858 | |
1859 coeff.w[0] = px->w[0]; | |
1860 coeff.w[1] = (px->w[1]) & SMALL_COEFF_MASK128; | |
1861 | |
1862 // 10^34 | |
1863 T34 = power10_table_128[34]; | |
1864 // check for non-canonical values | |
1865 if (__unsigned_compare_ge_128 (coeff, T34)) | |
1866 coeff.w[0] = coeff.w[1] = 0; | |
1867 | |
1868 pcoefficient_x->w[0] = coeff.w[0]; | |
1869 pcoefficient_x->w[1] = coeff.w[1]; | |
1870 | |
1871 ex = (px->w[1]) >> 49; | |
1872 *pexponent_x = ((int) ex) & EXPONENT_MASK128; | |
1873 | |
1874 return coeff.w[0] | coeff.w[1]; | |
1875 } | |
1876 | |
1877 // | |
1878 // Pack macro checks for overflow, but not underflow | |
1879 // | |
1880 __BID_INLINE__ UINT128 * | |
1881 get_BID128_very_fast_OF (UINT128 * pres, UINT64 sgn, int expon, | |
1882 UINT128 coeff, unsigned *prounding_mode, | |
1883 unsigned *fpsc) { | |
1884 UINT128 T; | |
1885 UINT64 tmp, tmp2; | |
1886 | |
1887 if ((unsigned) expon > DECIMAL_MAX_EXPON_128) { | |
1888 | |
1889 if (expon - MAX_FORMAT_DIGITS_128 <= DECIMAL_MAX_EXPON_128) { | |
1890 T = power10_table_128[MAX_FORMAT_DIGITS_128 - 1]; | |
1891 while (__unsigned_compare_gt_128 (T, coeff) | |
1892 && expon > DECIMAL_MAX_EXPON_128) { | |
1893 coeff.w[1] = | |
1894 (coeff.w[1] << 3) + (coeff.w[1] << 1) + (coeff.w[0] >> 61) + | |
1895 (coeff.w[0] >> 63); | |
1896 tmp2 = coeff.w[0] << 3; | |
1897 coeff.w[0] = (coeff.w[0] << 1) + tmp2; | |
1898 if (coeff.w[0] < tmp2) | |
1899 coeff.w[1]++; | |
1900 | |
1901 expon--; | |
1902 } | |
1903 } | |
1904 if ((unsigned) expon > DECIMAL_MAX_EXPON_128) { | |
1905 // OF | |
1906 #ifdef SET_STATUS_FLAGS | |
1907 __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
1908 #endif | |
1909 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
1910 #ifndef IEEE_ROUND_NEAREST | |
1911 if (*prounding_mode == ROUNDING_TO_ZERO | |
1912 || (sgn && *prounding_mode == ROUNDING_UP) || (!sgn | |
1913 && | |
1914 *prounding_mode | |
1915 == | |
1916 ROUNDING_DOWN)) | |
1917 { | |
1918 pres->w[1] = sgn | LARGEST_BID128_HIGH; | |
1919 pres->w[0] = LARGEST_BID128_LOW; | |
1920 } else | |
1921 #endif | |
1922 #endif | |
1923 { | |
1924 pres->w[1] = sgn | INFINITY_MASK64; | |
1925 pres->w[0] = 0; | |
1926 } | |
1927 return pres; | |
1928 } | |
1929 } | |
1930 | |
1931 pres->w[0] = coeff.w[0]; | |
1932 tmp = expon; | |
1933 tmp <<= 49; | |
1934 pres->w[1] = sgn | tmp | coeff.w[1]; | |
1935 | |
1936 return pres; | |
1937 } | |
1938 | |
1939 | |
1940 // | |
1941 // No overflow/underflow checks | |
1942 // No checking for coefficient == 10^34 (rounding artifact) | |
1943 // | |
1944 __BID_INLINE__ UINT128 * | |
1945 get_BID128_very_fast (UINT128 * pres, UINT64 sgn, int expon, | |
1946 UINT128 coeff) { | |
1947 UINT64 tmp; | |
1948 | |
1949 pres->w[0] = coeff.w[0]; | |
1950 tmp = expon; | |
1951 tmp <<= 49; | |
1952 pres->w[1] = sgn | tmp | coeff.w[1]; | |
1953 | |
1954 return pres; | |
1955 } | |
1956 | |
1957 // | |
1958 // No overflow/underflow checks | |
1959 // | |
1960 __BID_INLINE__ UINT128 * | |
1961 get_BID128_fast (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff) { | |
1962 UINT64 tmp; | |
1963 | |
1964 // coeff==10^34? | |
1965 if (coeff.w[1] == 0x0001ed09bead87c0ull | |
1966 && coeff.w[0] == 0x378d8e6400000000ull) { | |
1967 expon++; | |
1968 // set coefficient to 10^33 | |
1969 coeff.w[1] = 0x0000314dc6448d93ull; | |
1970 coeff.w[0] = 0x38c15b0a00000000ull; | |
1971 } | |
1972 | |
1973 pres->w[0] = coeff.w[0]; | |
1974 tmp = expon; | |
1975 tmp <<= 49; | |
1976 pres->w[1] = sgn | tmp | coeff.w[1]; | |
1977 | |
1978 return pres; | |
1979 } | |
1980 | |
1981 // | |
1982 // General BID128 pack macro | |
1983 // | |
1984 __BID_INLINE__ UINT128 * | |
1985 get_BID128 (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff, | |
1986 unsigned *prounding_mode, unsigned *fpsc) { | |
1987 UINT128 T; | |
1988 UINT64 tmp, tmp2; | |
1989 | |
1990 // coeff==10^34? | |
1991 if (coeff.w[1] == 0x0001ed09bead87c0ull | |
1992 && coeff.w[0] == 0x378d8e6400000000ull) { | |
1993 expon++; | |
1994 // set coefficient to 10^33 | |
1995 coeff.w[1] = 0x0000314dc6448d93ull; | |
1996 coeff.w[0] = 0x38c15b0a00000000ull; | |
1997 } | |
1998 // check OF, UF | |
1999 if (expon < 0 || expon > DECIMAL_MAX_EXPON_128) { | |
2000 // check UF | |
2001 if (expon < 0) { | |
2002 return handle_UF_128 (pres, sgn, expon, coeff, prounding_mode, | |
2003 fpsc); | |
2004 } | |
2005 | |
2006 if (expon - MAX_FORMAT_DIGITS_128 <= DECIMAL_MAX_EXPON_128) { | |
2007 T = power10_table_128[MAX_FORMAT_DIGITS_128 - 1]; | |
2008 while (__unsigned_compare_gt_128 (T, coeff) | |
2009 && expon > DECIMAL_MAX_EXPON_128) { | |
2010 coeff.w[1] = | |
2011 (coeff.w[1] << 3) + (coeff.w[1] << 1) + (coeff.w[0] >> 61) + | |
2012 (coeff.w[0] >> 63); | |
2013 tmp2 = coeff.w[0] << 3; | |
2014 coeff.w[0] = (coeff.w[0] << 1) + tmp2; | |
2015 if (coeff.w[0] < tmp2) | |
2016 coeff.w[1]++; | |
2017 | |
2018 expon--; | |
2019 } | |
2020 } | |
2021 if (expon > DECIMAL_MAX_EXPON_128) { | |
2022 if (!(coeff.w[1] | coeff.w[0])) { | |
2023 pres->w[1] = sgn | (((UINT64) DECIMAL_MAX_EXPON_128) << 49); | |
2024 pres->w[0] = 0; | |
2025 return pres; | |
2026 } | |
2027 // OF | |
2028 #ifdef SET_STATUS_FLAGS | |
2029 __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
2030 #endif | |
2031 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
2032 #ifndef IEEE_ROUND_NEAREST | |
2033 if (*prounding_mode == ROUNDING_TO_ZERO | |
2034 || (sgn && *prounding_mode == ROUNDING_UP) || (!sgn | |
2035 && | |
2036 *prounding_mode | |
2037 == | |
2038 ROUNDING_DOWN)) | |
2039 { | |
2040 pres->w[1] = sgn | LARGEST_BID128_HIGH; | |
2041 pres->w[0] = LARGEST_BID128_LOW; | |
2042 } else | |
2043 #endif | |
2044 #endif | |
2045 { | |
2046 pres->w[1] = sgn | INFINITY_MASK64; | |
2047 pres->w[0] = 0; | |
2048 } | |
2049 return pres; | |
2050 } | |
2051 } | |
2052 | |
2053 pres->w[0] = coeff.w[0]; | |
2054 tmp = expon; | |
2055 tmp <<= 49; | |
2056 pres->w[1] = sgn | tmp | coeff.w[1]; | |
2057 | |
2058 return pres; | |
2059 } | |
2060 | |
2061 | |
2062 // | |
2063 // Macro used for conversions from string | |
2064 // (no additional arguments given for rounding mode, status flags) | |
2065 // | |
2066 __BID_INLINE__ UINT128 * | |
2067 get_BID128_string (UINT128 * pres, UINT64 sgn, int expon, UINT128 coeff) { | |
2068 UINT128 D2, D8; | |
2069 UINT64 tmp; | |
2070 unsigned rmode = 0, status; | |
2071 | |
2072 // coeff==10^34? | |
2073 if (coeff.w[1] == 0x0001ed09bead87c0ull | |
2074 && coeff.w[0] == 0x378d8e6400000000ull) { | |
2075 expon++; | |
2076 // set coefficient to 10^33 | |
2077 coeff.w[1] = 0x0000314dc6448d93ull; | |
2078 coeff.w[0] = 0x38c15b0a00000000ull; | |
2079 } | |
2080 // check OF, UF | |
2081 if ((unsigned) expon > DECIMAL_MAX_EXPON_128) { | |
2082 // check UF | |
2083 if (expon < 0) | |
2084 return handle_UF_128 (pres, sgn, expon, coeff, &rmode, &status); | |
2085 | |
2086 // OF | |
2087 | |
2088 if (expon < DECIMAL_MAX_EXPON_128 + 34) { | |
2089 while (expon > DECIMAL_MAX_EXPON_128 && | |
2090 (coeff.w[1] < power10_table_128[33].w[1] || | |
2091 (coeff.w[1] == power10_table_128[33].w[1] | |
2092 && coeff.w[0] < power10_table_128[33].w[0]))) { | |
2093 D2.w[1] = (coeff.w[1] << 1) | (coeff.w[0] >> 63); | |
2094 D2.w[0] = coeff.w[0] << 1; | |
2095 D8.w[1] = (coeff.w[1] << 3) | (coeff.w[0] >> 61); | |
2096 D8.w[0] = coeff.w[0] << 3; | |
2097 | |
2098 __add_128_128 (coeff, D2, D8); | |
2099 expon--; | |
2100 } | |
2101 } else if (!(coeff.w[0] | coeff.w[1])) | |
2102 expon = DECIMAL_MAX_EXPON_128; | |
2103 | |
2104 if (expon > DECIMAL_MAX_EXPON_128) { | |
2105 pres->w[1] = sgn | INFINITY_MASK64; | |
2106 pres->w[0] = 0; | |
2107 switch (rmode) { | |
2108 case ROUNDING_DOWN: | |
2109 if (!sgn) { | |
2110 pres->w[1] = LARGEST_BID128_HIGH; | |
2111 pres->w[0] = LARGEST_BID128_LOW; | |
2112 } | |
2113 break; | |
2114 case ROUNDING_TO_ZERO: | |
2115 pres->w[1] = sgn | LARGEST_BID128_HIGH; | |
2116 pres->w[0] = LARGEST_BID128_LOW; | |
2117 break; | |
2118 case ROUNDING_UP: | |
2119 // round up | |
2120 if (sgn) { | |
2121 pres->w[1] = sgn | LARGEST_BID128_HIGH; | |
2122 pres->w[0] = LARGEST_BID128_LOW; | |
2123 } | |
2124 break; | |
2125 } | |
2126 | |
2127 return pres; | |
2128 } | |
2129 } | |
2130 | |
2131 pres->w[0] = coeff.w[0]; | |
2132 tmp = expon; | |
2133 tmp <<= 49; | |
2134 pres->w[1] = sgn | tmp | coeff.w[1]; | |
2135 | |
2136 return pres; | |
2137 } | |
2138 | |
2139 | |
2140 | |
2141 /***************************************************************************** | |
2142 * | |
2143 * BID32 pack/unpack macros | |
2144 * | |
2145 *****************************************************************************/ | |
2146 | |
2147 | |
2148 __BID_INLINE__ UINT32 | |
2149 unpack_BID32 (UINT32 * psign_x, int *pexponent_x, | |
2150 UINT32 * pcoefficient_x, UINT32 x) { | |
2151 UINT32 tmp; | |
2152 | |
2153 *psign_x = x & 0x80000000; | |
2154 | |
2155 if ((x & SPECIAL_ENCODING_MASK32) == SPECIAL_ENCODING_MASK32) { | |
2156 // special encodings | |
2157 if ((x & INFINITY_MASK32) == INFINITY_MASK32) { | |
2158 *pcoefficient_x = x & 0xfe0fffff; | |
2159 if ((x & 0x000fffff) >= 1000000) | |
2160 *pcoefficient_x = x & 0xfe000000; | |
2161 if ((x & NAN_MASK32) == INFINITY_MASK32) | |
2162 *pcoefficient_x = x & 0xf8000000; | |
2163 *pexponent_x = 0; | |
2164 return 0; // NaN or Infinity | |
2165 } | |
2166 // coefficient | |
2167 *pcoefficient_x = (x & SMALL_COEFF_MASK32) | LARGE_COEFF_HIGH_BIT32; | |
2168 // check for non-canonical value | |
2169 if (*pcoefficient_x >= 10000000) | |
2170 *pcoefficient_x = 0; | |
2171 // get exponent | |
2172 tmp = x >> 21; | |
2173 *pexponent_x = tmp & EXPONENT_MASK32; | |
2174 return 1; | |
2175 } | |
2176 // exponent | |
2177 tmp = x >> 23; | |
2178 *pexponent_x = tmp & EXPONENT_MASK32; | |
2179 // coefficient | |
2180 *pcoefficient_x = (x & LARGE_COEFF_MASK32); | |
2181 | |
2182 return *pcoefficient_x; | |
2183 } | |
2184 | |
2185 // | |
2186 // General pack macro for BID32 | |
2187 // | |
2188 __BID_INLINE__ UINT32 | |
2189 get_BID32 (UINT32 sgn, int expon, UINT64 coeff, int rmode, | |
2190 unsigned *fpsc) { | |
2191 UINT128 Q; | |
2192 UINT64 C64, remainder_h, carry, Stemp; | |
2193 UINT32 r, mask; | |
2194 int extra_digits, amount, amount2; | |
2195 unsigned status; | |
2196 | |
2197 if (coeff > 9999999ull) { | |
2198 expon++; | |
2199 coeff = 1000000ull; | |
2200 } | |
2201 // check for possible underflow/overflow | |
2202 if (((unsigned) expon) > DECIMAL_MAX_EXPON_32) { | |
2203 if (expon < 0) { | |
2204 // underflow | |
2205 if (expon + MAX_FORMAT_DIGITS_32 < 0) { | |
2206 #ifdef SET_STATUS_FLAGS | |
2207 __set_status_flags (fpsc, | |
2208 UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
2209 #endif | |
2210 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
2211 #ifndef IEEE_ROUND_NEAREST | |
2212 if (rmode == ROUNDING_DOWN && sgn) | |
2213 return 0x80000001; | |
2214 if (rmode == ROUNDING_UP && !sgn) | |
2215 return 1; | |
2216 #endif | |
2217 #endif | |
2218 // result is 0 | |
2219 return sgn; | |
2220 } | |
2221 // get digits to be shifted out | |
2222 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY | |
2223 rmode = 0; | |
2224 #endif | |
2225 #ifdef IEEE_ROUND_NEAREST | |
2226 rmode = 0; | |
2227 #endif | |
2228 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
2229 #ifndef IEEE_ROUND_NEAREST | |
2230 if (sgn && (unsigned) (rmode - 1) < 2) | |
2231 rmode = 3 - rmode; | |
2232 #endif | |
2233 #endif | |
2234 | |
2235 extra_digits = -expon; | |
2236 coeff += round_const_table[rmode][extra_digits]; | |
2237 | |
2238 // get coeff*(2^M[extra_digits])/10^extra_digits | |
2239 __mul_64x64_to_128 (Q, coeff, reciprocals10_64[extra_digits]); | |
2240 | |
2241 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 | |
2242 amount = short_recip_scale[extra_digits]; | |
2243 | |
2244 C64 = Q.w[1] >> amount; | |
2245 | |
2246 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY | |
2247 #ifndef IEEE_ROUND_NEAREST | |
2248 if (rmode == 0) //ROUNDING_TO_NEAREST | |
2249 #endif | |
2250 if (C64 & 1) { | |
2251 // check whether fractional part of initial_P/10^extra_digits is exactly .5 | |
2252 | |
2253 // get remainder | |
2254 amount2 = 64 - amount; | |
2255 remainder_h = 0; | |
2256 remainder_h--; | |
2257 remainder_h >>= amount2; | |
2258 remainder_h = remainder_h & Q.w[1]; | |
2259 | |
2260 if (!remainder_h && (Q.w[0] < reciprocals10_64[extra_digits])) { | |
2261 C64--; | |
2262 } | |
2263 } | |
2264 #endif | |
2265 | |
2266 #ifdef SET_STATUS_FLAGS | |
2267 | |
2268 if (is_inexact (fpsc)) | |
2269 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); | |
2270 else { | |
2271 status = INEXACT_EXCEPTION; | |
2272 // get remainder | |
2273 remainder_h = Q.w[1] << (64 - amount); | |
2274 | |
2275 switch (rmode) { | |
2276 case ROUNDING_TO_NEAREST: | |
2277 case ROUNDING_TIES_AWAY: | |
2278 // test whether fractional part is 0 | |
2279 if (remainder_h == 0x8000000000000000ull | |
2280 && (Q.w[0] < reciprocals10_64[extra_digits])) | |
2281 status = EXACT_STATUS; | |
2282 break; | |
2283 case ROUNDING_DOWN: | |
2284 case ROUNDING_TO_ZERO: | |
2285 if (!remainder_h && (Q.w[0] < reciprocals10_64[extra_digits])) | |
2286 status = EXACT_STATUS; | |
2287 break; | |
2288 default: | |
2289 // round up | |
2290 __add_carry_out (Stemp, carry, Q.w[0], | |
2291 reciprocals10_64[extra_digits]); | |
2292 if ((remainder_h >> (64 - amount)) + carry >= | |
2293 (((UINT64) 1) << amount)) | |
2294 status = EXACT_STATUS; | |
2295 } | |
2296 | |
2297 if (status != EXACT_STATUS) | |
2298 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION | status); | |
2299 } | |
2300 | |
2301 #endif | |
2302 | |
2303 return sgn | (UINT32) C64; | |
2304 } | |
2305 | |
2306 while (coeff < 1000000 && expon > DECIMAL_MAX_EXPON_32) { | |
2307 coeff = (coeff << 3) + (coeff << 1); | |
2308 expon--; | |
2309 } | |
2310 if (((unsigned) expon) > DECIMAL_MAX_EXPON_32) { | |
2311 __set_status_flags (fpsc, OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); | |
2312 // overflow | |
2313 r = sgn | INFINITY_MASK32; | |
2314 switch (rmode) { | |
2315 case ROUNDING_DOWN: | |
2316 if (!sgn) | |
2317 r = LARGEST_BID32; | |
2318 break; | |
2319 case ROUNDING_TO_ZERO: | |
2320 r = sgn | LARGEST_BID32; | |
2321 break; | |
2322 case ROUNDING_UP: | |
2323 // round up | |
2324 if (sgn) | |
2325 r = sgn | LARGEST_BID32; | |
2326 } | |
2327 return r; | |
2328 } | |
2329 } | |
2330 | |
2331 mask = 1 << 23; | |
2332 | |
2333 // check whether coefficient fits in DECIMAL_COEFF_FIT bits | |
2334 if (coeff < mask) { | |
2335 r = expon; | |
2336 r <<= 23; | |
2337 r |= ((UINT32) coeff | sgn); | |
2338 return r; | |
2339 } | |
2340 // special format | |
2341 | |
2342 r = expon; | |
2343 r <<= 21; | |
2344 r |= (sgn | SPECIAL_ENCODING_MASK32); | |
2345 // add coeff, without leading bits | |
2346 mask = (1 << 21) - 1; | |
2347 r |= (((UINT32) coeff) & mask); | |
2348 | |
2349 return r; | |
2350 } | |
2351 | |
2352 | |
2353 | |
2354 // | |
2355 // no overflow/underflow checks | |
2356 // | |
2357 __BID_INLINE__ UINT32 | |
2358 very_fast_get_BID32 (UINT32 sgn, int expon, UINT32 coeff) { | |
2359 UINT32 r, mask; | |
2360 | |
2361 mask = 1 << 23; | |
2362 | |
2363 // check whether coefficient fits in 10*2+3 bits | |
2364 if (coeff < mask) { | |
2365 r = expon; | |
2366 r <<= 23; | |
2367 r |= (coeff | sgn); | |
2368 return r; | |
2369 } | |
2370 // special format | |
2371 r = expon; | |
2372 r <<= 21; | |
2373 r |= (sgn | SPECIAL_ENCODING_MASK32); | |
2374 // add coeff, without leading bits | |
2375 mask = (1 << 21) - 1; | |
2376 coeff &= mask; | |
2377 r |= coeff; | |
2378 | |
2379 return r; | |
2380 } | |
2381 | |
2382 | |
2383 | |
2384 /************************************************************* | |
2385 * | |
2386 *************************************************************/ | |
2387 typedef | |
2388 ALIGN (16) | |
2389 struct { | |
2390 UINT64 w[6]; | |
2391 } UINT384; | |
2392 typedef ALIGN (16) | |
2393 struct { | |
2394 UINT64 w[8]; | |
2395 } UINT512; | |
2396 | |
2397 // #define P 34 | |
2398 #define MASK_STEERING_BITS 0x6000000000000000ull | |
2399 #define MASK_BINARY_EXPONENT1 0x7fe0000000000000ull | |
2400 #define MASK_BINARY_SIG1 0x001fffffffffffffull | |
2401 #define MASK_BINARY_EXPONENT2 0x1ff8000000000000ull | |
2402 //used to take G[2:w+3] (sec 3.3) | |
2403 #define MASK_BINARY_SIG2 0x0007ffffffffffffull | |
2404 //used to mask out G4:T0 (sec 3.3) | |
2405 #define MASK_BINARY_OR2 0x0020000000000000ull | |
2406 //used to prefix 8+G4 to T (sec 3.3) | |
2407 #define UPPER_EXPON_LIMIT 51 | |
2408 #define MASK_EXP 0x7ffe000000000000ull | |
2409 #define MASK_SPECIAL 0x7800000000000000ull | |
2410 #define MASK_NAN 0x7c00000000000000ull | |
2411 #define MASK_SNAN 0x7e00000000000000ull | |
2412 #define MASK_ANY_INF 0x7c00000000000000ull | |
2413 #define MASK_INF 0x7800000000000000ull | |
2414 #define MASK_SIGN 0x8000000000000000ull | |
2415 #define MASK_COEFF 0x0001ffffffffffffull | |
2416 #define BIN_EXP_BIAS (0x1820ull << 49) | |
2417 | |
2418 #define EXP_MIN 0x0000000000000000ull | |
2419 // EXP_MIN = (-6176 + 6176) << 49 | |
2420 #define EXP_MAX 0x5ffe000000000000ull | |
2421 // EXP_MAX = (6111 + 6176) << 49 | |
2422 #define EXP_MAX_P1 0x6000000000000000ull | |
2423 // EXP_MAX + 1 = (6111 + 6176 + 1) << 49 | |
2424 #define EXP_P1 0x0002000000000000ull | |
2425 // EXP_ P1= 1 << 49 | |
2426 #define expmin -6176 | |
2427 // min unbiased exponent | |
2428 #define expmax 6111 | |
2429 // max unbiased exponent | |
2430 #define expmin16 -398 | |
2431 // min unbiased exponent | |
2432 #define expmax16 369 | |
2433 // max unbiased exponent | |
2434 | |
2435 #define SIGNMASK32 0x80000000 | |
2436 #define BID64_SIG_MAX 0x002386F26FC0ffffull | |
2437 #define SIGNMASK64 0x8000000000000000ull | |
2438 | |
2439 // typedef unsigned int FPSC; // floating-point status and control | |
2440 // bit31: | |
2441 // bit30: | |
2442 // bit29: | |
2443 // bit28: | |
2444 // bit27: | |
2445 // bit26: | |
2446 // bit25: | |
2447 // bit24: | |
2448 // bit23: | |
2449 // bit22: | |
2450 // bit21: | |
2451 // bit20: | |
2452 // bit19: | |
2453 // bit18: | |
2454 // bit17: | |
2455 // bit16: | |
2456 // bit15: | |
2457 // bit14: RC:2 | |
2458 // bit13: RC:1 | |
2459 // bit12: RC:0 | |
2460 // bit11: PM | |
2461 // bit10: UM | |
2462 // bit9: OM | |
2463 // bit8: ZM | |
2464 // bit7: DM | |
2465 // bit6: IM | |
2466 // bit5: PE | |
2467 // bit4: UE | |
2468 // bit3: OE | |
2469 // bit2: ZE | |
2470 // bit1: DE | |
2471 // bit0: IE | |
2472 | |
2473 #define ROUNDING_MODE_MASK 0x00007000 | |
2474 | |
2475 typedef struct _DEC_DIGITS { | |
2476 unsigned int digits; | |
2477 UINT64 threshold_hi; | |
2478 UINT64 threshold_lo; | |
2479 unsigned int digits1; | |
2480 } DEC_DIGITS; | |
2481 | |
2482 extern DEC_DIGITS nr_digits[]; | |
2483 extern UINT64 midpoint64[]; | |
2484 extern UINT128 midpoint128[]; | |
2485 extern UINT192 midpoint192[]; | |
2486 extern UINT256 midpoint256[]; | |
2487 extern UINT64 ten2k64[]; | |
2488 extern UINT128 ten2k128[]; | |
2489 extern UINT256 ten2k256[]; | |
2490 extern UINT128 ten2mk128[]; | |
2491 extern UINT64 ten2mk64[]; | |
2492 extern UINT128 ten2mk128trunc[]; | |
2493 extern int shiftright128[]; | |
2494 extern UINT64 maskhigh128[]; | |
2495 extern UINT64 maskhigh128M[]; | |
2496 extern UINT64 maskhigh192M[]; | |
2497 extern UINT64 maskhigh256M[]; | |
2498 extern UINT64 onehalf128[]; | |
2499 extern UINT64 onehalf128M[]; | |
2500 extern UINT64 onehalf192M[]; | |
2501 extern UINT64 onehalf256M[]; | |
2502 extern UINT128 ten2mk128M[]; | |
2503 extern UINT128 ten2mk128truncM[]; | |
2504 extern UINT192 ten2mk192truncM[]; | |
2505 extern UINT256 ten2mk256truncM[]; | |
2506 extern int shiftright128M[]; | |
2507 extern int shiftright192M[]; | |
2508 extern int shiftright256M[]; | |
2509 extern UINT192 ten2mk192M[]; | |
2510 extern UINT256 ten2mk256M[]; | |
2511 extern unsigned char char_table2[]; | |
2512 extern unsigned char char_table3[]; | |
2513 | |
2514 extern UINT64 ten2m3k64[]; | |
2515 extern unsigned int shift_ten2m3k64[]; | |
2516 extern UINT128 ten2m3k128[]; | |
2517 extern unsigned int shift_ten2m3k128[]; | |
2518 | |
2519 | |
2520 | |
2521 /*************************************************************************** | |
2522 *************** TABLES FOR GENERAL ROUNDING FUNCTIONS ********************* | |
2523 ***************************************************************************/ | |
2524 | |
2525 extern UINT64 Kx64[]; | |
2526 extern unsigned int Ex64m64[]; | |
2527 extern UINT64 half64[]; | |
2528 extern UINT64 mask64[]; | |
2529 extern UINT64 ten2mxtrunc64[]; | |
2530 | |
2531 extern UINT128 Kx128[]; | |
2532 extern unsigned int Ex128m128[]; | |
2533 extern UINT64 half128[]; | |
2534 extern UINT64 mask128[]; | |
2535 extern UINT128 ten2mxtrunc128[]; | |
2536 | |
2537 extern UINT192 Kx192[]; | |
2538 extern unsigned int Ex192m192[]; | |
2539 extern UINT64 half192[]; | |
2540 extern UINT64 mask192[]; | |
2541 extern UINT192 ten2mxtrunc192[]; | |
2542 | |
2543 extern UINT256 Kx256[]; | |
2544 extern unsigned int Ex256m256[]; | |
2545 extern UINT64 half256[]; | |
2546 extern UINT64 mask256[]; | |
2547 extern UINT256 ten2mxtrunc256[]; | |
2548 | |
2549 typedef union __bid64_128 { | |
2550 UINT64 b64; | |
2551 UINT128 b128; | |
2552 } BID64_128; | |
2553 | |
2554 BID64_128 bid_fma (unsigned int P0, | |
2555 BID64_128 x1, unsigned int P1, | |
2556 BID64_128 y1, unsigned int P2, | |
2557 BID64_128 z1, unsigned int P3, | |
2558 unsigned int rnd_mode, FPSC * fpsc); | |
2559 | |
2560 #define P16 16 | |
2561 #define P34 34 | |
2562 | |
2563 union __int_double { | |
2564 UINT64 i; | |
2565 double d; | |
2566 }; | |
2567 typedef union __int_double int_double; | |
2568 | |
2569 | |
2570 union __int_float { | |
2571 UINT32 i; | |
2572 float d; | |
2573 }; | |
2574 typedef union __int_float int_float; | |
2575 | |
2576 #define SWAP(A,B,T) {\ | |
2577 T = A; \ | |
2578 A = B; \ | |
2579 B = T; \ | |
2580 } | |
2581 | |
2582 // this macro will find coefficient_x to be in [2^A, 2^(A+1) ) | |
2583 // ie it knows that it is A bits long | |
2584 #define NUMBITS(A, coefficient_x, tempx){\ | |
2585 temp_x.d=(float)coefficient_x;\ | |
2586 A=((tempx.i >>23) & EXPONENT_MASK32) - 0x7f;\ | |
2587 } | |
2588 | |
2589 enum class_types { | |
2590 signalingNaN, | |
2591 quietNaN, | |
2592 negativeInfinity, | |
2593 negativeNormal, | |
2594 negativeSubnormal, | |
2595 negativeZero, | |
2596 positiveZero, | |
2597 positiveSubnormal, | |
2598 positiveNormal, | |
2599 positiveInfinity | |
2600 }; | |
2601 | |
2602 typedef union { | |
2603 UINT64 ui64; | |
2604 double d; | |
2605 } BID_UI64DOUBLE; | |
2606 | |
2607 #endif |