]>
Commit | Line | Data |
---|---|---|
748086b7 | 1 | /* Copyright (C) 2007, 2009 Free Software Foundation, Inc. |
200359e8 L |
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 | |
748086b7 | 7 | Software Foundation; either version 3, or (at your option) any later |
200359e8 L |
8 | version. |
9 | ||
200359e8 L |
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 | ||
748086b7 JJ |
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/>. */ | |
200359e8 | 23 | |
b2a00c89 L |
24 | #include "bid_internal.h" |
25 | ||
26 | ||
27 | #if DECIMAL_CALL_BY_REFERENCE | |
28 | void | |
29 | bid64dq_add (UINT64 * pres, UINT64 * px, UINT128 * py | |
30 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
31 | _EXC_INFO_PARAM) { | |
32 | UINT64 x = *px; | |
33 | #if !DECIMAL_GLOBAL_ROUNDING | |
34 | unsigned int rnd_mode = *prnd_mode; | |
35 | #endif | |
36 | #else | |
37 | UINT64 | |
38 | bid64dq_add (UINT64 x, UINT128 y | |
39 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
40 | _EXC_INFO_PARAM) { | |
41 | #endif | |
42 | UINT64 res = 0xbaddbaddbaddbaddull; | |
43 | UINT128 x1; | |
44 | ||
45 | #if DECIMAL_CALL_BY_REFERENCE | |
46 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
47 | bid64qq_add (&res, &x1, py | |
48 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
49 | _EXC_INFO_ARG); | |
50 | #else | |
51 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
52 | res = bid64qq_add (x1, y | |
53 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
54 | _EXC_INFO_ARG); | |
55 | #endif | |
56 | BID_RETURN (res); | |
57 | } | |
58 | ||
59 | ||
60 | #if DECIMAL_CALL_BY_REFERENCE | |
61 | void | |
62 | bid64qd_add (UINT64 * pres, UINT128 * px, UINT64 * py | |
63 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
64 | _EXC_INFO_PARAM) { | |
65 | UINT64 y = *py; | |
66 | #if !DECIMAL_GLOBAL_ROUNDING | |
67 | unsigned int rnd_mode = *prnd_mode; | |
68 | #endif | |
69 | #else | |
70 | UINT64 | |
71 | bid64qd_add (UINT128 x, UINT64 y | |
72 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
73 | _EXC_INFO_PARAM) { | |
74 | #endif | |
75 | UINT64 res = 0xbaddbaddbaddbaddull; | |
76 | UINT128 y1; | |
77 | ||
78 | #if DECIMAL_CALL_BY_REFERENCE | |
79 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
80 | bid64qq_add (&res, px, &y1 | |
81 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
82 | _EXC_INFO_ARG); | |
83 | #else | |
84 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
85 | res = bid64qq_add (x, y1 | |
86 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
87 | _EXC_INFO_ARG); | |
88 | #endif | |
89 | BID_RETURN (res); | |
90 | } | |
91 | ||
92 | ||
93 | #if DECIMAL_CALL_BY_REFERENCE | |
94 | void | |
95 | bid64qq_add (UINT64 * pres, UINT128 * px, UINT128 * py | |
96 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
97 | _EXC_INFO_PARAM) { | |
98 | UINT128 x = *px, y = *py; | |
99 | #if !DECIMAL_GLOBAL_ROUNDING | |
100 | unsigned int rnd_mode = *prnd_mode; | |
101 | #endif | |
102 | #else | |
103 | UINT64 | |
104 | bid64qq_add (UINT128 x, UINT128 y | |
105 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
106 | _EXC_INFO_PARAM) { | |
107 | #endif | |
108 | ||
109 | UINT128 one = { {0x0000000000000001ull, 0x3040000000000000ull} | |
110 | }; | |
111 | UINT64 res = 0xbaddbaddbaddbaddull; | |
112 | ||
113 | BID_SWAP128 (one); | |
114 | #if DECIMAL_CALL_BY_REFERENCE | |
115 | bid64qqq_fma (&res, &one, &x, &y | |
116 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
117 | _EXC_INFO_ARG); | |
118 | #else | |
119 | res = bid64qqq_fma (one, x, y | |
120 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
121 | _EXC_INFO_ARG); | |
122 | #endif | |
123 | BID_RETURN (res); | |
124 | } | |
125 | ||
126 | ||
127 | #if DECIMAL_CALL_BY_REFERENCE | |
128 | void | |
129 | bid128dd_add (UINT128 * pres, UINT64 * px, UINT64 * py | |
130 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
131 | _EXC_INFO_PARAM) { | |
132 | UINT64 x = *px, y = *py; | |
133 | #if !DECIMAL_GLOBAL_ROUNDING | |
134 | unsigned int rnd_mode = *prnd_mode; | |
135 | #endif | |
136 | #else | |
137 | UINT128 | |
138 | bid128dd_add (UINT64 x, UINT64 y | |
139 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
140 | _EXC_INFO_PARAM) { | |
141 | #endif | |
142 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
143 | }; | |
144 | UINT128 x1, y1; | |
145 | ||
146 | #if DECIMAL_CALL_BY_REFERENCE | |
147 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
148 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
149 | bid128_add (&res, &x1, &y1 | |
150 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
151 | _EXC_INFO_ARG); | |
152 | #else | |
153 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
154 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
155 | res = bid128_add (x1, y1 | |
156 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
157 | _EXC_INFO_ARG); | |
158 | #endif | |
159 | BID_RETURN (res); | |
160 | } | |
161 | ||
162 | ||
163 | #if DECIMAL_CALL_BY_REFERENCE | |
164 | void | |
165 | bid128dq_add (UINT128 * pres, UINT64 * px, UINT128 * py | |
166 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
167 | _EXC_INFO_PARAM) { | |
168 | UINT64 x = *px; | |
169 | #if !DECIMAL_GLOBAL_ROUNDING | |
170 | unsigned int rnd_mode = *prnd_mode; | |
171 | #endif | |
172 | #else | |
173 | UINT128 | |
174 | bid128dq_add (UINT64 x, UINT128 y | |
175 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
176 | _EXC_INFO_PARAM) { | |
177 | #endif | |
178 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
179 | }; | |
180 | UINT128 x1; | |
181 | ||
182 | #if DECIMAL_CALL_BY_REFERENCE | |
183 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
184 | bid128_add (&res, &x1, py | |
185 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
186 | _EXC_INFO_ARG); | |
187 | #else | |
188 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
189 | res = bid128_add (x1, y | |
190 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
191 | _EXC_INFO_ARG); | |
192 | #endif | |
193 | BID_RETURN (res); | |
194 | } | |
195 | ||
196 | ||
197 | #if DECIMAL_CALL_BY_REFERENCE | |
198 | void | |
199 | bid128qd_add (UINT128 * pres, UINT128 * px, UINT64 * py | |
200 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
201 | _EXC_INFO_PARAM) { | |
202 | UINT64 y = *py; | |
203 | #if !DECIMAL_GLOBAL_ROUNDING | |
204 | unsigned int rnd_mode = *prnd_mode; | |
205 | #endif | |
206 | #else | |
207 | UINT128 | |
208 | bid128qd_add (UINT128 x, UINT64 y | |
209 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
210 | _EXC_INFO_PARAM) { | |
211 | #endif | |
212 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
213 | }; | |
214 | UINT128 y1; | |
215 | ||
216 | #if DECIMAL_CALL_BY_REFERENCE | |
217 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
218 | bid128_add (&res, px, &y1 | |
219 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
220 | _EXC_INFO_ARG); | |
221 | #else | |
222 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
223 | res = bid128_add (x, y1 | |
224 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
225 | _EXC_INFO_ARG); | |
226 | #endif | |
227 | BID_RETURN (res); | |
228 | } | |
229 | ||
230 | ||
231 | // bid128_add stands for bid128qq_add | |
232 | ||
233 | ||
200359e8 | 234 | /***************************************************************************** |
b2a00c89 | 235 | * BID64/BID128 sub |
200359e8 L |
236 | ****************************************************************************/ |
237 | ||
b2a00c89 L |
238 | #if DECIMAL_CALL_BY_REFERENCE |
239 | void | |
240 | bid64dq_sub (UINT64 * pres, UINT64 * px, UINT128 * py | |
241 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
242 | _EXC_INFO_PARAM) { | |
243 | UINT64 x = *px; | |
244 | #if !DECIMAL_GLOBAL_ROUNDING | |
245 | unsigned int rnd_mode = *prnd_mode; | |
246 | #endif | |
247 | #else | |
248 | UINT64 | |
249 | bid64dq_sub (UINT64 x, UINT128 y | |
250 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
251 | _EXC_INFO_PARAM) { | |
252 | #endif | |
253 | UINT64 res = 0xbaddbaddbaddbaddull; | |
254 | UINT128 x1; | |
255 | ||
256 | #if DECIMAL_CALL_BY_REFERENCE | |
257 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
258 | bid64qq_sub (&res, &x1, py | |
259 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
260 | _EXC_INFO_ARG); | |
261 | #else | |
262 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
263 | res = bid64qq_sub (x1, y | |
264 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
265 | _EXC_INFO_ARG); | |
266 | #endif | |
267 | BID_RETURN (res); | |
268 | } | |
269 | ||
270 | ||
271 | #if DECIMAL_CALL_BY_REFERENCE | |
272 | void | |
273 | bid64qd_sub (UINT64 * pres, UINT128 * px, UINT64 * py | |
274 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
275 | _EXC_INFO_PARAM) { | |
276 | UINT64 y = *py; | |
277 | #if !DECIMAL_GLOBAL_ROUNDING | |
278 | unsigned int rnd_mode = *prnd_mode; | |
279 | #endif | |
280 | #else | |
281 | UINT64 | |
282 | bid64qd_sub (UINT128 x, UINT64 y | |
283 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
284 | _EXC_INFO_PARAM) { | |
285 | #endif | |
286 | UINT64 res = 0xbaddbaddbaddbaddull; | |
287 | UINT128 y1; | |
288 | ||
289 | #if DECIMAL_CALL_BY_REFERENCE | |
290 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
291 | bid64qq_sub (&res, px, &y1 | |
292 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
293 | _EXC_INFO_ARG); | |
294 | #else | |
295 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
296 | res = bid64qq_sub (x, y1 | |
297 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
298 | _EXC_INFO_ARG); | |
299 | #endif | |
300 | BID_RETURN (res); | |
301 | } | |
302 | ||
200359e8 L |
303 | |
304 | #if DECIMAL_CALL_BY_REFERENCE | |
305 | void | |
b2a00c89 L |
306 | bid64qq_sub (UINT64 * pres, UINT128 * px, UINT128 * py |
307 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
308 | _EXC_INFO_PARAM) { | |
200359e8 L |
309 | UINT128 x = *px, y = *py; |
310 | #if !DECIMAL_GLOBAL_ROUNDING | |
311 | unsigned int rnd_mode = *prnd_mode; | |
312 | #endif | |
313 | #else | |
b2a00c89 L |
314 | UINT64 |
315 | bid64qq_sub (UINT128 x, UINT128 y | |
316 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
317 | _EXC_INFO_PARAM) { | |
318 | #endif | |
319 | ||
320 | UINT128 one = { {0x0000000000000001ull, 0x3040000000000000ull} | |
321 | }; | |
322 | UINT64 res = 0xbaddbaddbaddbaddull; | |
323 | UINT64 y_sign; | |
324 | ||
325 | BID_SWAP128 (one); | |
326 | if ((y.w[HIGH_128W] & MASK_NAN) != MASK_NAN) { // y is not NAN | |
327 | // change its sign | |
328 | y_sign = y.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative | |
329 | if (y_sign) | |
330 | y.w[HIGH_128W] = y.w[HIGH_128W] & 0x7fffffffffffffffull; | |
331 | else | |
332 | y.w[HIGH_128W] = y.w[HIGH_128W] | 0x8000000000000000ull; | |
333 | } | |
334 | #if DECIMAL_CALL_BY_REFERENCE | |
335 | bid64qqq_fma (&res, &one, &x, &y | |
336 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
337 | _EXC_INFO_ARG); | |
338 | #else | |
339 | res = bid64qqq_fma (one, x, y | |
340 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
341 | _EXC_INFO_ARG); | |
342 | #endif | |
343 | BID_RETURN (res); | |
344 | } | |
345 | ||
346 | ||
347 | #if DECIMAL_CALL_BY_REFERENCE | |
348 | void | |
349 | bid128dd_sub (UINT128 * pres, UINT64 * px, UINT64 * py | |
350 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
351 | _EXC_INFO_PARAM) { | |
352 | UINT64 x = *px, y = *py; | |
353 | #if !DECIMAL_GLOBAL_ROUNDING | |
354 | unsigned int rnd_mode = *prnd_mode; | |
355 | #endif | |
356 | #else | |
200359e8 | 357 | UINT128 |
b2a00c89 L |
358 | bid128dd_sub (UINT64 x, UINT64 y |
359 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
360 | _EXC_INFO_PARAM) { | |
200359e8 | 361 | #endif |
b2a00c89 L |
362 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} |
363 | }; | |
364 | UINT128 x1, y1; | |
200359e8 | 365 | |
b2a00c89 L |
366 | #if DECIMAL_CALL_BY_REFERENCE |
367 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
368 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
369 | bid128_sub (&res, &x1, &y1 | |
370 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
371 | _EXC_INFO_ARG); | |
372 | #else | |
373 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
374 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
375 | res = bid128_sub (x1, y1 | |
376 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
377 | _EXC_INFO_ARG); | |
378 | #endif | |
379 | BID_RETURN (res); | |
380 | } | |
381 | ||
382 | ||
383 | #if DECIMAL_CALL_BY_REFERENCE | |
384 | void | |
385 | bid128dq_sub (UINT128 * pres, UINT64 * px, UINT128 * py | |
386 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
387 | _EXC_INFO_PARAM) { | |
388 | UINT64 x = *px; | |
389 | #if !DECIMAL_GLOBAL_ROUNDING | |
390 | unsigned int rnd_mode = *prnd_mode; | |
391 | #endif | |
392 | #else | |
393 | UINT128 | |
394 | bid128dq_sub (UINT64 x, UINT128 y | |
395 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
396 | _EXC_INFO_PARAM) { | |
397 | #endif | |
398 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
399 | }; | |
400 | UINT128 x1; | |
401 | ||
402 | #if DECIMAL_CALL_BY_REFERENCE | |
403 | bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
404 | bid128_sub (&res, &x1, py | |
405 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
406 | _EXC_INFO_ARG); | |
407 | #else | |
408 | x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
409 | res = bid128_sub (x1, y | |
410 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
411 | _EXC_INFO_ARG); | |
412 | #endif | |
413 | BID_RETURN (res); | |
414 | } | |
415 | ||
416 | ||
417 | #if DECIMAL_CALL_BY_REFERENCE | |
418 | void | |
419 | bid128qd_sub (UINT128 * pres, UINT128 * px, UINT64 * py | |
420 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
421 | _EXC_INFO_PARAM) { | |
422 | UINT64 y = *py; | |
423 | #if !DECIMAL_GLOBAL_ROUNDING | |
424 | unsigned int rnd_mode = *prnd_mode; | |
425 | #endif | |
426 | #else | |
427 | UINT128 | |
428 | bid128qd_sub (UINT128 x, UINT64 y | |
429 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
430 | _EXC_INFO_PARAM) { | |
431 | #endif | |
432 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
433 | }; | |
434 | UINT128 y1; | |
435 | ||
436 | #if DECIMAL_CALL_BY_REFERENCE | |
437 | bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
438 | bid128_sub (&res, px, &y1 | |
439 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
440 | _EXC_INFO_ARG); | |
441 | #else | |
442 | y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); | |
443 | res = bid128_sub (x, y1 | |
444 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
445 | _EXC_INFO_ARG); | |
446 | #endif | |
447 | BID_RETURN (res); | |
448 | } | |
449 | ||
450 | #if DECIMAL_CALL_BY_REFERENCE | |
451 | void | |
452 | bid128_add (UINT128 * pres, UINT128 * px, UINT128 * py | |
453 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
454 | _EXC_INFO_PARAM) { | |
455 | UINT128 x = *px, y = *py; | |
456 | #if !DECIMAL_GLOBAL_ROUNDING | |
457 | unsigned int rnd_mode = *prnd_mode; | |
458 | #endif | |
459 | #else | |
460 | UINT128 | |
461 | bid128_add (UINT128 x, UINT128 y | |
462 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
463 | _EXC_INFO_PARAM) { | |
464 | #endif | |
465 | UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} | |
466 | }; | |
200359e8 | 467 | UINT64 x_sign, y_sign, tmp_sign; |
b2a00c89 | 468 | UINT64 x_exp, y_exp, tmp_exp; // e1 = x_exp, e2 = y_exp |
200359e8 L |
469 | UINT64 C1_hi, C2_hi, tmp_signif_hi; |
470 | UINT64 C1_lo, C2_lo, tmp_signif_lo; | |
b2a00c89 L |
471 | // Note: C1.w[1], C1.w[0] represent C1_hi, C1_lo (all UINT64) |
472 | // Note: C2.w[1], C2.w[0] represent C2_hi, C2_lo (all UINT64) | |
200359e8 L |
473 | UINT64 tmp64, tmp64A, tmp64B; |
474 | BID_UI64DOUBLE tmp1, tmp2; | |
475 | int x_nr_bits, y_nr_bits; | |
476 | int q1, q2, delta, scale, x1, ind, shift, tmp_inexact = 0; | |
477 | UINT64 halfulp64; | |
478 | UINT128 halfulp128; | |
479 | UINT128 C1, C2; | |
b2a00c89 L |
480 | UINT128 ten2m1; |
481 | UINT128 highf2star; // top 128 bits in f2*; low 128 bits in R256[1], R256[0] | |
200359e8 L |
482 | UINT256 P256, Q256, R256; |
483 | int is_inexact = 0, is_midpoint_lt_even = 0, is_midpoint_gt_even = 0; | |
484 | int is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; | |
485 | int second_pass = 0; | |
486 | ||
b2a00c89 L |
487 | BID_SWAP128 (x); |
488 | BID_SWAP128 (y); | |
489 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative | |
490 | y_sign = y.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative | |
491 | ||
200359e8 L |
492 | // check for NaN or Infinity |
493 | if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) | |
494 | || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) { | |
495 | // x is special or y is special | |
b2a00c89 L |
496 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
497 | // check first for non-canonical NaN payload | |
498 | if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || | |
499 | (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) | |
500 | && (x.w[0] > 0x38c15b09ffffffffull))) { | |
501 | x.w[1] = x.w[1] & 0xffffc00000000000ull; | |
502 | x.w[0] = 0x0ull; | |
200359e8 | 503 | } |
b2a00c89 L |
504 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN |
505 | // set invalid flag | |
506 | *pfpsf |= INVALID_EXCEPTION; | |
507 | // return quiet (x) | |
508 | res.w[1] = x.w[1] & 0xfc003fffffffffffull; | |
509 | // clear out also G[6]-G[16] | |
510 | res.w[0] = x.w[0]; | |
511 | } else { // x is QNaN | |
512 | // return x | |
513 | res.w[1] = x.w[1] & 0xfc003fffffffffffull; | |
514 | // clear out G[6]-G[16] | |
515 | res.w[0] = x.w[0]; | |
516 | // if y = SNaN signal invalid exception | |
517 | if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { | |
518 | // set invalid flag | |
519 | *pfpsf |= INVALID_EXCEPTION; | |
520 | } | |
521 | } | |
522 | BID_SWAP128 (res); | |
200359e8 | 523 | BID_RETURN (res); |
b2a00c89 L |
524 | } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN |
525 | // check first for non-canonical NaN payload | |
526 | if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || | |
527 | (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) | |
528 | && (y.w[0] > 0x38c15b09ffffffffull))) { | |
529 | y.w[1] = y.w[1] & 0xffffc00000000000ull; | |
530 | y.w[0] = 0x0ull; | |
531 | } | |
532 | if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN | |
533 | // set invalid flag | |
534 | *pfpsf |= INVALID_EXCEPTION; | |
535 | // return quiet (y) | |
536 | res.w[1] = y.w[1] & 0xfc003fffffffffffull; | |
537 | // clear out also G[6]-G[16] | |
538 | res.w[0] = y.w[0]; | |
539 | } else { // y is QNaN | |
540 | // return y | |
541 | res.w[1] = y.w[1] & 0xfc003fffffffffffull; | |
542 | // clear out G[6]-G[16] | |
543 | res.w[0] = y.w[0]; | |
200359e8 | 544 | } |
b2a00c89 | 545 | BID_SWAP128 (res); |
200359e8 | 546 | BID_RETURN (res); |
b2a00c89 L |
547 | } else { // neither x not y is NaN; at least one is infinity |
548 | if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x is infinity | |
549 | if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y is infinity | |
550 | // if same sign, return either of them | |
551 | if ((x.w[1] & MASK_SIGN) == (y.w[1] & MASK_SIGN)) { | |
552 | res.w[1] = x_sign | MASK_INF; | |
553 | res.w[0] = 0x0ull; | |
554 | } else { // x and y are infinities of opposite signs | |
555 | // set invalid flag | |
556 | *pfpsf |= INVALID_EXCEPTION; | |
557 | // return QNaN Indefinite | |
558 | res.w[1] = 0x7c00000000000000ull; | |
559 | res.w[0] = 0x0000000000000000ull; | |
560 | } | |
561 | } else { // y is 0 or finite | |
562 | // return x | |
563 | res.w[1] = x_sign | MASK_INF; | |
564 | res.w[0] = 0x0ull; | |
565 | } | |
566 | } else { // x is not NaN or infinity, so y must be infinity | |
567 | res.w[1] = y_sign | MASK_INF; | |
568 | res.w[0] = 0x0ull; | |
200359e8 | 569 | } |
b2a00c89 | 570 | BID_SWAP128 (res); |
200359e8 L |
571 | BID_RETURN (res); |
572 | } | |
573 | } | |
574 | // unpack the arguments | |
b2a00c89 | 575 | |
200359e8 | 576 | // unpack x |
200359e8 L |
577 | C1_hi = x.w[1] & MASK_COEFF; |
578 | C1_lo = x.w[0]; | |
200359e8 L |
579 | // test for non-canonical values: |
580 | // - values whose encoding begins with x00, x01, or x10 and whose | |
581 | // coefficient is larger than 10^34 -1, or | |
582 | // - values whose encoding begins with x1100, x1101, x1110 (if NaNs | |
583 | // and infinitis were eliminated already this test is reduced to | |
584 | // checking for x10x) | |
585 | ||
b2a00c89 L |
586 | // x is not infinity; check for non-canonical values - treated as zero |
587 | if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { | |
588 | // G0_G1=11; non-canonical | |
589 | x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits | |
590 | C1_hi = 0; // significand high | |
591 | C1_lo = 0; // significand low | |
592 | } else { // G0_G1 != 11 | |
593 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits | |
594 | if (C1_hi > 0x0001ed09bead87c0ull || | |
595 | (C1_hi == 0x0001ed09bead87c0ull | |
596 | && C1_lo > 0x378d8e63ffffffffull)) { | |
597 | // x is non-canonical if coefficient is larger than 10^34 -1 | |
598 | C1_hi = 0; | |
599 | C1_lo = 0; | |
600 | } else { // canonical | |
601 | ; | |
200359e8 | 602 | } |
200359e8 | 603 | } |
b2a00c89 L |
604 | |
605 | // unpack y | |
606 | C2_hi = y.w[1] & MASK_COEFF; | |
607 | C2_lo = y.w[0]; | |
608 | // y is not infinity; check for non-canonical values - treated as zero | |
609 | if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { | |
610 | // G0_G1=11; non-canonical | |
611 | y_exp = (y.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits | |
612 | C2_hi = 0; // significand high | |
613 | C2_lo = 0; // significand low | |
614 | } else { // G0_G1 != 11 | |
615 | y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bits | |
616 | if (C2_hi > 0x0001ed09bead87c0ull || | |
617 | (C2_hi == 0x0001ed09bead87c0ull | |
618 | && C2_lo > 0x378d8e63ffffffffull)) { | |
619 | // y is non-canonical if coefficient is larger than 10^34 -1 | |
620 | C2_hi = 0; | |
621 | C2_lo = 0; | |
622 | } else { // canonical | |
623 | ; | |
200359e8 | 624 | } |
200359e8 L |
625 | } |
626 | ||
627 | if ((C1_hi == 0x0ull) && (C1_lo == 0x0ull)) { | |
628 | // x is 0 and y is not special | |
629 | // if y is 0 return 0 with the smaller exponent | |
630 | if ((C2_hi == 0x0ull) && (C2_lo == 0x0ull)) { | |
631 | if (x_exp < y_exp) | |
b2a00c89 | 632 | res.w[1] = x_exp; |
200359e8 | 633 | else |
b2a00c89 | 634 | res.w[1] = y_exp; |
200359e8 | 635 | if (x_sign && y_sign) |
b2a00c89 | 636 | res.w[1] = res.w[1] | x_sign; // both negative |
200359e8 | 637 | else if (rnd_mode == ROUNDING_DOWN && x_sign != y_sign) |
b2a00c89 | 638 | res.w[1] = res.w[1] | 0x8000000000000000ull; // -0 |
200359e8 L |
639 | // else; // res = +0 |
640 | res.w[0] = 0; | |
641 | } else { | |
642 | // for 0 + y return y, with the preferred exponent | |
643 | if (y_exp <= x_exp) { | |
b2a00c89 L |
644 | res.w[1] = y.w[1]; |
645 | res.w[0] = y.w[0]; | |
646 | } else { // if y_exp > x_exp | |
647 | // return (C2 * 10^scale) * 10^(y_exp - scale) | |
648 | // where scale = min (P34-q2, y_exp-x_exp) | |
649 | // determine q2 = nr. of decimal digits in y | |
650 | // determine first the nr. of bits in y (y_nr_bits) | |
651 | ||
652 | if (C2_hi == 0) { // y_bits is the nr. of bits in C2_lo | |
653 | if (C2_lo >= 0x0020000000000000ull) { // y >= 2^53 | |
654 | // split the 64-bit value in two 32-bit halves to avoid | |
655 | // rounding errors | |
656 | if (C2_lo >= 0x0000000100000000ull) { // y >= 2^32 | |
657 | tmp2.d = (double) (C2_lo >> 32); // exact conversion | |
658 | y_nr_bits = | |
659 | 32 + | |
660 | ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
661 | } else { // y < 2^32 | |
662 | tmp2.d = (double) (C2_lo); // exact conversion | |
663 | y_nr_bits = | |
664 | ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
665 | } | |
666 | } else { // if y < 2^53 | |
667 | tmp2.d = (double) C2_lo; // exact conversion | |
668 | y_nr_bits = | |
669 | ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
670 | } | |
671 | } else { // C2_hi != 0 => nr. bits = 64 + nr_bits (C2_hi) | |
672 | tmp2.d = (double) C2_hi; // exact conversion | |
673 | y_nr_bits = | |
674 | 64 + ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
675 | } | |
676 | q2 = nr_digits[y_nr_bits].digits; | |
677 | if (q2 == 0) { | |
678 | q2 = nr_digits[y_nr_bits].digits1; | |
679 | if (C2_hi > nr_digits[y_nr_bits].threshold_hi || | |
680 | (C2_hi == nr_digits[y_nr_bits].threshold_hi && | |
681 | C2_lo >= nr_digits[y_nr_bits].threshold_lo)) | |
682 | q2++; | |
683 | } | |
684 | // return (C2 * 10^scale) * 10^(y_exp - scale) | |
685 | // where scale = min (P34-q2, y_exp-x_exp) | |
686 | scale = P34 - q2; | |
687 | ind = (y_exp - x_exp) >> 49; | |
688 | if (ind < scale) | |
689 | scale = ind; | |
690 | if (scale == 0) { | |
691 | res.w[1] = y.w[1]; | |
692 | res.w[0] = y.w[0]; | |
693 | } else if (q2 <= 19) { // y fits in 64 bits | |
694 | if (scale <= 19) { // 10^scale fits in 64 bits | |
695 | // 64 x 64 C2_lo * ten2k64[scale] | |
696 | __mul_64x64_to_128MACH (res, C2_lo, ten2k64[scale]); | |
697 | } else { // 10^scale fits in 128 bits | |
698 | // 64 x 128 C2_lo * ten2k128[scale - 20] | |
699 | __mul_128x64_to_128 (res, C2_lo, ten2k128[scale - 20]); | |
700 | } | |
701 | } else { // y fits in 128 bits, but 10^scale must fit in 64 bits | |
702 | // 64 x 128 ten2k64[scale] * C2 | |
703 | C2.w[1] = C2_hi; | |
704 | C2.w[0] = C2_lo; | |
705 | __mul_128x64_to_128 (res, ten2k64[scale], C2); | |
706 | } | |
707 | // subtract scale from the exponent | |
708 | y_exp = y_exp - ((UINT64) scale << 49); | |
709 | res.w[1] = res.w[1] | y_sign | y_exp; | |
200359e8 L |
710 | } |
711 | } | |
b2a00c89 | 712 | BID_SWAP128 (res); |
200359e8 L |
713 | BID_RETURN (res); |
714 | } else if ((C2_hi == 0x0ull) && (C2_lo == 0x0ull)) { | |
715 | // y is 0 and x is not special, and not zero | |
716 | // for x + 0 return x, with the preferred exponent | |
717 | if (x_exp <= y_exp) { | |
718 | res.w[1] = x.w[1]; | |
719 | res.w[0] = x.w[0]; | |
b2a00c89 | 720 | } else { // if x_exp > y_exp |
200359e8 L |
721 | // return (C1 * 10^scale) * 10^(x_exp - scale) |
722 | // where scale = min (P34-q1, x_exp-y_exp) | |
723 | // determine q1 = nr. of decimal digits in x | |
724 | // determine first the nr. of bits in x | |
b2a00c89 L |
725 | if (C1_hi == 0) { // x_bits is the nr. of bits in C1_lo |
726 | if (C1_lo >= 0x0020000000000000ull) { // x >= 2^53 | |
727 | // split the 64-bit value in two 32-bit halves to avoid | |
728 | // rounding errors | |
729 | if (C1_lo >= 0x0000000100000000ull) { // x >= 2^32 | |
730 | tmp1.d = (double) (C1_lo >> 32); // exact conversion | |
731 | x_nr_bits = | |
732 | 32 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - | |
733 | 0x3ff); | |
734 | } else { // x < 2^32 | |
735 | tmp1.d = (double) (C1_lo); // exact conversion | |
736 | x_nr_bits = | |
737 | ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
738 | } | |
739 | } else { // if x < 2^53 | |
740 | tmp1.d = (double) C1_lo; // exact conversion | |
741 | x_nr_bits = | |
742 | ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
743 | } | |
744 | } else { // C1_hi != 0 => nr. bits = 64 + nr_bits (C1_hi) | |
745 | tmp1.d = (double) C1_hi; // exact conversion | |
746 | x_nr_bits = | |
747 | 64 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
200359e8 | 748 | } |
b2a00c89 | 749 | q1 = nr_digits[x_nr_bits].digits; |
200359e8 | 750 | if (q1 == 0) { |
b2a00c89 L |
751 | q1 = nr_digits[x_nr_bits].digits1; |
752 | if (C1_hi > nr_digits[x_nr_bits].threshold_hi || | |
753 | (C1_hi == nr_digits[x_nr_bits].threshold_hi && | |
754 | C1_lo >= nr_digits[x_nr_bits].threshold_lo)) | |
755 | q1++; | |
200359e8 L |
756 | } |
757 | // return (C1 * 10^scale) * 10^(x_exp - scale) | |
758 | // where scale = min (P34-q1, x_exp-y_exp) | |
759 | scale = P34 - q1; | |
760 | ind = (x_exp - y_exp) >> 49; | |
761 | if (ind < scale) | |
b2a00c89 | 762 | scale = ind; |
200359e8 | 763 | if (scale == 0) { |
b2a00c89 L |
764 | res.w[1] = x.w[1]; |
765 | res.w[0] = x.w[0]; | |
766 | } else if (q1 <= 19) { // x fits in 64 bits | |
767 | if (scale <= 19) { // 10^scale fits in 64 bits | |
768 | // 64 x 64 C1_lo * ten2k64[scale] | |
769 | __mul_64x64_to_128MACH (res, C1_lo, ten2k64[scale]); | |
770 | } else { // 10^scale fits in 128 bits | |
771 | // 64 x 128 C1_lo * ten2k128[scale - 20] | |
772 | __mul_128x64_to_128 (res, C1_lo, ten2k128[scale - 20]); | |
773 | } | |
774 | } else { // x fits in 128 bits, but 10^scale must fit in 64 bits | |
775 | // 64 x 128 ten2k64[scale] * C1 | |
776 | C1.w[1] = C1_hi; | |
777 | C1.w[0] = C1_lo; | |
778 | __mul_128x64_to_128 (res, ten2k64[scale], C1); | |
200359e8 L |
779 | } |
780 | // subtract scale from the exponent | |
781 | x_exp = x_exp - ((UINT64) scale << 49); | |
782 | res.w[1] = res.w[1] | x_sign | x_exp; | |
783 | } | |
b2a00c89 | 784 | BID_SWAP128 (res); |
200359e8 | 785 | BID_RETURN (res); |
b2a00c89 | 786 | } else { // x and y are not canonical, not special, and are not zero |
200359e8 L |
787 | // note that the result may still be zero, and then it has to have the |
788 | // preferred exponent | |
b2a00c89 | 789 | if (x_exp < y_exp) { // if exp_x < exp_y then swap x and y |
200359e8 L |
790 | tmp_sign = x_sign; |
791 | tmp_exp = x_exp; | |
792 | tmp_signif_hi = C1_hi; | |
793 | tmp_signif_lo = C1_lo; | |
794 | x_sign = y_sign; | |
795 | x_exp = y_exp; | |
796 | C1_hi = C2_hi; | |
797 | C1_lo = C2_lo; | |
798 | y_sign = tmp_sign; | |
799 | y_exp = tmp_exp; | |
800 | C2_hi = tmp_signif_hi; | |
801 | C2_lo = tmp_signif_lo; | |
802 | } | |
803 | // q1 = nr. of decimal digits in x | |
804 | // determine first the nr. of bits in x | |
b2a00c89 L |
805 | if (C1_hi == 0) { // x_bits is the nr. of bits in C1_lo |
806 | if (C1_lo >= 0x0020000000000000ull) { // x >= 2^53 | |
807 | //split the 64-bit value in two 32-bit halves to avoid rounding errors | |
808 | if (C1_lo >= 0x0000000100000000ull) { // x >= 2^32 | |
809 | tmp1.d = (double) (C1_lo >> 32); // exact conversion | |
810 | x_nr_bits = | |
811 | 32 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
812 | } else { // x < 2^32 | |
813 | tmp1.d = (double) (C1_lo); // exact conversion | |
814 | x_nr_bits = | |
815 | ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
816 | } | |
817 | } else { // if x < 2^53 | |
818 | tmp1.d = (double) C1_lo; // exact conversion | |
819 | x_nr_bits = | |
820 | ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
200359e8 | 821 | } |
b2a00c89 L |
822 | } else { // C1_hi != 0 => nr. bits = 64 + nr_bits (C1_hi) |
823 | tmp1.d = (double) C1_hi; // exact conversion | |
200359e8 | 824 | x_nr_bits = |
b2a00c89 | 825 | 64 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 L |
826 | } |
827 | ||
b2a00c89 | 828 | q1 = nr_digits[x_nr_bits].digits; |
200359e8 | 829 | if (q1 == 0) { |
b2a00c89 L |
830 | q1 = nr_digits[x_nr_bits].digits1; |
831 | if (C1_hi > nr_digits[x_nr_bits].threshold_hi || | |
832 | (C1_hi == nr_digits[x_nr_bits].threshold_hi && | |
833 | C1_lo >= nr_digits[x_nr_bits].threshold_lo)) | |
834 | q1++; | |
200359e8 L |
835 | } |
836 | // q2 = nr. of decimal digits in y | |
837 | // determine first the nr. of bits in y (y_nr_bits) | |
b2a00c89 L |
838 | if (C2_hi == 0) { // y_bits is the nr. of bits in C2_lo |
839 | if (C2_lo >= 0x0020000000000000ull) { // y >= 2^53 | |
840 | //split the 64-bit value in two 32-bit halves to avoid rounding errors | |
841 | if (C2_lo >= 0x0000000100000000ull) { // y >= 2^32 | |
842 | tmp2.d = (double) (C2_lo >> 32); // exact conversion | |
843 | y_nr_bits = | |
844 | 32 + ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
845 | } else { // y < 2^32 | |
846 | tmp2.d = (double) (C2_lo); // exact conversion | |
847 | y_nr_bits = | |
848 | ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
849 | } | |
850 | } else { // if y < 2^53 | |
851 | tmp2.d = (double) C2_lo; // exact conversion | |
852 | y_nr_bits = | |
853 | ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
200359e8 | 854 | } |
b2a00c89 L |
855 | } else { // C2_hi != 0 => nr. bits = 64 + nr_bits (C2_hi) |
856 | tmp2.d = (double) C2_hi; // exact conversion | |
200359e8 | 857 | y_nr_bits = |
b2a00c89 | 858 | 64 + ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 L |
859 | } |
860 | ||
b2a00c89 | 861 | q2 = nr_digits[y_nr_bits].digits; |
200359e8 | 862 | if (q2 == 0) { |
b2a00c89 L |
863 | q2 = nr_digits[y_nr_bits].digits1; |
864 | if (C2_hi > nr_digits[y_nr_bits].threshold_hi || | |
865 | (C2_hi == nr_digits[y_nr_bits].threshold_hi && | |
866 | C2_lo >= nr_digits[y_nr_bits].threshold_lo)) | |
867 | q2++; | |
200359e8 L |
868 | } |
869 | ||
870 | delta = q1 + (int) (x_exp >> 49) - q2 - (int) (y_exp >> 49); | |
871 | ||
872 | if (delta >= P34) { | |
873 | // round the result directly because 0 < C2 < ulp (C1 * 10^(x_exp-e2)) | |
874 | // n = C1 * 10^e1 or n = C1 +/- 10^(q1-P34)) * 10^e1 | |
875 | // the result is inexact; the preferred exponent is the least possible | |
876 | ||
877 | if (delta >= P34 + 1) { | |
b2a00c89 L |
878 | // for RN the result is the operand with the larger magnitude, |
879 | // possibly scaled up by 10^(P34-q1) | |
880 | // an overflow cannot occur in this case (rounding to nearest) | |
881 | if (q1 < P34) { // scale C1 up by 10^(P34-q1) | |
882 | // Note: because delta >= P34+1 it is certain that | |
883 | // x_exp - ((UINT64)scale << 49) will stay above e_min | |
884 | scale = P34 - q1; | |
885 | if (q1 <= 19) { // C1 fits in 64 bits | |
886 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
887 | if (scale <= 19) { // 10^scale fits in 64 bits | |
888 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
889 | } else { // if 20 <= scale <= 33 | |
890 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
891 | // (C1 * 10^(scale-19)) fits in 64 bits | |
892 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
893 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
894 | } | |
895 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
896 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
897 | C1.w[1] = C1_hi; | |
898 | C1.w[0] = C1_lo; | |
899 | // C1 = ten2k64[P34 - q1] * C1 | |
900 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
901 | } | |
902 | x_exp = x_exp - ((UINT64) scale << 49); | |
903 | C1_hi = C1.w[1]; | |
904 | C1_lo = C1.w[0]; | |
905 | } | |
906 | // some special cases arise: if delta = P34 + 1 and C1 = 10^(P34-1) | |
907 | // (after scaling) and x_sign != y_sign and C2 > 5*10^(q2-1) => | |
908 | // subtract 1 ulp | |
909 | // Note: do this only for rounding to nearest; for other rounding | |
910 | // modes the correction will be applied next | |
911 | if ((rnd_mode == ROUNDING_TO_NEAREST | |
912 | || rnd_mode == ROUNDING_TIES_AWAY) && delta == (P34 + 1) | |
913 | && C1_hi == 0x0000314dc6448d93ull | |
914 | && C1_lo == 0x38c15b0a00000000ull && x_sign != y_sign | |
915 | && ((q2 <= 19 && C2_lo > midpoint64[q2 - 1]) || (q2 >= 20 | |
916 | && (C2_hi > | |
917 | midpoint128 | |
918 | [q2 - | |
919 | 20]. | |
920 | w[1] | |
921 | || | |
922 | (C2_hi | |
923 | == | |
924 | midpoint128 | |
925 | [q2 - | |
926 | 20]. | |
927 | w[1] | |
928 | && | |
929 | C2_lo | |
930 | > | |
931 | midpoint128 | |
932 | [q2 - | |
933 | 20]. | |
934 | w | |
935 | [0]))))) | |
936 | { | |
937 | // C1 = 10^34 - 1 and decrement x_exp by 1 (no underflow possible) | |
938 | C1_hi = 0x0001ed09bead87c0ull; | |
939 | C1_lo = 0x378d8e63ffffffffull; | |
940 | x_exp = x_exp - EXP_P1; | |
941 | } | |
942 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
943 | if ((rnd_mode == ROUNDING_DOWN && x_sign && y_sign) || | |
944 | (rnd_mode == ROUNDING_UP && !x_sign && !y_sign)) { | |
945 | // add 1 ulp and then check for overflow | |
946 | C1_lo = C1_lo + 1; | |
947 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
948 | C1_hi = C1_hi + 1; | |
949 | } | |
950 | if (C1_hi == 0x0001ed09bead87c0ull | |
951 | && C1_lo == 0x378d8e6400000000ull) { | |
952 | // C1 = 10^34 => rounding overflow | |
953 | C1_hi = 0x0000314dc6448d93ull; | |
954 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
955 | x_exp = x_exp + EXP_P1; | |
956 | if (x_exp == EXP_MAX_P1) { // overflow | |
957 | C1_hi = 0x7800000000000000ull; // +inf | |
958 | C1_lo = 0x0ull; | |
959 | x_exp = 0; // x_sign is preserved | |
960 | // set overflow flag (the inexact flag was set too) | |
961 | *pfpsf |= OVERFLOW_EXCEPTION; | |
962 | } | |
963 | } | |
964 | } else if ((rnd_mode == ROUNDING_DOWN && !x_sign && y_sign) || | |
965 | (rnd_mode == ROUNDING_UP && x_sign && !y_sign) || | |
966 | (rnd_mode == ROUNDING_TO_ZERO | |
967 | && x_sign != y_sign)) { | |
968 | // subtract 1 ulp from C1 | |
969 | // Note: because delta >= P34 + 1 the result cannot be zero | |
970 | C1_lo = C1_lo - 1; | |
971 | if (C1_lo == 0xffffffffffffffffull) | |
972 | C1_hi = C1_hi - 1; | |
973 | // if the coefficient is 10^33 - 1 then make it 10^34 - 1 and | |
974 | // decrease the exponent by 1 (because delta >= P34 + 1 the | |
975 | // exponent will not become less than e_min) | |
976 | // 10^33 - 1 = 0x0000314dc6448d9338c15b09ffffffff | |
977 | // 10^34 - 1 = 0x0001ed09bead87c0378d8e63ffffffff | |
978 | if (C1_hi == 0x0000314dc6448d93ull | |
979 | && C1_lo == 0x38c15b09ffffffffull) { | |
980 | // make C1 = 10^34 - 1 | |
981 | C1_hi = 0x0001ed09bead87c0ull; | |
982 | C1_lo = 0x378d8e63ffffffffull; | |
983 | x_exp = x_exp - EXP_P1; | |
984 | } | |
985 | } else { | |
986 | ; // the result is already correct | |
987 | } | |
988 | } | |
989 | // set the inexact flag | |
990 | *pfpsf |= INEXACT_EXCEPTION; | |
991 | // assemble the result | |
992 | res.w[1] = x_sign | x_exp | C1_hi; | |
993 | res.w[0] = C1_lo; | |
994 | } else { // delta = P34 | |
995 | // in most cases, the smaller operand may be < or = or > 1/2 ulp of the | |
996 | // larger operand | |
997 | // however, the case C1 = 10^(q1-1) and x_sign != y_sign is special due | |
998 | // to accuracy loss after subtraction, and will be treated separately | |
999 | if (x_sign == y_sign || (q1 <= 20 | |
1000 | && (C1_hi != 0 | |
1001 | || C1_lo != ten2k64[q1 - 1])) | |
1002 | || (q1 >= 21 && (C1_hi != ten2k128[q1 - 21].w[1] | |
1003 | || C1_lo != ten2k128[q1 - 21].w[0]))) { | |
1004 | // if x_sign == y_sign or C1 != 10^(q1-1) | |
1005 | // compare C2 with 1/2 ulp = 5 * 10^(q2-1), the latter read from table | |
1006 | // Note: cases q1<=19 and q1>=20 can be coalesced at some latency cost | |
1007 | if (q2 <= 19) { // C2 and 5*10^(q2-1) both fit in 64 bits | |
1008 | halfulp64 = midpoint64[q2 - 1]; // 5 * 10^(q2-1) | |
1009 | if (C2_lo < halfulp64) { // n2 < 1/2 ulp (n1) | |
1010 | // for RN the result is the operand with the larger magnitude, | |
1011 | // possibly scaled up by 10^(P34-q1) | |
1012 | // an overflow cannot occur in this case (rounding to nearest) | |
1013 | if (q1 < P34) { // scale C1 up by 10^(P34-q1) | |
1014 | // Note: because delta = P34 it is certain that | |
1015 | // x_exp - ((UINT64)scale << 49) will stay above e_min | |
1016 | scale = P34 - q1; | |
1017 | if (q1 <= 19) { // C1 fits in 64 bits | |
1018 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1019 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1020 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1021 | } else { // if 20 <= scale <= 33 | |
1022 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1023 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1024 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1025 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1026 | } | |
1027 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1028 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1029 | C1.w[1] = C1_hi; | |
1030 | C1.w[0] = C1_lo; | |
1031 | // C1 = ten2k64[P34 - q1] * C1 | |
1032 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1033 | } | |
1034 | x_exp = x_exp - ((UINT64) scale << 49); | |
1035 | C1_hi = C1.w[1]; | |
1036 | C1_lo = C1.w[0]; | |
1037 | } | |
1038 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
1039 | if ((rnd_mode == ROUNDING_DOWN && x_sign && y_sign) || | |
1040 | (rnd_mode == ROUNDING_UP && !x_sign && !y_sign)) { | |
1041 | // add 1 ulp and then check for overflow | |
1042 | C1_lo = C1_lo + 1; | |
1043 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1044 | C1_hi = C1_hi + 1; | |
1045 | } | |
1046 | if (C1_hi == 0x0001ed09bead87c0ull | |
1047 | && C1_lo == 0x378d8e6400000000ull) { | |
1048 | // C1 = 10^34 => rounding overflow | |
1049 | C1_hi = 0x0000314dc6448d93ull; | |
1050 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1051 | x_exp = x_exp + EXP_P1; | |
1052 | if (x_exp == EXP_MAX_P1) { // overflow | |
1053 | C1_hi = 0x7800000000000000ull; // +inf | |
1054 | C1_lo = 0x0ull; | |
1055 | x_exp = 0; // x_sign is preserved | |
1056 | // set overflow flag (the inexact flag was set too) | |
1057 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1058 | } | |
1059 | } | |
1060 | } else | |
1061 | if ((rnd_mode == ROUNDING_DOWN && !x_sign && y_sign) | |
1062 | || (rnd_mode == ROUNDING_UP && x_sign && !y_sign) | |
1063 | || (rnd_mode == ROUNDING_TO_ZERO | |
1064 | && x_sign != y_sign)) { | |
1065 | // subtract 1 ulp from C1 | |
1066 | // Note: because delta >= P34 + 1 the result cannot be zero | |
1067 | C1_lo = C1_lo - 1; | |
1068 | if (C1_lo == 0xffffffffffffffffull) | |
1069 | C1_hi = C1_hi - 1; | |
1070 | // if the coefficient is 10^33-1 then make it 10^34-1 and | |
1071 | // decrease the exponent by 1 (because delta >= P34 + 1 the | |
1072 | // exponent will not become less than e_min) | |
1073 | // 10^33 - 1 = 0x0000314dc6448d9338c15b09ffffffff | |
1074 | // 10^34 - 1 = 0x0001ed09bead87c0378d8e63ffffffff | |
1075 | if (C1_hi == 0x0000314dc6448d93ull | |
1076 | && C1_lo == 0x38c15b09ffffffffull) { | |
1077 | // make C1 = 10^34 - 1 | |
1078 | C1_hi = 0x0001ed09bead87c0ull; | |
1079 | C1_lo = 0x378d8e63ffffffffull; | |
1080 | x_exp = x_exp - EXP_P1; | |
1081 | } | |
1082 | } else { | |
1083 | ; // the result is already correct | |
1084 | } | |
1085 | } | |
1086 | // set the inexact flag | |
1087 | *pfpsf |= INEXACT_EXCEPTION; | |
1088 | // assemble the result | |
1089 | res.w[1] = x_sign | x_exp | C1_hi; | |
1090 | res.w[0] = C1_lo; | |
1091 | } else if ((C2_lo == halfulp64) | |
1092 | && (q1 < P34 || ((C1_lo & 0x1) == 0))) { | |
1093 | // n2 = 1/2 ulp (n1) and C1 is even | |
1094 | // the result is the operand with the larger magnitude, | |
1095 | // possibly scaled up by 10^(P34-q1) | |
1096 | // an overflow cannot occur in this case (rounding to nearest) | |
1097 | if (q1 < P34) { // scale C1 up by 10^(P34-q1) | |
1098 | // Note: because delta = P34 it is certain that | |
1099 | // x_exp - ((UINT64)scale << 49) will stay above e_min | |
1100 | scale = P34 - q1; | |
1101 | if (q1 <= 19) { // C1 fits in 64 bits | |
1102 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1103 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1104 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1105 | } else { // if 20 <= scale <= 33 | |
1106 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1107 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1108 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1109 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1110 | } | |
1111 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1112 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1113 | C1.w[1] = C1_hi; | |
1114 | C1.w[0] = C1_lo; | |
1115 | // C1 = ten2k64[P34 - q1] * C1 | |
1116 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1117 | } | |
1118 | x_exp = x_exp - ((UINT64) scale << 49); | |
1119 | C1_hi = C1.w[1]; | |
1120 | C1_lo = C1.w[0]; | |
1121 | } | |
1122 | if ((rnd_mode == ROUNDING_TO_NEAREST && x_sign == y_sign | |
1123 | && (C1_lo & 0x01)) || (rnd_mode == ROUNDING_TIES_AWAY | |
1124 | && x_sign == y_sign) | |
1125 | || (rnd_mode == ROUNDING_UP && !x_sign && !y_sign) | |
1126 | || (rnd_mode == ROUNDING_DOWN && x_sign && y_sign)) { | |
1127 | // add 1 ulp and then check for overflow | |
1128 | C1_lo = C1_lo + 1; | |
1129 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1130 | C1_hi = C1_hi + 1; | |
1131 | } | |
1132 | if (C1_hi == 0x0001ed09bead87c0ull | |
1133 | && C1_lo == 0x378d8e6400000000ull) { | |
1134 | // C1 = 10^34 => rounding overflow | |
1135 | C1_hi = 0x0000314dc6448d93ull; | |
1136 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1137 | x_exp = x_exp + EXP_P1; | |
1138 | if (x_exp == EXP_MAX_P1) { // overflow | |
1139 | C1_hi = 0x7800000000000000ull; // +inf | |
1140 | C1_lo = 0x0ull; | |
1141 | x_exp = 0; // x_sign is preserved | |
1142 | // set overflow flag (the inexact flag was set too) | |
1143 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1144 | } | |
1145 | } | |
1146 | } else | |
1147 | if ((rnd_mode == ROUNDING_TO_NEAREST && x_sign != y_sign | |
1148 | && (C1_lo & 0x01)) || (rnd_mode == ROUNDING_DOWN | |
1149 | && !x_sign && y_sign) | |
1150 | || (rnd_mode == ROUNDING_UP && x_sign && !y_sign) | |
1151 | || (rnd_mode == ROUNDING_TO_ZERO | |
1152 | && x_sign != y_sign)) { | |
1153 | // subtract 1 ulp from C1 | |
1154 | // Note: because delta >= P34 + 1 the result cannot be zero | |
1155 | C1_lo = C1_lo - 1; | |
1156 | if (C1_lo == 0xffffffffffffffffull) | |
1157 | C1_hi = C1_hi - 1; | |
1158 | // if the coefficient is 10^33 - 1 then make it 10^34 - 1 | |
1159 | // and decrease the exponent by 1 (because delta >= P34 + 1 | |
1160 | // the exponent will not become less than e_min) | |
1161 | // 10^33 - 1 = 0x0000314dc6448d9338c15b09ffffffff | |
1162 | // 10^34 - 1 = 0x0001ed09bead87c0378d8e63ffffffff | |
1163 | if (C1_hi == 0x0000314dc6448d93ull | |
1164 | && C1_lo == 0x38c15b09ffffffffull) { | |
1165 | // make C1 = 10^34 - 1 | |
1166 | C1_hi = 0x0001ed09bead87c0ull; | |
1167 | C1_lo = 0x378d8e63ffffffffull; | |
1168 | x_exp = x_exp - EXP_P1; | |
1169 | } | |
1170 | } else { | |
1171 | ; // the result is already correct | |
1172 | } | |
1173 | // set the inexact flag | |
1174 | *pfpsf |= INEXACT_EXCEPTION; | |
1175 | // assemble the result | |
1176 | res.w[1] = x_sign | x_exp | C1_hi; | |
1177 | res.w[0] = C1_lo; | |
1178 | } else { // if C2_lo > halfulp64 || | |
1179 | // (C2_lo == halfulp64 && q1 == P34 && ((C1_lo & 0x1) == 1)), i.e. | |
1180 | // 1/2 ulp(n1) < n2 < 1 ulp(n1) or n2 = 1/2 ulp(n1) and C1 odd | |
1181 | // res = x+1 ulp if n1*n2 > 0 and res = x-1 ulp if n1*n2 < 0 | |
1182 | if (q1 < P34) { // then 1 ulp = 10^(e1+q1-P34) < 10^e1 | |
1183 | // Note: if (q1 == P34) then 1 ulp = 10^(e1+q1-P34) = 10^e1 | |
1184 | // because q1 < P34 we must first replace C1 by | |
1185 | // C1 * 10^(P34-q1), and must decrease the exponent by | |
1186 | // (P34-q1) (it will still be at least e_min) | |
1187 | scale = P34 - q1; | |
1188 | if (q1 <= 19) { // C1 fits in 64 bits | |
1189 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1190 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1191 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1192 | } else { // if 20 <= scale <= 33 | |
1193 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1194 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1195 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1196 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1197 | } | |
1198 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1199 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1200 | C1.w[1] = C1_hi; | |
1201 | C1.w[0] = C1_lo; | |
1202 | // C1 = ten2k64[P34 - q1] * C1 | |
1203 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1204 | } | |
1205 | x_exp = x_exp - ((UINT64) scale << 49); | |
1206 | C1_hi = C1.w[1]; | |
1207 | C1_lo = C1.w[0]; | |
1208 | // check for rounding overflow | |
1209 | if (C1_hi == 0x0001ed09bead87c0ull | |
1210 | && C1_lo == 0x378d8e6400000000ull) { | |
1211 | // C1 = 10^34 => rounding overflow | |
1212 | C1_hi = 0x0000314dc6448d93ull; | |
1213 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1214 | x_exp = x_exp + EXP_P1; | |
1215 | } | |
1216 | } | |
1217 | if ((rnd_mode == ROUNDING_TO_NEAREST && x_sign != y_sign) | |
1218 | || (rnd_mode == ROUNDING_TIES_AWAY && x_sign != y_sign | |
1219 | && C2_lo != halfulp64) | |
1220 | || (rnd_mode == ROUNDING_DOWN && !x_sign && y_sign) | |
1221 | || (rnd_mode == ROUNDING_UP && x_sign && !y_sign) | |
1222 | || (rnd_mode == ROUNDING_TO_ZERO | |
1223 | && x_sign != y_sign)) { | |
1224 | // the result is x - 1 | |
1225 | // for RN n1 * n2 < 0; underflow not possible | |
1226 | C1_lo = C1_lo - 1; | |
1227 | if (C1_lo == 0xffffffffffffffffull) | |
1228 | C1_hi--; | |
1229 | // check if we crossed into the lower decade | |
1230 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
1231 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
1232 | C1_lo = 0x378d8e63ffffffffull; | |
1233 | x_exp = x_exp - EXP_P1; // no underflow, because n1 >> n2 | |
1234 | } | |
1235 | } else | |
1236 | if ((rnd_mode == ROUNDING_TO_NEAREST | |
1237 | && x_sign == y_sign) | |
1238 | || (rnd_mode == ROUNDING_TIES_AWAY | |
1239 | && x_sign == y_sign) | |
1240 | || (rnd_mode == ROUNDING_DOWN && x_sign && y_sign) | |
1241 | || (rnd_mode == ROUNDING_UP && !x_sign | |
1242 | && !y_sign)) { | |
1243 | // the result is x + 1 | |
1244 | // for RN x_sign = y_sign, i.e. n1*n2 > 0 | |
1245 | C1_lo = C1_lo + 1; | |
1246 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1247 | C1_hi = C1_hi + 1; | |
1248 | } | |
1249 | if (C1_hi == 0x0001ed09bead87c0ull | |
1250 | && C1_lo == 0x378d8e6400000000ull) { | |
1251 | // C1 = 10^34 => rounding overflow | |
1252 | C1_hi = 0x0000314dc6448d93ull; | |
1253 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1254 | x_exp = x_exp + EXP_P1; | |
1255 | if (x_exp == EXP_MAX_P1) { // overflow | |
1256 | C1_hi = 0x7800000000000000ull; // +inf | |
1257 | C1_lo = 0x0ull; | |
1258 | x_exp = 0; // x_sign is preserved | |
1259 | // set the overflow flag | |
1260 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1261 | } | |
1262 | } | |
1263 | } else { | |
1264 | ; // the result is x | |
1265 | } | |
1266 | // set the inexact flag | |
1267 | *pfpsf |= INEXACT_EXCEPTION; | |
1268 | // assemble the result | |
1269 | res.w[1] = x_sign | x_exp | C1_hi; | |
1270 | res.w[0] = C1_lo; | |
1271 | } | |
1272 | } else { // if q2 >= 20 then 5*10^(q2-1) and C2 (the latter in | |
1273 | // most cases) fit only in more than 64 bits | |
1274 | halfulp128 = midpoint128[q2 - 20]; // 5 * 10^(q2-1) | |
1275 | if ((C2_hi < halfulp128.w[1]) | |
1276 | || (C2_hi == halfulp128.w[1] | |
1277 | && C2_lo < halfulp128.w[0])) { | |
1278 | // n2 < 1/2 ulp (n1) | |
1279 | // the result is the operand with the larger magnitude, | |
1280 | // possibly scaled up by 10^(P34-q1) | |
1281 | // an overflow cannot occur in this case (rounding to nearest) | |
1282 | if (q1 < P34) { // scale C1 up by 10^(P34-q1) | |
1283 | // Note: because delta = P34 it is certain that | |
1284 | // x_exp - ((UINT64)scale << 49) will stay above e_min | |
1285 | scale = P34 - q1; | |
1286 | if (q1 <= 19) { // C1 fits in 64 bits | |
1287 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1288 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1289 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1290 | } else { // if 20 <= scale <= 33 | |
1291 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1292 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1293 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1294 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1295 | } | |
1296 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1297 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1298 | C1.w[1] = C1_hi; | |
1299 | C1.w[0] = C1_lo; | |
1300 | // C1 = ten2k64[P34 - q1] * C1 | |
1301 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1302 | } | |
1303 | C1_hi = C1.w[1]; | |
1304 | C1_lo = C1.w[0]; | |
1305 | x_exp = x_exp - ((UINT64) scale << 49); | |
1306 | } | |
1307 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
1308 | if ((rnd_mode == ROUNDING_DOWN && x_sign && y_sign) || | |
1309 | (rnd_mode == ROUNDING_UP && !x_sign && !y_sign)) { | |
1310 | // add 1 ulp and then check for overflow | |
1311 | C1_lo = C1_lo + 1; | |
1312 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1313 | C1_hi = C1_hi + 1; | |
1314 | } | |
1315 | if (C1_hi == 0x0001ed09bead87c0ull | |
1316 | && C1_lo == 0x378d8e6400000000ull) { | |
1317 | // C1 = 10^34 => rounding overflow | |
1318 | C1_hi = 0x0000314dc6448d93ull; | |
1319 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1320 | x_exp = x_exp + EXP_P1; | |
1321 | if (x_exp == EXP_MAX_P1) { // overflow | |
1322 | C1_hi = 0x7800000000000000ull; // +inf | |
1323 | C1_lo = 0x0ull; | |
1324 | x_exp = 0; // x_sign is preserved | |
1325 | // set overflow flag (the inexact flag was set too) | |
1326 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1327 | } | |
1328 | } | |
1329 | } else | |
1330 | if ((rnd_mode == ROUNDING_DOWN && !x_sign && y_sign) | |
1331 | || (rnd_mode == ROUNDING_UP && x_sign && !y_sign) | |
1332 | || (rnd_mode == ROUNDING_TO_ZERO | |
1333 | && x_sign != y_sign)) { | |
1334 | // subtract 1 ulp from C1 | |
1335 | // Note: because delta >= P34 + 1 the result cannot be zero | |
1336 | C1_lo = C1_lo - 1; | |
1337 | if (C1_lo == 0xffffffffffffffffull) | |
1338 | C1_hi = C1_hi - 1; | |
1339 | // if the coefficient is 10^33-1 then make it 10^34-1 and | |
1340 | // decrease the exponent by 1 (because delta >= P34 + 1 the | |
1341 | // exponent will not become less than e_min) | |
1342 | // 10^33 - 1 = 0x0000314dc6448d9338c15b09ffffffff | |
1343 | // 10^34 - 1 = 0x0001ed09bead87c0378d8e63ffffffff | |
1344 | if (C1_hi == 0x0000314dc6448d93ull | |
1345 | && C1_lo == 0x38c15b09ffffffffull) { | |
1346 | // make C1 = 10^34 - 1 | |
1347 | C1_hi = 0x0001ed09bead87c0ull; | |
1348 | C1_lo = 0x378d8e63ffffffffull; | |
1349 | x_exp = x_exp - EXP_P1; | |
1350 | } | |
1351 | } else { | |
1352 | ; // the result is already correct | |
1353 | } | |
1354 | } | |
1355 | // set the inexact flag | |
1356 | *pfpsf |= INEXACT_EXCEPTION; | |
1357 | // assemble the result | |
1358 | res.w[1] = x_sign | x_exp | C1_hi; | |
1359 | res.w[0] = C1_lo; | |
1360 | } else if ((C2_hi == halfulp128.w[1] | |
1361 | && C2_lo == halfulp128.w[0]) | |
1362 | && (q1 < P34 || ((C1_lo & 0x1) == 0))) { | |
1363 | // midpoint & lsb in C1 is 0 | |
1364 | // n2 = 1/2 ulp (n1) and C1 is even | |
1365 | // the result is the operand with the larger magnitude, | |
1366 | // possibly scaled up by 10^(P34-q1) | |
1367 | // an overflow cannot occur in this case (rounding to nearest) | |
1368 | if (q1 < P34) { // scale C1 up by 10^(P34-q1) | |
1369 | // Note: because delta = P34 it is certain that | |
1370 | // x_exp - ((UINT64)scale << 49) will stay above e_min | |
1371 | scale = P34 - q1; | |
1372 | if (q1 <= 19) { // C1 fits in 64 bits | |
1373 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1374 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1375 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1376 | } else { // if 20 <= scale <= 33 | |
1377 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1378 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1379 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1380 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1381 | } | |
1382 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1383 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1384 | C1.w[1] = C1_hi; | |
1385 | C1.w[0] = C1_lo; | |
1386 | // C1 = ten2k64[P34 - q1] * C1 | |
1387 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1388 | } | |
1389 | x_exp = x_exp - ((UINT64) scale << 49); | |
1390 | C1_hi = C1.w[1]; | |
1391 | C1_lo = C1.w[0]; | |
1392 | } | |
1393 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
1394 | if ((rnd_mode == ROUNDING_TIES_AWAY && x_sign == y_sign) | |
1395 | || (rnd_mode == ROUNDING_UP && !y_sign)) { | |
1396 | // add 1 ulp and then check for overflow | |
1397 | C1_lo = C1_lo + 1; | |
1398 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1399 | C1_hi = C1_hi + 1; | |
1400 | } | |
1401 | if (C1_hi == 0x0001ed09bead87c0ull | |
1402 | && C1_lo == 0x378d8e6400000000ull) { | |
1403 | // C1 = 10^34 => rounding overflow | |
1404 | C1_hi = 0x0000314dc6448d93ull; | |
1405 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1406 | x_exp = x_exp + EXP_P1; | |
1407 | if (x_exp == EXP_MAX_P1) { // overflow | |
1408 | C1_hi = 0x7800000000000000ull; // +inf | |
1409 | C1_lo = 0x0ull; | |
1410 | x_exp = 0; // x_sign is preserved | |
1411 | // set overflow flag (the inexact flag was set too) | |
1412 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1413 | } | |
1414 | } | |
1415 | } else if ((rnd_mode == ROUNDING_DOWN && y_sign) | |
1416 | || (rnd_mode == ROUNDING_TO_ZERO | |
1417 | && x_sign != y_sign)) { | |
1418 | // subtract 1 ulp from C1 | |
1419 | // Note: because delta >= P34 + 1 the result cannot be zero | |
1420 | C1_lo = C1_lo - 1; | |
1421 | if (C1_lo == 0xffffffffffffffffull) | |
1422 | C1_hi = C1_hi - 1; | |
1423 | // if the coefficient is 10^33 - 1 then make it 10^34 - 1 | |
1424 | // and decrease the exponent by 1 (because delta >= P34 + 1 | |
1425 | // the exponent will not become less than e_min) | |
1426 | // 10^33 - 1 = 0x0000314dc6448d9338c15b09ffffffff | |
1427 | // 10^34 - 1 = 0x0001ed09bead87c0378d8e63ffffffff | |
1428 | if (C1_hi == 0x0000314dc6448d93ull | |
1429 | && C1_lo == 0x38c15b09ffffffffull) { | |
1430 | // make C1 = 10^34 - 1 | |
1431 | C1_hi = 0x0001ed09bead87c0ull; | |
1432 | C1_lo = 0x378d8e63ffffffffull; | |
1433 | x_exp = x_exp - EXP_P1; | |
1434 | } | |
1435 | } else { | |
1436 | ; // the result is already correct | |
1437 | } | |
1438 | } | |
1439 | // set the inexact flag | |
1440 | *pfpsf |= INEXACT_EXCEPTION; | |
1441 | // assemble the result | |
1442 | res.w[1] = x_sign | x_exp | C1_hi; | |
1443 | res.w[0] = C1_lo; | |
1444 | } else { // if C2 > halfulp128 || | |
1445 | // (C2 == halfulp128 && q1 == P34 && ((C1 & 0x1) == 1)), i.e. | |
1446 | // 1/2 ulp(n1) < n2 < 1 ulp(n1) or n2 = 1/2 ulp(n1) and C1 odd | |
1447 | // res = x+1 ulp if n1*n2 > 0 and res = x-1 ulp if n1*n2 < 0 | |
1448 | if (q1 < P34) { // then 1 ulp = 10^(e1+q1-P34) < 10^e1 | |
1449 | // Note: if (q1 == P34) then 1 ulp = 10^(e1+q1-P34) = 10^e1 | |
1450 | // because q1 < P34 we must first replace C1 by C1*10^(P34-q1), | |
1451 | // and must decrease the exponent by (P34-q1) (it will still be | |
1452 | // at least e_min) | |
1453 | scale = P34 - q1; | |
1454 | if (q1 <= 19) { // C1 fits in 64 bits | |
1455 | // 1 <= q1 <= 19 => 15 <= scale <= 33 | |
1456 | if (scale <= 19) { // 10^scale fits in 64 bits | |
1457 | __mul_64x64_to_128MACH (C1, ten2k64[scale], C1_lo); | |
1458 | } else { // if 20 <= scale <= 33 | |
1459 | // C1 * 10^scale = (C1 * 10^(scale-19)) * 10^19 where | |
1460 | // (C1 * 10^(scale-19)) fits in 64 bits | |
1461 | C1_lo = C1_lo * ten2k64[scale - 19]; | |
1462 | __mul_64x64_to_128MACH (C1, ten2k64[19], C1_lo); | |
1463 | } | |
1464 | } else { //if 20 <= q1 <= 33=P34-1 then C1 fits only in 128 bits | |
1465 | // => 1 <= P34 - q1 <= 14 so 10^(P34-q1) fits in 64 bits | |
1466 | C1.w[1] = C1_hi; | |
1467 | C1.w[0] = C1_lo; | |
1468 | // C1 = ten2k64[P34 - q1] * C1 | |
1469 | __mul_128x64_to_128 (C1, ten2k64[P34 - q1], C1); | |
1470 | } | |
1471 | C1_hi = C1.w[1]; | |
1472 | C1_lo = C1.w[0]; | |
1473 | x_exp = x_exp - ((UINT64) scale << 49); | |
1474 | } | |
1475 | if ((rnd_mode == ROUNDING_TO_NEAREST && x_sign != y_sign) | |
1476 | || (rnd_mode == ROUNDING_TIES_AWAY && x_sign != y_sign | |
1477 | && (C2_hi != halfulp128.w[1] | |
1478 | || C2_lo != halfulp128.w[0])) | |
1479 | || (rnd_mode == ROUNDING_DOWN && !x_sign && y_sign) | |
1480 | || (rnd_mode == ROUNDING_UP && x_sign && !y_sign) | |
1481 | || (rnd_mode == ROUNDING_TO_ZERO | |
1482 | && x_sign != y_sign)) { | |
1483 | // the result is x - 1 | |
1484 | // for RN n1 * n2 < 0; underflow not possible | |
1485 | C1_lo = C1_lo - 1; | |
1486 | if (C1_lo == 0xffffffffffffffffull) | |
1487 | C1_hi--; | |
1488 | // check if we crossed into the lower decade | |
1489 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
1490 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
1491 | C1_lo = 0x378d8e63ffffffffull; | |
1492 | x_exp = x_exp - EXP_P1; // no underflow, because n1 >> n2 | |
1493 | } | |
1494 | } else | |
1495 | if ((rnd_mode == ROUNDING_TO_NEAREST | |
1496 | && x_sign == y_sign) | |
1497 | || (rnd_mode == ROUNDING_TIES_AWAY | |
1498 | && x_sign == y_sign) | |
1499 | || (rnd_mode == ROUNDING_DOWN && x_sign && y_sign) | |
1500 | || (rnd_mode == ROUNDING_UP && !x_sign | |
1501 | && !y_sign)) { | |
1502 | // the result is x + 1 | |
1503 | // for RN x_sign = y_sign, i.e. n1*n2 > 0 | |
1504 | C1_lo = C1_lo + 1; | |
1505 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1506 | C1_hi = C1_hi + 1; | |
1507 | } | |
1508 | if (C1_hi == 0x0001ed09bead87c0ull | |
1509 | && C1_lo == 0x378d8e6400000000ull) { | |
1510 | // C1 = 10^34 => rounding overflow | |
1511 | C1_hi = 0x0000314dc6448d93ull; | |
1512 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1513 | x_exp = x_exp + EXP_P1; | |
1514 | if (x_exp == EXP_MAX_P1) { // overflow | |
1515 | C1_hi = 0x7800000000000000ull; // +inf | |
1516 | C1_lo = 0x0ull; | |
1517 | x_exp = 0; // x_sign is preserved | |
1518 | // set the overflow flag | |
1519 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1520 | } | |
1521 | } | |
1522 | } else { | |
1523 | ; // the result is x | |
1524 | } | |
1525 | // set the inexact flag | |
1526 | *pfpsf |= INEXACT_EXCEPTION; | |
1527 | // assemble the result | |
1528 | res.w[1] = x_sign | x_exp | C1_hi; | |
1529 | res.w[0] = C1_lo; | |
1530 | } | |
1531 | } // end q1 >= 20 | |
1532 | // end case where C1 != 10^(q1-1) | |
1533 | } else { // C1 = 10^(q1-1) and x_sign != y_sign | |
1534 | // instead of C' = (C1 * 10^(e1-e2) + C2)rnd,P34 | |
1535 | // calculate C' = C1 * 10^(e1-e2-x1) + (C2 * 10^(-x1))rnd,P34 | |
1536 | // where x1 = q2 - 1, 0 <= x1 <= P34 - 1 | |
1537 | // Because C1 = 10^(q1-1) and x_sign != y_sign, C' will have P34 | |
1538 | // digits and n = C' * 10^(e2+x1) | |
1539 | // If the result has P34+1 digits, redo the steps above with x1+1 | |
1540 | // If the result has P34-1 digits or less, redo the steps above with | |
1541 | // x1-1 but only if initially x1 >= 1 | |
1542 | // NOTE: these two steps can be improved, e.g we could guess if | |
1543 | // P34+1 or P34-1 digits will be obtained by adding/subtracting | |
1544 | // just the top 64 bits of the two operands | |
1545 | // The result cannot be zero, and it cannot overflow | |
1546 | x1 = q2 - 1; // 0 <= x1 <= P34-1 | |
1547 | // Calculate C1 * 10^(e1-e2-x1) where 1 <= e1-e2-x1 <= P34 | |
1548 | // scale = (int)(e1 >> 49) - (int)(e2 >> 49) - x1; 0 <= scale <= P34-1 | |
1549 | scale = P34 - q1 + 1; // scale=e1-e2-x1 = P34+1-q1; 1<=scale<=P34 | |
1550 | // either C1 or 10^(e1-e2-x1) may not fit is 64 bits, | |
1551 | // but their product fits with certainty in 128 bits | |
1552 | if (scale >= 20) { //10^(e1-e2-x1) doesn't fit in 64 bits, but C1 does | |
1553 | __mul_128x64_to_128 (C1, C1_lo, ten2k128[scale - 20]); | |
1554 | } else { // if (scale >= 1 | |
1555 | // if 1 <= scale <= 19 then 10^(e1-e2-x1) fits in 64 bits | |
1556 | if (q1 <= 19) { // C1 fits in 64 bits | |
1557 | __mul_64x64_to_128MACH (C1, C1_lo, ten2k64[scale]); | |
1558 | } else { // q1 >= 20 | |
1559 | C1.w[1] = C1_hi; | |
1560 | C1.w[0] = C1_lo; | |
1561 | __mul_128x64_to_128 (C1, ten2k64[scale], C1); | |
1562 | } | |
1563 | } | |
1564 | tmp64 = C1.w[0]; // C1.w[1], C1.w[0] contains C1 * 10^(e1-e2-x1) | |
1565 | ||
1566 | // now round C2 to q2-x1 = 1 decimal digit | |
1567 | // C2' = C2 + 1/2 * 10^x1 = C2 + 5 * 10^(x1-1) | |
1568 | ind = x1 - 1; // -1 <= ind <= P34 - 2 | |
1569 | if (ind >= 0) { // if (x1 >= 1) | |
1570 | C2.w[0] = C2_lo; | |
1571 | C2.w[1] = C2_hi; | |
1572 | if (ind <= 18) { | |
1573 | C2.w[0] = C2.w[0] + midpoint64[ind]; | |
1574 | if (C2.w[0] < C2_lo) | |
1575 | C2.w[1]++; | |
1576 | } else { // 19 <= ind <= 32 | |
1577 | C2.w[0] = C2.w[0] + midpoint128[ind - 19].w[0]; | |
1578 | C2.w[1] = C2.w[1] + midpoint128[ind - 19].w[1]; | |
1579 | if (C2.w[0] < C2_lo) | |
1580 | C2.w[1]++; | |
1581 | } | |
1582 | // the approximation of 10^(-x1) was rounded up to 118 bits | |
1583 | __mul_128x128_to_256 (R256, C2, ten2mk128[ind]); // R256 = C2*, f2* | |
1584 | // calculate C2* and f2* | |
1585 | // C2* is actually floor(C2*) in this case | |
1586 | // C2* and f2* need shifting and masking, as shown by | |
1587 | // shiftright128[] and maskhigh128[] | |
1588 | // the top Ex bits of 10^(-x1) are T* = ten2mk128trunc[ind], e.g. | |
1589 | // if x1=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
1590 | // if (0 < f2* < 10^(-x1)) then | |
1591 | // if floor(C1+C2*) is even then C2* = floor(C2*) - logical right | |
1592 | // shift; C2* has p decimal digits, correct by Prop. 1) | |
1593 | // else if floor(C1+C2*) is odd C2* = floor(C2*)-1 (logical right | |
1594 | // shift; C2* has p decimal digits, correct by Pr. 1) | |
1595 | // else | |
1596 | // C2* = floor(C2*) (logical right shift; C has p decimal digits, | |
1597 | // correct by Property 1) | |
1598 | // n = C2* * 10^(e2+x1) | |
1599 | ||
1600 | if (ind <= 2) { | |
1601 | highf2star.w[1] = 0x0; | |
1602 | highf2star.w[0] = 0x0; // low f2* ok | |
1603 | } else if (ind <= 21) { | |
1604 | highf2star.w[1] = 0x0; | |
1605 | highf2star.w[0] = R256.w[2] & maskhigh128[ind]; // low f2* ok | |
1606 | } else { | |
1607 | highf2star.w[1] = R256.w[3] & maskhigh128[ind]; | |
1608 | highf2star.w[0] = R256.w[2]; // low f2* is ok | |
1609 | } | |
1610 | // shift right C2* by Ex-128 = shiftright128[ind] | |
1611 | if (ind >= 3) { | |
1612 | shift = shiftright128[ind]; | |
1613 | if (shift < 64) { // 3 <= shift <= 63 | |
1614 | R256.w[2] = | |
1615 | (R256.w[2] >> shift) | (R256.w[3] << (64 - shift)); | |
1616 | R256.w[3] = (R256.w[3] >> shift); | |
1617 | } else { // 66 <= shift <= 102 | |
1618 | R256.w[2] = (R256.w[3] >> (shift - 64)); | |
1619 | R256.w[3] = 0x0ULL; | |
1620 | } | |
1621 | } | |
1622 | // redundant | |
1623 | is_inexact_lt_midpoint = 0; | |
1624 | is_inexact_gt_midpoint = 0; | |
1625 | is_midpoint_lt_even = 0; | |
1626 | is_midpoint_gt_even = 0; | |
1627 | // determine inexactness of the rounding of C2* | |
1628 | // (cannot be followed by a second rounding) | |
1629 | // if (0 < f2* - 1/2 < 10^(-x1)) then | |
1630 | // the result is exact | |
1631 | // else (if f2* - 1/2 > T* then) | |
1632 | // the result of is inexact | |
1633 | if (ind <= 2) { | |
1634 | if (R256.w[1] > 0x8000000000000000ull || | |
1635 | (R256.w[1] == 0x8000000000000000ull | |
1636 | && R256.w[0] > 0x0ull)) { | |
1637 | // f2* > 1/2 and the result may be exact | |
1638 | tmp64A = R256.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
1639 | if ((tmp64A > ten2mk128trunc[ind].w[1] | |
1640 | || (tmp64A == ten2mk128trunc[ind].w[1] | |
1641 | && R256.w[0] >= ten2mk128trunc[ind].w[0]))) { | |
1642 | // set the inexact flag | |
1643 | *pfpsf |= INEXACT_EXCEPTION; | |
1644 | // this rounding is applied to C2 only! | |
1645 | // x_sign != y_sign | |
1646 | is_inexact_gt_midpoint = 1; | |
1647 | } // else the result is exact | |
1648 | // rounding down, unless a midpoint in [ODD, EVEN] | |
1649 | } else { // the result is inexact; f2* <= 1/2 | |
1650 | // set the inexact flag | |
1651 | *pfpsf |= INEXACT_EXCEPTION; | |
1652 | // this rounding is applied to C2 only! | |
1653 | // x_sign != y_sign | |
1654 | is_inexact_lt_midpoint = 1; | |
1655 | } | |
1656 | } else if (ind <= 21) { // if 3 <= ind <= 21 | |
1657 | if (highf2star.w[1] > 0x0 || (highf2star.w[1] == 0x0 | |
1658 | && highf2star.w[0] > | |
1659 | onehalf128[ind]) | |
1660 | || (highf2star.w[1] == 0x0 | |
1661 | && highf2star.w[0] == onehalf128[ind] | |
1662 | && (R256.w[1] || R256.w[0]))) { | |
1663 | // f2* > 1/2 and the result may be exact | |
1664 | // Calculate f2* - 1/2 | |
1665 | tmp64A = highf2star.w[0] - onehalf128[ind]; | |
1666 | tmp64B = highf2star.w[1]; | |
1667 | if (tmp64A > highf2star.w[0]) | |
1668 | tmp64B--; | |
1669 | if (tmp64B || tmp64A | |
1670 | || R256.w[1] > ten2mk128trunc[ind].w[1] | |
1671 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
1672 | && R256.w[0] > ten2mk128trunc[ind].w[0])) { | |
1673 | // set the inexact flag | |
1674 | *pfpsf |= INEXACT_EXCEPTION; | |
1675 | // this rounding is applied to C2 only! | |
1676 | // x_sign != y_sign | |
1677 | is_inexact_gt_midpoint = 1; | |
1678 | } // else the result is exact | |
1679 | } else { // the result is inexact; f2* <= 1/2 | |
1680 | // set the inexact flag | |
1681 | *pfpsf |= INEXACT_EXCEPTION; | |
1682 | // this rounding is applied to C2 only! | |
1683 | // x_sign != y_sign | |
1684 | is_inexact_lt_midpoint = 1; | |
1685 | } | |
1686 | } else { // if 22 <= ind <= 33 | |
1687 | if (highf2star.w[1] > onehalf128[ind] | |
1688 | || (highf2star.w[1] == onehalf128[ind] | |
1689 | && (highf2star.w[0] || R256.w[1] | |
1690 | || R256.w[0]))) { | |
1691 | // f2* > 1/2 and the result may be exact | |
1692 | // Calculate f2* - 1/2 | |
1693 | // tmp64A = highf2star.w[0]; | |
1694 | tmp64B = highf2star.w[1] - onehalf128[ind]; | |
1695 | if (tmp64B || highf2star.w[0] | |
1696 | || R256.w[1] > ten2mk128trunc[ind].w[1] | |
1697 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
1698 | && R256.w[0] > ten2mk128trunc[ind].w[0])) { | |
1699 | // set the inexact flag | |
1700 | *pfpsf |= INEXACT_EXCEPTION; | |
1701 | // this rounding is applied to C2 only! | |
1702 | // x_sign != y_sign | |
1703 | is_inexact_gt_midpoint = 1; | |
1704 | } // else the result is exact | |
1705 | } else { // the result is inexact; f2* <= 1/2 | |
1706 | // set the inexact flag | |
1707 | *pfpsf |= INEXACT_EXCEPTION; | |
1708 | // this rounding is applied to C2 only! | |
1709 | // x_sign != y_sign | |
1710 | is_inexact_lt_midpoint = 1; | |
1711 | } | |
1712 | } | |
1713 | // check for midpoints after determining inexactness | |
1714 | if ((R256.w[1] || R256.w[0]) && (highf2star.w[1] == 0) | |
1715 | && (highf2star.w[0] == 0) | |
1716 | && (R256.w[1] < ten2mk128trunc[ind].w[1] | |
1717 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
1718 | && R256.w[0] <= ten2mk128trunc[ind].w[0]))) { | |
1719 | // the result is a midpoint | |
1720 | if ((tmp64 + R256.w[2]) & 0x01) { // MP in [EVEN, ODD] | |
1721 | // if floor(C2*) is odd C = floor(C2*) - 1; the result may be 0 | |
1722 | R256.w[2]--; | |
1723 | if (R256.w[2] == 0xffffffffffffffffull) | |
1724 | R256.w[3]--; | |
1725 | // this rounding is applied to C2 only! | |
1726 | // x_sign != y_sign | |
1727 | is_midpoint_lt_even = 1; | |
1728 | is_inexact_lt_midpoint = 0; | |
1729 | is_inexact_gt_midpoint = 0; | |
1730 | } else { | |
1731 | // else MP in [ODD, EVEN] | |
1732 | // this rounding is applied to C2 only! | |
1733 | // x_sign != y_sign | |
1734 | is_midpoint_gt_even = 1; | |
1735 | is_inexact_lt_midpoint = 0; | |
1736 | is_inexact_gt_midpoint = 0; | |
1737 | } | |
1738 | } | |
1739 | } else { // if (ind == -1) only when x1 = 0 | |
1740 | R256.w[2] = C2_lo; | |
1741 | R256.w[3] = C2_hi; | |
1742 | is_midpoint_lt_even = 0; | |
1743 | is_midpoint_gt_even = 0; | |
1744 | is_inexact_lt_midpoint = 0; | |
1745 | is_inexact_gt_midpoint = 0; | |
1746 | } | |
1747 | // and now subtract C1 * 10^(e1-e2-x1) - (C2 * 10^(-x1))rnd,P34 | |
1748 | // because x_sign != y_sign this last operation is exact | |
1749 | C1.w[0] = C1.w[0] - R256.w[2]; | |
1750 | C1.w[1] = C1.w[1] - R256.w[3]; | |
1751 | if (C1.w[0] > tmp64) | |
1752 | C1.w[1]--; // borrow | |
1753 | if (C1.w[1] >= 0x8000000000000000ull) { // negative coefficient! | |
1754 | C1.w[0] = ~C1.w[0]; | |
1755 | C1.w[0]++; | |
1756 | C1.w[1] = ~C1.w[1]; | |
1757 | if (C1.w[0] == 0x0) | |
1758 | C1.w[1]++; | |
1759 | tmp_sign = y_sign; // the result will have the sign of y | |
1760 | } else { | |
1761 | tmp_sign = x_sign; | |
1762 | } | |
1763 | // the difference has exactly P34 digits | |
1764 | x_sign = tmp_sign; | |
1765 | if (x1 >= 1) | |
1766 | y_exp = y_exp + ((UINT64) x1 << 49); | |
1767 | C1_hi = C1.w[1]; | |
1768 | C1_lo = C1.w[0]; | |
1769 | // general correction from RN to RA, RM, RP, RZ; result uses y_exp | |
1770 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
1771 | if ((!x_sign | |
1772 | && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) | |
1773 | || | |
1774 | ((rnd_mode == ROUNDING_TIES_AWAY | |
1775 | || rnd_mode == ROUNDING_UP) | |
1776 | && is_midpoint_gt_even))) || (x_sign | |
1777 | && | |
1778 | ((rnd_mode == | |
1779 | ROUNDING_DOWN | |
1780 | && | |
1781 | is_inexact_lt_midpoint) | |
1782 | || | |
1783 | ((rnd_mode == | |
1784 | ROUNDING_TIES_AWAY | |
1785 | || rnd_mode == | |
1786 | ROUNDING_DOWN) | |
1787 | && | |
1788 | is_midpoint_gt_even)))) | |
1789 | { | |
1790 | // C1 = C1 + 1 | |
1791 | C1_lo = C1_lo + 1; | |
1792 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
1793 | C1_hi = C1_hi + 1; | |
1794 | } | |
1795 | if (C1_hi == 0x0001ed09bead87c0ull | |
1796 | && C1_lo == 0x378d8e6400000000ull) { | |
1797 | // C1 = 10^34 => rounding overflow | |
1798 | C1_hi = 0x0000314dc6448d93ull; | |
1799 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
1800 | y_exp = y_exp + EXP_P1; | |
1801 | } | |
1802 | } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) | |
1803 | && | |
1804 | ((x_sign | |
1805 | && (rnd_mode == ROUNDING_UP | |
1806 | || rnd_mode == ROUNDING_TO_ZERO)) | |
1807 | || (!x_sign | |
1808 | && (rnd_mode == ROUNDING_DOWN | |
1809 | || rnd_mode == ROUNDING_TO_ZERO)))) { | |
1810 | // C1 = C1 - 1 | |
1811 | C1_lo = C1_lo - 1; | |
1812 | if (C1_lo == 0xffffffffffffffffull) | |
1813 | C1_hi--; | |
1814 | // check if we crossed into the lower decade | |
1815 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
1816 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
1817 | C1_lo = 0x378d8e63ffffffffull; | |
1818 | y_exp = y_exp - EXP_P1; | |
1819 | // no underflow, because delta + q2 >= P34 + 1 | |
1820 | } | |
1821 | } else { | |
1822 | ; // exact, the result is already correct | |
1823 | } | |
1824 | } | |
1825 | // assemble the result | |
1826 | res.w[1] = x_sign | y_exp | C1_hi; | |
1827 | res.w[0] = C1_lo; | |
1828 | } | |
1829 | } // end delta = P34 | |
1830 | } else { // if (|delta| <= P34 - 1) | |
1831 | if (delta >= 0) { // if (0 <= delta <= P34 - 1) | |
1832 | if (delta <= P34 - 1 - q2) { | |
1833 | // calculate C' directly; the result is exact | |
1834 | // in this case 1<=q1<=P34-1, 1<=q2<=P34-1 and 0 <= e1-e2 <= P34-2 | |
1835 | // The coefficient of the result is C1 * 10^(e1-e2) + C2 and the | |
1836 | // exponent is e2; either C1 or 10^(e1-e2) may not fit is 64 bits, | |
1837 | // but their product fits with certainty in 128 bits (actually in 113) | |
1838 | scale = delta - q1 + q2; // scale = (int)(e1 >> 49) - (int)(e2 >> 49) | |
1839 | ||
1840 | if (scale >= 20) { // 10^(e1-e2) does not fit in 64 bits, but C1 does | |
1841 | __mul_128x64_to_128 (C1, C1_lo, ten2k128[scale - 20]); | |
1842 | C1_hi = C1.w[1]; | |
1843 | C1_lo = C1.w[0]; | |
1844 | } else if (scale >= 1) { | |
1845 | // if 1 <= scale <= 19 then 10^(e1-e2) fits in 64 bits | |
1846 | if (q1 <= 19) { // C1 fits in 64 bits | |
1847 | __mul_64x64_to_128MACH (C1, C1_lo, ten2k64[scale]); | |
1848 | } else { // q1 >= 20 | |
1849 | C1.w[1] = C1_hi; | |
1850 | C1.w[0] = C1_lo; | |
1851 | __mul_128x64_to_128 (C1, ten2k64[scale], C1); | |
1852 | } | |
1853 | C1_hi = C1.w[1]; | |
1854 | C1_lo = C1.w[0]; | |
1855 | } else { // if (scale == 0) C1 is unchanged | |
1856 | C1.w[0] = C1_lo; // C1.w[1] = C1_hi; | |
1857 | } | |
1858 | // now add C2 | |
1859 | if (x_sign == y_sign) { | |
1860 | // the result cannot overflow | |
1861 | C1_lo = C1_lo + C2_lo; | |
1862 | C1_hi = C1_hi + C2_hi; | |
1863 | if (C1_lo < C1.w[0]) | |
1864 | C1_hi++; | |
1865 | } else { // if x_sign != y_sign | |
1866 | C1_lo = C1_lo - C2_lo; | |
1867 | C1_hi = C1_hi - C2_hi; | |
1868 | if (C1_lo > C1.w[0]) | |
1869 | C1_hi--; | |
1870 | // the result can be zero, but it cannot overflow | |
1871 | if (C1_lo == 0 && C1_hi == 0) { | |
1872 | // assemble the result | |
1873 | if (x_exp < y_exp) | |
1874 | res.w[1] = x_exp; | |
1875 | else | |
1876 | res.w[1] = y_exp; | |
1877 | res.w[0] = 0; | |
1878 | if (rnd_mode == ROUNDING_DOWN) { | |
1879 | res.w[1] |= 0x8000000000000000ull; | |
1880 | } | |
1881 | BID_SWAP128 (res); | |
1882 | BID_RETURN (res); | |
1883 | } | |
1884 | if (C1_hi >= 0x8000000000000000ull) { // negative coefficient! | |
1885 | C1_lo = ~C1_lo; | |
1886 | C1_lo++; | |
1887 | C1_hi = ~C1_hi; | |
1888 | if (C1_lo == 0x0) | |
1889 | C1_hi++; | |
1890 | x_sign = y_sign; // the result will have the sign of y | |
1891 | } | |
1892 | } | |
1893 | // assemble the result | |
1894 | res.w[1] = x_sign | y_exp | C1_hi; | |
1895 | res.w[0] = C1_lo; | |
1896 | } else if (delta == P34 - q2) { | |
1897 | // calculate C' directly; the result may be inexact if it requires | |
1898 | // P34+1 decimal digits; in this case the 'cutoff' point for addition | |
1899 | // is at the position of the lsb of C2, so 0 <= e1-e2 <= P34-1 | |
1900 | // The coefficient of the result is C1 * 10^(e1-e2) + C2 and the | |
1901 | // exponent is e2; either C1 or 10^(e1-e2) may not fit is 64 bits, | |
1902 | // but their product fits with certainty in 128 bits (actually in 113) | |
1903 | scale = delta - q1 + q2; // scale = (int)(e1 >> 49) - (int)(e2 >> 49) | |
1904 | if (scale >= 20) { // 10^(e1-e2) does not fit in 64 bits, but C1 does | |
1905 | __mul_128x64_to_128 (C1, C1_lo, ten2k128[scale - 20]); | |
1906 | } else if (scale >= 1) { | |
1907 | // if 1 <= scale <= 19 then 10^(e1-e2) fits in 64 bits | |
1908 | if (q1 <= 19) { // C1 fits in 64 bits | |
1909 | __mul_64x64_to_128MACH (C1, C1_lo, ten2k64[scale]); | |
1910 | } else { // q1 >= 20 | |
1911 | C1.w[1] = C1_hi; | |
1912 | C1.w[0] = C1_lo; | |
1913 | __mul_128x64_to_128 (C1, ten2k64[scale], C1); | |
1914 | } | |
1915 | } else { // if (scale == 0) C1 is unchanged | |
1916 | C1.w[1] = C1_hi; | |
1917 | C1.w[0] = C1_lo; // only the low part is necessary | |
1918 | } | |
1919 | C1_hi = C1.w[1]; | |
1920 | C1_lo = C1.w[0]; | |
1921 | // now add C2 | |
1922 | if (x_sign == y_sign) { | |
1923 | // the result can overflow! | |
1924 | C1_lo = C1_lo + C2_lo; | |
1925 | C1_hi = C1_hi + C2_hi; | |
1926 | if (C1_lo < C1.w[0]) | |
1927 | C1_hi++; | |
1928 | // test for overflow, possible only when C1 >= 10^34 | |
1929 | if (C1_hi > 0x0001ed09bead87c0ull || (C1_hi == 0x0001ed09bead87c0ull && C1_lo >= 0x378d8e6400000000ull)) { // C1 >= 10^34 | |
1930 | // in this case q = P34 + 1 and x = q - P34 = 1, so multiply | |
1931 | // C'' = C'+ 5 = C1 + 5 by k1 ~ 10^(-1) calculated for P34 + 1 | |
1932 | // decimal digits | |
1933 | // Calculate C'' = C' + 1/2 * 10^x | |
1934 | if (C1_lo >= 0xfffffffffffffffbull) { // low half add has carry | |
1935 | C1_lo = C1_lo + 5; | |
1936 | C1_hi = C1_hi + 1; | |
1937 | } else { | |
1938 | C1_lo = C1_lo + 5; | |
1939 | } | |
1940 | // the approximation of 10^(-1) was rounded up to 118 bits | |
1941 | // 10^(-1) =~ 33333333333333333333333333333400 * 2^-129 | |
1942 | // 10^(-1) =~ 19999999999999999999999999999a00 * 2^-128 | |
1943 | C1.w[1] = C1_hi; | |
1944 | C1.w[0] = C1_lo; // C'' | |
1945 | ten2m1.w[1] = 0x1999999999999999ull; | |
1946 | ten2m1.w[0] = 0x9999999999999a00ull; | |
1947 | __mul_128x128_to_256 (P256, C1, ten2m1); // P256 = C*, f* | |
1948 | // C* is actually floor(C*) in this case | |
1949 | // the top Ex = 128 bits of 10^(-1) are | |
1950 | // T* = 0x00199999999999999999999999999999 | |
1951 | // if (0 < f* < 10^(-x)) then | |
1952 | // if floor(C*) is even then C = floor(C*) - logical right | |
1953 | // shift; C has p decimal digits, correct by Prop. 1) | |
1954 | // else if floor(C*) is odd C = floor(C*) - 1 (logical right | |
1955 | // shift; C has p decimal digits, correct by Pr. 1) | |
1956 | // else | |
1957 | // C = floor(C*) (logical right shift; C has p decimal digits, | |
1958 | // correct by Property 1) | |
1959 | // n = C * 10^(e2+x) | |
1960 | if ((P256.w[1] || P256.w[0]) | |
1961 | && (P256.w[1] < 0x1999999999999999ull | |
1962 | || (P256.w[1] == 0x1999999999999999ull | |
1963 | && P256.w[0] <= 0x9999999999999999ull))) { | |
1964 | // the result is a midpoint | |
1965 | if (P256.w[2] & 0x01) { | |
1966 | is_midpoint_gt_even = 1; | |
1967 | // if floor(C*) is odd C = floor(C*) - 1; the result is not 0 | |
1968 | P256.w[2]--; | |
1969 | if (P256.w[2] == 0xffffffffffffffffull) | |
1970 | P256.w[3]--; | |
1971 | } else { | |
1972 | is_midpoint_lt_even = 1; | |
1973 | } | |
1974 | } | |
1975 | // n = Cstar * 10^(e2+1) | |
1976 | y_exp = y_exp + EXP_P1; | |
1977 | // C* != 10^P because C* has P34 digits | |
1978 | // check for overflow | |
1979 | if (y_exp == EXP_MAX_P1 | |
1980 | && (rnd_mode == ROUNDING_TO_NEAREST | |
1981 | || rnd_mode == ROUNDING_TIES_AWAY)) { | |
1982 | // overflow for RN | |
1983 | res.w[1] = x_sign | 0x7800000000000000ull; // +/-inf | |
1984 | res.w[0] = 0x0ull; | |
1985 | // set the inexact flag | |
1986 | *pfpsf |= INEXACT_EXCEPTION; | |
1987 | // set the overflow flag | |
1988 | *pfpsf |= OVERFLOW_EXCEPTION; | |
1989 | BID_SWAP128 (res); | |
1990 | BID_RETURN (res); | |
1991 | } | |
1992 | // if (0 < f* - 1/2 < 10^(-x)) then | |
1993 | // the result of the addition is exact | |
1994 | // else | |
1995 | // the result of the addition is inexact | |
1996 | if (P256.w[1] > 0x8000000000000000ull || (P256.w[1] == 0x8000000000000000ull && P256.w[0] > 0x0ull)) { // the result may be exact | |
1997 | tmp64 = P256.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
1998 | if ((tmp64 > 0x1999999999999999ull | |
1999 | || (tmp64 == 0x1999999999999999ull | |
2000 | && P256.w[0] >= 0x9999999999999999ull))) { | |
2001 | // set the inexact flag | |
2002 | *pfpsf |= INEXACT_EXCEPTION; | |
2003 | is_inexact = 1; | |
2004 | } // else the result is exact | |
2005 | } else { // the result is inexact | |
2006 | // set the inexact flag | |
2007 | *pfpsf |= INEXACT_EXCEPTION; | |
2008 | is_inexact = 1; | |
2009 | } | |
2010 | C1_hi = P256.w[3]; | |
2011 | C1_lo = P256.w[2]; | |
2012 | if (!is_midpoint_gt_even && !is_midpoint_lt_even) { | |
2013 | is_inexact_lt_midpoint = is_inexact | |
2014 | && (P256.w[1] & 0x8000000000000000ull); | |
2015 | is_inexact_gt_midpoint = is_inexact | |
2016 | && !(P256.w[1] & 0x8000000000000000ull); | |
2017 | } | |
2018 | // general correction from RN to RA, RM, RP, RZ; | |
2019 | // result uses y_exp | |
2020 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
2021 | if ((!x_sign | |
2022 | && | |
2023 | ((rnd_mode == ROUNDING_UP | |
2024 | && is_inexact_lt_midpoint) | |
2025 | || | |
2026 | ((rnd_mode == ROUNDING_TIES_AWAY | |
2027 | || rnd_mode == ROUNDING_UP) | |
2028 | && is_midpoint_gt_even))) || (x_sign | |
2029 | && | |
2030 | ((rnd_mode == | |
2031 | ROUNDING_DOWN | |
2032 | && | |
2033 | is_inexact_lt_midpoint) | |
2034 | || | |
2035 | ((rnd_mode == | |
2036 | ROUNDING_TIES_AWAY | |
2037 | || rnd_mode == | |
2038 | ROUNDING_DOWN) | |
2039 | && | |
2040 | is_midpoint_gt_even)))) | |
2041 | { | |
2042 | // C1 = C1 + 1 | |
2043 | C1_lo = C1_lo + 1; | |
2044 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
2045 | C1_hi = C1_hi + 1; | |
2046 | } | |
2047 | if (C1_hi == 0x0001ed09bead87c0ull | |
2048 | && C1_lo == 0x378d8e6400000000ull) { | |
2049 | // C1 = 10^34 => rounding overflow | |
2050 | C1_hi = 0x0000314dc6448d93ull; | |
2051 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
2052 | y_exp = y_exp + EXP_P1; | |
2053 | } | |
2054 | } else | |
2055 | if ((is_midpoint_lt_even || is_inexact_gt_midpoint) | |
2056 | && | |
2057 | ((x_sign | |
2058 | && (rnd_mode == ROUNDING_UP | |
2059 | || rnd_mode == ROUNDING_TO_ZERO)) | |
2060 | || (!x_sign | |
2061 | && (rnd_mode == ROUNDING_DOWN | |
2062 | || rnd_mode == ROUNDING_TO_ZERO)))) { | |
2063 | // C1 = C1 - 1 | |
2064 | C1_lo = C1_lo - 1; | |
2065 | if (C1_lo == 0xffffffffffffffffull) | |
2066 | C1_hi--; | |
2067 | // check if we crossed into the lower decade | |
2068 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
2069 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
2070 | C1_lo = 0x378d8e63ffffffffull; | |
2071 | y_exp = y_exp - EXP_P1; | |
2072 | // no underflow, because delta + q2 >= P34 + 1 | |
2073 | } | |
2074 | } else { | |
2075 | ; // exact, the result is already correct | |
2076 | } | |
2077 | // in all cases check for overflow (RN and RA solved already) | |
2078 | if (y_exp == EXP_MAX_P1) { // overflow | |
2079 | if ((rnd_mode == ROUNDING_DOWN && x_sign) || // RM and res < 0 | |
2080 | (rnd_mode == ROUNDING_UP && !x_sign)) { // RP and res > 0 | |
2081 | C1_hi = 0x7800000000000000ull; // +inf | |
2082 | C1_lo = 0x0ull; | |
2083 | } else { // RM and res > 0, RP and res < 0, or RZ | |
2084 | C1_hi = 0x5fffed09bead87c0ull; | |
2085 | C1_lo = 0x378d8e63ffffffffull; | |
2086 | } | |
2087 | y_exp = 0; // x_sign is preserved | |
2088 | // set the inexact flag (in case the exact addition was exact) | |
2089 | *pfpsf |= INEXACT_EXCEPTION; | |
2090 | // set the overflow flag | |
2091 | *pfpsf |= OVERFLOW_EXCEPTION; | |
2092 | } | |
2093 | } | |
2094 | } // else if (C1 < 10^34) then C1 is the coeff.; the result is exact | |
2095 | } else { // if x_sign != y_sign the result is exact | |
2096 | C1_lo = C1_lo - C2_lo; | |
2097 | C1_hi = C1_hi - C2_hi; | |
2098 | if (C1_lo > C1.w[0]) | |
2099 | C1_hi--; | |
2100 | // the result can be zero, but it cannot overflow | |
2101 | if (C1_lo == 0 && C1_hi == 0) { | |
2102 | // assemble the result | |
2103 | if (x_exp < y_exp) | |
2104 | res.w[1] = x_exp; | |
2105 | else | |
2106 | res.w[1] = y_exp; | |
2107 | res.w[0] = 0; | |
2108 | if (rnd_mode == ROUNDING_DOWN) { | |
2109 | res.w[1] |= 0x8000000000000000ull; | |
2110 | } | |
2111 | BID_SWAP128 (res); | |
2112 | BID_RETURN (res); | |
2113 | } | |
2114 | if (C1_hi >= 0x8000000000000000ull) { // negative coefficient! | |
2115 | C1_lo = ~C1_lo; | |
2116 | C1_lo++; | |
2117 | C1_hi = ~C1_hi; | |
2118 | if (C1_lo == 0x0) | |
2119 | C1_hi++; | |
2120 | x_sign = y_sign; // the result will have the sign of y | |
2121 | } | |
2122 | } | |
2123 | // assemble the result | |
2124 | res.w[1] = x_sign | y_exp | C1_hi; | |
2125 | res.w[0] = C1_lo; | |
2126 | } else { // if (delta >= P34 + 1 - q2) | |
2127 | // instead of C' = (C1 * 10^(e1-e2) + C2)rnd,P34 | |
2128 | // calculate C' = C1 * 10^(e1-e2-x1) + (C2 * 10^(-x1))rnd,P34 | |
2129 | // where x1 = q1 + e1 - e2 - P34, 1 <= x1 <= P34 - 1 | |
2130 | // In most cases C' will have P34 digits, and n = C' * 10^(e2+x1) | |
2131 | // If the result has P34+1 digits, redo the steps above with x1+1 | |
2132 | // If the result has P34-1 digits or less, redo the steps above with | |
2133 | // x1-1 but only if initially x1 >= 1 | |
2134 | // NOTE: these two steps can be improved, e.g we could guess if | |
2135 | // P34+1 or P34-1 digits will be obtained by adding/subtracting just | |
2136 | // the top 64 bits of the two operands | |
2137 | // The result cannot be zero, but it can overflow | |
2138 | x1 = delta + q2 - P34; // 1 <= x1 <= P34-1 | |
2139 | roundC2: | |
2140 | // Calculate C1 * 10^(e1-e2-x1) where 0 <= e1-e2-x1 <= P34 - 1 | |
2141 | // scale = (int)(e1 >> 49) - (int)(e2 >> 49) - x1; 0 <= scale <= P34-1 | |
2142 | scale = delta - q1 + q2 - x1; // scale = e1 - e2 - x1 = P34 - q1 | |
2143 | // either C1 or 10^(e1-e2-x1) may not fit is 64 bits, | |
2144 | // but their product fits with certainty in 128 bits (actually in 113) | |
2145 | if (scale >= 20) { //10^(e1-e2-x1) doesn't fit in 64 bits, but C1 does | |
2146 | __mul_128x64_to_128 (C1, C1_lo, ten2k128[scale - 20]); | |
2147 | } else if (scale >= 1) { | |
2148 | // if 1 <= scale <= 19 then 10^(e1-e2-x1) fits in 64 bits | |
2149 | if (q1 <= 19) { // C1 fits in 64 bits | |
2150 | __mul_64x64_to_128MACH (C1, C1_lo, ten2k64[scale]); | |
2151 | } else { // q1 >= 20 | |
2152 | C1.w[1] = C1_hi; | |
2153 | C1.w[0] = C1_lo; | |
2154 | __mul_128x64_to_128 (C1, ten2k64[scale], C1); | |
2155 | } | |
2156 | } else { // if (scale == 0) C1 is unchanged | |
2157 | C1.w[1] = C1_hi; | |
2158 | C1.w[0] = C1_lo; | |
2159 | } | |
2160 | tmp64 = C1.w[0]; // C1.w[1], C1.w[0] contains C1 * 10^(e1-e2-x1) | |
2161 | ||
2162 | // now round C2 to q2-x1 decimal digits, where 1<=x1<=q2-1<=P34-1 | |
2163 | // (but if we got here a second time after x1 = x1 - 1, then | |
2164 | // x1 >= 0; note that for x1 = 0 C2 is unchanged) | |
2165 | // C2' = C2 + 1/2 * 10^x1 = C2 + 5 * 10^(x1-1) | |
2166 | ind = x1 - 1; // 0 <= ind <= q2-2<=P34-2=32; but note that if x1 = 0 | |
2167 | // during a second pass, then ind = -1 | |
2168 | if (ind >= 0) { // if (x1 >= 1) | |
2169 | C2.w[0] = C2_lo; | |
2170 | C2.w[1] = C2_hi; | |
2171 | if (ind <= 18) { | |
2172 | C2.w[0] = C2.w[0] + midpoint64[ind]; | |
2173 | if (C2.w[0] < C2_lo) | |
2174 | C2.w[1]++; | |
2175 | } else { // 19 <= ind <= 32 | |
2176 | C2.w[0] = C2.w[0] + midpoint128[ind - 19].w[0]; | |
2177 | C2.w[1] = C2.w[1] + midpoint128[ind - 19].w[1]; | |
2178 | if (C2.w[0] < C2_lo) | |
2179 | C2.w[1]++; | |
2180 | } | |
2181 | // the approximation of 10^(-x1) was rounded up to 118 bits | |
2182 | __mul_128x128_to_256 (R256, C2, ten2mk128[ind]); // R256 = C2*, f2* | |
2183 | // calculate C2* and f2* | |
2184 | // C2* is actually floor(C2*) in this case | |
2185 | // C2* and f2* need shifting and masking, as shown by | |
2186 | // shiftright128[] and maskhigh128[] | |
2187 | // the top Ex bits of 10^(-x1) are T* = ten2mk128trunc[ind], e.g. | |
2188 | // if x1=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
2189 | // if (0 < f2* < 10^(-x1)) then | |
2190 | // if floor(C1+C2*) is even then C2* = floor(C2*) - logical right | |
2191 | // shift; C2* has p decimal digits, correct by Prop. 1) | |
2192 | // else if floor(C1+C2*) is odd C2* = floor(C2*)-1 (logical right | |
2193 | // shift; C2* has p decimal digits, correct by Pr. 1) | |
2194 | // else | |
2195 | // C2* = floor(C2*) (logical right shift; C has p decimal digits, | |
2196 | // correct by Property 1) | |
2197 | // n = C2* * 10^(e2+x1) | |
2198 | ||
2199 | if (ind <= 2) { | |
2200 | highf2star.w[1] = 0x0; | |
2201 | highf2star.w[0] = 0x0; // low f2* ok | |
2202 | } else if (ind <= 21) { | |
2203 | highf2star.w[1] = 0x0; | |
2204 | highf2star.w[0] = R256.w[2] & maskhigh128[ind]; // low f2* ok | |
2205 | } else { | |
2206 | highf2star.w[1] = R256.w[3] & maskhigh128[ind]; | |
2207 | highf2star.w[0] = R256.w[2]; // low f2* is ok | |
2208 | } | |
2209 | // shift right C2* by Ex-128 = shiftright128[ind] | |
2210 | if (ind >= 3) { | |
2211 | shift = shiftright128[ind]; | |
2212 | if (shift < 64) { // 3 <= shift <= 63 | |
2213 | R256.w[2] = | |
2214 | (R256.w[2] >> shift) | (R256.w[3] << (64 - shift)); | |
2215 | R256.w[3] = (R256.w[3] >> shift); | |
2216 | } else { // 66 <= shift <= 102 | |
2217 | R256.w[2] = (R256.w[3] >> (shift - 64)); | |
2218 | R256.w[3] = 0x0ULL; | |
2219 | } | |
2220 | } | |
2221 | if (second_pass) { | |
2222 | is_inexact_lt_midpoint = 0; | |
2223 | is_inexact_gt_midpoint = 0; | |
2224 | is_midpoint_lt_even = 0; | |
2225 | is_midpoint_gt_even = 0; | |
2226 | } | |
2227 | // determine inexactness of the rounding of C2* (this may be | |
2228 | // followed by a second rounding only if we get P34+1 | |
2229 | // decimal digits) | |
2230 | // if (0 < f2* - 1/2 < 10^(-x1)) then | |
2231 | // the result is exact | |
2232 | // else (if f2* - 1/2 > T* then) | |
2233 | // the result of is inexact | |
2234 | if (ind <= 2) { | |
2235 | if (R256.w[1] > 0x8000000000000000ull || | |
2236 | (R256.w[1] == 0x8000000000000000ull | |
2237 | && R256.w[0] > 0x0ull)) { | |
2238 | // f2* > 1/2 and the result may be exact | |
2239 | tmp64A = R256.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
2240 | if ((tmp64A > ten2mk128trunc[ind].w[1] | |
2241 | || (tmp64A == ten2mk128trunc[ind].w[1] | |
2242 | && R256.w[0] >= ten2mk128trunc[ind].w[0]))) { | |
2243 | // set the inexact flag | |
2244 | // *pfpsf |= INEXACT_EXCEPTION; | |
2245 | tmp_inexact = 1; // may be set again during a second pass | |
2246 | // this rounding is applied to C2 only! | |
2247 | if (x_sign == y_sign) | |
2248 | is_inexact_lt_midpoint = 1; | |
2249 | else // if (x_sign != y_sign) | |
2250 | is_inexact_gt_midpoint = 1; | |
2251 | } // else the result is exact | |
2252 | // rounding down, unless a midpoint in [ODD, EVEN] | |
2253 | } else { // the result is inexact; f2* <= 1/2 | |
2254 | // set the inexact flag | |
2255 | // *pfpsf |= INEXACT_EXCEPTION; | |
2256 | tmp_inexact = 1; // just in case we will round a second time | |
2257 | // rounding up, unless a midpoint in [EVEN, ODD] | |
2258 | // this rounding is applied to C2 only! | |
2259 | if (x_sign == y_sign) | |
2260 | is_inexact_gt_midpoint = 1; | |
2261 | else // if (x_sign != y_sign) | |
2262 | is_inexact_lt_midpoint = 1; | |
2263 | } | |
2264 | } else if (ind <= 21) { // if 3 <= ind <= 21 | |
2265 | if (highf2star.w[1] > 0x0 || (highf2star.w[1] == 0x0 | |
2266 | && highf2star.w[0] > | |
2267 | onehalf128[ind]) | |
2268 | || (highf2star.w[1] == 0x0 | |
2269 | && highf2star.w[0] == onehalf128[ind] | |
2270 | && (R256.w[1] || R256.w[0]))) { | |
2271 | // f2* > 1/2 and the result may be exact | |
2272 | // Calculate f2* - 1/2 | |
2273 | tmp64A = highf2star.w[0] - onehalf128[ind]; | |
2274 | tmp64B = highf2star.w[1]; | |
2275 | if (tmp64A > highf2star.w[0]) | |
2276 | tmp64B--; | |
2277 | if (tmp64B || tmp64A | |
2278 | || R256.w[1] > ten2mk128trunc[ind].w[1] | |
2279 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
2280 | && R256.w[0] > ten2mk128trunc[ind].w[0])) { | |
2281 | // set the inexact flag | |
2282 | // *pfpsf |= INEXACT_EXCEPTION; | |
2283 | tmp_inexact = 1; // may be set again during a second pass | |
2284 | // this rounding is applied to C2 only! | |
2285 | if (x_sign == y_sign) | |
2286 | is_inexact_lt_midpoint = 1; | |
2287 | else // if (x_sign != y_sign) | |
2288 | is_inexact_gt_midpoint = 1; | |
2289 | } // else the result is exact | |
2290 | } else { // the result is inexact; f2* <= 1/2 | |
2291 | // set the inexact flag | |
2292 | // *pfpsf |= INEXACT_EXCEPTION; | |
2293 | tmp_inexact = 1; // may be set again during a second pass | |
2294 | // rounding up, unless a midpoint in [EVEN, ODD] | |
2295 | // this rounding is applied to C2 only! | |
2296 | if (x_sign == y_sign) | |
2297 | is_inexact_gt_midpoint = 1; | |
2298 | else // if (x_sign != y_sign) | |
2299 | is_inexact_lt_midpoint = 1; | |
2300 | } | |
2301 | } else { // if 22 <= ind <= 33 | |
2302 | if (highf2star.w[1] > onehalf128[ind] | |
2303 | || (highf2star.w[1] == onehalf128[ind] | |
2304 | && (highf2star.w[0] || R256.w[1] | |
2305 | || R256.w[0]))) { | |
2306 | // f2* > 1/2 and the result may be exact | |
2307 | // Calculate f2* - 1/2 | |
2308 | // tmp64A = highf2star.w[0]; | |
2309 | tmp64B = highf2star.w[1] - onehalf128[ind]; | |
2310 | if (tmp64B || highf2star.w[0] | |
2311 | || R256.w[1] > ten2mk128trunc[ind].w[1] | |
2312 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
2313 | && R256.w[0] > ten2mk128trunc[ind].w[0])) { | |
2314 | // set the inexact flag | |
2315 | // *pfpsf |= INEXACT_EXCEPTION; | |
2316 | tmp_inexact = 1; // may be set again during a second pass | |
2317 | // this rounding is applied to C2 only! | |
2318 | if (x_sign == y_sign) | |
2319 | is_inexact_lt_midpoint = 1; | |
2320 | else // if (x_sign != y_sign) | |
2321 | is_inexact_gt_midpoint = 1; | |
2322 | } // else the result is exact | |
2323 | } else { // the result is inexact; f2* <= 1/2 | |
2324 | // set the inexact flag | |
2325 | // *pfpsf |= INEXACT_EXCEPTION; | |
2326 | tmp_inexact = 1; // may be set again during a second pass | |
2327 | // rounding up, unless a midpoint in [EVEN, ODD] | |
2328 | // this rounding is applied to C2 only! | |
2329 | if (x_sign == y_sign) | |
2330 | is_inexact_gt_midpoint = 1; | |
2331 | else // if (x_sign != y_sign) | |
2332 | is_inexact_lt_midpoint = 1; | |
2333 | } | |
2334 | } | |
2335 | // check for midpoints | |
2336 | if ((R256.w[1] || R256.w[0]) && (highf2star.w[1] == 0) | |
2337 | && (highf2star.w[0] == 0) | |
2338 | && (R256.w[1] < ten2mk128trunc[ind].w[1] | |
2339 | || (R256.w[1] == ten2mk128trunc[ind].w[1] | |
2340 | && R256.w[0] <= ten2mk128trunc[ind].w[0]))) { | |
2341 | // the result is a midpoint | |
2342 | if ((tmp64 + R256.w[2]) & 0x01) { // MP in [EVEN, ODD] | |
2343 | // if floor(C2*) is odd C = floor(C2*) - 1; the result may be 0 | |
2344 | R256.w[2]--; | |
2345 | if (R256.w[2] == 0xffffffffffffffffull) | |
2346 | R256.w[3]--; | |
2347 | // this rounding is applied to C2 only! | |
2348 | if (x_sign == y_sign) | |
2349 | is_midpoint_gt_even = 1; | |
2350 | else // if (x_sign != y_sign) | |
2351 | is_midpoint_lt_even = 1; | |
2352 | is_inexact_lt_midpoint = 0; | |
2353 | is_inexact_gt_midpoint = 0; | |
2354 | } else { | |
2355 | // else MP in [ODD, EVEN] | |
2356 | // this rounding is applied to C2 only! | |
2357 | if (x_sign == y_sign) | |
2358 | is_midpoint_lt_even = 1; | |
2359 | else // if (x_sign != y_sign) | |
2360 | is_midpoint_gt_even = 1; | |
2361 | is_inexact_lt_midpoint = 0; | |
2362 | is_inexact_gt_midpoint = 0; | |
2363 | } | |
2364 | } | |
2365 | // end if (ind >= 0) | |
2366 | } else { // if (ind == -1); only during a 2nd pass, and when x1 = 0 | |
2367 | R256.w[2] = C2_lo; | |
2368 | R256.w[3] = C2_hi; | |
2369 | tmp_inexact = 0; | |
2370 | // to correct a possible setting to 1 from 1st pass | |
2371 | if (second_pass) { | |
2372 | is_midpoint_lt_even = 0; | |
2373 | is_midpoint_gt_even = 0; | |
2374 | is_inexact_lt_midpoint = 0; | |
2375 | is_inexact_gt_midpoint = 0; | |
2376 | } | |
2377 | } | |
2378 | // and now add/subtract C1 * 10^(e1-e2-x1) +/- (C2 * 10^(-x1))rnd,P34 | |
2379 | if (x_sign == y_sign) { // addition; could overflow | |
2380 | // no second pass is possible this way (only for x_sign != y_sign) | |
2381 | C1.w[0] = C1.w[0] + R256.w[2]; | |
2382 | C1.w[1] = C1.w[1] + R256.w[3]; | |
2383 | if (C1.w[0] < tmp64) | |
2384 | C1.w[1]++; // carry | |
2385 | // if the sum has P34+1 digits, i.e. C1>=10^34 redo the calculation | |
2386 | // with x1=x1+1 | |
2387 | if (C1.w[1] > 0x0001ed09bead87c0ull || (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] >= 0x378d8e6400000000ull)) { // C1 >= 10^34 | |
2388 | // chop off one more digit from the sum, but make sure there is | |
2389 | // no double-rounding error (see table - double rounding logic) | |
2390 | // now round C1 from P34+1 to P34 decimal digits | |
2391 | // C1' = C1 + 1/2 * 10 = C1 + 5 | |
2392 | if (C1.w[0] >= 0xfffffffffffffffbull) { // low half add has carry | |
2393 | C1.w[0] = C1.w[0] + 5; | |
2394 | C1.w[1] = C1.w[1] + 1; | |
2395 | } else { | |
2396 | C1.w[0] = C1.w[0] + 5; | |
2397 | } | |
2398 | // the approximation of 10^(-1) was rounded up to 118 bits | |
2399 | __mul_128x128_to_256 (Q256, C1, ten2mk128[0]); // Q256 = C1*, f1* | |
2400 | // C1* is actually floor(C1*) in this case | |
2401 | // the top 128 bits of 10^(-1) are | |
2402 | // T* = ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
2403 | // if (0 < f1* < 10^(-1)) then | |
2404 | // if floor(C1*) is even then C1* = floor(C1*) - logical right | |
2405 | // shift; C1* has p decimal digits, correct by Prop. 1) | |
2406 | // else if floor(C1*) is odd C1* = floor(C1*) - 1 (logical right | |
2407 | // shift; C1* has p decimal digits, correct by Pr. 1) | |
2408 | // else | |
2409 | // C1* = floor(C1*) (logical right shift; C has p decimal digits | |
2410 | // correct by Property 1) | |
2411 | // n = C1* * 10^(e2+x1+1) | |
2412 | if ((Q256.w[1] || Q256.w[0]) | |
2413 | && (Q256.w[1] < ten2mk128trunc[0].w[1] | |
2414 | || (Q256.w[1] == ten2mk128trunc[0].w[1] | |
2415 | && Q256.w[0] <= ten2mk128trunc[0].w[0]))) { | |
2416 | // the result is a midpoint | |
2417 | if (is_inexact_lt_midpoint) { // for the 1st rounding | |
2418 | is_inexact_gt_midpoint = 1; | |
2419 | is_inexact_lt_midpoint = 0; | |
2420 | is_midpoint_gt_even = 0; | |
2421 | is_midpoint_lt_even = 0; | |
2422 | } else if (is_inexact_gt_midpoint) { // for the 1st rounding | |
2423 | Q256.w[2]--; | |
2424 | if (Q256.w[2] == 0xffffffffffffffffull) | |
2425 | Q256.w[3]--; | |
2426 | is_inexact_gt_midpoint = 0; | |
2427 | is_inexact_lt_midpoint = 1; | |
2428 | is_midpoint_gt_even = 0; | |
2429 | is_midpoint_lt_even = 0; | |
2430 | } else if (is_midpoint_gt_even) { // for the 1st rounding | |
2431 | // Note: cannot have is_midpoint_lt_even | |
2432 | is_inexact_gt_midpoint = 0; | |
2433 | is_inexact_lt_midpoint = 1; | |
2434 | is_midpoint_gt_even = 0; | |
2435 | is_midpoint_lt_even = 0; | |
2436 | } else { // the first rounding must have been exact | |
2437 | if (Q256.w[2] & 0x01) { // MP in [EVEN, ODD] | |
2438 | // the truncated result is correct | |
2439 | Q256.w[2]--; | |
2440 | if (Q256.w[2] == 0xffffffffffffffffull) | |
2441 | Q256.w[3]--; | |
2442 | is_inexact_gt_midpoint = 0; | |
2443 | is_inexact_lt_midpoint = 0; | |
2444 | is_midpoint_gt_even = 1; | |
2445 | is_midpoint_lt_even = 0; | |
2446 | } else { // MP in [ODD, EVEN] | |
2447 | is_inexact_gt_midpoint = 0; | |
2448 | is_inexact_lt_midpoint = 0; | |
2449 | is_midpoint_gt_even = 0; | |
2450 | is_midpoint_lt_even = 1; | |
2451 | } | |
2452 | } | |
2453 | tmp_inexact = 1; // in all cases | |
2454 | } else { // the result is not a midpoint | |
2455 | // determine inexactness of the rounding of C1 (the sum C1+C2*) | |
2456 | // if (0 < f1* - 1/2 < 10^(-1)) then | |
2457 | // the result is exact | |
2458 | // else (if f1* - 1/2 > T* then) | |
2459 | // the result of is inexact | |
2460 | // ind = 0 | |
2461 | if (Q256.w[1] > 0x8000000000000000ull | |
2462 | || (Q256.w[1] == 0x8000000000000000ull | |
2463 | && Q256.w[0] > 0x0ull)) { | |
2464 | // f1* > 1/2 and the result may be exact | |
2465 | Q256.w[1] = Q256.w[1] - 0x8000000000000000ull; // f1* - 1/2 | |
2466 | if ((Q256.w[1] > ten2mk128trunc[0].w[1] | |
2467 | || (Q256.w[1] == ten2mk128trunc[0].w[1] | |
2468 | && Q256.w[0] > ten2mk128trunc[0].w[0]))) { | |
2469 | is_inexact_gt_midpoint = 0; | |
2470 | is_inexact_lt_midpoint = 1; | |
2471 | is_midpoint_gt_even = 0; | |
2472 | is_midpoint_lt_even = 0; | |
2473 | // set the inexact flag | |
2474 | tmp_inexact = 1; | |
2475 | // *pfpsf |= INEXACT_EXCEPTION; | |
2476 | } else { // else the result is exact for the 2nd rounding | |
2477 | if (tmp_inexact) { // if the previous rounding was inexact | |
2478 | if (is_midpoint_lt_even) { | |
2479 | is_inexact_gt_midpoint = 1; | |
2480 | is_midpoint_lt_even = 0; | |
2481 | } else if (is_midpoint_gt_even) { | |
2482 | is_inexact_lt_midpoint = 1; | |
2483 | is_midpoint_gt_even = 0; | |
2484 | } else { | |
2485 | ; // no change | |
2486 | } | |
2487 | } | |
2488 | } | |
2489 | // rounding down, unless a midpoint in [ODD, EVEN] | |
2490 | } else { // the result is inexact; f1* <= 1/2 | |
2491 | is_inexact_gt_midpoint = 1; | |
2492 | is_inexact_lt_midpoint = 0; | |
2493 | is_midpoint_gt_even = 0; | |
2494 | is_midpoint_lt_even = 0; | |
2495 | // set the inexact flag | |
2496 | tmp_inexact = 1; | |
2497 | // *pfpsf |= INEXACT_EXCEPTION; | |
2498 | } | |
2499 | } // end 'the result is not a midpoint' | |
2500 | // n = C1 * 10^(e2+x1) | |
2501 | C1.w[1] = Q256.w[3]; | |
2502 | C1.w[0] = Q256.w[2]; | |
2503 | y_exp = y_exp + ((UINT64) (x1 + 1) << 49); | |
2504 | } else { // C1 < 10^34 | |
2505 | // C1.w[1] and C1.w[0] already set | |
2506 | // n = C1 * 10^(e2+x1) | |
2507 | y_exp = y_exp + ((UINT64) x1 << 49); | |
2508 | } | |
2509 | // check for overflow | |
2510 | if (y_exp == EXP_MAX_P1 | |
2511 | && (rnd_mode == ROUNDING_TO_NEAREST | |
2512 | || rnd_mode == ROUNDING_TIES_AWAY)) { | |
2513 | res.w[1] = 0x7800000000000000ull | x_sign; // +/-inf | |
2514 | res.w[0] = 0x0ull; | |
2515 | // set the inexact flag | |
2516 | *pfpsf |= INEXACT_EXCEPTION; | |
2517 | // set the overflow flag | |
2518 | *pfpsf |= OVERFLOW_EXCEPTION; | |
2519 | BID_SWAP128 (res); | |
2520 | BID_RETURN (res); | |
2521 | } // else no overflow | |
2522 | } else { // if x_sign != y_sign the result of this subtract. is exact | |
2523 | C1.w[0] = C1.w[0] - R256.w[2]; | |
2524 | C1.w[1] = C1.w[1] - R256.w[3]; | |
2525 | if (C1.w[0] > tmp64) | |
2526 | C1.w[1]--; // borrow | |
2527 | if (C1.w[1] >= 0x8000000000000000ull) { // negative coefficient! | |
2528 | C1.w[0] = ~C1.w[0]; | |
2529 | C1.w[0]++; | |
2530 | C1.w[1] = ~C1.w[1]; | |
2531 | if (C1.w[0] == 0x0) | |
2532 | C1.w[1]++; | |
2533 | tmp_sign = y_sign; | |
2534 | // the result will have the sign of y if last rnd | |
2535 | } else { | |
2536 | tmp_sign = x_sign; | |
2537 | } | |
2538 | // if the difference has P34-1 digits or less, i.e. C1 < 10^33 then | |
2539 | // redo the calculation with x1=x1-1; | |
2540 | // redo the calculation also if C1 = 10^33 and | |
2541 | // (is_inexact_gt_midpoint or is_midpoint_lt_even); | |
2542 | // (the last part should have really been | |
2543 | // (is_inexact_lt_midpoint or is_midpoint_gt_even) from | |
2544 | // the rounding of C2, but the position flags have been reversed) | |
2545 | // 10^33 = 0x0000314dc6448d93 0x38c15b0a00000000 | |
2546 | if ((C1.w[1] < 0x0000314dc6448d93ull || (C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] < 0x38c15b0a00000000ull)) || (C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b0a00000000ull && (is_inexact_gt_midpoint || is_midpoint_lt_even))) { // C1=10^33 | |
2547 | x1 = x1 - 1; // x1 >= 0 | |
2548 | if (x1 >= 0) { | |
2549 | // clear position flags and tmp_inexact | |
2550 | is_midpoint_lt_even = 0; | |
2551 | is_midpoint_gt_even = 0; | |
2552 | is_inexact_lt_midpoint = 0; | |
2553 | is_inexact_gt_midpoint = 0; | |
2554 | tmp_inexact = 0; | |
2555 | second_pass = 1; | |
2556 | goto roundC2; // else result has less than P34 digits | |
2557 | } | |
2558 | } | |
2559 | // if the coefficient of the result is 10^34 it means that this | |
2560 | // must be the second pass, and we are done | |
2561 | if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34 | |
2562 | C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33 | |
2563 | C1.w[0] = 0x38c15b0a00000000ull; | |
2564 | y_exp = y_exp + ((UINT64) 1 << 49); | |
2565 | } | |
2566 | x_sign = tmp_sign; | |
2567 | if (x1 >= 1) | |
2568 | y_exp = y_exp + ((UINT64) x1 << 49); | |
2569 | // x1 = -1 is possible at the end of a second pass when the | |
2570 | // first pass started with x1 = 1 | |
2571 | } | |
2572 | C1_hi = C1.w[1]; | |
2573 | C1_lo = C1.w[0]; | |
2574 | // general correction from RN to RA, RM, RP, RZ; result uses y_exp | |
2575 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
2576 | if ((!x_sign | |
2577 | && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) | |
2578 | || | |
2579 | ((rnd_mode == ROUNDING_TIES_AWAY | |
2580 | || rnd_mode == ROUNDING_UP) | |
2581 | && is_midpoint_gt_even))) || (x_sign | |
2582 | && | |
2583 | ((rnd_mode == | |
2584 | ROUNDING_DOWN | |
2585 | && | |
2586 | is_inexact_lt_midpoint) | |
2587 | || | |
2588 | ((rnd_mode == | |
2589 | ROUNDING_TIES_AWAY | |
2590 | || rnd_mode == | |
2591 | ROUNDING_DOWN) | |
2592 | && | |
2593 | is_midpoint_gt_even)))) | |
2594 | { | |
2595 | // C1 = C1 + 1 | |
2596 | C1_lo = C1_lo + 1; | |
2597 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
2598 | C1_hi = C1_hi + 1; | |
2599 | } | |
2600 | if (C1_hi == 0x0001ed09bead87c0ull | |
2601 | && C1_lo == 0x378d8e6400000000ull) { | |
2602 | // C1 = 10^34 => rounding overflow | |
2603 | C1_hi = 0x0000314dc6448d93ull; | |
2604 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
2605 | y_exp = y_exp + EXP_P1; | |
2606 | } | |
2607 | } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) | |
2608 | && | |
2609 | ((x_sign | |
2610 | && (rnd_mode == ROUNDING_UP | |
2611 | || rnd_mode == ROUNDING_TO_ZERO)) | |
2612 | || (!x_sign | |
2613 | && (rnd_mode == ROUNDING_DOWN | |
2614 | || rnd_mode == ROUNDING_TO_ZERO)))) { | |
2615 | // C1 = C1 - 1 | |
2616 | C1_lo = C1_lo - 1; | |
2617 | if (C1_lo == 0xffffffffffffffffull) | |
2618 | C1_hi--; | |
2619 | // check if we crossed into the lower decade | |
2620 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
2621 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
2622 | C1_lo = 0x378d8e63ffffffffull; | |
2623 | y_exp = y_exp - EXP_P1; | |
2624 | // no underflow, because delta + q2 >= P34 + 1 | |
2625 | } | |
2626 | } else { | |
2627 | ; // exact, the result is already correct | |
2628 | } | |
2629 | // in all cases check for overflow (RN and RA solved already) | |
2630 | if (y_exp == EXP_MAX_P1) { // overflow | |
2631 | if ((rnd_mode == ROUNDING_DOWN && x_sign) || // RM and res < 0 | |
2632 | (rnd_mode == ROUNDING_UP && !x_sign)) { // RP and res > 0 | |
2633 | C1_hi = 0x7800000000000000ull; // +inf | |
2634 | C1_lo = 0x0ull; | |
2635 | } else { // RM and res > 0, RP and res < 0, or RZ | |
2636 | C1_hi = 0x5fffed09bead87c0ull; | |
2637 | C1_lo = 0x378d8e63ffffffffull; | |
2638 | } | |
2639 | y_exp = 0; // x_sign is preserved | |
2640 | // set the inexact flag (in case the exact addition was exact) | |
2641 | *pfpsf |= INEXACT_EXCEPTION; | |
2642 | // set the overflow flag | |
2643 | *pfpsf |= OVERFLOW_EXCEPTION; | |
2644 | } | |
2645 | } | |
2646 | // assemble the result | |
2647 | res.w[1] = x_sign | y_exp | C1_hi; | |
2648 | res.w[0] = C1_lo; | |
2649 | if (tmp_inexact) | |
2650 | *pfpsf |= INEXACT_EXCEPTION; | |
2651 | } | |
2652 | } else { // if (-P34 + 1 <= delta <= -1) <=> 1 <= -delta <= P34 - 1 | |
2653 | // NOTE: the following, up to "} else { // if x_sign != y_sign | |
2654 | // the result is exact" is identical to "else if (delta == P34 - q2) {" | |
2655 | // from above; also, the code is not symmetric: a+b and b+a may take | |
2656 | // different paths (need to unify eventually!) | |
2657 | // calculate C' = C2 + C1 * 10^(e1-e2) directly; the result may be | |
2658 | // inexact if it requires P34 + 1 decimal digits; in either case the | |
2659 | // 'cutoff' point for addition is at the position of the lsb of C2 | |
2660 | // The coefficient of the result is C1 * 10^(e1-e2) + C2 and the | |
2661 | // exponent is e2; either C1 or 10^(e1-e2) may not fit is 64 bits, | |
2662 | // but their product fits with certainty in 128 bits (actually in 113) | |
2663 | // Note that 0 <= e1 - e2 <= P34 - 2 | |
2664 | // -P34 + 1 <= delta <= -1 <=> -P34 + 1 <= delta <= -1 <=> | |
2665 | // -P34 + 1 <= q1 + e1 - q2 - e2 <= -1 <=> | |
2666 | // q2 - q1 - P34 + 1 <= e1 - e2 <= q2 - q1 - 1 <=> | |
2667 | // 1 - P34 - P34 + 1 <= e1-e2 <= P34 - 1 - 1 => 0 <= e1-e2 <= P34 - 2 | |
2668 | scale = delta - q1 + q2; // scale = (int)(e1 >> 49) - (int)(e2 >> 49) | |
2669 | if (scale >= 20) { // 10^(e1-e2) does not fit in 64 bits, but C1 does | |
2670 | __mul_128x64_to_128 (C1, C1_lo, ten2k128[scale - 20]); | |
2671 | } else if (scale >= 1) { | |
2672 | // if 1 <= scale <= 19 then 10^(e1-e2) fits in 64 bits | |
2673 | if (q1 <= 19) { // C1 fits in 64 bits | |
2674 | __mul_64x64_to_128MACH (C1, C1_lo, ten2k64[scale]); | |
2675 | } else { // q1 >= 20 | |
2676 | C1.w[1] = C1_hi; | |
2677 | C1.w[0] = C1_lo; | |
2678 | __mul_128x64_to_128 (C1, ten2k64[scale], C1); | |
2679 | } | |
2680 | } else { // if (scale == 0) C1 is unchanged | |
2681 | C1.w[1] = C1_hi; | |
2682 | C1.w[0] = C1_lo; // only the low part is necessary | |
2683 | } | |
2684 | C1_hi = C1.w[1]; | |
2685 | C1_lo = C1.w[0]; | |
2686 | // now add C2 | |
2687 | if (x_sign == y_sign) { | |
2688 | // the result can overflow! | |
2689 | C1_lo = C1_lo + C2_lo; | |
2690 | C1_hi = C1_hi + C2_hi; | |
2691 | if (C1_lo < C1.w[0]) | |
2692 | C1_hi++; | |
2693 | // test for overflow, possible only when C1 >= 10^34 | |
2694 | if (C1_hi > 0x0001ed09bead87c0ull || (C1_hi == 0x0001ed09bead87c0ull && C1_lo >= 0x378d8e6400000000ull)) { // C1 >= 10^34 | |
2695 | // in this case q = P34 + 1 and x = q - P34 = 1, so multiply | |
2696 | // C'' = C'+ 5 = C1 + 5 by k1 ~ 10^(-1) calculated for P34 + 1 | |
2697 | // decimal digits | |
2698 | // Calculate C'' = C' + 1/2 * 10^x | |
2699 | if (C1_lo >= 0xfffffffffffffffbull) { // low half add has carry | |
2700 | C1_lo = C1_lo + 5; | |
2701 | C1_hi = C1_hi + 1; | |
2702 | } else { | |
2703 | C1_lo = C1_lo + 5; | |
2704 | } | |
2705 | // the approximation of 10^(-1) was rounded up to 118 bits | |
2706 | // 10^(-1) =~ 33333333333333333333333333333400 * 2^-129 | |
2707 | // 10^(-1) =~ 19999999999999999999999999999a00 * 2^-128 | |
2708 | C1.w[1] = C1_hi; | |
2709 | C1.w[0] = C1_lo; // C'' | |
2710 | ten2m1.w[1] = 0x1999999999999999ull; | |
2711 | ten2m1.w[0] = 0x9999999999999a00ull; | |
2712 | __mul_128x128_to_256 (P256, C1, ten2m1); // P256 = C*, f* | |
2713 | // C* is actually floor(C*) in this case | |
2714 | // the top Ex = 128 bits of 10^(-1) are | |
2715 | // T* = 0x00199999999999999999999999999999 | |
2716 | // if (0 < f* < 10^(-x)) then | |
2717 | // if floor(C*) is even then C = floor(C*) - logical right | |
2718 | // shift; C has p decimal digits, correct by Prop. 1) | |
2719 | // else if floor(C*) is odd C = floor(C*) - 1 (logical right | |
2720 | // shift; C has p decimal digits, correct by Pr. 1) | |
2721 | // else | |
2722 | // C = floor(C*) (logical right shift; C has p decimal digits, | |
2723 | // correct by Property 1) | |
2724 | // n = C * 10^(e2+x) | |
2725 | if ((P256.w[1] || P256.w[0]) | |
2726 | && (P256.w[1] < 0x1999999999999999ull | |
2727 | || (P256.w[1] == 0x1999999999999999ull | |
2728 | && P256.w[0] <= 0x9999999999999999ull))) { | |
2729 | // the result is a midpoint | |
2730 | if (P256.w[2] & 0x01) { | |
2731 | is_midpoint_gt_even = 1; | |
2732 | // if floor(C*) is odd C = floor(C*) - 1; the result is not 0 | |
2733 | P256.w[2]--; | |
2734 | if (P256.w[2] == 0xffffffffffffffffull) | |
2735 | P256.w[3]--; | |
2736 | } else { | |
2737 | is_midpoint_lt_even = 1; | |
2738 | } | |
2739 | } | |
2740 | // n = Cstar * 10^(e2+1) | |
2741 | y_exp = y_exp + EXP_P1; | |
2742 | // C* != 10^P34 because C* has P34 digits | |
2743 | // check for overflow | |
2744 | if (y_exp == EXP_MAX_P1 | |
2745 | && (rnd_mode == ROUNDING_TO_NEAREST | |
2746 | || rnd_mode == ROUNDING_TIES_AWAY)) { | |
2747 | // overflow for RN | |
2748 | res.w[1] = x_sign | 0x7800000000000000ull; // +/-inf | |
2749 | res.w[0] = 0x0ull; | |
2750 | // set the inexact flag | |
2751 | *pfpsf |= INEXACT_EXCEPTION; | |
2752 | // set the overflow flag | |
2753 | *pfpsf |= OVERFLOW_EXCEPTION; | |
2754 | BID_SWAP128 (res); | |
2755 | BID_RETURN (res); | |
2756 | } | |
2757 | // if (0 < f* - 1/2 < 10^(-x)) then | |
2758 | // the result of the addition is exact | |
2759 | // else | |
2760 | // the result of the addition is inexact | |
2761 | if (P256.w[1] > 0x8000000000000000ull || (P256.w[1] == 0x8000000000000000ull && P256.w[0] > 0x0ull)) { // the result may be exact | |
2762 | tmp64 = P256.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
2763 | if ((tmp64 > 0x1999999999999999ull | |
2764 | || (tmp64 == 0x1999999999999999ull | |
2765 | && P256.w[0] >= 0x9999999999999999ull))) { | |
2766 | // set the inexact flag | |
2767 | *pfpsf |= INEXACT_EXCEPTION; | |
2768 | is_inexact = 1; | |
2769 | } // else the result is exact | |
2770 | } else { // the result is inexact | |
2771 | // set the inexact flag | |
2772 | *pfpsf |= INEXACT_EXCEPTION; | |
2773 | is_inexact = 1; | |
2774 | } | |
2775 | C1_hi = P256.w[3]; | |
2776 | C1_lo = P256.w[2]; | |
2777 | if (!is_midpoint_gt_even && !is_midpoint_lt_even) { | |
2778 | is_inexact_lt_midpoint = is_inexact | |
2779 | && (P256.w[1] & 0x8000000000000000ull); | |
2780 | is_inexact_gt_midpoint = is_inexact | |
2781 | && !(P256.w[1] & 0x8000000000000000ull); | |
2782 | } | |
2783 | // general correction from RN to RA, RM, RP, RZ; result uses y_exp | |
2784 | if (rnd_mode != ROUNDING_TO_NEAREST) { | |
2785 | if ((!x_sign | |
2786 | && ((rnd_mode == ROUNDING_UP | |
2787 | && is_inexact_lt_midpoint) | |
2788 | || ((rnd_mode == ROUNDING_TIES_AWAY | |
2789 | || rnd_mode == ROUNDING_UP) | |
2790 | && is_midpoint_gt_even))) || (x_sign | |
2791 | && | |
2792 | ((rnd_mode == | |
2793 | ROUNDING_DOWN | |
2794 | && | |
2795 | is_inexact_lt_midpoint) | |
2796 | || | |
2797 | ((rnd_mode == | |
2798 | ROUNDING_TIES_AWAY | |
2799 | || rnd_mode | |
2800 | == | |
2801 | ROUNDING_DOWN) | |
2802 | && | |
2803 | is_midpoint_gt_even)))) | |
2804 | { | |
2805 | // C1 = C1 + 1 | |
2806 | C1_lo = C1_lo + 1; | |
2807 | if (C1_lo == 0) { // rounding overflow in the low 64 bits | |
2808 | C1_hi = C1_hi + 1; | |
2809 | } | |
2810 | if (C1_hi == 0x0001ed09bead87c0ull | |
2811 | && C1_lo == 0x378d8e6400000000ull) { | |
2812 | // C1 = 10^34 => rounding overflow | |
2813 | C1_hi = 0x0000314dc6448d93ull; | |
2814 | C1_lo = 0x38c15b0a00000000ull; // 10^33 | |
2815 | y_exp = y_exp + EXP_P1; | |
2816 | } | |
2817 | } else | |
2818 | if ((is_midpoint_lt_even || is_inexact_gt_midpoint) && | |
2819 | ((x_sign && (rnd_mode == ROUNDING_UP || | |
2820 | rnd_mode == ROUNDING_TO_ZERO)) || | |
2821 | (!x_sign && (rnd_mode == ROUNDING_DOWN || | |
2822 | rnd_mode == ROUNDING_TO_ZERO)))) { | |
2823 | // C1 = C1 - 1 | |
2824 | C1_lo = C1_lo - 1; | |
2825 | if (C1_lo == 0xffffffffffffffffull) | |
2826 | C1_hi--; | |
2827 | // check if we crossed into the lower decade | |
2828 | if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) { // 10^33 - 1 | |
2829 | C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 | |
2830 | C1_lo = 0x378d8e63ffffffffull; | |
2831 | y_exp = y_exp - EXP_P1; | |
2832 | // no underflow, because delta + q2 >= P34 + 1 | |
2833 | } | |
2834 | } else { | |
2835 | ; // exact, the result is already correct | |
2836 | } | |
2837 | // in all cases check for overflow (RN and RA solved already) | |
2838 | if (y_exp == EXP_MAX_P1) { // overflow | |
2839 | if ((rnd_mode == ROUNDING_DOWN && x_sign) || // RM and res < 0 | |
2840 | (rnd_mode == ROUNDING_UP && !x_sign)) { // RP and res > 0 | |
2841 | C1_hi = 0x7800000000000000ull; // +inf | |
2842 | C1_lo = 0x0ull; | |
2843 | } else { // RM and res > 0, RP and res < 0, or RZ | |
2844 | C1_hi = 0x5fffed09bead87c0ull; | |
2845 | C1_lo = 0x378d8e63ffffffffull; | |
2846 | } | |
2847 | y_exp = 0; // x_sign is preserved | |
2848 | // set the inexact flag (in case the exact addition was exact) | |
2849 | *pfpsf |= INEXACT_EXCEPTION; | |
2850 | // set the overflow flag | |
2851 | *pfpsf |= OVERFLOW_EXCEPTION; | |
2852 | } | |
2853 | } | |
2854 | } // else if (C1 < 10^34) then C1 is the coeff.; the result is exact | |
2855 | // assemble the result | |
2856 | res.w[1] = x_sign | y_exp | C1_hi; | |
2857 | res.w[0] = C1_lo; | |
2858 | } else { // if x_sign != y_sign the result is exact | |
2859 | C1_lo = C2_lo - C1_lo; | |
2860 | C1_hi = C2_hi - C1_hi; | |
2861 | if (C1_lo > C2_lo) | |
2862 | C1_hi--; | |
2863 | if (C1_hi >= 0x8000000000000000ull) { // negative coefficient! | |
2864 | C1_lo = ~C1_lo; | |
2865 | C1_lo++; | |
2866 | C1_hi = ~C1_hi; | |
2867 | if (C1_lo == 0x0) | |
2868 | C1_hi++; | |
2869 | x_sign = y_sign; // the result will have the sign of y | |
2870 | } | |
2871 | // the result can be zero, but it cannot overflow | |
2872 | if (C1_lo == 0 && C1_hi == 0) { | |
2873 | // assemble the result | |
2874 | if (x_exp < y_exp) | |
2875 | res.w[1] = x_exp; | |
2876 | else | |
2877 | res.w[1] = y_exp; | |
2878 | res.w[0] = 0; | |
2879 | if (rnd_mode == ROUNDING_DOWN) { | |
2880 | res.w[1] |= 0x8000000000000000ull; | |
2881 | } | |
2882 | BID_SWAP128 (res); | |
2883 | BID_RETURN (res); | |
2884 | } | |
2885 | // assemble the result | |
2886 | res.w[1] = y_sign | y_exp | C1_hi; | |
2887 | res.w[0] = C1_lo; | |
2888 | } | |
200359e8 L |
2889 | } |
2890 | } | |
b2a00c89 | 2891 | BID_SWAP128 (res); |
200359e8 L |
2892 | BID_RETURN (res) |
2893 | } | |
2894 | } | |
2895 | ||
b2a00c89 L |
2896 | |
2897 | ||
2898 | // bid128_sub stands for bid128qq_sub | |
2899 | ||
200359e8 L |
2900 | /***************************************************************************** |
2901 | * BID128 sub | |
2902 | ****************************************************************************/ | |
2903 | ||
2904 | #if DECIMAL_CALL_BY_REFERENCE | |
2905 | void | |
b2a00c89 L |
2906 | bid128_sub (UINT128 * pres, UINT128 * px, UINT128 * py |
2907 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
2908 | _EXC_INFO_PARAM) { | |
200359e8 L |
2909 | UINT128 x = *px, y = *py; |
2910 | #if !DECIMAL_GLOBAL_ROUNDING | |
2911 | unsigned int rnd_mode = *prnd_mode; | |
2912 | #endif | |
2913 | #else | |
2914 | UINT128 | |
b2a00c89 L |
2915 | bid128_sub (UINT128 x, UINT128 y |
2916 | _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
2917 | _EXC_INFO_PARAM) { | |
200359e8 L |
2918 | #endif |
2919 | ||
2920 | UINT128 res; | |
2921 | UINT64 y_sign; | |
2922 | ||
b2a00c89 | 2923 | if ((y.w[HIGH_128W] & MASK_NAN) != MASK_NAN) { // y is not NAN |
200359e8 | 2924 | // change its sign |
b2a00c89 | 2925 | y_sign = y.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 | 2926 | if (y_sign) |
b2a00c89 | 2927 | y.w[HIGH_128W] = y.w[HIGH_128W] & 0x7fffffffffffffffull; |
200359e8 | 2928 | else |
b2a00c89 | 2929 | y.w[HIGH_128W] = y.w[HIGH_128W] | 0x8000000000000000ull; |
200359e8 L |
2930 | } |
2931 | #if DECIMAL_CALL_BY_REFERENCE | |
b2a00c89 L |
2932 | bid128_add (&res, &x, &y |
2933 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
2934 | _EXC_INFO_ARG); | |
200359e8 | 2935 | #else |
b2a00c89 L |
2936 | res = bid128_add (x, y |
2937 | _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG | |
2938 | _EXC_INFO_ARG); | |
200359e8 L |
2939 | #endif |
2940 | BID_RETURN (res); | |
2941 | } |